Curved gmsh meshes in Firedrake

[ethrelfall]

Dear community,

I'd like to be able to use higher-order geometries in Firedrake - eventually 3D, but I'm starting with 2D. I'm using meshes in gmsh .msh format.

There is a very simple example at https://github.com/ethrelfall/simple-curved-meshes/tree/main. It doesn't work; it gives the error

ValueError: cannot reshape array of size 168 into shape (12,2)

when trying to set up a FunctionSpace. Note there are 14 elements, quadratic order, 168=14*12 wh. 12 is probably 6 nodes * 2 dimensions.

Is there a way to make this work, please?

Search for matching discussions done - nada.

Thanks,

Ed

[colinjcotter]

Hi Ed,
I don't think that reading higher order meshes from Gmsh is currently supported (but I could be wrong).
The best way to use higher order meshes in Firedrake is via Netgen. It is a much more sophisticated and reliable meshing tool than Gmsh.
There is a tutorial here with a section on higher order meshes:
https://www.firedrakeproject.org/demos/netgen_mesh.py.html

[ethrelfall]

Many thanks for the answer.

I have .msh files generated by some of the Nektar++ guys from their own format (xml) - Nektar's high-order geom capability is pretty good. Can these be converted to Netgen do you know?

I'll take a good look at the Netgen demo anyhow.

Cheers,

Ed

[colinjcotter]

It looks like Netgen has a gmsh reader. @UZerbinati might be able to advise on how comprehensive it is.

[ethrelfall]

V good - I've got some meshing experts via Nektar collab so I can ask them too ...

[UZerbinati]

A quick look at the implementation of the ReadGmsh function in NGSolve (https://github.com/NGSolve/netgen/blob/master/python/read_gmsh.py) makes me think that we don't support reading curved mesh from msh format. Maybe worth directly exporting the high order mesh from Nektar++ to NGSolve ?

[colinjcotter]

By the way, if it is specifically a high order disk mesh that you want, you can get one from a low order disk mesh by interpolating the distance R to the centre of the disk to the CG1 space. Then you build a Vector CGk field by interpolating
‘’’
(x,y)*R/sqrt(x^2 + y^2)
‘’’
after dealing appropriately with a point at the centre, if there is one.
Then you make a new mesh with that vector field as the coordinates, as per the Changing Mesh Coordinates section of the manual.

Similar tricks are possible whenever an analytic equation defining the boundary (or part of it) is known.

[ethrelfall]

Many thanks for the answers. I've found that I can make / save / load simple meshes inc. curved elements using the Netgen integration linked by Colin. (Uploaded simple script with things I tried to my repo linked above.) I haven't yet found a way to use my gmsh-format curved meshes in Firedrake (tried converting to Netgen .vol format using meshio, which works for straight meshes but not curved ones). I'm talking to the Nektar++ guys about how to resolve this - the meshes are fairly involved as they come from a CAD rep of a tokamak divertor.

[UZerbinati]

As mentioned in the previous comment, I don't think Netgen imports the curved information. Thus even importing a curved mesh you will end up having only it's linear bit, no ?

[ethrelfall]

I can save a Netgen mesh (.vol) using Firedrake's Netgen integration, re-load it and make a curved mesh using Firedrake's curve_field method. So it seems to me Netgen is successfully loading the info needed for a curved mesh in this case. (See test_netgen.py in my repo.)
This doesn't work with curved mesh from gmsh ported to .vol using meshio - it looks like some necessary info (edgesegmentsgi2) is missing from the .vol in this case - I assume the problem is meshio not being able to port the info needed for curved mesh.
I get that read_gmsh doesn't currently support reading higher-order info.
I'm absolutely happy to be corrected ...

[colinjcotter]

What's meshio?

[ethrelfall]

mesh convertor https://pypi.org/project/meshio/

[ethrelfall]

Further to this - I've been trying to generate curved meshes in Netgen format. I want to be able to load these from file, but I'm not sure how. I can save a mesh and re-load it, but the call to the curve_field method seems to need the geo object that I used to create the mesh i.e. the curved mesh info is not actually saved in the mesh file. See https://github.com/ethrelfall/simple-curved-meshes/tree/main/further_example_1.
Is there a way to save / load a curved mesh? And also can the .vol file handle higher than order-2 geometry?

[ethrelfall]

This is great - all your scripts work for me. I wasn't calling the ngmesh.Curve method and so my saved mesh was linear (my error was a div zero in computation of barycentres).
I can see this generates curved_mesh.vol as a second-order geometry. How can I generate a higher-order one?
I think I now need to repeat these steps in 3D and that should enable us to work out how to convert a higher-order Gmsh mesh to .vol.
Again, many thanks!

[UZerbinati]

I'm looking into converting high-order meshes from gmsh to .vol, as well.
It is interesting because a .vol file generated by meshio seems to contain some curving information, see image attached (where I import the .vol file in Netgen). Yet the curving information is lost as soon as we wrap in Firedrake (I'm looking in Netgen codebase to figure out why now).

[ethrelfall]

Yeah, the triangles I see from meshio are triangle6s, with 6 point indices, but Firedrake won't read this (KeyError for me). It will read the format in your generated curved_mesh file, which has the first 3 points as integer indices then 3 additional (x,y) pairs describing points (I think).

[ChrLackner]

Netgen has 2 ways of curving:

1. secondorder where edge midpoints are really stored in the mesh and trigs are stored with 6 points, this can be done by

geo = OCCGeometry("mystepfile.step")
ngmesh = geo.GenerateMesh(...)
ngmesh.SecondOrder()

2. High (arbitrary) order curving where curvature coefficients are computed by projecting edges onto the geometry, this is what you use with firedrake bindings afaik

geo = OCCGeometry("mystepfile.step")
mesh = ngsolve.Mesh(geo.GenerateMesh(...))
mesh.Curve(3)

If you save both of these mesh files you should see that in case 1 6 points per trig are stored and in case 2 only 3 points, but additionally the geometry to recalculate the curvature

For 2 to be able to save and load the geometry needs to be stored with the mesh (netgen should do this automatically when the geometry is attached to the mesh as ngmesh.SetGeometry or when creating the mesh from the geometry).
And then on mesh load you can just call the Curve function to recompute the curvature factors.

Another option is to pickle the mesh, this should also save the curvature coefficients.

[ethrelfall]

Many thanks for this! Is there documentation please?