smooth sets (#207)

* adding smoothness factor to set ops

Co-authored-by: Keith Roberts <krober@usp.edu>

README.md

@@ -819,6 +819,10 @@ Changelog
 
 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
+### Added
+- Smoothed sets (e.g., intersections, differences, and unions)
+
 ## [3.5.0]-2021-03-09
 ### Added
 - Rotations for all geometric primitives

SeismicMesh/geometry/signed_distance_functions.py

@@ -240,8 +240,15 @@ def show(self, filename=None, samples=10000):
         _show(self, filename=None, samples=samples)
 
 
+def _loop_call(func, d):
+    tmp = d[0]
+    for i in range(0, len(d) - 1):
+        tmp = func(tmp, d[i + 1])
+    return tmp
+
+
 class Union:
-    def __init__(self, domains, smoothness=None):
+    def __init__(self, domains, smoothness=0.0):
         geom_dim = [d.dim for d in domains]
         assert np.all(geom_dim != 2) or np.all(geom_dim != 3)
         self.dim = geom_dim[0]
@@ -265,23 +272,27 @@ def __init__(self, domains, smoothness=None):
         self.corners = _gather_corners(domains)
         self.domains = domains
 
+    def _smooth_union(self, d1, d2):
+        h = np.maximum(self.k - np.abs(d1 - d2), 0.0)
+        return np.minimum(d1, d2) - np.divide(h * h * 0.25, self.k)
+
     def eval(self, x):
-        if self.k is None:
-            return np.minimum.reduce([d.eval(x) for d in self.domains])
+        d = [d.eval(x) for d in self.domains]
+        if self.k == 0.0:
+            return np.minimum.reduce(d)
         else:
-            d = [d.eval(x) for d in self.domains]
-            h = np.maximum(self.k - np.abs(d[0] - d[1]), 0.0)
-            return np.minimum(d[0], d[1]) - np.divide(h * h * 0.25, self.k)
+            return _loop_call(self._smooth_union, d)
 
     def show(self, filename=None, samples=10000):
         _show(self, filename=None, samples=samples)
 
 
 class Intersection:
-    def __init__(self, domains):
+    def __init__(self, domains, smoothness=0.0):
         geom_dim = [d.dim for d in domains]
         assert np.all(geom_dim != 2) or np.all(geom_dim != 3)
         self.dim = geom_dim[0]
+        self.k = smoothness
         if self.dim == 2:
             self.bbox = (
                 min(d.bbox[0] for d in domains),
@@ -301,18 +312,27 @@ def __init__(self, domains):
         self.corners = _gather_corners(domains)
         self.domains = domains
 
+    def _smooth_intersection(self, d1, d2):
+        h = np.maximum(self.k - np.abs(d1 - d2), 0.0)
+        return np.maximum(d1, d2) + h * h * 0.25 / self.k
+
     def eval(self, x):
-        return np.maximum.reduce([d.eval(x) for d in self.domains])
+        d = [d.eval(x) for d in self.domains]
+        if self.k == 0.0:
+            return np.maximum.reduce(d)
+        else:
+            return _loop_call(self._smooth_intersection, d)
 
     def show(self, filename=None, samples=10000):
         _show(self, filename=None, samples=samples)
 
 
 class Difference:
-    def __init__(self, domains):
+    def __init__(self, domains, smoothness=0.0):
         geom_dim = [d.dim for d in domains]
         assert np.all(geom_dim != 2) or np.all(geom_dim != 3)
         self.dim = geom_dim[0]
+        self.k = smoothness
         if self.dim == 2:
             self.bbox = (
                 min(d.bbox[0] for d in domains),
@@ -332,10 +352,18 @@ def __init__(self, domains):
         self.corners = _gather_corners(domains)
         self.domains = domains
 
+    def _smooth_difference(self, d1, d2):
+        h = np.maximum(self.k - np.abs(-d1 - d2), 0.0)
+        return np.maximum(-d1, d2) + np.divide(h * h * 0.25, self.k)
+
     def eval(self, x):
-        return np.maximum.reduce(
-            [-d.eval(x) if n > 0 else d.eval(x) for n, d in enumerate(self.domains)]
-        )
+        if self.k == 0.0:
+            return np.maximum.reduce(
+                [-d.eval(x) if n > 0 else d.eval(x) for n, d in enumerate(self.domains)]
+            )
+        else:
+            d = [d.eval(x) for d in self.domains]
+            return _loop_call(self._smooth_difference, d[::-1])
 
     def show(self, filename=None, samples=10000):
         _show(self, filename=None, samples=samples)

tests/test_smooth_sets.py

@@ -0,0 +1,26 @@
+import numpy as np
+import pytest
+
+import SeismicMesh as sm
+
+
+@pytest.mark.serial
+def test_smooth_diff():
+    cube1 = sm.Cube((-0.5, 0.5, -0.5, 0.5, -0.5, 0.5))
+    ball1 = sm.Ball((0.0, 0.0, 0.5), 0.85)
+
+    domain = sm.Difference([ball1, cube1], smoothness=0.20)
+    points, cells = sm.generate_mesh(
+        domain=domain,
+        edge_length=0.10,
+    )
+    points, cells = sm.sliver_removal(
+        points=points,
+        domain=domain,
+        edge_length=0.10,
+    )
+    assert np.abs(cells.shape[0] - 9004) < 100
+
+
+if __name__ == "__main__":
+    test_smooth_diff()