global variables (#105)

* removing global dictionaries and placing them inside their functions to avoid mutable default arguments.

README.md

@@ -249,7 +249,12 @@ points, cells = generate_mesh(
 )
 
 points, cells = sliver_removal(
-    points=points, domain=cylinder, edge_length=fh, h0=hmin, min_dh_angle_bound=5.0
+    points=points,
+    domain=cylinder,
+    edge_length=fh,
+    h0=hmin,
+    min_dh_angle_bound=5.0,
+    bbox=bbox,
 )
 
 
@@ -425,7 +430,7 @@ and this project (tries to) adhere to [Semantic Versioning](https://semver.org/s
 
 ### Unreleased
 - Added more examples on README
-### Fixed 
+### Fixed
 - Silence messages about pfix when verbose=0
 
 ### [3.0.5] - 2020-10-18

SeismicMesh/generation/mesh_generator.py

@@ -28,20 +28,6 @@
 
 __all__ = ["sliver_removal", "generate_mesh"]
 
-opts = {
-    "verbose": 1,
-    "max_iter": 50,
-    "seed": 0,
-    "perform_checks": False,
-    "pfix": None,
-    "axis": 1,
-    "min_dh_angle_bound": 10.0,
-    "max_dh_angle_bound": 180.0,
-    "points": None,
-    "delta_t": 0.30,
-    "geps_mult": 0.1,
-}
-
 
 def silence(func):
     def wrapper(*args, **kwargs):
@@ -117,10 +103,22 @@ def sliver_removal(points, domain, edge_length, comm=None, **kwargs):  # noqa: C
             warnings.warn("Sliver removal only works in serial for now")
         return True, True
 
-    opts.update(kwargs)
+    sliver_opts = {
+        "verbose": 1,
+        "max_iter": 50,
+        "perform_checks": False,
+        "axis": 1,
+        "min_dh_angle_bound": 10.0,
+        "max_dh_angle_bound": 180.0,
+        "points": None,
+        "delta_t": 0.30,
+        "geps_mult": 0.1,
+    }
+
+    sliver_opts.update(kwargs)
     _parse_kwargs(kwargs)
 
-    verbosity1, verbosity2 = _select_verbosity(opts)
+    verbosity1, verbosity2 = _select_verbosity(sliver_opts)
 
     @verbosity1
     def print_msg1(msg):
@@ -135,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)
+    fd, bbox0 = _unpack_domain(domain, sliver_opts)
 
     fh, bbox1, hmin = _unpack_sizing(edge_length)
 
@@ -156,25 +154,25 @@ def print_msg2(msg):
     if hmin is not None:
         h0 = hmin
     else:
-        h0 = opts["h0"]
+        h0 = sliver_opts["h0"]
 
     # check h0
     if h0 < 0:
         raise ValueError("`h0` must be > 0")
 
     deps = np.sqrt(np.finfo(np.double).eps) * h0
 
-    if opts["max_iter"] < 0:
+    if sliver_opts["max_iter"] < 0:
         raise ValueError("`max_iter` must be > 0")
-    max_iter = opts["max_iter"]
+    max_iter = sliver_opts["max_iter"]
 
     print_msg1(
         "Will attempt " + str(max_iter) + " iterations to bound the dihedral angles..."
     )
 
-    geps = opts["geps_mult"] * h0
-    min_dh_bound = opts["min_dh_angle_bound"] * math.pi / 180
-    max_dh_bound = opts["max_dh_angle_bound"] * math.pi / 180
+    geps = sliver_opts["geps_mult"] * 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
 
     print_msg1(
         "Enforcing a min. dihedral bound of: "
@@ -228,7 +226,7 @@ def print_msg2(msg):
             print_msg1(
                 "FAILURE: Termination...maximum number of iterations reached.",
             )
-            p, t = _termination(p, t, opts, comm, sliver=True)
+            p, t = _termination(p, t, sliver_opts, comm, sliver=True)
             break
 
         print_msg1(
@@ -329,13 +327,23 @@ def generate_mesh(domain, edge_length, comm=None, **kwargs):  # noqa: C901
 
     """
     comm = comm or MPI.COMM_WORLD
-
+    gen_opts = {
+        "verbose": 1,
+        "max_iter": 50,
+        "seed": 0,
+        "perform_checks": False,
+        "pfix": None,
+        "axis": 1,
+        "points": None,
+        "delta_t": 0.30,
+        "geps_mult": 0.1,
+    }
     # check call was correct
-    opts.update(kwargs)
+    gen_opts.update(kwargs)
     _parse_kwargs(kwargs)
 
     # verbosity decorators
-    verbosity1, verbosity2 = _select_verbosity(opts)
+    verbosity1, verbosity2 = _select_verbosity(gen_opts)
 
     @verbosity1
     def print_msg1(msg):
@@ -346,7 +354,7 @@ def print_msg2(msg):
         print(msg, flush=True)
 
     # unpack domain
-    fd, bbox0 = _unpack_domain(domain)
+    fd, bbox0 = _unpack_domain(domain, gen_opts)
 
     fh, bbox1, hmin = _unpack_sizing(edge_length)
 
@@ -376,29 +384,31 @@ def print_msg2(msg):
     if hmin is not None:
         h0 = hmin
     else:
-        h0 = opts["h0"]
+        h0 = gen_opts["h0"]
 
     # check h0
     if h0 < 0:
         raise ValueError("`h0` must be > 0")
 
-    if opts["max_iter"] < 0:
+    if gen_opts["max_iter"] < 0:
         raise ValueError("`max_iter` must be > 0")
-    max_iter = opts["max_iter"]
+    max_iter = gen_opts["max_iter"]
 
     # these parameters originate from the original DistMesh
     L0mult = 1 + 0.4 / 2 ** (dim - 1)
-    delta_t = opts["delta_t"]
-    geps = opts["geps_mult"] * h0
+    delta_t = gen_opts["delta_t"]
+    geps = gen_opts["geps_mult"] * h0
     deps = np.sqrt(np.finfo(np.double).eps) * h0
 
     DT = _select_cgal_dim(dim)
 
-    pfix, nfix = _unpack_pfix(dim, opts, comm)
+    pfix, nfix = _unpack_pfix(dim, gen_opts, comm)
     if comm.rank == 0:
         print_msg1("Constraining " + str(nfix) + " fixed points..")
 
-    fh, p, extents = _initialize_points(dim, geps, bbox, fh, fd, h0, opts, pfix, comm)
+    fh, p, extents = _initialize_points(
+        dim, geps, bbox, fh, fd, h0, gen_opts, pfix, comm
+    )
 
     N = p.shape[0]
 
@@ -440,7 +450,7 @@ def print_msg2(msg):
                 print_msg1(
                     "Termination reached...maximum number of iterations reached.",
                 )
-            p, t = _termination(p, t, opts, comm)
+            p, t = _termination(p, t, gen_opts, comm)
             if comm.rank == 0:
                 p = _improve_level_set_newton(p, t, fd, deps, deps * 1000)
             break
@@ -534,7 +544,7 @@ def func(x):
     return fh, bbox, hmin
 
 
-def _unpack_domain(domain):
+def _unpack_domain(domain, opts):
     if isinstance(domain, geometry.Rectangle):
         bbox = (domain.x1, domain.x2, domain.y1, domain.y2)
         fd = domain.eval
@@ -746,15 +756,14 @@ def _generate_initial_points(h0, geps, dim, bbox, fh, fd, pfix, comm, opts):
             p[np.random.rand(p.shape[0]) < r0m ** dim / r0 ** dim],
         )
     )
-    extents = _form_extents(p, h0, comm)
+    extents = _form_extents(p, h0, comm, opts)
     return fh, p, extents
 
 
 def _initialize_points(dim, geps, bbox, fh, fd, h0, opts, pfix, comm):
     """Form initial point set to mesh with"""
     points = opts["points"]
     if points is None:
-        # def _generate_initial_points(h0, geps, dim, bbox, fh, fd, pfix, comm, opts):
         fh, p, extents = _generate_initial_points(
             h0, geps, dim, bbox, fh, fd, pfix, comm, opts
         )
@@ -763,7 +772,7 @@ def _initialize_points(dim, geps, bbox, fh, fd, h0, opts, pfix, comm):
     return fh, p, extents
 
 
-def _form_extents(p, h0, comm):
+def _form_extents(p, h0, comm, opts):
     dim = p.shape[1]
     _axis = opts["axis"]
     if comm.size > 1:

tests/test_3dmesher_SDF.py

@@ -39,7 +39,7 @@ def EF(p):
     )
 
     points, cells = sliver_removal(
-        points=points, domain=cylinder, edge_length=EF, h0=hmin
+        points=points, domain=cylinder, edge_length=EF, h0=hmin, bbox=bbox
     )
 
     assert allclose(sum(geometry.simp_vol(points, cells)), 6.28, atol=hmin)