diff --git a/README.md b/README.md
index e38ffc60..22630ef7 100644
--- a/README.md
+++ b/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.
diff --git a/SeismicMesh/__init__.py b/SeismicMesh/__init__.py
index 5c35a207..5067a8c4 100644
--- a/SeismicMesh/__init__.py
+++ b/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",
diff --git a/SeismicMesh/generation/mesh_generator.py b/SeismicMesh/generation/mesh_generator.py
index a8629aec..7bdd6040 100644
--- a/SeismicMesh/generation/mesh_generator.py
+++ b/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,15 +548,21 @@ 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"]
@@ -561,7 +570,7 @@ def _unpack_domain(domain, opts):
     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
diff --git a/SeismicMesh/geometry/__init__.py b/SeismicMesh/geometry/__init__.py
index d0e89f4e..2afee92b 100644
--- a/SeismicMesh/geometry/__init__.py
+++ b/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",
 ]
diff --git a/SeismicMesh/geometry/signed_distance_functions.py b/SeismicMesh/geometry/signed_distance_functions.py
index 15f29549..1ae7d3af 100644
--- a/SeismicMesh/geometry/signed_distance_functions.py
+++ b/SeismicMesh/geometry/signed_distance_functions.py
@@ -1,42 +1,161 @@
 import numpy as np
+import itertools
+
+
+def corners(bbox):
+    """Get the corners of a box in N-dim"""
+    mins = bbox[::2]
+    maxs = bbox[1::2]
+    return list(itertools.product(*zip(mins, maxs)))
+
+
+def _gather_corners(domains):
+    corners = [d.corners for d in domains if d.corners is not None]
+    corners = np.concatenate(corners)
+    if len(corners) == 0:
+        corners = None
+    return corners
+
+
+class Union:
+    def __init__(self, domains):
+        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]
+        if self.dim == 2:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+            )
+        elif self.dim == 3:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+                min(d.bbox[4] for d in domains),
+                max(d.bbox[5] for d in domains),
+            )
+        self.corners = _gather_corners(domains)
+        self.domains = domains
+
+    def eval(self, x):
+        return np.minimum.reduce([d.eval(x) for d in self.domains])
+
+
+class Intersection:
+    def __init__(self, domains):
+        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]
+        if self.dim == 2:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+            )
+        elif self.dim == 3:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+                min(d.bbox[4] for d in domains),
+                max(d.bbox[5] for d in domains),
+            )
+        self.corners = _gather_corners(domains)
+        self.domains = domains
+
+    def eval(self, x):
+        return np.maximum.reduce([d.eval(x) for d in self.domains])
+
+
+class Difference:
+    def __init__(self, domains):
+        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]
+        if self.dim == 2:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+            )
+        elif self.dim == 3:
+            self.bbox = (
+                min(d.bbox[0] for d in domains),
+                max(d.bbox[1] for d in domains),
+                min(d.bbox[2] for d in domains),
+                max(d.bbox[3] for d in domains),
+                min(d.bbox[4] for d in domains),
+                max(d.bbox[5] for d in domains),
+            )
+        self.corners = _gather_corners(domains)
+        self.domains = domains
+
+    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)]
+        )
 
 
 class Disk:
     def __init__(self, x0, r):
+        self.dim = 2
+        self.corners = None
         self.xc = x0[0]
         self.yc = x0[1]
         self.r = r
-        self.x1 = x0[0] - r
-        self.x2 = x0[0] + r
-        self.y1 = x0[1] - r
-        self.y2 = x0[1] + r
+        self.bbox = (x0[0] - r, x0[0] + r, x0[1] - r, x0[1] + r)
 
     def eval(self, x):
         return _ddisk(x, self.xc, self.yc, self.r)
 
 
+class Ball:
+    def __init__(self, x0, r):
+        self.dim = 3
+        self.corners = None
+        self.xc = x0[0]
+        self.yc = x0[1]
+        self.zc = x0[2]
+        self.r = r
+        self.bbox = (x0[0] - r, x0[0] + r, x0[1] - r, x0[1] + r, x0[2] - r, x0[2] + r)
+
+    def eval(self, x):
+        return dball(x, self.xc, self.yc, self.zc, self.r)
+
+
 class Rectangle:
     def __init__(self, bbox):
-        self.x1 = bbox[0]
-        self.x2 = bbox[1]
-        self.y1 = bbox[2]
-        self.y2 = bbox[3]
+        self.dim = 2
+        self.corners = corners(bbox)
+        self.bbox = bbox
 
     def eval(self, x):
