Skip to content

Changes to C0 polyhedron tetrahedralization loops #4359

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
GiudGiud wants to merge 9 commits into
base: devel
Choose a base branch
from
+422 −69
Open

Changes to C0 polyhedron tetrahedralization loops #4359

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d47670e
Add generation of a tile of polygons+tris+quad for the polygon option
GiudGiud Dec 20, 2025
72c23c5
Add a unit test for hexagonal tiling
GiudGiud Dec 21, 2025
4a805f6
Use only C0Polygons in the tile for now to avoid
GiudGiud Dec 22, 2025
dd2654c
Fixup on logic for retriangulation
GiudGiud Dec 22, 2025
f38d74c
Add heuristics to tetrahedralize hexahedrons with high valence polygo…
GiudGiud Dec 22, 2025
e6d4ae0
Add another heuristic:
GiudGiud Dec 25, 2025
7be614e
Add a routine to get the poly sides the polyhedron tets are part of
GiudGiud Dec 26, 2025
80b8cef
Fix the metric:
GiudGiud Dec 27, 2025
313097a
Don't resort nodes by geometry for the final tet
GiudGiud Dec 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions include/geom/cell_c0polyhedron.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,6 +170,14 @@ class C0Polyhedron : public Polyhedron

Point side_vertex_average_normal(const unsigned int s) const override final;

/**
* \returns the local side-indices of subelement sides of the
* polyhedron
*
* Each subelement here is a tetrahedron
*/
virtual std::array<int, 4> subelement_sides_to_poly_sides(unsigned int i) const;

