MeshGenerator for Breaking Mesh Along Embedded Side Sets

[Kavan-Khaledi]

Dear MOOSE Team,

I’m currently working on developing a MeshGenerator to break a mesh along embedded side sets (breakMeshByEmbeddSidesetGenerator). This will be particularly useful for applications such as cohesive zone modeling of arbitrarily oriented fractures embedded within a 3D subdomain.

To implement this, I’ve been following the logic used in BreakMeshByBlockGenerator:

1. Define a nodeToElementMap that maps the nodes on the side set to the elements that need to be renumbered.
2. Duplicate the nodes on the side set and renumber the elements on one side of the interface accordingly.
3. Create a new interface while preserving the element-neighbor connectivity information.

However, I’m encountering a couple of issues and would appreciate your guidance:

---

Question 1:
A key step in BreakMeshByBlockGenerator is disconnecting elements on the two sides of the interface (e.g., positive and negative sides). In that case, this separation is straightforward because the elements are distinguished by their subdomain IDs.

In my case, subdomain IDs are not applicable. My idea was to identify the positive and negative sides of the interface using the normal vector of the embedded side set. Specifically, for each element connected to the embedded side set, we can compute the face normal using the cross product of two in-plane vectors:

n = (p2 - p0) × (p1 - p0)

where p0, p1, and p2 are three nodal points of the face.

Here’s the issue: it seems that libMesh reorders the nodes of a face shared between two elements in a way that the computed normal vector is always inward-facing relative to the element. Is this understanding correct?

If so, is there a way to obtain a consistent face normal that is independent of node ordering? Alternatively, is there a recommended approach in libMesh for determining which side of an arbitrary interface an element lies on?

---

Question 2:
I'm looking for a robust way to exclude terminal nodes from the nodeToElementMap. These are the nodes at the ends of an embedded side set — for example, crack tips — where the interface terminates inside the domain. Although these nodes lie on the interface, they should not be treated as shared (and must not be duplicate

I didn’t find an existing method in Boundary_Info or the Elem class to detect such nodes. Do you have any recommendations on how best to identify and exclude them?

Best regards,
Kavan

[GiudGiud]

Hello

That mesh generator is a great idea. It would fit well into MOOSE

Here’s the issue: it seems that libMesh reorders the nodes of a face shared between two elements in a way that the computed normal vector is always inward-facing relative to the element. Is this understanding correct?

Alternatively, is there a recommended approach in libMesh for determining which side of an arbitrary interface an element lies on?

The normal vector for a sideset is always facing outwards (from primary to secondary side). This is the vector you want to use to distinguish the two sides of the mesh wrt to the sideset. You get that from the (elem, side, bid) tuple.

I didn’t find an existing method in Boundary_Info or the Elem class to detect such nodes. Do you have any recommendations on how best to identify and exclude them?

To find these nodes, you look at the nodeToElem map in mooseMesh (or build one if you only have a MeshBase). These nodes will have at least one of the elements in the nodeToElemMap that has no side that is part of the sideset

[Kavan-Khaledi]

Thank you, @GiudGiud, for your response. I forgot to mention in my initial post that I’m looking for a general solution that works for both TET4 and HEX8 elements.

With TET4 elements—unlike HEX8—there can be many elements on both sides of the sideset that only share a single node or edge with the sideset, rather than a full face. These elements are part of the sideset and still appear in my nodeToElemMap, which makes it difficult to distinguish them from the actual terminal elements (even when checking their neighbors).

I’ve managed to solve the first issue related to normal vector calculation, particularly for elements that do not have a full face on the sideset. However, I still haven’t found a good solution for the second issue.

A third issue I’ve encountered involves elements located at the intersection of two sidesets (e.g., two intersecting fractures). These elements may be considered primary for one sideset and secondary for the other, which leads to a contradiction that needs to be resolved. I’m skipping this issue for now, but I may return to it and reach out again if I can't figure it out.

Best regards,
Kavan

[GiudGiud]

With TET4 elements—unlike HEX8—there can be many elements on both sides of the sideset that only share a single node or edge with the sideset, rather than a full face. These elements are part of the sideset and still appear in my nodeToElemMap, which makes it difficult to distinguish them from the actual terminal elements (even when checking their neighbors).

You are looking for terminal nodes though, rather than terminal elements. You could sweep along the terminal edge. Find a node that is only connected to a single face in the sideset, that's a terminal node. Then find all other nodes that are:

• not on an external boundary of the mesh
• connected to the previous node by the edge of a face in the sideset, and connected to only one other face in the sideset

and repeat

A third issue I’ve encountered involves elements located at the intersection of two sidesets (e.g., two intersecting fractures). These elements may be considered primary for one sideset and secondary for the other, which leads to a contradiction that needs to be resolved. I’m skipping this issue for now, but I may return to it and reach out again if I can't figure it out.

Where is the problem in that? Is it for creating lowerD meshes on their faces?

[Kavan-Khaledi]

Dear GiudGiud, thank you for the solution! it was a very helpful suggestion to sweep along the edges rather than faces or nodes.
I’m happy that the two initial issues described in my first post are now solved.

