improvements (#171)

* performance improvements in parallel.

Author: Keith Roberts

README.md

@@ -593,7 +593,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## Unreleased
-- Adding basic periodic domains with the Repeat SDF.
+
+## [3.2.0] -2020-12-14
+- Adding basic periodic domains with the `Repeat` SDF.
+- Reworking CPP code and bottlenecks...20-30% faster `generate_mesh` in parallel for 2D/3D from previous versions.
+- `sliver_removal` has optional variable step size when perturbing vertices. Helps to remove the "last sliver".
+- Faster rectangle and cube primitives. 
 
 ## [3.1.7] - 2020-11-27
 ### Improved

SeismicMesh/generation/mesh_generator.py

@@ -114,6 +114,7 @@ def sliver_removal(points, domain, edge_length, comm=None, **kwargs):  # noqa: C
         "delta_t": 0.30,
         "geps_mult": 0.1,
         "subdomains": [],
+        "gamma": 1.0,
     }
 
     sliver_opts.update(kwargs)
@@ -170,6 +171,7 @@ def print_msg2(msg):
     )
 
     geps = sliver_opts["geps_mult"] * h0
+
     # deps = np.sqrt(np.finfo(np.double).eps) * h0
     min_dh_bound = sliver_opts["min_dh_angle_bound"] * math.pi / 180
     max_dh_bound = sliver_opts["max_dh_angle_bound"] * math.pi / 180
@@ -197,10 +199,12 @@ def print_msg2(msg):
 
     count = 0
     pold = None
-    push = 0.10
+    step = 0.10
+    gamma = sliver_opts["gamma"]
+    num_old_bad = np.inf
 
     dt = DT()
-    dt.insert(p.flatten().tolist())
+    dt.insert(p.ravel().tolist())
     while True:
 
         start = time.time()
@@ -209,7 +213,7 @@ def print_msg2(msg):
         # Using CGAL's incremental Delaunay triangulation capabilities.
         if count > 0:
             to_move = np.where(_dist(p, pold) > 0)[0]
-            dt.move(to_move.flatten().tolist(), p[to_move].flatten().tolist())
+            dt.move(to_move.ravel().tolist(), p[to_move].ravel().tolist())
 
         # Get the current topology of the triangulation
         p, t = _get_topology(dt)
@@ -224,7 +228,7 @@ def print_msg2(msg):
         # Number of iterations reached, stop.
         if count == (max_iter - 1):
             print_msg1(
-                "FAILURE: Termination...maximum number of iterations reached.",
+                "FAILURE: Termination...maximum number of iterations reached. Try increasing max_iter when generating the mesh",
             )
             p, t = _termination(p, t, sliver_opts, comm, sliver=True)
             break
@@ -266,11 +270,26 @@ def print_msg2(msg):
         perturb_norm = np.sum(np.abs(perturb) ** 2, axis=-1) ** (1.0 / 2)
         perturb /= perturb_norm[:, None]
 
-        # perturb push % of minimum mesh size
-        p[move] += push * h0 * perturb
+        num_bad = len(ele_nums)
+
+        if num_bad < num_old_bad:
+            # increase step
+            step /= gamma
+        elif num_bad > num_old_bad:
+            # decrease step
+            step *= gamma
+        elif num_bad == num_old_bad:
+            # algorithm appears stuck
+            # increase step
+            step /= 0.8
+
+        # perturb step % of minimum mesh size
+        p[move] += step * h0 * perturb
 
         count += 1
 
+        num_old_bad = len(ele_nums)
+
         end = time.time()
         if comm.rank == 0:
             print_msg2("     Elapsed wall-clock time %f : " % (end - start))
@@ -444,11 +463,11 @@ def print_msg2(msg):
         start = time.time()
 
         # Remove non-unique points
-        p = np.array(list(set(tuple(p) for p in p)))
+        # p = np.array(list(set(tuple(p) for p in p)))
 
         # (Re)-triangulation by the Delaunay algorithm
         dt = DT()
-        dt.insert(p.flatten().tolist())
+        dt.insert(p.ravel().tolist())
 
         # Get the current topology of the triangulation
         p, t = _get_topology(dt)
@@ -621,6 +640,7 @@ def _parse_kwargs(kwargs):
             "h0",
             "geps_mult",
             "subdomains",
+            "gamma",
         }:
             pass
         else:
@@ -696,7 +716,7 @@ def _add_ghost_vertices(p, t, dt, extents, comm):
     exports = migration.enqueue(extents, p, t, comm.rank, comm.size, dim=dim)
     recv = migration.exchange(comm, comm.rank, comm.size, exports, dim=dim)
     recv_ix = len(recv)
