Implement insert points

examples/test_insert_points.py

@@ -0,0 +1,42 @@
+import logging
+import numpy as np
+
+from meshpy.tet import MeshInfo, build
+
+if __name__ == "__main__":
+    logger = logging.getLogger("test_insert_points.py")
+
+    points = [(0, 0, 0),
+              (0, 0, 1),
+              (1, 0, 0),
+              (1, 0, 1),
+              (1, 1, 0),
+              (1, 1, 1),
+              (0, 1, 0),
+              (0, 1, 1)]
+
+    facets = [(0, 1, 3, 2),
+              (2, 3, 5, 4),
+              (4, 5, 7, 6),
+              (6, 7, 1, 0),
+              (0, 2, 4, 6),
+              (1, 3, 5, 7)]
+
+    mesh_info = MeshInfo()
+    mesh_info.set_points(points)
+    mesh_info.set_facets(facets)
+
+    # Insert an interior point of the cube as a constrained point
+    interior_point = (0.33, 0.7, 0.91)
+
+    insert_points_mesh_info = MeshInfo()
+    insert_points_mesh_info.set_points([interior_point])
+
+    mesh = build(mesh_info, max_volume=0.25,
+                 insert_points=insert_points_mesh_info)
+
+    mesh_points = np.array(mesh.points)
+    min_dist = np.sqrt(np.sum((interior_point - mesh_points), axis=1)**2).min()
+
+    if min_dist > 0:
+        logger.error("tetrahedron mesh does not contain contrained point")

meshpy/tet.py

@@ -141,7 +141,7 @@ def __init__(self, switches, **kwargs):
                 setattr(self, k, v)
 
 
-def tetrahedralize(mesh_info, options):
+def tetrahedralize(mesh_info, options, insert_points=None):
     mesh = MeshInfo()
 
     # restore "C" locale--otherwise tetgen might mis-parse stuff like "a0.01"
@@ -155,7 +155,7 @@ def tetrahedralize(mesh_info, options):
         locale.setlocale(locale.LC_NUMERIC, "C")
 
     try:
-        internals.tetrahedralize(options, mesh_info, mesh)
+        internals.tetrahedralize(options, mesh_info, mesh, insert_points)
     finally:
         # restore previous locale if we've changed it
         if have_locale:
@@ -171,6 +171,8 @@ def build(mesh_info, options=Options("pq"), verbose=False,
         options.quiet = 1
 
     if insert_points is not None:
+        if not isinstance(insert_points, MeshInfo):
+            raise ValueError('insert_points should be an instance of MeshInfo')
         options.insertaddpoints = 1
 
     if attributes:
@@ -183,4 +185,4 @@ def build(mesh_info, options=Options("pq"), verbose=False,
     if diagnose:
         options.diagnose = 1
 
-    return tetrahedralize(mesh_info, options)
+    return tetrahedralize(mesh_info, options, insert_points)

test/test_meshpy.py

@@ -162,6 +162,56 @@ def test_tetgen_points():
 
     #mesh.write_vtk("test.vtk")
 
+
+def test_tetgen_insert_points():
+    import numpy as np
+
+    from meshpy.tet import MeshInfo, build
+
+    points = [(0, 0, 0),
+              (0, 0, 1),
+              (1, 0, 0),
+              (1, 0, 1),
+              (1, 1, 0),
+              (1, 1, 1),
+              (0, 1, 0),
+              (0, 1, 1)]
+
+    facets = [(0, 1, 3, 2),
+              (2, 3, 5, 4),
+              (4, 5, 7, 6),
+              (6, 7, 1, 0),
+              (0, 2, 4, 6),
+              (1, 3, 5, 7)]
+
+    mesh_info = MeshInfo()
+    mesh_info.set_points(points)
+    mesh_info.set_facets(facets)
+
+    interior_point = (0.33, 0.7, 0.91)
+
+    insert_points_mesh_info = MeshInfo()
+    insert_points_mesh_info.set_points([interior_point])
+
+    mesh = build(mesh_info, max_volume=0.25,
+                 insert_points=insert_points_mesh_info)
+
+    mesh_points = np.array(mesh.points)
+
+    # because of the max_volume==0.25 constraint, TetGen should have
+    # inserted some points in addition to the 8 cube vertices and the
+    # constranied interior point
+    assert mesh_points.shape[0] > 9
+
+    interior_point_index = \
+        np.where((mesh_points == np.array(interior_point)).all(1))[0][0]
+
+    mesh_elements = np.array(mesh.elements, dtype=int)
+
+    # make sure that there cells including the constrained point have
+    # been added to the mesh
+    assert (mesh_elements == interior_point_index).any()
+
 # }}}