reading 3D hexmesh from GMSH file

[florian98765]

Dear Firedrakers, I would like to solve some 3D Maxwell problems on a structured hexmesh. I achieved to create a mesh using gmsh, but it seems that firedrake currently only supports tetrahedra (except for extruded meshes). So here my question:

• is reading the meshfile the only limitation or are there other missing pieces (do e.g. edge-elements work on the hex elements?)
• if I understand the code correctly firedrake uses PETSc to read gmsh-files. But I think PETSc supports hexmeshes. So what are the missing pieces and how can they be implemented?

thanks for any input best wishes, Florian

[florian98765]

After some reading and coding I am able to create a plex object for some hexahedral mesh (either read from gmsh or manually created by a modified version of BoxMesh.

Unfortunately I end up with the following error during creation of a FunctionSpace on the created mesh:

  File "/firedrake/src/firedrake/firedrake/functionspacedata.py", line 53, in cached
    result = f(mesh, key, *args, **kwargs)
  File "/firedrake/src/firedrake/firedrake/functionspacedata.py", line 179, in get_map_cache
    mesh.interior_facets.set: None,
  File "/firedrake/src/PyOP2/pyop2/utils.py", line 62, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/firedrake/src/firedrake/firedrake/mesh.py", line 873, in interior_facets
    return self._facets("interior")
  File "/firedrake/src/firedrake/firedrake/mesh.py", line 857, in _facets
    self.cell_closure)
  File "/firedrake/src/PyOP2/pyop2/utils.py", line 62, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/firedrake/src/firedrake/firedrake/mesh.py", line 827, in cell_closure
    raise NotImplementedError("Cell type '%s' not supported." % cell)
NotImplementedError: Cell type 'hexahedron' not supported.

Inspecting the implementation of quadrilateral meshes in 2D showed that edges (and in 3D maybe faces) need to be re-oriented. Is this a hard requirement within Firedrake, or could one work around somehow?

Details on the 2D reorder algorithm can be found e.g. here. It is mentioned that the algorithm can be extended to 3D, but there are some cases which make the method fail (e.g. something like a Möbius stripe).

Why is this reordering necessary at all? Does not PETSc also use a canonical library ordering which could be used? (e.g. 3D simplices -> first 3 vertices give an outward normal, hexahedra -> first 4 and last 4 vertices give an outward normal). What would happen if I just provide one arbitrary closure_ordering?

best wishes Florian

[wence-]

We rely on a canonical ordering so that when we go from global dofs to cell dofs for assembly, we end up with the "same" view from both sides of a face. There are two ways this could be achieved. Either we apply per entity permutations of dofs when we do global to local, or we do that "up front" [this is the one we currently do]. Getting the closure ordering right is necessary so that the second approach works.

[florian98765]

Ok, so this means that one needs to generalize the closure ordering (as partly described in the publication above) from 2D to 3D. But it also implies that one could not assemble on certain pathologic meshes (like the mentioned Möbius stripe), right?

By the way, how could a best recompile firedrake if a change something within the cython code. Right now I am using firedrake-install inside of a docker container.

best wishes Florian

[wence-]

In the top-level firedrake source directory, type make after activating the virtualenv.

[wence-]

Ok, so this means that one needs to generalize the closure ordering (as partly described in the publication above) from 2D to 3D. But it also implies that one could not assemble on certain pathologic meshes (like the mentioned Möbius stripe), right?

Yes, that's right. A broader question is if we should do option 1 instead (that would be the right thing if one also went in the direction of nonconforming mesh refinement).

[rckirby]

As a stop gap, aren’t once-refined meshes orientable? Do I recall that Wolfgang proved that? This costs computer resources if orientation fails but saves developer time to figure out orientations in the code.

On Oct 14, 2020, at 8:46 AM, Lawrence Mitchell notifications@github.com wrote:

 Ok, so this means that one needs to generalize the closure ordering (as partly described in the publication above) from 2D to 3D. But it also implies that one could not assemble on certain pathologic meshes (like the mentioned Möbius stripe), right?

Yes, that's right. A broader question is if we should do option 1 instead (that would be the right thing if one also went in the direction of nonconforming mesh refinement).

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[wence-]

As a stop gap, aren’t once-refined meshes orientable? Do I recall that Wolfgang proved that?

Yep.

[rckirby]

We have a global refinement feature, so if Florian can read the mesh, he can just try orienting after reading and handle the exception by refinement. Done for now, feature added until we get around to handling orientations locally. Running on a refinement of a general hex mesh is better than not having a general hex mesh?

On Oct 14, 2020, at 9:02 AM, Lawrence Mitchell notifications@github.com wrote:

 As a stop gap, aren’t once-refined meshes orientable? Do I recall that Wolfgang proved that?

Yep.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[wence-]

Even orientable hex meshes require writing a bunch of code to implement the orientation.

[vabhi6]

Hi, I would like to solve a structural dynamics problem using a higher order quadrature element based FEM scheme. I generated a hexahedral mesh in GMSH and tried to import it. However I get a KeyError from UFL.Cell since the hexahedron cell is not listed in the _cells dict in mesh.py. Florian earlier suggested that he was able to create a plex object but ran into trouble with FunctionSpace. So is there a solution to this as of now? Is my only option to work with hexahedral meshes created by extruding quadrilateral meshes? hexerror.txt

[florian98765]

As far as I know, there is no implemented solution so far. However the orientation-method used right now for quadrilaterals in 2D could be extended to 3D. For me this would be the way to go, but I don't know how complicated it will be to generalize the orientation method.

[EstebanFO]

Hello everyone,

Thanks to your post I understand why the traditional Hexahedral mesh created in GMSH cannot be imported in Firedrake. I have run some tests with extruded meshes created directly in Firedrake and they work well. But I was wondering, is there a way to configure GMSH to export an extruded mesh that is compatible with Firedrake?. When I extrude a hexahedral mesh in GMSH and import it into Firedrake, I get the same error (KeyError 6) than the one for non extruded meshes.

Thanks!

[vabhi6]

Dear Esteban, Sorry, I couldnt get back to you earlier. I suggest you create the 2D mesh in GMSH and then extrude it inside firedrake like so:

Create the 2D mesh in GMSH and save it as mesh2d.msh :

mesh2d=Mesh('mesh2d.msh') mesh3d=ExtrudedMesh(mesh2d,layer_height=1.0,extrusion_type='uniform')

You can use mesh3d after this

On Tue, Jun 29, 2021 at 4:10 PM EstebanFO @.***> wrote:

Hello everyone,

Thanks to your post I understand why the regular Hexahedral mesh created in GMSH cannot be imported in Firedrake. I have run some tests with extruded meshes created directly in Firedrake and they work well. But I was wondering, is there a way to configure GMSH to export an extruded mesh that is compatible with Firedrake?. When I extrude a hexahedral mesh in GMSH and import it into Firedrake, I get the same error (KeyError 6) than the one for non extruded meshes.

Thanks!

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/firedrakeproject/firedrake/issues/1859#issuecomment-870882562, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFMERSMBNUPBKAOYWOVPPN3TVISCTANCNFSM4SF24NAQ .

-- Abhishek Venketeswaran, PhD

[EstebanFO]

Hello Abhishek,

Thanks for your answer, it's still completely relevant.

That's actually the method I have been using so far and it works correctly.

Let's hope in the future we can import the 3D HEXmesh directly from Gmsh and maybe be able to use more complex geometries. For example, I've been working with heat sinks with changing cross section so I would need to make three extrusions (in different planes) to build the entire geometry.

Best regards