support for immersed domains

[krober10nd]

• Let us consider a domain ω totally included in the larger one, Ω. It will be said that ω is immersed into Ω. More precisely, if Γ = δω is its boundary, immersion means that ω ∩ Ω = ω and Γ ∩Ω = Γ .

• A 0-level set can be approximately immersed inside a larger domain through the changes in this PR.

• A new name value pair for generate_mesh is used called subdomain which is a signed distance function that is completely immersed inside the parent domain.

• A multi-Newton method is called to reproject interior points onto the 0-level set of the subdomain each meshing iteration.

• Along the surface of the immersed domain quality appears to suffer. Not sure why yet.

• MWE

 import SeismicMesh
 import meshio

 ball = SeismicMesh.Ball([0.0, 0.0, 0.0], 1.0)
 ball2 = SeismicMesh.Ball([0.0, 0.0, 0.0], 0.5)

 points, cells = SeismicMesh.generate_mesh(
     domain=ball,
     bbox=(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0),
     edge_length=0.10,
     subdomain=[ball2],
 )
 meshio.write_points_cells(
     "ball_inside_ball.vtk",
     points,
     [("tetra", cells)],
     file_format="vtk",
 )

[krober10nd]

 import SeismicMesh
 import meshio

 ball = SeismicMesh.Ball([0.0, 0.0, 0.0], 1.0)
 ball2 = SeismicMesh.Ball([0.0, 0.0, 0.0], 0.5)

 points, cells = SeismicMesh.generate_mesh(
     domain=ball,
     edge_length=0.10,
     subdomain=[ball2],
 )

 meshio.write_points_cells(
     "ball_inside_ball.vtk",
     points,
     [("tetra", cells)],
     point_data={"sd": ball2.eval(points)},
     file_format="vtk",
 )

[krober10nd]

 import SeismicMesh
 import meshio

 ball = SeismicMesh.Ball([0.0, 0.0, 0.0], 1.0)
 cube2 = SeismicMesh.Cube((-0.5, 0.5, -0.5, 0.5, -0.5, 0.5))
 torus = SeismicMesh.Torus(r1=0.25, r2=0.10)

 hmin = 0.01

 points, cells = SeismicMesh.generate_mesh(
     domain=ball,
     edge_length=lambda x: hmin + 0.1 * (x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2),
     h0=hmin,
     subdomains=[torus],
 )

 meshio.write_points_cells(
     "torus_inside_ball.vtk",
     points,
     [("tetra", cells)],
     point_data={"sd": torus.eval(points)},
     file_format="vtk",
 )

[krober10nd]

• improving boundary intersections by removing "outside" simplices at the end after improving the boundary, but some problems still remain.

 import SeismicMesh
 import meshio

 cube1 = SeismicMesh.Cube((-1.0, 1.0, -1.0, 1.0, -1.0, 1.0))
 cube2 = SeismicMesh.Cube((-2.0, 0.5, -2.5, 0.5, -2.5, 0.5))

 hmin = 0.10

 union = SeismicMesh.Union([cube1, cube2])

 points, cells = SeismicMesh.generate_mesh(
     domain=union,
     edge_length=hmin,  # lambda x: hmin + 0.05 * (x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2),
     h0=hmin,
     subdomains=[cube2],
     verbose=2,
     max_iter=100,
 )

 meshio.write_points_cells(
     "out.vtk",
     points,
     [("tetra", cells)],
     point_data={"sd": cube2.eval(points)},
     file_format="vtk",
 )

[krober10nd]

In 2D

 import meshio
 import SeismicMesh

 box0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0))
 disk0 = SeismicMesh.Disk([0.5, 0.5], 0.25)

 points, cells = SeismicMesh.generate_mesh(
     domain=box0, edge_length=0.10, subdomains=[disk0]
 )
 meshio.write_points_cells(
     "immersed_disk.vtk",
     points,
     [("triangle", cells)],
     file_format="vtk",
 )

[krober10nd]

@jorgensd perhaps this feature would be useful for topology optimization? maybe you have some suggestions?

[lgtm-com[bot]]

This pull request fixes 1 alert when merging 4c8953c339c2fb552b8e1539b2106e2f29c54ecd into 6bbd5e576460d2771cdb750f0c62d8df283822fe - view on LGTM.com

fixed alerts:

• 1 for Unhashable object hashed

[jorgensd]

@jorgensd perhaps this feature would be useful for topology optimization? maybe you have some suggestions?

This is quite interesting, not only for topology optimization, but for multiphysics problems as fluid structure interaction, and multi material problems in for instance elasticity.

It can of course also be used in topology/shape optimization to track domain changes. Will it be possible to have multiple level sets intersecting with a background level set (this would be something similar to gmsh's boolean fragments feature)?

[krober10nd]

Will it be possible to have multiple level sets intersecting with a background level set (this would be something similar to gmsh's boolean fragments feature)?

Maybe. I'm not quite sure I understand what that Gmsh boolean feature does.

In my case, the subdomain to be immersed is just simply a SDF which avoids the need to have explicit geometry information regarding its immersed boundary.

So in the example above I just create a disk but I could also immerse the union/intersection/difference of any number of SDFs using the SeismicMesh.geometry stuff:

 import meshio

 import SeismicMesh

 box0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0))
 disk0 = SeismicMesh.Disk([0.5, 0.5], 0.25)
 sqr0 = SeismicMesh.Rectangle((0.5, 0.75, 0.5, 0.75))

 domain = SeismicMesh.Union([disk0, sqr0])

 hmin = 0.0025

 fh = lambda p: 0.05 * np.abs(domain.eval(p)) + hmin

 points, cells = SeismicMesh.generate_mesh(
     domain=box0,
     edge_length=fh,
     h0=hmin,
     subdomains=[domain],
 )

 pmid = points[cells].sum(1) / 3  # Compute centroids
 sd = domain.eval(pmid)

 tags = np.zeros((len(cells))) + 10
 tags[sd < 0] = 11

 meshio.write_points_cells(
     "immersed_disk.msh",
     points,
     cells=[("triangle", cells)],
     field_data={"OutsideDisk": np.array([10, 2]), "InsideDisk": np.array([11, 2])},
     cell_data={
         "gmsh:physical": np.array([tags]),  # np.array([np.repeat(11, len(cells))]),
         "gmsh:geometrical": np.array([tags]),  # np.array([np.repeat(1, len(cells))]),
     },
     file_format="gmsh22",
     binary=False,
 )

[lgtm-com[bot]]

This pull request fixes 1 alert when merging 08044d4744962445ab0b3f576b7c3d479b433602 into 6bbd5e576460d2771cdb750f0c62d8df283822fe - view on LGTM.com

fixed alerts:

• 1 for Unhashable object hashed