Bubble elements for a mesh with orphan nodes

[hgpeterson]

Hi all,

I came across this issue with the newly-added bubble elements for a mesh that contains an "orphan" node. I am using Julia 1.12.2 and Gridap 0.19.6. Here is a MWE:

using Gridap
using GridapGmsh
using Gmsh

# make simple 2D mesh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("mesh")
gmsh.model.geo.addPoint(0, 0, 0)
gmsh.model.geo.addPoint(1, 2, 0)  # control point
gmsh.model.geo.addPoint(2, 0, 0) 
gmsh.model.geo.addBezier([1, 2, 3])
gmsh.model.geo.addLine(3, 1)
gmsh.model.geo.addCurveLoop(1:2, 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
# gmsh.option.setNumber("Mesh.SaveWithoutOrphans", 1)
gmsh.write("mesh.msh")
gmsh.finalize()

# load mesh in Gridap
model = GmshDiscreteModel("mesh.msh")

# ✅ define lagrangian element space on mesh
reffe_l = ReferenceFE(lagrangian, Float64, 1)
V = TestFESpace(model, reffe_l)
U = TrialFESpace(V)

# ❌ define bubble element space on mesh
reffe_b = ReferenceFE(bubble, Float64)
R = TestFESpace(model, reffe_b)
B = TrialFESpace(R)

The first part makes a simple 2D mesh that contains a Bezier curve. Making this curve creates an "orphan" node that is not actually part of the mesh but is used as a control point. The lagrangian elements in Gridap don't seem to mind this node, but when I try to define bubble elements, I get the following error:

ERROR: LoadError: BoundsError: attempt to access 3-element Vector{Int64} at index [0]
Stacktrace:
  [1] throw_boundserror(A::Vector{Int64}, I::Tuple{Int64})
    @ Base ./essentials.jl:15
  [2] getindex
    @ ./essentials.jl:919 [inlined]
  [3] getindex
    @ ./abstractarray.jl:1345 [inlined]
  [4] _generate_face_to_own_dofs_count_d!(face_to_own_dofs_ptrs::Vector{…}, offset::Int64, cell_to_ctype::Vector{…}, dface_to_cell_owner::Gridap.Arrays.LocalItemFromTable{…}, dface_to_ldface::Gridap.Arrays.LocalIndexFromTable{…}, ctype_to_ldface_to_num_own_ldofs::Vector{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/ConformingFESpaces.jl:374
  [5] _generate_face_to_own_dofs(n_faces::Int64, cell_to_ctype::Vector{…}, d_to_cell_to_dfaces::Vector{…}, d_to_dface_to_cells::Vector{…}, d_to_offset::Vector{…}, d_to_ctype_to_ldface_to_own_ldofs::Vector{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/ConformingFESpaces.jl:342
  [6] compute_conforming_cell_dofs(cell_fe::Gridap.FESpaces.CellFE{…}, cell_conformity::Gridap.FESpaces.CellConformity{…}, grid_topology::Gridap.Geometry.UnstructuredGridTopology{…}, face_labeling::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{…}, dirichlet_components::Nothing)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/ConformingFESpaces.jl:263
  [7] _ConformingFESpace(vector_type::Type, model::Gridap.Geometry.UnstructuredDiscreteModel{…}, face_labeling::Gridap.Geometry.FaceLabeling, cell_fe::Gridap.FESpaces.CellFE{…}, dirichlet_tags::Vector{…}, dirichlet_components::Nothing, trian::Gridap.Geometry.BodyFittedTriangulation{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/ConformingFESpaces.jl:175
  [8] #FESpace#33
    @ ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:23 [inlined]
  [9] FESpace
    @ ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:8 [inlined]
 [10] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{…}, cell_reffe::Gridap.Arrays.CompressedArray{…}; conformity::Nothing, trian::Gridap.Geometry.BodyFittedTriangulation{…}, labels::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{…}, dirichlet_masks::Nothing, constraint::Nothing, vector_type::Nothing)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:99
 [11] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{…}, cell_reffe::Gridap.Arrays.CompressedArray{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:65
 [12] #FESpace#35
    @ ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:114 [inlined]
 [13] FESpace
    @ ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:110 [inlined]
 [14] #TestFESpace#37
    @ ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:126 [inlined]
 [15] TestFESpace(::Gridap.Geometry.UnstructuredDiscreteModel{…}, ::Tuple{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/dk0DA/src/FESpaces/FESpaceFactories.jl:125
...

The issue is resolved if the line

# gmsh.option.setNumber("Mesh.SaveWithoutOrphans", 1)

is uncommented in the MWE. It would be nice if this wasn't necessary, as with the lagrangian elements.

[DimhamT]

Hi, @hgpeterson .
I think the deeper root cause is related to the construction of the conforming space based on the mesh you provided. For example, you can reproduce the same error with:

reffe_l = ReferenceFE(lagrangian, Float64, 0)
V = TestFESpace(model, reffe_l)

This probably needs further consideration. For now, the workaround is to use L2, i.e.

using Gridap.ReferenceFEs
ReferenceFEs.Conformity(::GenericRefFE{Bubble}, ::Symbol) = L2Conformity()

model = GmshDiscreteModel("mesh.msh")
reffe_b = ReferenceFE(bubble, Float64)
R = TestFESpace(model, reffe_b, conformity=:L2)

[DimhamT]

Is this control point required? It can also cause other issues. For example, in writevtk, this point does not belong to any cell, so its type cannot be updated, which eventually triggers an invalid index (0) access.

[JordiManyer]

Why are bubbles not discontinuous by default? Isnt that the definition of bubbles?

[JordiManyer]

Why are bubbles not discontinuous by default? Isnt that the definition of bubbles?

[DimhamT]

The Bubble space is indeed $C^0$, but from the programming perspective, its degrees of freedom are of the PointValue type located at the barycenter of each element. Therefore, I believe using an :L2 space should not be an issue. The real problem here lies in the presence of nodes in the mesh (*.msh) file that do not belong to any element. If _ConformingFESpace is used, it will cause errors. As a temporary workaround, _DiscontinuousFESpace is employed instead.

[JordiManyer]

Therefore, I believe using an :L2 space should not be an issue.

Yes, I agree. What I was saying is: why is this not the default? I don't think bubbles can be anything but L2. So I would just change this in the code.

[DimhamT]

The default is :L2 conformity, because

julia> reffe = BubbleRefFE(Float64, TRI)
GenericRefFE()

julia> ReferenceFEs.Conformity(reffe, nothing)
L2Conformity()

However, when constructing the FESpace, if :L2 is not explicitly specified, it will internally call FESpace->_ConformingFESpace. This implementation runs into problems when handling the .msh file provided in the current issue (I’m not sure whether this is a problem with the code implementation or with the .msh file itself).
The above approach simply bypasses the problematic implementation by explicitly specifying :L2 to force a call to _DiscontinuousFESpace. This is only a temporary workaround, not the final solution.

[JordiManyer]

Ok I see. I will think about this.... there is something about that FESpace constructor I don't like. Anyway, I guess just use the workaround.

[hgpeterson]

Thank you for the reply, @DimhamT. Your workaround

using Gridap.ReferenceFEs
ReferenceFEs.Conformity(::GenericRefFE{Bubble}, ::Symbol) = L2Conformity()

model = GmshDiscreteModel("mesh.msh")
reffe_b = ReferenceFE(bubble, Float64)
R = TestFESpace(model, reffe_b, conformity=:L2)

does seem to do the trick. I will also note that this workaround is also needed for meshes generated with the RecombineAll flag turned on in Gmsh to make a quadrilateral mesh. This might suggest that the problem goes beyond just orphan nodes (unless Gmsh is still using those under-the-hood for quad meshes?). I can post an example if there's interest.

Is this control point required? It can also cause other issues. For example, in writevtk, this point does not belong to any cell, so its type cannot be updated, which eventually triggers an invalid index (0) access.

I've never run into issues with writevtk for meshes generated using control points for bezier curves. It's is true that they aren't required for the actual solution of PDEs, though; they're just needed for the mesh generation step. In that sense, it's probably best to use the SaveWithoutOrphans flag anyways.

[JordiManyer]

the RecombineAll flag turned on

what does that do, exactly? I am not that familiar with gmsh

[hgpeterson]

By default, Gmsh meshes with triangles. To use quadrangles instead, you can Recombine. Here's an example:

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("mesh")
gmsh.model.geo.addPoint(0, 0, 0)
gmsh.model.geo.addPoint(3, 0, 0)
gmsh.model.geo.addPoint(3, 3, 0) 
gmsh.model.geo.addPoint(0, 3, 0) 
gmsh.model.geo.addLine(1, 2)
gmsh.model.geo.addLine(2, 3)
gmsh.model.geo.addLine(3, 4)
gmsh.model.geo.addLine(4, 1)
gmsh.model.geo.addCurveLoop(1:4, 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(0, [1], 1, "tag_1")
gmsh.model.addPhysicalGroup(0, [2], 2, "dirichlet")
gmsh.model.addPhysicalGroup(0, 3:4, 3, "neumann")
gmsh.model.addPhysicalGroup(1, [1], 2, "dirichlet")
gmsh.model.addPhysicalGroup(1, 2:4, 3, "neumann")
gmsh.model.addPhysicalGroup(2, [1], 4, "interior")
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 1)  # resolution
gmsh.option.setNumber("Mesh.RecombineAll", 1)  # rectangles
gmsh.option.setNumber("Mesh.Algorithm", 8)  # rectangles
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()

I am using this mesh to attempt to reproduce the test in tests/GridapTests/StokesMiniTests.jl using a mesh generated by Gmsh instead of Gridap (although so far it is failing...)

[DimhamT]

Hi, @hgpeterson .

I've never run into issues with writevtk for meshes generated using control points for bezier curves.

Have you tried the following code?

model = DiscreteModelFromFile("mesh.msh")
writevtk(model, "model")

Oh my side, this throws an error

ERROR: BoundsError: attempt to access 1-element Vector{Int64} at index [0]
Stacktrace:
  [1] getindex
    @ ./essentials.jl:13 [inlined]
  [2] getindex
    @ ./abstractarray.jl:1293 [inlined]
  [3] _generate_vtk_cells(cache::Gridap.Arrays.CachedVector{…}, cell_to_nodes::Gridap.Arrays.Table{…}, cell_to_type::Vector{…}, type_to_vtkid::Vector{…}, type_to_vtknodes::Vector{…})
...

[DimhamT]

I am using this mesh to attempt to reproduce the test in tests/GridapTests/StokesMiniTests.jl using a mesh generated by Gmsh instead of Gridap (although so far it is failing...)

test/GridapTests/StokesMiniTests.jl is only intended to test whether the mini element can be solved exactly when $u \in B^3 + P_1$. It is a specially constructed example. The bubble function is defined based on the domain $(0,1)^2$. Therefore, if you use the domain $(0,3)^2$, you need to construct a different b((x,y)). If you set the domain to $(0,1)^2$, the test should pass.

[hgpeterson]

Oh my side, this throws an error

Same here. I think I see what's going on (and why I haven't had this problem in the past). Once you define "Physical Groups" in Gmsh, the errors (both for defining bubble elements and using writevtk) go away. So the following mesh generation code fixes the MWE (without the need for the conformity workaround):

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("mesh")
gmsh.model.geo.addPoint(0, 0, 0)
gmsh.model.geo.addPoint(1, 2, 0)  # control point
gmsh.model.geo.addPoint(2, 0, 0) 
gmsh.model.geo.addBezier([1, 2, 3])
gmsh.model.geo.addLine(3, 1)
gmsh.model.geo.addCurveLoop(1:2, 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(0, [1,3], 1, "boundary")
gmsh.model.addPhysicalGroup(1, [1,2], 1, "boundary")
gmsh.model.addPhysicalGroup(2, [1], 1, "boundary")
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()

However, the conformity workaround is still needed for the quad mesh using Recombine.

test/GridapTests/StokesMiniTests.jl is only intended to test whether the mini element can be solved exactly when u ∈ B 3 + P 1 . It is a specially constructed example. The bubble function is defined based on the domain ( 0 , 1 ) 2 . Therefore, if you use the domain ( 0 , 3 ) 2 , you need to construct a different b((x,y)). If you set the domain to ( 0 , 1 ) 2 , the test should pass.

From what I can tell, the line

domain = (0, 3, 0, 3)

here would imply that the test is constructed for the domain $(0, 3)^2$. My version with a Gmsh-generated mesh still fails either way, but it might just be user error.