/**
* Create a triangulation (tetrahedralization) based on the current
* sides' triangulations.
Expand Down
244 changes: 183 additions & 61 deletions src/geom/cell_c0polyhedron.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,39 @@ C0Polyhedron::side_vertex_average_normal(const unsigned int s) const
}


std::array<int, 4>
C0Polyhedron::subelement_sides_to_poly_sides(unsigned int tri_i) const
{
std::array<int, 4> sides = {invalid_int, invalid_int, invalid_int, invalid_int};
const auto & tet_nodes = _triangulation[tri_i];
for (const auto poly_side : make_range(n_sides()))
{
const auto sd_nodes = nodes_on_side(poly_side);
bool zero_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[0]) != sd_nodes.end();
bool one_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[1]) != sd_nodes.end();
bool two_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[2]) != sd_nodes.end();
bool three_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[3]) != sd_nodes.end();

// Take advantage of the fact that the poly side can only contain one tet side
// Note: we could optimize by replacing the booleans with the searches above
if (zero_in)
{
if (one_in)
{
if (two_in)
sides[0] = poly_side;
else if (three_in)
sides[1] = poly_side;
}
else if (two_in && three_in)
sides[3] = poly_side;
}
else if (three_in && two_in && one_in)
sides[2] = poly_side;
}
return sides;
}


void C0Polyhedron::retriangulate()
{
Expand Down Expand Up @@ -373,7 +406,7 @@ void C0Polyhedron::retriangulate()
// sometimes with a similarly-sorted vector of surrounding elements
// to go with it.
auto surroundings_of =
[&nodes_to_elem_map, & surface]
[& nodes_to_elem_map, & surface]
(const Node & node,
std::vector<Elem *> * surrounding_elems)
{
Expand Down Expand Up @@ -439,7 +472,11 @@ void C0Polyhedron::retriangulate()
v02 = static_cast<Point>(*surrounding_nodes[n+1]) - node,
v03 = static_cast<Point>(*surrounding_nodes[n+2]) - node;

total_solid_angle += solid_angle(v01, v02, v03);
// Solid angles are all negative because the surface triangulation triangles were pointing
// outwards and the 4th node is inward because the polyhedron is convex
// Some are zero when all the nodes selected are from the same (planar) polygon
// Make it positive to store smallest angle first
total_solid_angle += std::abs(solid_angle(v01, v02, v03));
}

return std::make_pair(n_surrounding, total_solid_angle);
Expand All @@ -456,8 +493,11 @@ void C0Polyhedron::retriangulate()
// in parallel.
typedef std::multimap<std::pair<int, Real>, Node*> node_map_type;
node_map_type nodes_by_geometry;
// Keep track of index in map
std::map<Node *, node_map_type::iterator> node_index;

// The map is first sorted by the valence
// then by the angle. Max angle, if they are negative (positive right now), first
for (auto node : surface.node_ptr_range())
node_index[node] =
nodes_by_geometry.emplace(geometry_at(*node), node);
Expand All @@ -466,6 +506,9 @@ void C0Polyhedron::retriangulate()
// each vertex, and an inner loop to remove multiple tetrahedra in
// cases where the vertex has more than 3 neighboring triangles.

// One of the checks does not work if we skip adding a 0 volume tet
bool has_skipped_adding_tets_and_retriangulating = false;

// We'll be done when there are only three "unremoved" nodes left,
// so they don't actually enclose any volume.
for (auto i : make_range(nodes_by_geometry.size()-3))
Expand Down Expand Up @@ -500,19 +543,56 @@ void C0Polyhedron::retriangulate()
// predecessor and successor nodes in order afterward.
auto find_valid_nodes_around =
[n_surrounding, & surrounding_nodes]
(unsigned int j)
(unsigned int j, unsigned int skip=1e6)
{
unsigned int jnext = (j+1)%n_surrounding;
while (!surrounding_nodes[jnext])
while (!surrounding_nodes[jnext] || jnext == skip)
jnext = (jnext+1)%n_surrounding;

unsigned int jprev = (j+n_surrounding-1)%n_surrounding;
while (!surrounding_nodes[jprev])
while (!surrounding_nodes[jprev] || jprev == skip)
jprev = (jprev+n_surrounding-1)%n_surrounding;

return std::make_pair(jprev, jnext);
};


// Vectors from the center node to each of its surrounding
// nodes are helpful for calculating prospective tet
// quality.
std::vector<Point> v0s(n_surrounding);
for (auto j : make_range(n_surrounding))
v0s[j] = *(Point *)surrounding_nodes[j] - *node;

// Find the tet quality we'd potentially get from each
// possible choice of tet
auto local_tet_quality_of =
[& surrounding_nodes, & v0s, & find_valid_nodes_around]
(unsigned int j, unsigned int skip=1e6)
{
auto [jminus, jplus] = find_valid_nodes_around(j, skip);

// Anything proportional to the ratio of volume to
// total-edge-length-cubed should peak for perfect tets
// but hit 0 for pancakes and slivers.

const Real total_len =
v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
(*(Point *)surrounding_nodes[jplus] -
*(Point *)surrounding_nodes[j]).norm() +
(*(Point *)surrounding_nodes[j] -
*(Point *)surrounding_nodes[jminus]).norm() +
(*(Point *)surrounding_nodes[jminus] -
*(Point *)surrounding_nodes[jplus]).norm();

// Orientation here is tricky. Think of the triple
// product as (v1 cross v2) dot v3, with right hand rule.
const Real six_vol =
triple_product(v0s[jminus], v0s[jplus], v0s[j]);

return six_vol / (total_len * total_len * total_len);
};

// We may have too many surrounding nodes to handle with
// just one tet. In that case we'll keep a cache of the
// element qualities that we'd get by making a tet with the
Expand All @@ -525,8 +605,9 @@ void C0Polyhedron::retriangulate()
std::vector<Real> local_tet_quality(n_surrounding, 1);

// From our center node with N surrounding nodes we can make N-2
// tetrahedra. The first N-3 each replace two surface tets with
// two new surface tets.
// tetrahedra.
// For the first N-3 tets we build, we may have to replace two surface
// triangles with two new surface triangles.
//
// My first idea was to greedily pick nodes with the smallest
// local (solid) angles to get the best quality. This works in
Expand All @@ -539,15 +620,39 @@ void C0Polyhedron::retriangulate()
// My third idea is to try to fix the lowest quality tets first,
// by picking cases where they have higher quality neighbors,
// and creating those neighbors so as to change them.
//
// A fourth idea was to make sure to fix as many neighbor tets as
// possible.

auto find_new_tet_nodes =
[& local_tet_quality, & find_valid_nodes_around]
[& local_tet_quality, & local_tet_quality_of, & find_valid_nodes_around, & surrounding_nodes]
()
{
unsigned int jbest = 0;
auto [jminus, jplus] = find_valid_nodes_around(jbest);
Real qneighbest = std::min(local_tet_quality[jminus],
local_tet_quality[jplus]);
// Count number of bad neighbors
auto num_bad_neigh_best = (local_tet_quality[jminus] <= 0) + (local_tet_quality[jplus] <= 0);
// Skip our heuristics if we are considering a bad tet
if (local_tet_quality[jbest] == -1e6)
{
num_bad_neigh_best = 0;
qneighbest = 1e6;
}

// Count the number of valid surrounding nodes
unsigned int n_valid_surrounding = 0;
for (const auto surr_node : surrounding_nodes)
if (surr_node)
n_valid_surrounding++;

bool bad_future_neighbors = false;
if ((n_valid_surrounding > 3) && local_tet_quality_of(jminus, /*without*/jbest) <= 0 &&
local_tet_quality_of(jplus, /*without*/jbest) <= 0)
bad_future_neighbors = true;

