Check local numeration of quads and hexs

[fverdugo]

To be fixed if one wants to import meshes of quads / hexs

[AndiMD]

Here is a simple med file with 16 linear hexahedrons exported from Salome Platform. It seems to open fine in gmsh, but import fails

model = GmshDiscreteModel("Mesh_SimpleHex.med")

Info    : Reading '/home/andi/projects/juliaModules/gridapTest/bugreport/Mesh_SimpleHex.med'...
Info    : Reading MED file V4.0.0 using MED library V4.0.0
Info    : Reading 3-D unstructured mesh 'Mesh_1'
Info    : Done reading '/home/andi/projects/juliaModules/gridapTest/bugreport/Mesh_SimpleHex.med'
ERROR: LoadError: BoundsError: attempt to access 114-element Array{Int64,1} at index [0]
Stacktrace:
 [1] setindex! at ./array.jl:847 [inlined]
 [2] _apply_offset_for_faces!(::Array{Int64,1}, ::Array{Int64,1}, ::Array{Int64,1}, ::Int64) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:464
 [3] _fill_dim_to_face_to_label!(::Array{Array{Int64,1},1}, ::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:436
 [4] _setup_face_labels(::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Int64,1}, ::Array{Array{Array{Int64,1},1},1}, ::Array{Array{String,1},1}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:394
 [5] _setup_labeling(::Module, ::Gridap.Geometry.UnstructuredGrid{3,3,Float64,false}, ::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:284
 [6] GmshDiscreteModel(::String) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:18
 [7] top-level scope at /home/andi/projects/juliaModules/gridapTest/main.jl:9
 [8] include(::String) at ./client.jl:457
 [9] top-level scope at REPL[1]:1
in expression starting at /home/andi/projects/juliaModules/gridapTest/main.jl:9

EDIT: Simpler .msh Mesh containing 2 Hexahedra:
twoHexa.msh.txt

Mesh_SimpleHex.zip

[fverdugo]

We have always used msh format instead of med. Are they equivalent?

[AndiMD]

No, the med format is a HDF5 based format used by Salome Platform (for Code_Aster). Gmsh converts it to msh. Opening tetrahedra-based med files in GridapGmsh works just fine, so I assume the format is not the issue here.

[fverdugo]

OK. Thanks for reporting!

[fverdugo]

@jesusbonilla if your fork fixes this issue for general quad and hex meshes, we definitively want to merge it. Does it fix the MWE of @AndiMD ?

[jesusbonilla]

Yes. I have just tried, and I am able to create the model using file Mesh_SimpleHex.med than @AndiMD attached.

[AndiMD]

@fverdugo @jesusbonilla
I experimented with a few different meshes:

• Strangely, it looks like first saving the .med file with gmsh to .msh before using in GridapGmsh fixes some cases where import works, but the .vtu export looks visually wrong (see Reorder local numbering for quad and hex #25 (comment)).
• The option renumber=true is usually required. With renumber=false, the mesh might open but again, the .vtu export may show something is wrong.
• Some big hex meshes (>500k cells) created by Salome Platform seem to work flawlessly with this technique
• Something is wrong with Mesh_SimpleHex.med/msh: It imports only with renumber=false, but the mesh export is visually wrong (Reorder local numbering for quad and hex #25 (comment)).

Here is an additional mesh that has 1d,2d,3d entities and could make a nice test:
Cylinder5Hex.msh.txt

[jesusbonilla]

@AndiMD
I have been trying the examples that you sent and if I use the .med file I can reproduce the error you show in your comment. However, if I use the .msh files with renumber=true then it works (the vtu looks good).

[fverdugo]

@AndiMD @jesusbonilla Do you think we want renumber=true by default?

[jesusbonilla]

I would suggest so.

[AndiMD]

It is the better default, but possibly breaks backwards compatibility in some corner cases.

[AndiMD]

@jesusbonilla All other hex meshes open fine for me too, just this one does not (see below). My gmsh version is gmsh-4.6.0-Linux64-sdk from the website. I am willing to blame this on my setup, and merge your branch.

The med variant:

julia> a = GridapGmsh.GmshDiscreteModel("Mesh_SimpleHex.med",renumber=true)
Info    : No current model available: creating one
Info    : Reading 'Mesh_SimpleHex.med'...
Info    : Reading MED file V4.0.0 using MED library V4.0.0
Info    : Reading 3-D unstructured mesh 'Mesh_1'
Info    : Done reading 'Mesh_SimpleHex.med'
ERROR: Only consecutive elem tags allowed.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] _check_cell_tags(::Array{Array{UInt64,1},1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:136
 [3] _setup_connectivity(::Module, ::Int64) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:105
 [4] _setup_grid(::Module) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:34
 [5] GmshDiscreteModel(::String; renumber::Bool) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:19
 [6] top-level scope at REPL[13]:1

The msh variant:

julia> a = GridapGmsh.GmshDiscreteModel("Mesh_SimpleHex.msh",renumber=true)
Info    : No current model available: creating one
Info    : Reading 'Mesh_SimpleHex.msh'...
Info    : 13 entities
Info    : 44 nodes
Info    : 34 elements
Info    : Done reading 'Mesh_SimpleHex.msh'
ERROR: BoundsError: attempt to access 54-element Array{Int64,1} at index [0]
Stacktrace:
 [1] setindex! at ./array.jl:847 [inlined]
 [2] _apply_offset_for_faces!(::Array{Int64,1}, ::Array{Int64,1}, ::Array{Int64,1}, ::Int64) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:491
 [3] _fill_dim_to_face_to_label!(::Array{Array{Int64,1},1}, ::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:463
 [4] _setup_face_labels(::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Int64,1}, ::Array{Array{Array{Int64,1},1},1}, ::Array{Array{String,1},1}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:421
 [5] _setup_labeling(::Module, ::Gridap.Geometry.UnstructuredGrid{3,3,Float64,false}, ::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,false}, ::Array{Int64,1}) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:311
 [6] GmshDiscreteModel(::String; renumber::Bool) at /home/andi/.julia/dev/GridapGmsh/src/GmshDiscreteModels.jl:21
 [7] top-level scope at REPL[14]:1

[AndiMD]

Indeed when looking at both meshes in Gmsh, the converted .msh file (right) is missing something:

EDIT: Adding the ungrouped hexahedra into a volume group fixes this issue.

[jesusbonilla]

@AndiMD I use the same version of Gmsh that you have. Can you try if this (cube.zip) med mesh works? It works properly in my case. Or can you send an even simpler mesh that does not work (just with a couple of elements)?

[AndiMD]

@jesusbonilla Your cube.med mesh works fine.

Here is a simpler mesh that reproduces the issue:
Mesh_SimpleHex2.zip

This seems to be an issue with how gmsh deals with the .med format. gmsh seems to omit entities that are not found in a group when saving the .med file as .msh:

When I export the same mesh from Salome to the .unv mesh format, it works flawlessly. (unfortuately, .unv is not a very good replacement in general)

julia> model = GridapGmsh.GmshDiscreteModel("Mesh_SimpleHex2.med",renumber=true)
...
ERROR: Only consecutive elem tags allowed.

[AndiMD]

See issue for Gmsh: https://gitlab.onelab.info/gmsh/gmsh/-/issues/1004

[AndiMD]

Thanks to Max Orok's kind reference to the manual at onelab.info, this issue seems to be completely solved for me.
See my PR #27

@jesusbonilla Could you please confirm that the (non-default) option Mesh.SaveAll = 1; explains the difference between our systems? Thank you.

[jesusbonilla]

@AndiMD yes that option explains the different behaviors in my system too. I think that PR #27 fixes this issue :)