improvements to geometry API (#106)

* adding Union/Intersection/Difference classes for geometric primitives.
* adding a ball SDF
* Automatic corner constraints in serial

README.md

@@ -291,8 +291,7 @@ import SeismicMesh
 
 bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
 cube = SeismicMesh.Cube(bbox)
-corners = SeismicMesh.geometry.corners(bbox)
-points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05, pfix=corners)
+points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05)
 points, cells = SeismicMesh.sliver_removal(points=points, domain=cube, edge_length=0.05)
 meshio.write_points_cells(
     "cube.vtk",
@@ -311,8 +310,7 @@ import SeismicMesh
 
 bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
 cube = SeismicMesh.Cube(bbox)
-corners = SeismicMesh.geometry.corners(bbox)
-points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05, pfix=corners)
+points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05)
 meshio.write_points_cells(
     "cube.vtk",
     points,
@@ -401,6 +399,94 @@ meshio.write_points_cells(
 )
 ```
 
+![Union](https://user-images.githubusercontent.com/18619644/96755280-e3aaff00-13a8-11eb-9f88-95a6684e928b.png)
+
+
+```python
+# Compute the union of several SDFs to create more complex geometries
+import meshio
+import SeismicMesh
+
+h = 0.10
+rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 0.5))
+rect1 = SeismicMesh.Rectangle((0.0, 0.5, 0.0, 1.0))
+disk0 = SeismicMesh.Disk([0.5, 0.5], 0.5)
+union = SeismicMesh.Union([rect0, rect1, disk0])
+points, cells = SeismicMesh.generate_mesh(domain=union, edge_length=h)
+meshio.write_points_cells(
+    "Lshape_wDisk.vtk",
+    points,
+    [("triangle", cells)],
+    file_format="vtk",
+)
+```
+![Leaf](https://user-images.githubusercontent.com/18619644/96755336-f6bdcf00-13a8-11eb-99ec-bd7e7d9cad1d.png)
+
+```python
+# Compute the intersection of several SDFs to create more complex geometries
+import meshio
+import SeismicMesh
+
+h = 0.05
+rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0))
+disk0 = SeismicMesh.Disk([0.25, 0.25], 0.5)
+disk1 = SeismicMesh.Disk([0.75, 0.75], 0.5)
+intersection = SeismicMesh.Intersection([rect0, disk0, disk1])
+points, cells = SeismicMesh.generate_mesh(domain=intersection, edge_length=h)
+meshio.write_points_cells(
+    "Leaf.vtk",
+    points,
+    [("triangle", cells)],
+    file_format="vtk",
+)
+```
+
+![Hole](https://user-images.githubusercontent.com/18619644/96766828-0ab9fe80-13b2-11eb-8bca-6306934008d4.png)
+
+```python
+# Compute the difference of two SDFs to create more complex geometries.
+import meshio
+import SeismicMesh
+
+h = 0.05
+rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0))
+disk0 = SeismicMesh.Disk([0.5, 0.5], 0.1)
+disk1 = SeismicMesh.Disk([0.75, 0.75], 0.20)
+difference = SeismicMesh.Difference([rect0, disk0, disk1])
+points, cells = SeismicMesh.generate_mesh(domain=difference, edge_length=h)
+meshio.write_points_cells(
+    "Hole.vtk",
+    points,
+    [("triangle", cells)],
+    file_format="vtk",
+)
+```
+
+![Cube_wHoles](https://user-images.githubusercontent.com/18619644/96785337-0a772e80-13c5-11eb-88fb-311b5bfdfed4.png)
+
+```python
+import meshio
+import SeismicMesh
+
+h = 0.10
+cube0 = SeismicMesh.Cube((0.0, 1.0, 0.0, 1.0, 0.0, 1.0))
+ball1 = SeismicMesh.Ball([0.5, 0.0, 0.5], 0.30)
+ball2 = SeismicMesh.Ball([0.5, 0.5, 0.0], 0.30)
+ball3 = SeismicMesh.Ball([0.0, 0.5, 0.5], 0.30)
+ball4 = SeismicMesh.Ball([0.5, 0.5, 0.5], 0.45)
+difference = SeismicMesh.Difference([cube0, ball1, ball2, ball3, ball4])
+points, cells = SeismicMesh.generate_mesh(domain=difference, edge_length=h, verbose=1)
+points, cells = SeismicMesh.sliver_removal(
+    points=points, domain=difference, edge_length=h, verbose=1
+)
+meshio.write_points_cells(
+    "Ball_wHoles.vtk",
+    points,
+    [("tetra", cells)],
+    file_format="vtk",
+)
+```
+
 How does performance and cell quality compare to `gmsh` and `cgal` mesh generators?
 ===================================================================================
 
@@ -428,12 +514,16 @@ Changelog
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project (tries to) adhere to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-### Unreleased
-- Added more examples on README
+## Unreleased
+
 ### Fixed
 - Silence messages about pfix when verbose=0
+### Added
+- Added more examples on README
+- New unions/intersections/differences with several SDF primivitives
+- Automatic corner constraints in serial
 
-### [3.0.5] - 2020-10-18
+## [3.0.5] - 2020-10-18
 
 ### Fixed
 - Preserving fixed points in serial.
@@ -446,7 +536,7 @@ and this project (tries to) adhere to [Semantic Versioning](https://semver.org/s
 - More support for reading binary files packed in a binary format.
 - Check to make sure bbox is composed of all floats.
 
-### [3.0.4] - 2020-10-12
+## [3.0.4] - 2020-10-12
 
 ### Added
 - Improve conformity of level-set in final mesh through additional set of Newton boundary projection iterations.

SeismicMesh/__init__.py

@@ -7,7 +7,7 @@
 #  see <http://www.gnu.org/licenses/>.
 
 from . import decomp, geometry, migration
-from .geometry import Disk, Cube, Rectangle
+from .geometry import Ball, Disk, Cube, Rectangle, Union, Intersection, Difference
 from .generation import generate_mesh, sliver_removal
 from .sizing import (
     get_sizing_function_from_segy,
@@ -23,6 +23,10 @@
     "Rectangle",
     "Cube",
     "Disk",
+    "Union",
+    "Ball",
+    "Intersection",
+    "Difference",
     "get_sizing_function_from_segy",
     "write_velocity_model",
     "plot_sizing_function",

SeismicMesh/generation/mesh_generator.py

@@ -133,7 +133,7 @@ def print_msg2(msg):
         if comm.rank == 0:
             raise Exception("Mesh improvement currently on works in 3D")
 
-    fd, bbox0 = _unpack_domain(domain, sliver_opts)
+    fd, bbox0, _ = _unpack_domain(domain, sliver_opts)
 
     fh, bbox1, hmin = _unpack_sizing(edge_length)
 
@@ -354,7 +354,7 @@ def print_msg2(msg):
         print(msg, flush=True)
 
     # unpack domain
-    fd, bbox0 = _unpack_domain(domain, gen_opts)
+    fd, bbox0, corners = _unpack_domain(domain, gen_opts)
 
     fh, bbox1, hmin = _unpack_sizing(edge_length)
 
@@ -403,6 +403,9 @@ def print_msg2(msg):
     DT = _select_cgal_dim(dim)
 
     pfix, nfix = _unpack_pfix(dim, gen_opts, comm)
+    if corners is not None and comm.size == 1:
+        pfix = np.append(pfix, corners, axis=0)
+        nfix = len(pfix)
     if comm.rank == 0:
         print_msg1("Constraining " + str(nfix) + " fixed points..")
 
@@ -545,23 +548,29 @@ def func(x):
 
 
 def _unpack_domain(domain, opts):
-    if isinstance(domain, geometry.Rectangle):
-        bbox = (domain.x1, domain.x2, domain.y1, domain.y2)
-        fd = domain.eval
-    elif isinstance(domain, geometry.Cube):
-        bbox = (domain.x1, domain.x2, domain.y1, domain.y2, domain.z1, domain.z2)
-        fd = domain.eval
-    elif isinstance(domain, geometry.Disk):
-        bbox = (domain.x1, domain.x2, domain.y1, domain.y2)
+    corners = None
+    domains = [
+        geometry.Ball,
+        geometry.Rectangle,
+        geometry.Disk,
+        geometry.Cube,
+        geometry.Cube,
+        geometry.Union,
+        geometry.Intersection,
+        geometry.Difference,
+    ]
+    if np.any([isinstance(domain, d) for d in domains]):
+        bbox = domain.bbox
         fd = domain.eval
+        corners = domain.corners
     elif callable(domain):
         # get the bbox from the name value pairs or quit
         bbox = opts["bbox"]
         fd = domain
     else:
         raise ValueError("`domain` must be a function or a :class:`geometry` object")
     _check_bbox(bbox)
-    return fd, bbox
+    return fd, bbox, corners
 
 
 def _parse_kwargs(kwargs):
@@ -615,7 +624,6 @@ def _get_edges(t):
     return geometry.unique_edges(edges)
 
 
-# @profile
 def _compute_forces(p, t, fh, h0, L0mult):
     """Compute the forces on each edge based on the sizing function"""
     dim = p.shape[1]
@@ -792,12 +800,9 @@ def _dist(p1, p2):
 
 def _unpack_pfix(dim, opts, comm):
     """Unpack fixed points"""
-    if opts["pfix"] is not None:
-        if comm.size > 1:
-            raise Exception("Fixed points are not yet supported in parallel.")
-        else:
-            pfix = np.array(opts["pfix"], dtype="d")
-            nfix = len(pfix)
+    if opts["pfix"] is not None and comm.size == 1:
+        pfix = np.array(opts["pfix"], dtype="d")
+        nfix = len(pfix)
     else:
         pfix = np.empty((0, dim))
         nfix = 0

SeismicMesh/geometry/__init__.py

@@ -6,7 +6,18 @@
     calc_volume_grad,
     unique_edges,
 )
-from .signed_distance_functions import Rectangle, Cube, Disk, drectangle, dblock
+from .signed_distance_functions import (
+    Rectangle,
+    Cube,
+    Disk,
+    Ball,
+    drectangle,
+    dblock,
+    corners,
+    Union,
+    Intersection,
+    Difference,
+)
 from .utils import (
     calc_re_ratios,
     delete_boundary_entities,
@@ -28,7 +39,6 @@
     unique_rows,
     vertex_to_entities,
     vertex_in_entity3,
-    corners,
 )
 
 __all__ = [
@@ -39,6 +49,7 @@
     "calc_4x4determinant",
     "corners",
     "Rectangle",
+    "Ball",
     "Cube",
     "Disk",
     "drectangle",
@@ -64,4 +75,7 @@
     "get_winded_boundary_edges",
     "vertex_in_entity3",
     "unique_edges",
+    "Union",
+    "Intersection",
+    "Difference",
 ]