// Expected qualities post removing jbest from the surrounding nodes
for (auto j : make_range(std::size_t(1),
local_tet_quality.size()))
{
Expand All @@ -565,71 +670,58 @@ void C0Polyhedron::retriangulate()
qneighj > 0)
continue;

// Avoid chosing a tet that when constructed would leave both neighbors in bad shape
// TODO: simply pre-compute that for every vertex and exclude those nodes
// Note: It does not matter if we are down to 3 surrounding nodes, as we are building the final tet
if ((n_valid_surrounding > 3) && local_tet_quality_of(jminus, /*without*/j) <= 0 &&
local_tet_quality_of(jplus, /*without*/j) <= 0)
continue;

const auto num_bad_neigh = (local_tet_quality[jminus] <= 0) +
(local_tet_quality[jplus] <= 0);
// We want to try for the best possible fix.
if ((local_tet_quality[j] - qneighj) >
(local_tet_quality[jbest] - qneighj))
// - The more bad (volume <= 0) neighbors we fix the better
// - If same number of bad neighbors, pick best quality one
// - If the current pick for best tet will make the two neighbor tets bad,
// don't keep it and always try another
if (((num_bad_neigh > num_bad_neigh_best) || ((num_bad_neigh == num_bad_neigh_best)
&& ((local_tet_quality[j] - qneighj) >
(local_tet_quality[jbest] - qneighj))))
|| bad_future_neighbors)
{
jbest = j;
qneighbest = qneighj;
num_bad_neigh_best = num_bad_neigh;
bad_future_neighbors = false;
}
}

// If there are only 3 neighors and the tet is bad, we will skip it downstream
libmesh_error_msg_if
(local_tet_quality[jbest] <= 0,
(n_valid_surrounding > 3 && local_tet_quality[jbest] <= 0,
"Cannot build non-singular non-inverted tet");

std::tie(jminus, jplus) = find_valid_nodes_around(jbest);

return std::make_tuple(jbest, jminus, jplus);
};