below you can see a simple example (a box with two embedded sidesets under tensile loading).
With that, I would like to come back to the third issue, which I already mentioned earlier:
the problem of intersection points between embedded sidesets (or between different branches of one embedded sideset).

As shown in Fig. 3, a group of “patchy” elements appears along the intersection line.
This happens because the node duplication and element renumbering are incorrect for these elements.

These intersection points clearly need to be treated differently from regular interface points, because:

• During splitting of the interface, all other nodes must be duplicated only once, but intersection points require multiple copies (since they belong to different branches of the interface — see the sketch in the figure).

• Determining the primary/secondary side based on the normal vector becomes more complicated for the elements connected to the intersection line, especially in the situation where all embedded elements belong to the same sideset ID (pair option is not applicable).

So my questions are:

1. How can we reliably detect such intersection points?
   I guess any edge that is connected to more than 2 sideset faces could be considered an intersection edge — is that correct?

2. Once detected, what would be the correct strategy for duplicating the node at that intersection?
   In other words: how many copies need to be created in order to allow the individual interface branches to open/slide independently?

I would really appreciate your thoughts on this.

Best regards,
kavan

[GiudGiud]

How can we reliably detect such intersection points?

In your current setup, are you giving both boundaries to the "BreakMeshOnSideset(s)Generator"
For this next sentence (the cited one) below, I think we need the mesh generator to be aware of ALL boundaries

Have you tried the other way, chaining two "BreakMeshOnSideset(s)Generator" that are not aware of each other, and do you get the wrong result?

I guess any edge that is connected to more than 2 sideset faces could be considered an intersection edge — is that correct?

Connected to 2 sideset faces, one from each boundary, then yes I agree. If both faces from the same sideset, then no it's probably internal rather than an external edge. Just being pedantic, I m sure that's what you meant.

Once detected, what would be the correct strategy for duplicating the node at that intersection?
In other words: how many copies need to be created in order to allow the individual interface branches to open/slide independently?

Can you just do one copy per intersection sideset?

[GiudGiud]

note that the MeshDiagnosticsGenerator can be used to see the number and location of overlapping nodes, to see if things are being duplicated as expected

[Kavan-Khaledi]

Dear @GiudGiud, Following your suggestion, I tested two setups:

1. General Case: Single ID for All Sets with One Generator

In this setup, the user assigns a single ID to all embedded sidesets that need to be split, and the mesh is processed by one generator. The sideset can include multiple surfaces, branches, and intersections. A typical application is in Discrete Fracture Network (DFN) modeling, where numerous randomly distributed fractures intersect. The user assigns the same physical name (ID) to all fractures. An example input looks like this:

[mesh_Generator]
type = breakMeshOnSidesetsGenerators
input = mesh_file
Sideset = fracture_network
new_interface_name = interface
[]

Pros: Simplifies the mesh definition when many sidesets need to be split.

Cons: Managing intersections is challenging. Each intersection point requires three copies of the node. Since all sideset faces share the same ID, it is difficult to implement a general and robust solution that works across different element types and shapes.

A potential solution (not yet tested) is to internally group faces based on their normal vectors, essentially creating temporary sideset IDs and then iterating over these groups to split them individually (similar to Setup 2 below).

2. Single Generator with Multiple Sidesets (Each with a Unique ID)

In this configuration, i used a single generator to split the mesh. The sidesets are listed in a vector, each with a unique ID. Each sideset represents a flat or curved surface without branches or intersections. I tested this setup, and it works well: both intersecting and terminal nodes are handled correctly.

Because each embedded surface has its own ID, it is straightforward to identify intersection points and duplicate them properly by iterating over the sidesets and copying both the original and newly generated nodes.

[mesh_Generator]
type = breakMeshOnSidesetsGenerators
input = mesh_file
Sideset_list = 'fracture1 fracture2 fracture3 ...'
new_interfaces_list = 'interface1 interface2 interface3 ...'
[]

Pros: Easy to implement. Full control over node duplication and renumbering, especially at intersections, since everything is handled within a single generator.

Cons: For large models, the vector list of sidesets and interfaces can become very long.

For now, I will continue with Setup 2, although it would also be useful to have Setup 1 since it is more general. However, I am encountering a new issue! After splitting the mesh, if I use any other MeshGenerator that calls prepared_for_use() (e.g., in LowerDBlockFromSidesetGenerator), the element-neighbor information is deleted, and the cohesive interface no longer works.
When I remove the call to prepared_for_use(), the CZ interface becomes active again, and the lower-dimensional block is generated as expected. Do you know the reason? and how to avoid this?

Best regards,
Kavan

[GiudGiud]

Cons: For large models, the vector list of sidesets and interfaces can become very long.

dont worry about this. With input file includes, you can put this mesh generator in a separate file.

When I remove the call to prepared_for_use(), the CZ interface becomes active again, and the lower-dimensional block is generated as expected. Do you know the reason? and how to avoid this?

I think this is the same problem as the BreakMeshByBlockGenerator.
There is some special logic in MooseMesh to prevent that. Look for calls to has_bbmg which is checking for the use of that mesh generator