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

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

+      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 1510b5f wanted to post the following:

Coverage

	be322b	#4359 1510b5
	Total	Total	+/-	New
Rate	65.38%	65.46%	+0.08%	100.00%
Hits	78029	78165	+136	69
Misses	41320	41242	-78	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

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

+  // 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

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-17T16:23:03Z

@@ -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

… anymore

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

… instead in the generic projector

roystgnr

This looks practically ready.

roystgnr · 2026-03-23T16:13:51Z

+      // Treat polyhedron midnode as a vertex
+      // NOTE: if we start having edge or side nodes on polyhedra, we need to use that +1 offset
+      // in the edge and side nodes code as well!
+      for (unsigned int v=0; v != n_vertices + has_poly_midnode; ++v)


Let's switch to a make_range() here? I think compilers are smart enough to realize we're not doing any const_cast tricks on n_vertices or has_poly_midnode and so we don't need to do that addition on every loop iteration, but nicer to be certain.

Hmm... more importantly: treating the midnode as a vertex means we're going to be interpolating onto it regardless of FEType. That's actually what we want to do anyway for LAGRANGE variables, and the test of FEInterface::n_dofs_at_node(fe_type, elem, 0, add_p_level)) means we should be fine for MONOMIAL and XYZ too, but I could see this biting us in a very subtle way when we try to add higher order or more exotic continuous FE options. Could we add something like libmesh_assert(!has_poly_midnode || fe_type.family == LAGRANGE) above this loop?

Done and added

roystgnr · 2026-03-23T16:20:54Z

-      for (unsigned int v=0; v != n_vertices; ++v)
+      // Treat polyhedron midnode as a vertex
+      // NOTE: if we start having edge or side nodes on polyhedra, we need to use that +1 offset
+      // in the edge and side nodes code as well!


IMHO we'll want to number the midnode last (as we do with e.g. Quad9 and Hex27), but that's speculative; this comment is fine as-is and we definitely want some reminder of the issue.

roystgnr · 2026-03-23T16:31:45Z

 #include "libmesh/libmesh_common.h"
 #include "libmesh/cell_polyhedron.h"
 #include "libmesh/enum_order.h"
+#include "libmesh/node.h"


Can we get away with a forward declaration of Node? No big deal if not.

yeah we can. changed

roystgnr · 2026-03-23T16:33:44Z

+bool C0Polyhedron::is_vertex(const unsigned int i) const
 {
-  libmesh_assert (i < this->n_nodes());
+  libmesh_assert_less (i, this->n_vertices() + _has_mid_elem_node);


libmesh_assert_less (i, this->n_nodes()) should still work here, right?

roystgnr · 2026-03-23T16:55:01Z

+   *
+   * Each subelement here is a tetrahedron
+   */
+  virtual std::array<int, 4> subelement_sides_to_poly_sides(unsigned int i) const;


A little more doxygen here? I assume returnval[j] of subelement_sides_to_poly_sides(i) is the polyhedron-local index of the subtet-local index j of subtet i, but we should spell that out.

roystgnr · 2026-03-23T16:56:45Z

+  std::array<int, 4> sides = {invalid_int, invalid_int, invalid_int, invalid_int};
+  libmesh_assert(tri_i < this->n_subelements());
+  const auto & tet_nodes = _triangulation[tri_i];
+  for (const auto poly_side : make_range(n_sides()))


At some point we should document which of the "you'd think it'd be O(1)" polyhedron APIs is actually O(n_sides()), but I'm pretty sure I committed this sin before you did and we don't need to fix it in this PR.

added a note

roystgnr · 2026-03-23T17:01:27Z

+  // This might not succeed, not every surface triangulation gives a tetrahedralization
+  // with no additional interior point
+  // But if it succeeds, it uses less tetrahedra to cut the polyhedron
+  try


We probably need something to fix --disable-exceptions builds; I'll add that recipe.

We could just do libmesh_try/libmesh_catch ... which basically ifdef's out the catch block, and gives you code that still works, but only in the cases where no exceptions are thrown.

Perhaps here we should be (manually) ifdef'ing out the try block instead? If a user doesn't have exceptions disabled for some reason, we just put a midpoint node in every poly so we always get a valid (albeit often less efficient) element.

ok used that last solution

- add an assert in generic projector to avoid getting surprised with the midnode later - when disabling exceptions, always create the midnode as midnode need detection relies on exceptions - try a forward declare for Node

roystgnr · 2026-03-23T20:20:58Z

I can't tell for sure whether the SIGABRT in vectorpostprocessors/extra_id_integral.functor/missing_vector_name_error is unrelated, but that's pretty much always from something like the test harness, isn't it?

roystgnr · 2026-03-23T20:27:34Z

Double-free in the --disable-exceptions build. I'll look at that one.

GiudGiud · 2026-03-23T21:45:55Z

I think the tests might be hardcoded for a given configuration. I ll amend that

GiudGiud commented Dec 23, 2025

View reviewed changes

Comment thread src/geom/cell_polyhedron.C Outdated

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

Comment thread src/geom/cell_c0polyhedron.C Outdated

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

Comment thread tests/fe/fe_test.h

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 added 11 commits March 17, 2026 22:35

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

roystgnr requested changes Mar 23, 2026

View reviewed changes

GiudGiud force-pushed the PR_poly_fixup branch from a318f84 to b13e250 Compare March 23, 2026 18:40

Address review:

7c6b620

- add an assert in generic projector to avoid getting surprised with the midnode later - when disabling exceptions, always create the midnode as midnode need detection relies on exceptions - try a forward declare for Node

GiudGiud force-pushed the PR_poly_fixup branch from b13e250 to 7c6b620 Compare March 23, 2026 18:42

roystgnr approved these changes Mar 23, 2026

View reviewed changes

GiudGiud force-pushed the PR_poly_fixup branch 2 times, most recently from 47ea3a9 to 8262bb7 Compare March 24, 2026 01:10

GiudGiud added 3 commits March 23, 2026 19:25

Adapt test to forcing midnode creation

63ed020

Add an additional override to account for the midnode

eda9e8c

Dont delete manually the clone midnode, it's held by a unique ptr

1510b5f

GiudGiud force-pushed the PR_poly_fixup branch from 8262bb7 to 1510b5f Compare March 24, 2026 01:44

roystgnr approved these changes Mar 24, 2026

View reviewed changes

roystgnr enabled auto-merge March 24, 2026 02:46

roystgnr merged commit 2ecfa0b into libMesh:devel Mar 24, 2026
27 of 28 checks passed

github-project-automation Bot moved this from In progress to Done in NEAMS MP Backend 26 Mar 24, 2026

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!

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!

roystgnr left a comment

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!

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!

Choose a reason for hiding this comment

Uh oh!

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 commented Mar 17, 2026 •

edited

Loading