Hex to Tet mesh generator

[aprilnovak]

Converts from a HEX8 mesh to a TET4 mesh by splitting each HEX8 element into 6 tets. Currently does not preserve sideset indexing, but wanted to get this here so that others could take it from here because I believe this is of larger interest as well.

Reason

Automate DAGMC model generation (which requires TET meshes) via MOOSE.

Design

Split each HEX8 into 6 tets.

Impact

New feature.

Refs #23822

[moosebuild]

Job Precheck on 18cd05c wanted to post the following:

Your code requires style changes.

A patch was auto generated and copied here
You can directly apply the patch by running, in the top level of your repository:

curl -s https://mooseframework.inl.gov/docs/PRs/25433/clang_format/style.patch | git apply -v

Alternatively, with your repository up to date and in the top level of your repository:

git clang-format d00bbe771ad86c50421707edcf7c889ccb955254

[moosebuild]

Job Documentation on 4ee7570 wanted to post the following:

View the site here

This comment will be updated on new commits.

[roystgnr]

So, this is probably my fault for not getting the Doxygen up to date ... but have you looked at libMesh::MeshModification::all_tri()?

[roystgnr]

all_tri() splits a hex into 24 tets, IIRC; not as efficient if your mesh is already fine enough, but it avoids conflicts in diagonal selection when you have a hex mesh that doesn't have a simple uniform orientation.

[aprilnovak]

Was not aware of that, thanks! But I am a bit concerned with restricting people to have 24 tets for every starting hex (for our use case, our hex meshes will be sufficiently fine already, and this converter is purely to get around a compatibility restriction in DAGMC).

What if we had multiple options of how many tets to create per hex? That way if users are concerned about the ambiguity in tet orientation going from 1 -> 6, they could have the 1 -> 24 option as well.

[roystgnr]

I do like the idea of having an option, but I wouldn't want to leave it up to users to figure out when the option is necessary. Either the 6-split works or it doesn't. Maybe we could either test the orientations as the 6-split is performed and make sure there are no inconsistencies, or test the neighbors afterward and make sure there are none missing?

[YaqiWang]

24 is a pretty big factor. I think it is appealing to have a minimum factor for converting HEX to TET.

[aprilnovak]

@roystgnr can you elaborate more on this? I'm not sure I follow what types of checks we should add:

make sure there are no inconsistencies, or test the neighbors afterward and make sure there are none missing?

[roystgnr]

I haven't looked deeply into the code, but on every side of every hex you ought to check elem->neighbor_ptr(s), and if it's not null then do things like check that it's also a hex and check that on its which_side_am_i side it's subdivision will be splitting the same diagonal.

Imagine a hex mesh on which your code gives a consistent tet mesh. Now rotate one of those hexes (by swapping node indices) 90 degrees about any of xi/eta/zeta, and give that to your code instead. You're hard-coding _tet4_nodes here, right? So the diagonal edges added on the 2 rotated faces of the hex will also be rotated, and four of the tets on that hex won't match their neighbors on those sides ... and libMesh find_neighbors won't find any matches, so it won't set any new neighbor_ptr values, so you'll have two quadrilateral slits in your domain there. This is worse than a segfault; at least with Undefined Behavior the compiler is only allowed to corrupt your output; if you set up a well-defined mesh with unintended slits it then the result is required to be wrong.

[moosebuild]

Job Coverage on b18b747 wanted to post the following:

Framework coverage

	848bef	#25433 b18b74
	Total	Total	+/-	New
Rate	85.70%	85.71%	+0.01%	98.18%
Hits	91485	91542	+57	54
Misses	15268	15266	-2	1

Diff coverage report

Full coverage report

Modules coverage

Coverage did not change

Full coverage reports

Reports

• framework
• chemical_reactions
• combined
• contact
• electromagnetics
• external_petsc_solver
• fluid_properties
• fsi
• functional_expansion_tools
• geochemistry
• heat_conduction
• level_set
• misc
• navier_stokes
• optimization
• peridynamics
• phase_field
• porous_flow
• ray_tracing
• rdg
• reactor
• richards
• scalar_transport
• solid_properties
• stochastic_tools
• tensor_mechanics
• thermal_hydraulics
• xfem

This comment will be updated on new commits.

[aprilnovak]

Thanks @roystgnr, this all makes sense to me now. I'll see what types of error checks can be included to protect against that possibility. But I am not sure how I would go about creating a mesh where 1 element is rotated relative to its original construction. Is there some way I could do that with MOOSE mesh generators?

[GiudGiud]

ElementGwnerator and TransformGenerator

[aprilnovak]

But to replicate the problem, wouldn't I need to have a rotated element within a larger mesh of other elements? I think ElementGenerator just makes one isolated element? Can it be used to build a contiguous mesh of elements one-by-one?

[GiudGiud]

You’d have to stitch them together.

Depending on what you want it might be easier to block restrict one element and use a transform (or a parsedNodeTransoform ? Any support for block restriction there?) on that block

[roystgnr]

I don't think there yet exists a MOOSE MeshGenerator using it, because there's no practical use for it other than unit testing, but libMesh::MeshTools::MeshModification::permute_elements() will pseudorandomly rotate and flip every element in a mesh - while leaving the mesh domain unchanged, because it's just local permutations, the same geometry but with different xi/eta/zeta->x/y/z maps. IIRC the underlying Elem::permute() also consistently permutes boundary conditions.

