GMSH curved meshes support

[Arpit-Babbar]

Thanks to Daniel and everyone involved for bringing the wonderful GMSH mesh support to Trixi in #1799!!

I am unable to get the current parser to handle a GMSH generated mesh with curved boundaries. Here is how I generate a curved mesh with GMSH (I only know how to do it with GUI at this moment). I tried this with NACA6412.geo available here and another one here

1. Open the .geo file in GMSH
2. Click Modules -> Mesh -> 2D (Generates the straight sided mesh)
3. Click Modules -> Mesh -> Set order 2 (curves the mesh boundary to be order 2. Higher is also possible)
4. File -> Export -> mesh_NACA6412.inp, making sure to tick "Save all elements" and "Save groups of nodes" (Exporting with these options ticked is the step I cannot do via command line)

If I skip step 3, the parser can handle it fine and I am able to run the available elixir https://github.com/trixi-framework/Trixi.jl/blob/main/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl. However, if I try to run same the elixir when I did the above step 3 to curve the mesh, I get a segmentation fault in the line that reads the mesh, i.e., mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols)

It starts with

[p4est] No elements or nodes found in mesh file.
[p4est] Failed to read /home2/arpit/repositories/cloned/Trixi.jl/examples/p4est_2d_dgsem/NACA6412_curved.inp: pass 1

Full error message

P4est segmentation fault

[p4est] No elements or nodes found in mesh file.
[p4est] Failed to read /home2/arpit/repositories/cloned/Trixi.jl/examples/p4est_2d_dgsem/NACA6412_curved.inp: pass 1

[595759] signal (11.1): Segmentation fault
in expression starting at /home2/arpit/repositories/cloned/Trixi.jl/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl:78
unsafe_load at ./pointer.jl:119 [inlined]
unsafe_load at /home2/arpit/.julia/packages/P4est/0rTeV/src/pointerwrappers.jl:86 [inlined]
getindex at /home2/arpit/.julia/packages/P4est/0rTeV/src/pointerwrappers.jl:82 [inlined]
getindex at /home2/arpit/.julia/packages/P4est/0rTeV/src/pointerwrappers.jl:82 [inlined]
p4est_mesh_from_standard_abaqus at /home2/arpit/repositories/cloned/Trixi.jl/src/meshes/p4est_mesh.jl:460
#_#279 at /home2/arpit/repositories/cloned/Trixi.jl/src/meshes/p4est_mesh.jl:376
P4estMesh at /home2/arpit/repositories/cloned/Trixi.jl/src/meshes/p4est_mesh.jl:353
unknown function (ip: 0x7f420e7c66a0)
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/julia.h:1982 [inlined]
do_call at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:126
eval_value at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:223
eval_stmt_value at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:174 [inlined]
eval_body at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:617
jl_interpret_toplevel_thunk at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:775
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:934
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:877
ijl_toplevel_eval_in at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:985
eval at ./boot.jl:385 [inlined]
include_string at ./loading.jl:2070
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
_include at ./loading.jl:2130
include at ./client.jl:489
unknown function (ip: 0x7f420e7b5b55)
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/julia.h:1982 [inlined]
do_call at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:126
eval_value at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:223
eval_stmt_value at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:174 [inlined]
eval_body at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:617
jl_interpret_toplevel_thunk at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/interpreter.c:775
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:934
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:877
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:877
jl_toplevel_eval_flex at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:877
ijl_toplevel_eval_in at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/toplevel.c:985
eval at ./boot.jl:385 [inlined]
eval_user_input at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:150
repl_backend_loop at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:246
#start_repl_backend#46 at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:231
start_repl_backend at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:228
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
#run_repl#59 at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:389
run_repl at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:375
jfptr_run_repl_91689.1 at /home2/arpit/.local/julia-1.10.0/lib/julia/sys.so (unknown line)
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
#1013 at ./client.jl:432
jfptr_YY.1013_82677.1 at /home2/arpit/.local/julia-1.10.0/lib/julia/sys.so (unknown line)
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/julia.h:1982 [inlined]
jl_f__call_latest at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/builtins.c:812
#invokelatest#2 at ./essentials.jl:887 [inlined]
invokelatest at ./essentials.jl:884 [inlined]
run_main_repl at ./client.jl:416
exec_options at ./client.jl:333
_start at ./client.jl:552
jfptr__start_82703.1 at /home2/arpit/.local/julia-1.10.0/lib/julia/sys.so (unknown line)
_jl_invoke at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:2894 [inlined]
ijl_apply_generic at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/julia.h:1982 [inlined]
true_main at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/jlapi.c:582
jl_repl_entrypoint at /cache/build/builder-amdci4-6/julialang/julia-release-1-dot-10/src/jlapi.c:731
main at julia (unknown line)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
unknown function (ip: 0x4010b8)
Allocations: 25176877 (Pool: 25169453; Big: 7424); GC: 90
Segmentation fault

