Error with quadrilateral mesh generated from Gmsh + Lowest-order Nedelec FEs

[amartinhuertas]

The error is reproduced with MWE below (mesh file attached).

Please see my explanation below on the cause of the problem.

Versions of packages:

• GridapGmsh v0.6.0
• Gridap v0.17.13 https://github.com/gridap/Gridap.jl#master (46b5388)

refined_lshape.zip

using Gridap
using GridapGmsh
mesh_file="refined_lshape.msh"
order = 0
reffe = ReferenceFE(nedelec, Float64, order)
model = GmshDiscreteModel(mesh_file)
writevtk(model,"m")
V0 = FESpace(model, reffe; dirichlet_tags=["dirichlet_surfaces", "dirichlet_points"])

And the error is:

ERROR: BoundsError: attempt to access 1-element Vector{Vector{Int64}} at index [2]
Stacktrace:
  [1] getindex
    @ ./array.jl:861 [inlined]
  [2] getindex
    @ ./abstractarray.jl:1221 [inlined]
  [3] getindex!(cache::Gridap.Arrays.CachedVector{Int32, Vector{Int32}}, a::Gridap.FESpaces.CellDofsNonOriented, cell::Int64)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/ConformingFESpaces.jl:555
  [4] _generate_data_and_ptrs_fill_ptrs!(ptrs::Vector{Int32}, vv::Gridap.FESpaces.CellDofsNonOriented)
    @ Gridap.Arrays ~/.julia/packages/Gridap/JCqMt/src/Arrays/Tables.jl:151
  [5] generate_data_and_ptrs(vv::Gridap.FESpaces.CellDofsNonOriented)
    @ Gridap.Arrays ~/.julia/packages/Gridap/JCqMt/src/Arrays/Tables.jl:139
  [6] Table
    @ ~/.julia/packages/Gridap/JCqMt/src/Arrays/Tables.jl:27 [inlined]
  [7] compute_conforming_cell_dofs(cell_fe::Gridap.FESpaces.CellFE{Vector{Int8}}, cell_conformity::Gridap.FESpaces.CellConformity{Vector{Int8}}, grid_topology::Gridap.Geometry.UnstructuredGridTopology{2, 2, Float64, Gridap.Geometry.NonOriented}, face_labeling::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{String}, dirichlet_components::Nothing)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/ConformingFESpaces.jl:269
  [8] _ConformingFESpace(vector_type::Type, model::Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, face_labeling::Gridap.Geometry.FaceLabeling, cell_fe::Gridap.FESpaces.CellFE{Vector{Int8}}, dirichlet_tags::Vector{String}, dirichlet_components::Nothing, trian::Gridap.Geometry.BodyFittedTriangulation{2, 2, Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, Gridap.Geometry.UnstructuredGrid{2, 2, Float64, Gridap.Geometry.NonOriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/ConformingFESpaces.jl:148
  [9] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, cell_fe::Gridap.FESpaces.CellFE{Vector{Int8}}; trian::Gridap.Geometry.BodyFittedTriangulation{2, 2, Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, Gridap.Geometry.UnstructuredGrid{2, 2, Float64, Gridap.Geometry.NonOriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}, labels::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{String}, dirichlet_masks::Nothing, constraint::Nothing, vector_type::Type)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/FESpaceFactories.jl:23
 [10] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, cell_reffe::Gridap.Arrays.CompressedArray{Gridap.ReferenceFEs.GenericRefFE{Nedelec, 2}, 1, Vector{Gridap.ReferenceFEs.GenericRefFE{Nedelec, 2}}, Vector{Int8}}; conformity::Nothing, trian::Gridap.Geometry.BodyFittedTriangulation{2, 2, Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, Gridap.Geometry.UnstructuredGrid{2, 2, Float64, Gridap.Geometry.NonOriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}, labels::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{String}, dirichlet_masks::Nothing, constraint::Nothing, vector_type::Nothing)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/FESpaceFactories.jl:111
 [11] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{2, 2, Float64, Gridap.Geometry.NonOriented}, reffe::Tuple{Nedelec, Tuple{DataType, Int64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; kwargs::Base.Pairs{Symbol, Vector{String}, Tuple{Symbol}, NamedTuple{(:dirichlet_tags,), Tuple{Vector{String}}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JCqMt/src/FESpaces/FESpaceFactories.jl:126
 [12] top-level scope
    @ ~/git-repos/Gridap.jl/mwe.jl:10

Explanation

As far as I understand GridapGmsh is returning a quadrilateral mesh which is non-oriented

And the current implementation of Nedelec FEs on quadrilaterals does not seem to support non-oriented quadrilateral meshes

Let me explain how I have arrived to that conclusion.

If you execute Table(get_cell_permutations(grid_topology)) with grid_topology being the grid_topology object related to the model

you will get a cell-wise array of arrays, where the array associated to each cell is of size 9 (1..4 vertices, 5..8 edges, 9 cell interior)

the relevant portion of these arrays is the one corresponding to the edges, i.e., positions 5..8 of these arrays

the two possible orientations of an edge within a cell are coded with 1 and 2 resp.

if there is at least one edge for which the orientations of the edge from the perspective of both cells around do not match, i.e., if 1 versus 2 or 2 versus 1, then the mesh is a non-oriented mesh

If you look at the array returned by Table(get_cell_permutations(grid_topology)) you will see 2s and 1s there, so that I guess that the mesh is non-oriented

the indices in the array returned by Table(get_cell_permutations(grid_topology)) are used to indirectly address another array within the CellConformity object, which is called ctype_lface_pindex_pdofs

the third nested index from left to right is the one which is feed from Table(get_cell_permutations(grid_topology))

ctype_lface_pindex_pdofs[ctype][lface][pindex]

the cause of the problem is that the length of ctype_lface_pindex_pdofs[ctype][lface] is 1

so that when we use the 2s in Table(get_cell_permutations(grid_topology)) to feed ctype_lface_pindex_pdofs[ctype][lface][pindex] we get an out-of-bounds error

[Gccxeon]

I am facing the same error with a mesh contains quadrilateral and hexahedrons, is there any plan to fix this?

[ovanvincq]

The same error occurs when using an higher-order mesh generated by Gmsh.
Is a fix expected?

[amartinhuertas]

What do you exactly mean by higher-order mesh? (e.g., curved elements?)
Do you use simplices or bricks?

[amartinhuertas]

Can you write a MWE so that we can reproduce the error?

[ovanvincq]

@amartinhuertas I mean curved 2D elements generated by Gmsh. Below a MWE.

Here are the files used for the mesh:

• geo file opened with Gmsh: https://github.com/ovanvincq/test_cudss/raw/main/test.geo
• msh file containing 1st order elements (just click on Mesh->2D in Gmsh): https://github.com/ovanvincq/test_cudss/raw/main/test1.msh
• msh file containing 2nd order elements (click on Mesh->2D and Mesh->Set order2 in Gmsh): https://github.com/ovanvincq/test_cudss/raw/main/test2.msh

using Gridap, GridapGmsh
model1 = GmshDiscreteModel("test1.msh"); #1st order elements
model2 = GmshDiscreteModel("test2.msh"); #2nd order elements
Ω1 = Triangulation(model1);
Ω2 = Triangulation(model2);
reffe_nedelec = ReferenceFE(nedelec,Float64,2)
reffe_lagrangian = ReferenceFE(lagrangian,Float64,2)
V1_lagrangian = TestFESpace(Ω1,reffe_lagrangian,conformity=:H1) #ok
V1_nedelec = TestFESpace(Ω1,reffe_nedelec,conformity=:Hcurl) #ok
V2_lagrangian = TestFESpace(Ω2,reffe_lagrangian,conformity=:H1) #ok
V2_nedelec = TestFESpace(Ω2,reffe_nedelec,conformity=:Hcurl) #error

The last line generates an error with Gridap v0.18.9, GridapGmsh v0.7.2 and julia v1.10.8:

ERROR: BoundsError: attempt to access 1-element Vector{Vector{Int64}} at index [2]
Stacktrace:
  [1] getindex
    @ .\essentials.jl:13 [inlined]
  [2] getindex
    @ .\abstractarray.jl:1293 [inlined]
  [3] getindex!(cache::Gridap.Arrays.CachedVector{Int32, Vector{Int32}}, a::Gridap.FESpaces.CellDofsNonOriented, cell::Int64)
    @ Gridap.FESpaces C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\ConformingFESpaces.jl:562
  [4] _generate_data_and_ptrs_fill_ptrs!(ptrs::Vector{Int32}, vv::Gridap.FESpaces.CellDofsNonOriented)
    @ Gridap.Arrays C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\Arrays\Tables.jl:165
  [5] generate_data_and_ptrs(vv::Gridap.FESpaces.CellDofsNonOriented)
    @ Gridap.Arrays C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\Arrays\Tables.jl:153
  [6] Table
    @ C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\Arrays\Tables.jl:27 [inlined]
  [7] 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 C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\ConformingFESpaces.jl:276
  [8] _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 C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\ConformingFESpaces.jl:149
  [9] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{…}, cell_fe::Gridap.FESpaces.CellFE{…}; trian::Gridap.Geometry.BodyFittedTriangulation{…}, labels::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{…}, dirichlet_masks::Nothing, constraint::Nothing, vector_type::Type)
    @ Gridap.FESpaces C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\FESpaceFactories.jl:23
 [10] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{…}, cell_reffe::Gridap.Arrays.CompressedArray{…}; conformity::Symbol, trian::Gridap.Geometry.BodyFittedTriangulation{…}, labels::Gridap.Geometry.FaceLabeling, dirichlet_tags::Vector{…}, dirichlet_masks::Nothing, constraint::Nothing, vector_type::Nothing)
    @ Gridap.FESpaces C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\FESpaceFactories.jl:111
 [11] FESpace(model::Gridap.Geometry.UnstructuredDiscreteModel{…}, reffe::Tuple{…}; kwargs::@Kwargs{…})
    @ Gridap.FESpaces C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\FESpaceFactories.jl:126
 [12] FESpace(t::Gridap.Geometry.BodyFittedTriangulation{…}, reffes::Tuple{…}; trian::Nothing, kwargs::@Kwargs{…})
    @ Gridap.FESpaces C:\Users\Olivier\.julia\packages\Gridap\Ce4E8\src\FESpaces\FESpaceFactories.jl:5
 [13] top-level scope
    @ REPL[28]:1