209

See code below....Final assert fails when pfix is not passed to the laplacian2 routine (which is correct). Make sure to deactivate mesh improvement in generate_mesh if you care about this stuff remaining fixed.

 import numpy as np
 import meshio

 import SeismicMesh as sm

 bbox = (0.0, 1.0, 0.0, 1.0)
 channel = sm.Rectangle(bbox)
 suspension = sm.Disk([0.5, 0.5], 0.25)

 pfix = np.array(
     [
         [0.4, 0.4],
         [0.45, 0.45],
         [0.5, 0.5],
         [0.55, 0.55],
         [0.6, 0.60],
         [0.65, 0.65],
         [0.70, 0.70],
     ]
 )
 hmin = 0.10
 fh = lambda p: 0.15 * np.abs(suspension.eval(p)) + hmin
 points, cells = sm.generate_mesh(
     domain=channel, edge_length=fh, h0=hmin, mesh_improvement=False, pfix=pfix
 )
 meshio.write_points_cells("before_mesh.vtk", points, [("triangle", cells)])

 points_w_fix, cells_w_fix = sm.geometry.laplacian2(points, cells, pfix=pfix)
 meshio.write_points_cells("mesh_w_fix.vtk", points_w_fix, [("triangle", cells_w_fix)])

 points_wo_fix, cells_wo_fix = sm.geometry.laplacian2(points, cells)
 meshio.write_points_cells(
     "mesh_wo_fix.vtk", points_wo_fix, [("triangle", cells_wo_fix)]
 )

 def _closest_node(node, nodes):
     nodes = np.asarray(nodes)
     deltas = nodes - node
     dist_2 = np.einsum("ij,ij->i", deltas, deltas)
     return dist_2

 for p in pfix:
     d = _closest_node(p, points)
     assert np.isclose(np.min(d), 0.0)

 for p in pfix:
     d = _closest_node(p, points_w_fix)
     assert np.isclose(np.min(d), 0.0)

 for p in pfix:
     d = _closest_node(p, points_wo_fix)
     assert np.isclose(np.min(d), 0.0)

krober10nd / SeismicMesh

fixing points in iterative laplacian2 #210

209