[Feature request] Mixed dimensional embedment meshing

[BinWang0213]

Is possible to perform mixed-dimensional meshing by embedment?
Such as: https://gitlab.onelab.info/gmsh/gmsh/blob/gmsh_4_7_1/tutorial/t15.geo

If so, we can mesh a geology model with fault and fractures, such as

[krober10nd]

Hey Bin, thanks for your interest! This is a intriguing idea. Firstly no, SM does not currently support mixed-dimensional embedment. I had thought about this at one point because faults are a big thing in seismology but never came up with a great solution.

It's a bit tricky. Theoretically, first thing to get this mixed-dimension embedment, we would need to immerse a 2D signed distance function in a 3D one. If that's possible, then it should also be possible to compute intersections of these mixed-dimensional geometry. And then if that's possible, it would be possible to immerse the set like I can do already in SM.

However, the problem with the 2D SDF in a 3D geometry is that the notion of "in" or "out" the domain isn't well defined anymore. In other words, implicit representation can only be used for closed surfaces.

One solution: it might be possible to mesh an open subset of a surface geometry using a second implicit function that defines the new boundaries (i.e., one SDF is used for the x-y plane and a second for the y-z plane).

As a current possible hack, what you could do is supply the points of the 2d immersed object as pfix to the generator function. This would give you the overall shape of the structure but some edges wouldn't be aligned with it.

I'll think about it more.

[krober10nd]

Hey @BinWang0213 do you have an example of a fault you'd like to constrain in a 3d domain? It'd be helpful to have a dataset to experiment with here.

[BinWang0213]

Thanks for your detailed explanations.

Sure. Here is an example model you can test with:

Box domain: [-0.6, -0.6, -0.6] -> [0.6, 0.6, 0.6]
Falults/Fractures: 4 fractures intersected together

It includes the paraview vtp files and gmsh input file:
input_files.zip

[krober10nd]

I haven't forgot about this! I'm nearly finished with the geometric primitives that are necessary to embed structures like this (translation, rotation, and stretching).

[BinWang0213]

Nice to hear from you! Looking forward to your updates!

[krober10nd]

@BinWang0213 I've added translation, rotation, and stretching to all the geometric primitives in the latest release SM 3.9 which now should allow you to get creative with these faults.

For example to create a similar example to the one you showed, I computed the union of several relatively thin cubes (i.e., fault-like) each one rotated and translated to resemble the image you gave above. You can visualize your creation with the domain.show() functionality. Note that all rotations and translations are relative to the origin so best build the mesh around the origin and then translate it later on.

Using the immerse option in generate_mesh I'm able to respect these structures more or less. It's not perfect but it's a start! To smooth intersections between faults, SDFs can be intersected with a smoothing parameter, which is not implemented yet.

 import meshio
 
 import SeismicMesh as sm
 
 parent = sm.Cube((-0.6, 0.6, -0.6, 0.6, -0.6, 0.6))
 
 cubes = []
 cubes.append(sm.Cube((-0.5, 0.5, -0.5, 0.5, 0.0, 0.1)))
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1), rotate=0.50 * 3.14, translate=[0.0, 0.0, 0.25]
     )
 )
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1), rotate=0.50 * 3.14, translate=[0, 0.25, -0.25]
     )
 )
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1),
         rotate=0.50 * 3.14,
         translate=[0, -0.25, -0.25],
     )
 )
 
 faults = sm.Union(cubes)
 # faults.show(samples=1e6)
 
 points, cells = sm.generate_mesh(
     domain=parent, edge_length=0.05, subdomains=[faults], verbose=2, max_iter=5
 )
 
 meshio.write_points_cells(
     "cube_with_faults.vtk",
     points,
     [("tetra", cells)],
     point_data={"sd": faults.eval(points)},
     file_format="vtk",
 )

[krober10nd]

I've added some more complexity to show how you can add mesh refinement near the immersed structures and ensure the mesh is sliver-free.

 import meshio
 
 import SeismicMesh as sm
 
 parent = sm.Cube((-0.6, 0.6, -0.6, 0.6, -0.6, 0.6))
 
 cubes = []
 cubes.append(sm.Cube((-0.5, 0.5, -0.5, 0.5, 0.0, 0.1)))
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1), rotate=0.50 * 3.14, translate=[0.0, 0.0, 0.25]
     )
 )
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1), rotate=0.50 * 3.14, translate=[0, 0.25, -0.25]
     )
 )
 cubes.append(
     sm.Cube(
         (-0.5, 0.5, -0.5, 0.5, 0.0, 0.1),
         rotate=0.50 * 3.14,
         translate=[0, -0.25, -0.25],
     )
 )
 
 faults = sm.Union(cubes)
 # faults.show(samples=1e6)
 
 
 def edge_length(x):
     return 0.05 + 0.15 * faults.eval(x)
 
 
 points, cells = sm.generate_mesh(
     domain=parent,
     h0=0.05,
     edge_length=edge_length,
     subdomains=[faults],
     verbose=2,
 )
 
 points, cells = sm.sliver_removal(
     points=points, domain=parent, h0=0.05, edge_length=edge_length, preserve=True
 )
 
 meshio.write_points_cells(
     "cube_with_faults.vtk",
     points,
     [("tetra", cells)],
     point_data={"sd": faults.eval(points)},
     file_format="vtk",
 )

[krober10nd]

Hey @BinWang0213 did you try my workaround? Using signed distance functions to define geometry constraints, this is the best we can get I think.

[BinWang0213]

@krober10nd Thanks a lot for your updates! I didn't have a chance to try your code as I'm working hard on my dissertation recently.

Does seismicMesh support define polygon faults/fractures below? It is not necessary a rectangle. In 3D it looks like a polyhedron slab.

[krober10nd]

yes, this is very similar to the case above! Wow talk about complex. We'd have to somehow describe each structure as a signed distance function, but in theory it's possible.