if (n_surrounding > 3)
{
// We'll be searching local_tet_quality even after tet
// extraction disconnects us from some nodes; when we do we
// don't want to get one.
constexpr Real far_node = -1e6;

// Vectors from the center node to each of its surrounding
// nodes are helpful for calculating prospective tet
// quality.
std::vector<Point> v0s(n_surrounding);
for (auto j : make_range(n_surrounding))
v0s[j] = *(Point *)surrounding_nodes[j] - *node;

// Find the tet quality we'd potentially get from each
// possible choice of tet
auto local_tet_quality_of =
[& surrounding_nodes, & v0s, & find_valid_nodes_around]
(unsigned int j)
{
auto [jminus, jplus] = find_valid_nodes_around(j);

// Anything proportional to the ratio of volume to
// total-edge-length-cubed should peak for perfect tets
// but hit 0 for pancakes and slivers.

const Real total_len =
v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
(*(Point *)surrounding_nodes[jplus] -
*(Point *)surrounding_nodes[j]).norm() +
(*(Point *)surrounding_nodes[j] -
*(Point *)surrounding_nodes[jminus]).norm() +
(*(Point *)surrounding_nodes[jminus] -
*(Point *)surrounding_nodes[jplus]).norm();

// Orientation here is tricky. Think of the triple
// product as (v1 cross v2) dot v3, with right hand rule.
const Real six_vol =
triple_product(v0s[jminus], v0s[jplus], v0s[j]);
// Get the quality of the tets for every node around the node:
// we have chosen 1 node. If we chose another node for the tet, there
// are only two other nodes that are also neighbors, hence 1 tet per node
// Note: we always need to know the quality, because as we delete vertices
// we could end up with a a node with 3 surrounding nodes, but all coplanar!
for (auto j : make_range(n_surrounding))
local_tet_quality[j] = local_tet_quality_of(j);

return six_vol / (total_len * total_len * total_len);
};
// We'll be searching local_tet_quality even after tet
// extraction disconnects us from some nodes; when we do we
// don't want to re-use it so we mark the quality as -1e6
constexpr Real far_node = -1e6;

for (auto j : make_range(n_surrounding))
local_tet_quality[j] = local_tet_quality_of(j);

// If we have N surrounding nodes, we can make N tets and
if (n_surrounding > 3)
{
// If we have N surrounding nodes, we can make N-3 tets and
// that'll bring us back to the 3-surrounding-node case to
// finish.
for (auto t : make_range(n_surrounding-3))
Expand Down Expand Up @@ -697,6 +789,7 @@ void C0Polyhedron::retriangulate()
// unchanged until we're out of this inner loop; let's
// recalculate it here, and then we'll be done with it.
Node * & nbestref = surrounding_nodes[jbest];
libmesh_assert(nbestref);
nodes_by_geometry.erase(node_index[nbestref]);
node_index[nbestref] =
nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref);
Expand Down Expand Up @@ -724,6 +817,18 @@ void C0Polyhedron::retriangulate()
// counterclockwise" we get from the lambda.
auto [j2, j1, j3] = find_new_tet_nodes();

// If the final tet has zero volume, skip, we have covered all we need for this node
// This occurs when considering a polygonal prism.
// The opposite node on the other extruded face was deleted first
// and depending on the triangulation, this leaves the node only attached
// to 3 other coplanar nodes on one of the side
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;
Copy link
Contributor Author

@GiudGiud GiudGiud Dec 28, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

@roystgnr roystgnr Dec 31, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor Author

@GiudGiud GiudGiud Jan 1, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

@GiudGiud GiudGiud Jan 1, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

@roystgnr roystgnr Jan 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

@jwpeterson jwpeterson Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Fine with C0PolyhedronP1Plus

Copy link
Contributor Author

@GiudGiud GiudGiud Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

@GiudGiud GiudGiud Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

@roystgnr roystgnr Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

@lindsayad lindsayad Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

continue;
}

// Turn these last four nodes into a tet
Node * n1 = surrounding_nodes[j1],
* n2 = surrounding_nodes[j2],
Expand Down Expand Up @@ -767,6 +872,26 @@ void C0Polyhedron::retriangulate()
surface.delete_elem(oldtri2);
surface.delete_elem(oldtri3);

// We've used up our center node, so it's not something we can
// eliminate again.
nodes_by_geometry.erase(geometry_it);

// Recompute the valence and angles of the nodes we used
// The idea is that if one node is being isolated on one side,
// we need to treat it asap
if (nodes_by_geometry.size() > 4)
{
Node * & node_j2 = surrounding_nodes[j2];
nodes_by_geometry.erase(node_index[node_j2]);
node_index[node_j2] = nodes_by_geometry.emplace(geometry_at(*node_j2), node_j2);
Node * & node_j1 = surrounding_nodes[j1];
nodes_by_geometry.erase(node_index[node_j1]);
node_index[node_j1] = nodes_by_geometry.emplace(geometry_at(*node_j1), node_j1);
Node * & node_j3 = surrounding_nodes[j3];
nodes_by_geometry.erase(node_index[node_j3]);
node_index[node_j3] = nodes_by_geometry.emplace(geometry_at(*node_j3), node_j3);
}

// We should have used up all our surrounding nodes now, and we
// shouldn't have messed up our surface in the process, and our
// center node should no longer be part of the surface.
Expand All @@ -786,14 +911,11 @@ void C0Polyhedron::retriangulate()
libmesh_assert_not_equal_to
(elem->node_ptr(p), node);
#endif

// We've used up our center node, so it's not something we can
// eliminate again.
nodes_by_geometry.erase(geometry_it);
}

// At this point our surface should just have two triangles left.
libmesh_assert_equal_to(surface.n_elem(), 2);
if (!has_skipped_adding_tets_and_retriangulating)
libmesh_assert_equal_to(surface.n_elem(), 2);
}


Expand Down
Loading