Polygonal/polyhedral element support

[dschwen]

I would like to explore the addition of a general element class to support arbitrary polygonal and polyhedral elements (i.e. two classes, one for 2d one for 3d, the face elements for polyhedral elements would be polygonal elements)

Polygonal elements would be useful for FV fluid modeling, and FE mechanics (among other things) and they would largely solve the meshing problem. Assuming tri/tet meshing is a solved problem, any tri/tet mesh can be transformed into a dual polygonal/polyhedral mesh using Voronoi tessellation. The Voro++ library (BSD licensed) could do almost all heavy lifting here!

A generic polyhedral element class design could look like this :

• std::vector<Node *> _local_nodes - a list of node pointers
• auto n_nodes() { return _local_nodes.size(); } - the variable number of nodes that make up the element
• std::vector<std::vector<unsigned int>> _side_nodes - a dynamically sized list of local node IDs for each element side (note that not every node needs to be a corner node, we could easily permit "free floating nodes" for enrichment DOFs - like Abaqus does)
• auto n_sides() { return _side_nodes.size(); } - the number of sides making up the element

In 2d each side would have 2 nodes for first order elements. Higher order 3d polyhedral element support would require more info to distinguish corner nodes from midpoint nodes etc.

This element would support neighbor pointers, but not adaptivity for example. It would not be hard to implement point locator support, but a reference element probably doesn't make sense here (this might be a issue for FE...).

In a first stab we would not support shape functions etc., only volume and face area calculations for FV. However even the minimal support for such a super flexible element type would make it possible to support user defined elements (like Abaqus does).

[jwpeterson]

this might be a issue for FE...

Do you have a link to the approach that Abaqus is using for polyhedral elements? I found this but it's pretty recent and looks limited to 2D. My initial reaction is to wonder why people would have spent so much effort and research on the problem of hexahedral mesh generation if polyhedral elements could have just been used instead...

[dschwen]

Abaqus has no internal support for those elements. They completely rely on UELs. UEL is a plugin type, in which users can implement all calculations to compute a residual and a Jacobian contribution. What they get from Abaqus are just the DOF values at the nodes belonging to the current element. Abaqus has no GUI support for creating or even visualizing these elements, the users will have to build their own .inp files somehow. A common pattern seems to be overlaying a mesh of dummy elements (with 1e-11 stiffness) of known element type and writing a custom UMAT plugin to transfer the solution. Bottom line, whatever we do won;t be much worse than what Abaqus does (and we could even do a lot better).

I believe OpenFOAM supports poly elements. Again, the biggest advantage is probably in the FV space, also poly elements have a much younger history. Shape functions for them are pretty complicated and they (according to some review article) only received real attention around 2000 (Wachspress functions are the simplest approach).

[roystgnr]

We have a third party working on polygonal element support, though I'm not sure how much progress they're making. I don't think anyone's done anything with libMesh towards polyhedra. I know there's user interest in a few specific polyhedra (hexagonal prisms) that would be easy to tack on to polygon support.

IIRC people spend so much effort on hex mesh generation because they're faster to compute with (not asymptotically, but a much nicer constant) and because not all formulations (and not all codes! doesn't deal.II still require tensor-product geometric elements?) support anything else.

The effort might go too far, though. I still remember a bunch of International Meshing Roundtable talks that were basically "last year we could only turn this into a 90%-hex mesh, but this year we have 92% hexes! hooray!"

[drwells]

(and not all codes! doesn't deal.II still require tensor-product geometric elements?) support anything else.

As of 9.3 (2021) we now support simplex and mixed meshes in deal.II. We need to advertise this better :)