-        return drectangle(x, self.x1, self.x2, self.y1, self.y2)
+        return drectangle(x, self.bbox[0], self.bbox[1], self.bbox[2], self.bbox[3])
 
 
 class Cube:
     def __init__(self, bbox):
-        self.x1 = bbox[0]
-        self.x2 = bbox[1]
-        self.y1 = bbox[2]
-        self.y2 = bbox[3]
-        self.z1 = bbox[4]
-        self.z2 = bbox[5]
+        self.dim = 3
+        self.corners = corners(bbox)
+        self.bbox = bbox
 
     def eval(self, x):
-        return dblock(x, self.x1, self.x2, self.y1, self.y2, self.z1, self.z2)
+        return dblock(
+            x,
+            self.bbox[0],
+            self.bbox[1],
+            self.bbox[2],
+            self.bbox[3],
+            self.bbox[4],
+            self.bbox[5],
+        )
 
 
 def _ddisk(p, xc, yc, r):
@@ -44,6 +163,11 @@ def _ddisk(p, xc, yc, r):
     return np.sqrt(((p - np.array([xc, yc])) ** 2).sum(-1)) - r
 
 
+def dball(p, xc, yc, zc, r):
+    """Signed distance function for a ball centered at xc,yc,zc with radius  r."""
+    return np.sqrt((p[:, 0] - xc) ** 2 + (p[:, 1] - yc) ** 2 + (p[:, 2] - zc) ** 2) - r
+
+
 def drectangle(p, x1, x2, y1, y2):
     min = np.minimum
     """Signed distance function for rectangle with corners (x1,y1), (x2,y1),
@@ -96,21 +220,3 @@ def dblock(p, x1, x2, y1, y2, z1, z2):
         ),
         x2 - p[:, 0],
     )
-
-
-def dintersect(d1, d2):
-    """Signed distance to set intersection of two regions described by signed
-    distance functions d1 and d2.
-    Not exact the true signed distance function for the difference,
-    for example around corners.
-    """
-    return np.maximum(d1, d2)
-
-
-def ddiff(d1, d2):
-    """Signed distance to set difference between two regions described by
-    signed distance functions d1 and d2.
-    Not exact the true signed distance function for the difference,
-    for example around corners.
-    """
-    return np.maximum(d1, -d2)
diff --git a/SeismicMesh/geometry/utils.py b/SeismicMesh/geometry/utils.py
index 5238ca84..dfad8a4a 100644
--- a/SeismicMesh/geometry/utils.py
+++ b/SeismicMesh/geometry/utils.py
@@ -1,4 +1,3 @@
-import itertools
 import numpy as np
 import scipy.sparse as spsparse
 
@@ -10,13 +9,6 @@
 dete = gutils.calc_4x4determinant
 
 
-def corners(bbox):
-    """Get the corners of a box in N-dim"""
-    mins = bbox[::2]
-    maxs = bbox[1::2]
-    return list(itertools.product(*zip(mins, maxs)))
-
-
 def calc_re_ratios(vertices, entities, dim=2):
     """Calculate radius edge ratios--mesh quality metric
 
diff --git a/tests/test_2dmesher.py b/tests/test_2dmesher.py
index c5bbd28d..fcbc1e1c 100644
--- a/tests/test_2dmesher.py
+++ b/tests/test_2dmesher.py
@@ -1,11 +1,11 @@
 import os
 
-from numpy import allclose
 import pytest
