Implement insert points

examples/test_insert_points.py

@@ -0,0 +1,42 @@
+import logging
+import numpy as np
+
+from meshpy.tet import MeshInfo, Options, 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)