-    dt.insert(recv.flatten().tolist())
+    dt.insert(recv.ravel().tolist())
     p, t = _get_topology(dt)
     p, t, inv = geometry.remove_external_entities(
         p,

SeismicMesh/geometry/__init__.py

@@ -5,6 +5,8 @@
     calc_dihedral_angles,
     calc_volume_grad,
     unique_edges,
+    drectangle_fast,
+    dblock_fast,
 )
 from .signed_distance_functions import (
     Rectangle,
@@ -56,6 +58,8 @@
     "Ball",
     "Cube",
     "Cylinder",
+    "drectangle_fast",
+    "dblock_fast",
     "Torus",
     "Prism",
     "Disk",

SeismicMesh/geometry/cpp/fast_geometry.cpp

@@ -101,6 +101,118 @@ double dot_product(const std::array<double, 3> &vect_A,
   return product;
 }
 
+std::vector<double> c_dblock(std::vector<double> &points, double x1, double x2,
+                             double y1, double y2, double z1, double z2) {
+  // Signed distance function for block with corners (x1,y1,z1), (x2,y2,z2)
+  int sz = points.size() / 3;
+  std::vector<double> dist;
+  dist.resize(sz);
+  // return -min(
+  //    min(
+  //        min(min(min(-z1 + p[:, 2], z2 - p[:, 2]), -y1 + p[:, 1]), y2 - p[:,
+  //        1]), -x1 + p[:, 0],
+  //    ),
+  //    x2 - p[:, 0],
+  //)
+  for (std::size_t iv = 0; iv < sz; ++iv) {
+
+    double sum1 = -z1 + points[3 * iv + 2];
+    double sum2 = z2 - points[3 * iv + 2];
+    double sum3 = -y1 + points[3 * iv + 1];
+    double sum4 = y2 - points[3 * iv + 1];
+    double sum5 = -x1 + points[3 * iv];
+    double sum6 = x2 - points[3 * iv];
+
+    double tmp1 = std::min(sum1, sum2);
+    double tmp2 = std::min(tmp1, sum3);
+    double tmp3 = std::min(tmp2, sum4);
+    double tmp4 = std::min(tmp3, sum5);
+
+    dist[iv] = -std::min(tmp4, sum6);
+  }
+  return dist;
+}
+
+py::array dblock_fast(
+    py::array_t<double, py::array::c_style | py::array::forcecast> points,
+    double x1, double x2, double y1, double y2, double z1, double z2) {
+
+  // check input dimensions
+  int num_points = points.size() / 3;
+  // allocate std::vector (to pass to the C++ function)
+  std::vector<double> cppPts(num_points * 3);
+
+  // copy py::array -> std::vector
+  std::memcpy(cppPts.data(), points.data(), 3 * num_points * sizeof(double));
+
+  std::vector<double> dist = c_dblock(cppPts, x1, x2, y1, y2, z1, z2);
+
+  ssize_t sodble = sizeof(double);
+  std::vector<ssize_t> shape = {num_points};
+  std::vector<ssize_t> strides = {sodble};
+
+  // return 1-D NumPy array
+  return py::array(
+      py::buffer_info(dist.data(),    /* data as contiguous array  */
+                      sizeof(double), /* size of one scalar        */
+                      py::format_descriptor<double>::format(), /* data type */
+                      1,      /* number of dimensions      */
+                      shape,  /* shape of the matrix       */
+                      strides /* strides for each axis     */
+                      ));
+}
+
+std::vector<double> c_drectangle(std::vector<double> &points, double x1,
+                                 double x2, double y1, double y2) {
+  // Signed distance function for rectangle with corners (x1,y1), (x2,y1)
+  int sz = points.size() / 2;
+  std::vector<double> dist;
+  dist.resize(sz);
+  // return -min(min(min(-y1 + p[:, 1], y2 - p[:, 1]), -x1 + p[:, 0]), x2 - p[:,
+  // 0])
+  for (std::size_t iv = 0; iv < sz; ++iv) {
+
+    double sum1 = -y1 + points[2 * iv + 1];
+    double sum2 = y2 - points[2 * iv + 1];
+    double sum3 = -x1 + points[2 * iv];
+    double sum4 = x2 - points[2 * iv];
+
+    double tmp1 = std::min(sum1, sum2);
+    double tmp2 = std::min(tmp1, sum3);
+    dist[iv] = -std::min(tmp2, sum4);
+  }
+  return dist;
+}
+
+py::array drectangle_fast(
+    py::array_t<double, py::array::c_style | py::array::forcecast> points,
+    double x1, double x2, double y1, double y2) {
+
+  // check input dimensions
+  int num_points = points.size() / 2;
+  // allocate std::vector (to pass to the C++ function)
+  std::vector<double> cppPts(num_points * 2);
+
+  // copy py::array -> std::vector
+  std::memcpy(cppPts.data(), points.data(), 2 * num_points * sizeof(double));
+
+  std::vector<double> dist = c_drectangle(cppPts, x1, x2, y1, y2);
+
+  ssize_t sodble = sizeof(double);
+  std::vector<ssize_t> shape = {num_points};
+  std::vector<ssize_t> strides = {sodble};
+
+  // return 1-D NumPy array
+  return py::array(
+      py::buffer_info(dist.data(),    /* data as contiguous array  */
+                      sizeof(double), /* size of one scalar        */
+                      py::format_descriptor<double>::format(), /* data type */
+                      1,      /* number of dimensions      */
+                      shape,  /* shape of the matrix       */
+                      strides /* strides for each axis     */
+                      ));
+}
+
 std::vector<double> c_calc_dihedral_angles(std::vector<double> &points,
                                            std::vector<int> &cells) {
   // compute the 6 dihedral angles of all tetrahedrons in the mesh
@@ -456,6 +568,8 @@ py::array calc_circumsphere_grad(
 }
 
 PYBIND11_MODULE(fast_geometry, m) {
+  m.def("drectangle_fast", &drectangle_fast);
+  m.def("dblock_fast", &dblock_fast);
   m.def("unique_edges", &unique_edges);
   m.def("calc_volume_grad", &calc_volume_grad);
   m.def("calc_circumsphere_grad", &calc_circumsphere_grad);

SeismicMesh/geometry/signed_distance_functions.py

@@ -1,6 +1,8 @@
 import numpy as np
 import itertools
 
+from .cpp.fast_geometry import drectangle_fast, dblock_fast
+
 
 def _length(x):
     return np.sum(np.abs(x) ** 2, axis=-1) ** (1.0 / 2)
@@ -154,7 +156,9 @@ def __init__(self, bbox):
         self.bbox = bbox
 
     def eval(self, x):
-        return drectangle(x, self.bbox[0], self.bbox[1], self.bbox[2], self.bbox[3])
+        return drectangle_fast(
+            x, self.bbox[0], self.bbox[1], self.bbox[2], self.bbox[3]
+        )
 
 
 class Cube:
@@ -164,7 +168,7 @@ def __init__(self, bbox):
         self.bbox = bbox
 
     def eval(self, x):
-        return dblock(
+        return dblock_fast(
             x,
             self.bbox[0],
             self.bbox[1],
@@ -244,40 +248,6 @@ def drectangle(p, x1, x2, y1, y2):
     return -min(min(min(-y1 + p[:, 1], y2 - p[:, 1]), -x1 + p[:, 0]), x2 - p[:, 0])
 
 
-def dblock0(p, x1, x2, y1, y2, z1, z2):
-    # adapted from:
-    # https://github.com/nschloe/dmsh/blob/3305c417d373d509c78491b24e77409411aa18c2/dmsh/geometry/rectangle.py#L31
-    # outside dist
-    # https://gamedev.stackexchange.com/a/44496
-    w = x2 - x1
-    h = y2 - y1
-    d = z2 - z1
-    cx = (x1 + x2) / 2
-    cy = (y1 + y2) / 2
-    cz = (z1 + z2) / 2
-    dx = np.abs(p[:, 0] - cx) - w / 2
-    dy = np.abs(p[:, 1] - cy) - h / 2
-    dz = np.abs(p[:, 2] - cz) - d / 2
-    is_inside = (dx <= 0) & (dy <= 0) & (dz <= 0)
-    dx[dx < 0.0] = 0.0
-    dy[dy < 0.0] = 0.0
-    dz[dz < 0.0] = 0.0
-    dist = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
-    # inside dist
-    a = np.array(
-        [
-            p[is_inside, 0] - x1,
-            x2 - p[is_inside, 0],
-            p[is_inside, 1] - y1,
-            y2 - p[is_inside, 1],
-            p[is_inside, 2] - z1,
-            z2 - p[is_inside, 2],
-        ]
-    )
-    dist[is_inside] = -np.min(a, axis=0)
-    return dist
-
-
 def dblock(p, x1, x2, y1, y2, z1, z2):
     min = np.minimum
     return -min(