+from numpy import allclose
 
 from SeismicMesh import (
-    generate_mesh,
     Rectangle,
+    generate_mesh,
     get_sizing_function_from_segy,
     plot_sizing_function,
     write_velocity_model,
diff --git a/tests/test_2dmesher_SDF.py b/tests/test_2dmesher_SDF.py
index 58397646..8e2f9a51 100644
--- a/tests/test_2dmesher_SDF.py
+++ b/tests/test_2dmesher_SDF.py
@@ -1,5 +1,5 @@
-from numpy import allclose, sum
 import pytest
+from numpy import allclose, sum
 
 from SeismicMesh import generate_mesh, geometry
 
diff --git a/tests/test_2dmesher_domain_extension.py b/tests/test_2dmesher_domain_extension.py
index b5ce9cc2..a00af067 100644
--- a/tests/test_2dmesher_domain_extension.py
+++ b/tests/test_2dmesher_domain_extension.py
@@ -1,11 +1,11 @@
 import os
 
-from numpy import allclose
 import pytest
+from numpy import allclose
 
 from SeismicMesh import (
-    generate_mesh,
     Rectangle,
+    generate_mesh,
     get_sizing_function_from_segy,
     plot_sizing_function,
     write_velocity_model,
diff --git a/tests/test_2dmesher_par.py b/tests/test_2dmesher_par.py
index 2ddfb927..4d96a19d 100644
--- a/tests/test_2dmesher_par.py
+++ b/tests/test_2dmesher_par.py
@@ -5,10 +5,10 @@
 from mpi4py import MPI
 
 from SeismicMesh import (
-    get_sizing_function_from_segy,
-    generate_mesh,
     Rectangle,
+    generate_mesh,
     geometry,
+    get_sizing_function_from_segy,
 )
 
 comm = MPI.COMM_WORLD
diff --git a/tests/test_2dmesher_par_adapt.py b/tests/test_2dmesher_par_adapt.py
index f7700635..24fb4475 100644
--- a/tests/test_2dmesher_par_adapt.py
+++ b/tests/test_2dmesher_par_adapt.py
@@ -5,10 +5,10 @@
 from mpi4py import MPI
 
 from SeismicMesh import (
-    get_sizing_function_from_segy,
+    Rectangle,
     generate_mesh,
     geometry,
-    Rectangle,
+    get_sizing_function_from_segy,
 )
 
 comm = MPI.COMM_WORLD
diff --git a/tests/test_3dmesher.py b/tests/test_3dmesher.py
index 6d852c45..6b28e623 100644
--- a/tests/test_3dmesher.py
+++ b/tests/test_3dmesher.py
@@ -1,13 +1,13 @@
 import os
 
-from numpy import allclose
 import pytest
+from numpy import allclose
 
 from SeismicMesh import (
-    generate_mesh,
-    sliver_removal,
     Cube,
+    generate_mesh,
     get_sizing_function_from_segy,
+    sliver_removal,
     write_velocity_model,
 )
 
diff --git a/tests/test_3dmesher_SDF.py b/tests/test_3dmesher_SDF.py
index 8e7a6b7d..a1bfd69a 100644
--- a/tests/test_3dmesher_SDF.py
+++ b/tests/test_3dmesher_SDF.py
@@ -1,7 +1,7 @@
-from numpy import allclose, sum, array, maximum, sqrt
 import pytest
+from numpy import allclose, array, maximum, sqrt, sum
 
-from SeismicMesh import generate_mesh, sliver_removal, geometry
+from SeismicMesh import generate_mesh, geometry, sliver_removal
 
 
 @pytest.mark.serial
diff --git a/tests/test_3dmesher_domain_extension.py b/tests/test_3dmesher_domain_extension.py
index bb70a9c8..2052b2a5 100644
--- a/tests/test_3dmesher_domain_extension.py
+++ b/tests/test_3dmesher_domain_extension.py
@@ -1,13 +1,13 @@
 import os
 
-from numpy import allclose
 import pytest
+from numpy import allclose
 
 from SeismicMesh import (
-    generate_mesh,
-    sliver_removal,
     Cube,
+    generate_mesh,
     get_sizing_function_from_segy,
+    sliver_removal,
     write_velocity_model,
 )
 
diff --git a/tests/test_3dmesher_par.py b/tests/test_3dmesher_par.py
index ee1f51fe..02ce0bc8 100644
--- a/tests/test_3dmesher_par.py
+++ b/tests/test_3dmesher_par.py
@@ -5,11 +5,11 @@
 from mpi4py import MPI
 
 from SeismicMesh import (
-    get_sizing_function_from_segy,
+    Cube,
     generate_mesh,
-    sliver_removal,
     geometry,
-    Cube,
+    get_sizing_function_from_segy,
+    sliver_removal,
 )
 
 comm = MPI.COMM_WORLD
diff --git a/tests/test_3dmesher_par_adapt.py b/tests/test_3dmesher_par_adapt.py
index 9ac504b6..e44dbb06 100644
--- a/tests/test_3dmesher_par_adapt.py
+++ b/tests/test_3dmesher_par_adapt.py
@@ -4,7 +4,7 @@
 import pytest
 from mpi4py import MPI
 
-from SeismicMesh import get_sizing_function_from_segy, generate_mesh, Cube, geometry
+from SeismicMesh import Cube, generate_mesh, geometry, get_sizing_function_from_segy
 
 comm = MPI.COMM_WORLD
 
diff --git a/tests/test_pfix.py b/tests/test_pfix.py
index c2e80618..314151aa 100644
--- a/tests/test_pfix.py
+++ b/tests/test_pfix.py
@@ -1,5 +1,6 @@
-import pytest
 import numpy as np
+import pytest
+
 import SeismicMesh