Gridap.jl icon indicating copy to clipboard operation
Gridap.jl copied to clipboard

changing map to lazy_map to fix issue #905

Open kishore-nori opened this issue 1 year ago • 7 comments

Changing map to lazy_map to fix issue #905. This is done to return identical CellPoint objects. Tests on my local machine have passed.

The issue was get_cell_point(dΩ), when Measure is built from a GmshDiscreteModel based on QUADs, returns CellPoint objects which are non-identical i.e. (get_cell_points(dΩ) === get_cell_points(dΩ) is not satisfied). This is due to non-identical cell_phys_point field in CellPoint. This results in the issue #905, where CellState constructed based on doesn't satisfy the compatibility test during integration.

Haven't included a test, as I ll then have to include the .msh file as well, or may be I can add a test to GridapGmsh.jl pakcage, which ever seems reasonable..

kishore-nori avatar May 30 '23 09:05 kishore-nori

Hav't included a test, as I ll then have to include the .msh file as well, or may be I can add a test to GridapGmsh.jl pakcage, which ever seems reasonable..

Can the issue be reproduced with a manually built UnstructuredDiscreteModel of quadrilaterals? see, e.g., for an example snippet https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384635611

amartinhuertas avatar May 30 '23 10:05 amartinhuertas

(get_cell_points(dΩ) === get_cell_points(dΩ) is not satisfied).

How can it be that this condition is not satistied? Isn't dΩ the same object in your comment?

amartinhuertas avatar May 30 '23 10:05 amartinhuertas

Hav't included a test, as I ll then have to include the .msh file as well, or may be I can add a test to GridapGmsh.jl pakcage, which ever seems reasonable..

Can the issue be reproduced with a manually built UnstructuredDiscreteModel of quadrilaterals? see, e.g., for an example snippet gridap/GridapP4est.jl#21 (comment)

I tried building UnstructuredDiscreteModel on top of CartesiandDiscreteModel but it doesn't fail in that case. I ll try out manual building and get back

kishore-nori avatar May 30 '23 10:05 kishore-nori

(get_cell_points(dΩ) === get_cell_points(dΩ) is not satisfied).

How can it be that this condition is not satistied? Isn't dΩ the same object in your comment?

No.. that is exactly the source of error for this issue, and this PR fixes it. Even if is the same, get_cell_points(dΩ) === get_cell_points(dΩ) don't end up same.. due to non identical cell_phys_point which somewhere deep inside uses map instead of lazy_map, this is only peculiar to Quad mesh with GmshDiscreteModel

kishore-nori avatar May 30 '23 10:05 kishore-nori

I tried building UnstructuredDiscreteModel on top of CartesiandDiscreteModel but it doesn't fail in that case. I ll try out manual building and get back

Don't use the example provided verbatim. In the example, I am generating a non-conforming mesh, which is currently not well supported by Gridap. You can, e.g., extract the submesh corresponding to the first four cells in the example

amartinhuertas avatar May 30 '23 10:05 amartinhuertas

Codecov Report

Merging #906 (f077bd9) into master (122da8b) will increase coverage by 0.06%. The diff coverage is 100.00%.

:exclamation: Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

@@            Coverage Diff             @@
##           master     #906      +/-   ##
==========================================
+ Coverage   88.51%   88.58%   +0.06%     
==========================================
  Files         173      173              
  Lines       20269    20540     +271     
==========================================
+ Hits        17942    18196     +254     
- Misses       2327     2344      +17     
Impacted Files Coverage Δ
src/Arrays/CompressedArrays.jl 94.00% <100.00%> (ø)

... and 2 files with indirect coverage changes

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

codecov-commenter avatar May 30 '23 10:05 codecov-commenter

see, e.g., for an example snippet https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384635611

Thank you for the pointer.. now the issue is reproducible without a .msh file :) and I also checked that the PR fixes it. See following, also including the issue #905 example. Where should I put the below test? GridapTests?

using Gridap
using GridapGmsh
using FillArrays

ptr  = [ 1, 5, 9, 13, 17 ]
data = [ 1,2,4,5, 2,3,5,6, 4,5,7,8, 5,6,8,9 ]
cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
node_coordinates = Vector{Point{2,Float64}}(undef,11)

for j in 1:3
  for i in 1:3
    node_coordinates[3*(j-1)+i]=Point{2,Float64}((i-1)/2.0,(j-1)/2.0)
  end
end

polytope = QUAD
scalar_reffe = Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
cell_types = collect(Fill(1,length(cell_vertex_lids)))
cell_reffes = [scalar_reffe]
grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                        cell_vertex_lids,
                                        cell_reffes,
                                        cell_types,
                                        Gridap.Geometry.NonOriented())
model = Gridap.Geometry.UnstructuredDiscreteModel(grid)

Ω = Triangulation(model)
dΩ = Measure(Ω,2)

# the problem is below:
dΩ.quad.cell_point === dΩ.quad.cell_point # true
dΩ.quad.cell_point === get_cell_points(dΩ) # false
get_cell_points(dΩ) === get_cell_points(dΩ) # false, the source of the issue #905

# issue #905 example

order = 1
function project(order,Ω)
  degree = 2*order
  dΩ = Measure(Ω,degree)
  # f(x) = sin(norm(x))
  f = CellState(1.0,dΩ)
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,reffe,conformity=:L2)
  a(u,v) = ∫( u*v )dΩ
  l(v) = ∫( v*f )*dΩ
  op = AffineFEOperator(a,l,V,V)
  fh = solve(op)
  fh
end

fh = project(order,Ω) # fails

ERROR:

It is not possible to evaluate the given CellState on the given CellPoint.

a CellState can only be evaluated at the CellPoint it was created from.
If you want to evaluate at another location, you would need first to project the CellState
to a FESpace (e.g. via a L2 projection).

kishore-nori avatar May 30 '23 10:05 kishore-nori

Where should I put the below test? GridapTests?

Yes. Let us put it in that directory. issue_905.jl file.

amartinhuertas avatar Apr 05 '24 01:04 amartinhuertas