For a more practical albeit less comprehensive use case, maybe try SphereMeshGenerator? There are 12 planes in a sphere mesh where hexes attach to each other with slightly-unusual orientations, and I'll bet at least one of those ends up not working well with a fixed tet subdivision pattern.

[aprilnovak]

I was not able to get any catastrophic behavior using SphereMeshGenerator out of the box. I next tried adding that permutation call when generating the sphere mesh, like so:

std::unique_ptr<MeshBase>
SphereMeshGenerator::generate()
{
  /// ...

  libMesh::MeshTools::Modification::permute_elements(*mesh);

  return dynamic_pointer_cast<MeshBase>(mesh);
}

and then (i) generated a sphere mesh and (ii) passed it in to HexToTetMeshGenerator. Again I did not see any catastrophic behavior, and the MMS solution gives the expected convergence rates.

Is my approach incorrect in some way? I've added my permuted (I think) mesh to this PR, as permute_in.e.

[GiudGiud]

I haven't looked deeply into the code, but on every side of every hex you ought to check elem->neighbor_ptr(s), and if it's not null then do things like check that it's also a hex and check that on its which_side_am_i side it's subdivision will be splitting the same diagonal.

I think the way forward here could be to write a loop over the mesh that checks that the other-diagonal-split did not happen.

[roystgnr]

Perhaps you could count the number of hexes (H) and the number of hexes' quad faces that have neighbors (N) before the splitting, then after the splitting you do a find_neighbors() and count the number of tets' tri faces that have neighbors (T), so you can assert that T=2N+12*H (I think I got that 12 right; 12 internal tet faces per hex?). Anything that passes that should be a valid split, because if you do ever select two different diagonals on the same quad that'll be 4 tri faces which should have neighbors but don't.

[aprilnovak]

How does this neighbor thing work with AMR? Adaptively refined meshes do allow nodes to sit on element "edges", right?

[roystgnr]

If an element e has a side s where there's no neighboring element at the same e.level(), then the e.neighbor_ptr(s) points to the most refined (albeit still coarser than e) element on that side. Usually that's one level "up", but the level-one smoothing rule is something users can disable.

I didn't realize this generator was intended to work with refined meshes, though; was it? There's a comment with the word "parent" referring to the hex/tet relationship, for instance, which is not vocabulary I'd expect someone thinking about Elem::parent() to have overloaded. Doing hex-to-tet conversion correctly on a refined (even a uniformly refined!) mesh doesn't require merely converting every hex into a half dozen tets, and converting all child hexes consistently with their parents (which this PR actually does, if only by accident, I think!), it would require setting all the child_ptr and refinement_flag values on tets-of-parents too, plus the parent values on tets-of-children.

Once that's done, though ... I think T=2N+12H would still be correct to assert, even in the AMR/hanging-nodes cases.

[aprilnovak]

No, this is not intended for refined meshes - sorry for any confusion there. I guess I'm still trying to get a mental image of what exactly is the potential problem at hand. I'm not fluent with all the libmesh terminology, so it's a knowledge gap on my end.

Suppose I have two adjacent hexes (A,B), which touch each other though a common face f. If I split those hexes into tets, suppose on the original face f I have two tets from A, and two tets from B, but they're not sharing that mid-f side. Is that the nature of the problem? Or is it instead that some data structure is not initialized which I still need to initialize, but perhaps there's nothing fundamentally wrong with the mesh?

I am confused because everything seemed to work fine with the element permutation.

[roystgnr]

suppose on the original face f I have two tets from A, and two tets from B, but they're not sharing that mid-f side. Is that the nature of the problem?

They're sharing the side, but in two different ways. Suppose the face f is a Quad4 with nodes (1,2,3,4). That might end up being two tri sides with nodes (1,2,3) and (1,3,4) (or permutations of those, which are equivalent), or it might end up being two tri sides with nodes (1,2,4) and (2,3,4) (or permutations). Either way is fine, as long as it's consistent - if hex A has (1,2,3) and (1,3,4) sides while B has (2,1,3) and (4,3,1), you're great. But if hex A gets split into tets with (1,2,3) and (1,3,4) sides while hex B gets split into tets with (2,1,4) and (3,2,4) sides, now you have a slit in your mesh: the two tets from A will not be neighbors of the two tets from B. Instead of sharing a common diagonal edge, they have two distinct new edges crossing each other.

I am confused because everything seemed to work fine with the element permutation.

Was that because you definitely didn't add any slits to your meshes, or was that because libMesh thinks a slit mesh is a perfectly valid thing to have and presumes you just intended it that way? Having elements which touch but aren't neighbors is a perfectly valid thing to have in a domain, if you have a flat plate or some other infinitely thin obstruction which is supposed to act as a boundary from both sides. If you don't have that but you introduce its discrete equivalent to the mesh, though, then the very same mesh is wrong.

[github-actions[bot]]

This pull request has been automatically marked as stale because it has not had recent activity in the last 100 days. It will be closed in 7 days if no further activity occurs. Thank you for your contributions.