Note that the .geo file does have enough information to curve the mesh, see images in spoilers below. Here is gmsh showing a straight sided mesh

Straight sided GMSH mesh

The curved airfoil already indicates that it is possible to curve the mesh. And, by setting a high order with Set order 3, I get

GMSH mesh of order 3

`HOHQ` generated meshes also give curved boundaries already.

Here is a gist with .inp file generated with curved boundaries - https://gist.githubusercontent.com/Arpit-Babbar/d9238c194b3e87574fdcbbd39fa59180/raw/6c985a65b58e5eda38630b6a8d3e54652293650a/NACA6412_curved.inp

Here is a gist with .inp file generated with straight sided boundaries - https://gist.githubusercontent.com/Arpit-Babbar/f4271dc9c57724bf62681c725aa16cd6/raw/e83e50d4f67ee1853f189d6b7269e8fb1cde5ab0/NACA6412.inp

[andrewwinters5000]

Can you post the .inp file that gmsh generates are you curve it? It would be helpful to see what the curved ABAQUS file format looks like.

The issue is with the creation of the connectivity versus the curvature. The connectivity is handled by p4est which can only handle straight-sided elements, as discussed briefly here. The curving of the mesh is handled internally within Trixi. We just need to work out how to read this boundary point data into appropriate arrays in order to create the boundary interpolants and compute the metric terms.

[Arpit-Babbar]

Can you post the .inp file that gmsh generates are you curve it? It would be helpful to see what the curved ABAQUS file format looks like.

Here a gist with .inp file generated with curved boundaries - https://gist.githubusercontent.com/Arpit-Babbar/d9238c194b3e87574fdcbbd39fa59180/raw/6c985a65b58e5eda38630b6a8d3e54652293650a/NACA6412_curved.inp

Here is a gist with .inp file generated with straight sided boundaries - https://gist.githubusercontent.com/Arpit-Babbar/f4271dc9c57724bf62681c725aa16cd6/raw/e83e50d4f67ee1853f189d6b7269e8fb1cde5ab0/NACA6412.inp

The issue is with the creation of the connectivity versus the curvature. The connectivity is handled by p4est which can only handle straight-sided elements, as discussed briefly here. The curving of the mesh is handled internally within Trixi. We just need to work out how to read this boundary point data into appropriate arrays in order to create the boundary interpolants and compute the metric terms.

I see, thanks for sharing that!

[andrewwinters5000]

This strategy of having p4est create and manage the connectivity and Trixi handle the curvature was the reason for the special .inp mesh file format for HOHQMesh, see here. Basically, all the curved boundary information is encoded as comments in the ABAQUS format such that the p4est parser ignores it and we can let Trixi handle everything curved boundary related.

[DanielDoehring]

We just need to work out how to read this boundary point data into appropriate arrays in order to create the boundary interpolants and compute the metric terms.

That would actually be a nice student project for which relatively little knowledge on math is required!

[JoshuaLampert]

Can this be closed after having merged #2217?

[DanielDoehring]

Hm so the 3D case is still not covered, so maybe keep this open?

[JoshuaLampert]

Ok, I wasn't sure if everything was resolved. So let's keep this open.