L-shaped in SeismicMesh

[krober10nd]

Adding L-shaped domains for SeismicMesh 3.0.5. The 3D L-shaped domain indeed has some trouble with the pockets between SDFs so that remains to be improved but it's a valid mesh with a bounded minimum dihedral angle.

Note: min() and max() of two SDFs do not always produce a valid SDF. The result is actually not a distance field, but just a lower bound to the real distance to the resulting surface. Here I use the union to approximate the domain which results in errors "deep" in the interior of the geometry but this okay for meshing in these examples...

• 2D

 import numpy as np
 import meshio
 from SeismicMesh import Rectangle, generate_mesh, geometry
 
 hmin = 0.05
 bbox = (0.0, 1.0, 0.0, 1.0)
 
 
 bbox0 = (0.0, 1.0, 0.0, 0.5)
 rect0 = Rectangle(bbox0)
 
 bbox1 = (0.0, 0.5, 0.0, 1.0)
 rect1 = Rectangle(bbox1)
 
 
 corners = geometry.corners
 
 
 def union(x):
     return np.minimum.reduce([rect0.eval(x), rect1.eval(x)])
 
 
 pfix = np.vstack((corners(bbox0), corners(bbox1)))
 
 points, cells = generate_mesh(
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     pfix=pfix,
     verbose=2,
 )
 
 meshio.write_points_cells(
     "L_shape_2d.vtk",
     points,
     [("triangle", cells)],
 )

• 3D

 from SeismicMesh import Cube, generate_mesh, sliver_removal, geometry
 
 hmin = 0.05
 bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
 
 tol = hmin/10
 
 
 bbox0 = (0.0, 1.0, 0.0, 1.0, 0.0, 0.50 + tol)
 cube0 = Cube(bbox0)
 
 bbox1 = (0.0, 1.0, 0.5, 1.0, 0.50 - tol, 1.0)
 cube1 = Cube(bbox1)
 
 bbox2 = (0.5 - tol, 1.0, 0.0, 1.0, 0.50 - tol, 1.0)
 cube2 = Cube(bbox2)
 
 corners = geometry.corners
 
 
 def union(x):
     return np.minimum.reduce([cube0.eval(x), cube1.eval(x), cube2.eval(x)])
 
 
 pfix = np.vstack((corners(bbox0), corners(bbox1), corners(bbox2)))
 
 points, cells = generate_mesh(
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     pfix=pfix,
     verbose=2,
 )
 
 points, cells = sliver_removal(
     points=points,
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     verbose=2,
 )
 
 meshio.write_points_cells(
     "L_shape.vtk",
     points,
     [("tetra", cells)],
 )

[nschloe]

Blocked by krober10nd/SeismicMesh#102.

[krober10nd]

Okay, now with SM 3.0.6 this is a lot simpler thanks to your helpful comments:

• 2D

 import meshio
 from SeismicMesh import Rectangle, generate_mesh, Union
 
 hmin = 0.01
 rect0 = Rectangle((0.0, 1.0, 0.0, 0.5))
 rect1 = Rectangle((0.0, 0.5, 0.0, 1.0))
 union = Union([rect0, rect1])
 points, cells = generate_mesh(
     domain=union,
     edge_length=hmin,
 )
 meshio.write_points_cells(
     "L_shape_2d.vtk",
     points,
     [("triangle", cells)],
 )

• 3D

 import meshio
 from SeismicMesh import Cube, generate_mesh, sliver_removal, Union
 
 hmin = 0.05
 tol = hmin / 10
 cube0 = Cube((0.0, 1.0, 0.0, 1.0, 0.0, 0.50 + tol))
 cube1 = Cube((0.0, 1.0, 0.5, 1.0, 0.50 - tol, 1.0))
 cube2 = Cube((0.5 - tol, 1.0, 0.0, 1.0, 0.50 - tol, 1.0))
 union = Union([cube0, cube1, cube2])
 points, cells = generate_mesh(
     domain=union,
     edge_length=hmin,
 )
 points, cells = sliver_removal(
     points=points,
     domain=union,
     edge_length=hmin,
 )
 meshio.write_points_cells(
     "L_shape.vtk",
     points,
     [("tetra", cells)],
 )

[nschloe]

Got it.