Changes to C0 polyhedron tetrahedralization loops by GiudGiud · Pull Request #4359 · libMesh/libmesh

GiudGiud · 2025-12-22T17:59:55Z

based on top of the hexagon tiling PR, ignore that content for now! It's in another PR

see companion PR in MOOSE creating a bunch of polys with exodus viz of the tets:
idaholab/moose#32141

src/geom/cell_polyhedron.C

GiudGiud · 2025-12-28T12:30:12Z

src/geom/cell_c0polyhedron.C

+      if (local_tet_quality[j2] == 0 || local_tet_quality[j2] == far_node)
+      {
+        nodes_by_geometry.erase(geometry_it);
+        has_skipped_adding_tets_and_retriangulating = true;


this is the saddest part of this PR. If we could still work with a valid surface mesh it would be a lot safer

unfortunately I seem to need this for every planar 4 node face, we always end up with an empty final tet

unfortunately I seem to need this for every planar 4 node face, we always end up with an empty final tet

Oh, man - are we getting hit again by the nastiest open secret of tet meshing?

It is not always possible to generate a mesh that has only non-degenerate tets and conforms to an existing triangulation, unless we add at least one "Steiner point" on the volume interior! See this horrifyingly simple example: https://github.com/libMesh/libmesh/blob/devel/tests/geom/elem_test.h#L224-L240 I think this is probably why it's so hard to find a tet mesher that doesn't habitually add interior points and screw up Yaqi's plans: if you know you're sometimes going to need to add points just to make a mesh possible, why hesitate to add points to make a mesh better?

I fear the only completely robust solution for our polyhedra code may be to add a new element (C0IPPolyhedron? I'd hope to come up with a better name...) that adds a single Interior Point, after which getting a nice tetrahedralization becomes trivial.

But in the short run, if you're allowed to swap diagonals then I think you can generally fix the problem with some diagonal swapping.

Is this PR ready for review now that it's no longer marked WIP?

Happy New Year!

But in the short run, if you're allowed to swap diagonals then I think you can generally fix the problem with some diagonal swapping.

I'm doing this to extrude polygons into prisms. On the sides I have quads. Are you saying I just need to do something to the diagonals of these (quad) polygons? There's no API for it right now right? It just comes from the ordering of the nodes in the polygon?

I fear the only completely robust solution for our polyhedra code may be to add a new element (C0IPPolyhedron? I'd hope to come up with a better name...) that adds a single Interior Point, after which getting a nice tetrahedralization becomes trivial.

Could we just use the centroid as the interior node? Seems like the tetrahedralization is trivial after doing that? (centroid to every triangle)
The thing is that I am doing this for finite volume sims, so I don't really need extra nodes. In fact I am not even sure I need a tetrahedralization (well except to compute the volume)

Is this PR ready for review now that it's no longer marked WIP?

Certainly ready for a look. if you use my MOOSE PR you can look at the current tetrahedralization of the hexagonal prism in exodus. You should see there is a non-matching edge inside, but is it a problem for the FE math? Certainly is a problem for me when converting the polyhedra to tets in some mesh generators

the only completely robust solution for our polyhedra code may be to add a new element (C0IPPolyhedron

What about a C0PolygonalPrism element? We could always triangulate it as the union of triangular prisms, and we could possibly have custom FE that would be worth including in libmesh? "Tensor product" of the 1D line and the C0Polygon

EDIT: well except if we do the FE math with the tets, the 1D line won't really appear in one direction. We'd have to do something different. Maybe we would not use tets at all? Just triangular prism sub-elements

Are you saying I just need to do something to the diagonals of these (quad) polygons? There's no API for it right now right? It just comes from the ordering of the nodes in the polygon?

Correct, correct, and correct (at construction time), I'm afraid. If you retriangulate() after construction time then we'll try to come up with a geometrically nicer triangulation, but you can't really control that without changing the geometry, a no-no in this context. The only way to assign the quad diagonals is to make use of knowing the default triangulation and then use that to choose which of the nodes of the quad you assign to node 0 after constructing it.

Could we just use the centroid as the interior node? Seems like the tetrahedralization is trivial after doing that?

In our mesh generators IMHO the vertex_average would be the best choice; the tetrahedralization for a convex polyhedron is trivial for any interior point choice.

In fact I am not even sure I need a tetrahedralization (well except to compute the volume)

Not directly, but we do use it in the Lagrange shape functions, so anything trying to compute a mapping will scream if it's not there. We might be able to get by without it thanks to having thrown up our hands at the idea of figuring out what a "master" element looks like for an arbitrary polyhedron, though...

What about a C0PolygonalPrism element?

Not crazy, but I think the IP polyhedron will end up being less new work.

Pinging @jwpeterson since he's had opinions on these names in the original PRs.

Fine with C0PolyhedronP1Plus

Ok I ll work on this. Don't think we have funding for this anymore though
I m planning to make it inherit from C0Polyedron. That makes sense right, adding interior points in the derived class

or do you prefer an intermediate base class? C0PolyhedronBase and two derived children, the regular one and the P1Plus

I have "Polyhedral/tetrahedral mesh enhancements" in my notes under "FY25-26 NEAMS Tasks"/"Named M2 Tasks", with "Polyhedra: very high priority" even, but I'll ping @lindsayad in case my meeting stenography skills failed me. The biggest goal there (IMHO - no notes on that) was I/O support, but if even a simple extrusion is giving us tetrahedralization problems then we probably can't hope to start reading arbitrary files until we've fixed that.

We probably want a C0PolyhedronBase here.

Your notes are spot-on and perfectly timed with my slack NEAMS-MP posts today!

moosebuild · 2025-12-28T20:58:42Z

Job Coverage, step Generate coverage on 26c1123 wanted to post the following:

Coverage

	3463d9	#4359 26c112
	Total	Total	+/-	New
Rate	65.38%	65.46%	+0.08%	100.00%
Hits	78029	78155	+126	64
Misses	41320	41247	-73	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

GiudGiud · 2026-03-16T05:27:19Z

src/geom/cell_c0polyhedron.C

+  // with no additional interior point
+  try
+  {
  // In 3D, this will require nested loops: an outer loop to remove


prefer as is or with an indentation (and a big diff)

comment goes line 400 now

src/geom/cell_c0polyhedron.C

GiudGiud · 2026-03-16T13:53:22Z

Going to try unit tests on the FE / quadratures.
My suspicion is that the element with the midnode as I have it has one more dof (as you'd expect) and is not equivalent. This would be the right reason to make it another class?

GiudGiud · 2026-03-17T01:08:37Z

tests/fe/fe_test.h

  public:                                                               \
-  FETest_##order##_##family##_##elemtype() :                            \
-    FETest<order,family,elemtype>() {                                   \
+  FETest_##order##_##family##_##elemtype##_##case_name() :               \


@roystgnr I m not sure how to "not append _ + something" when case_name isnt set.

Except if I copy paste the macro and duplicate it.
For now it appends _Default to all existing unit test names

OR
if I append "_" for the Default instead of "_Default"

GiudGiud · 2026-03-17T16:23:03Z

src/geom/cell_polyhedron.C

@@ -93,6 +93,10 @@ Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
          }
      }

+    // Plan room for an extra node so we don't invalidate this->_nodes when we add to
+    // nodelinks_data in the derived class
+    _nodelinks_data.reserve(_nodelinks_data.size() + 1);


could also modify the derived class but this is cheaper?

GiudGiud

This grew but I think it's ready for a look now!

GiudGiud · 2026-03-17T16:57:13Z

looks like we are not consistently getting the midnode with the cube poly. I ll try a few diagonal changes to see if it stabilizes on civet
on mac I get it consistently
EDIT: was opt vs devel not linux vs arm

include/geom/cell_c0polyhedron.h

src/geom/cell_c0polyhedron.C

GiudGiud · 2026-03-18T00:44:07Z

ok found the problem. Testing occurs in opt mode, and we did not detect the flat tets being added in opt mode.
Should be good now (except maybe dbg mode? we'll see)

GiudGiud · 2026-03-18T01:20:53Z

tests/fe/fe_lagrange_test.C

@@ -43,4 +43,7 @@ INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID14);
 INSTANTIATE_FETEST(THIRD, LAGRANGE, PYRAMID18);

 INSTANTIATE_FETEST(FIRST, LAGRANGE, C0POLYHEDRON);
+// We don't match the linear solution at this time
+// struct MidNode {};
+// INSTANTIATE_FETEST_CASE(FIRST, LAGRANGE, C0POLYHEDRON, MidNode);


am I interpreting this right. We are just not 1st order because of the additional midnode dof so we don't match the analytic first order solution
EDIT: does not seem right though

1) test: FETest_FIRST_LAGRANGE_C0POLYHEDRON_MidNode::testU (F) line: 814 ../../tests/fe/fe_test.h double equality assertion failed - Expected: 0.03125 - Actual : 0.00705710457070803 - Delta : 1e-09 2) test: FETest_FIRST_LAGRANGE_C0POLYHEDRON_MidNode::testGradU (F) line: 851 ../../tests/fe/fe_test.h double equality assertion failed - Expected: 0.955297792418862 - Actual : 1 - Delta : 1.5e-08 3) test: FETest_FIRST_LAGRANGE_C0POLYHEDRON_MidNode::testGradUComp (F) line: 886 ../../tests/fe/fe_test.h double equality assertion failed - Expected: 0.955297792418862 - Actual : 1 - Delta : 1.5e-08

or does this indicate a mistake with the mapping / the FE math?
The FE math for Lagrange looks like it would not change from having more sub-tets

ok changed the generic projector and it works now

…rst technique fails - adapt code for one point not being a node Clean up tetrahedralization routines when using the mid-node

Revert the master subelement back to the parent class

… bad poly meshes

- if created, add to the mesh after creating the element

Adapt other elem and FE tests to new constructor

… anymore

Add another error to catch for reverting to the mid-node creating algorithm

… instead in the generic projector

GiudGiud commented Dec 23, 2025

View reviewed changes

src/geom/cell_polyhedron.C Outdated Show resolved Hide resolved

GiudGiud force-pushed the PR_poly_fixup branch 2 times, most recently from 5990f8c to 45a4142 Compare December 25, 2025 17:39

GiudGiud changed the title ~~Poly fixups WIP~~ Changes to C0 polyhedron tetrahedralization loops Dec 28, 2025

GiudGiud commented Dec 28, 2025

View reviewed changes

GiudGiud force-pushed the PR_poly_fixup branch from 3f652d6 to 313097a Compare December 28, 2025 16:22

roystgnr added this to NEAMS MP Backend 26 Jan 7, 2026

roystgnr moved this to In progress in NEAMS MP Backend 26 Jan 7, 2026

GiudGiud force-pushed the PR_poly_fixup branch from 313097a to 5315df4 Compare March 16, 2026 05:25

GiudGiud commented Mar 16, 2026

View reviewed changes

src/geom/cell_c0polyhedron.C Outdated Show resolved Hide resolved

GiudGiud force-pushed the PR_poly_fixup branch 2 times, most recently from 156336b to 105e26c Compare March 17, 2026 01:07

GiudGiud commented Mar 17, 2026

View reviewed changes

GiudGiud force-pushed the PR_poly_fixup branch 10 times, most recently from c8b481b to 827a736 Compare March 17, 2026 16:21

GiudGiud commented Mar 17, 2026

View reviewed changes

GiudGiud force-pushed the PR_poly_fixup branch from 827a736 to 6f329cb Compare March 17, 2026 16:56

GiudGiud force-pushed the PR_poly_fixup branch 2 times, most recently from aa4f570 to 6c30f70 Compare March 17, 2026 22:39

GiudGiud commented Mar 17, 2026

View reviewed changes

include/geom/cell_c0polyhedron.h Outdated Show resolved Hide resolved

GiudGiud commented Mar 17, 2026

View reviewed changes

src/geom/cell_c0polyhedron.C Outdated Show resolved Hide resolved

GiudGiud force-pushed the PR_poly_fixup branch 3 times, most recently from 0588f58 to ed9f740 Compare March 18, 2026 00:42

GiudGiud force-pushed the PR_poly_fixup branch from ed9f740 to 9ae33a6 Compare March 18, 2026 01:11

GiudGiud commented Mar 18, 2026

View reviewed changes

GiudGiud added 19 commits March 17, 2026 22:35

Add a routine to get the poly sides the polyhedron tets are part of

ed95c14

Use try catch to tetrahedralize another way into the same class if fi…

961d849

…rst technique fails - adapt code for one point not being a node Clean up tetrahedralization routines when using the mid-node

Fix build error

2e69019

Re-implement master point for C0 polys to add the mid-element node

1ae7ef0

Revert the master subelement back to the parent class

Move up the try condition to hit various asserts within the 'try' for…

a751786

… bad poly meshes

Pass a reference to a unique ptr to the mid element node:

023a9b3

- if created, add to the mesh after creating the element

Add a volume test for polys made with new way

dd15dd1

Adapt other elem and FE tests to new constructor

Add an option for a case name when using the macro to create FE tests

705a74f

Add FE tests when using the mid node for XYZ monomial and lagrange

b721c0e

Fix re-allocation issue where _nodes does not point to nodelinks data…

1c6af23

… anymore

Fix number of vertices

7b0985a

Add a test for mapping the subtet sides to poly sides

8434a40

Try a diagonal swap to get a consistent midnode

ca1a7fd

Switch to a hexahedral prism to be sure to have a midnode

81a96e8

prepare the mesh as we added a node and an element

6979a41

Revert misc changes

66e6438

Remove a memory leak

bf95257

Add another error to catch for reverting to the mid-node creating algorithm

Detect the interior node in the generic projector

31a366c

Dont recognize the mid node as an interior node, but as a vertex node…

26c1123

… instead in the generic projector

GiudGiud force-pushed the PR_poly_fixup branch from 9ae33a6 to 26c1123 Compare March 18, 2026 15:52

Conversation

GiudGiud commented Dec 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

GiudGiud Jan 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

GiudGiud Jan 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

moosebuild commented Dec 28, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Coverage

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

GiudGiud commented Mar 16, 2026

Uh oh!

GiudGiud Mar 17, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

GiudGiud left a comment

Choose a reason for hiding this comment

Uh oh!

GiudGiud commented Mar 17, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

GiudGiud commented Mar 18, 2026

Uh oh!

GiudGiud Mar 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

GiudGiud Mar 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

GiudGiud commented Dec 22, 2025 •

edited

Loading

GiudGiud Jan 1, 2026 •

edited

Loading

GiudGiud Jan 1, 2026 •

edited

Loading

moosebuild commented Dec 28, 2025 •

edited

Loading

GiudGiud Mar 17, 2026 •

edited

Loading

GiudGiud commented Mar 17, 2026 •

edited

Loading

GiudGiud Mar 18, 2026 •

edited

Loading

GiudGiud Mar 18, 2026 •

edited

Loading