dampening in Newton optmi. (#157)

* dampening in Newton optmi.
* boundary projection occurs after `sliver_removal` too

Author: Keith Roberts

README.md

@@ -511,9 +511,13 @@ 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
+## Unreleased 
+- None
+
+## [3.1.5] - 2020-11-24
 - Support for constraining/immersing subdomains represented as signed distance functions.
 - Faster cell manipulation operations for ~5-10% better speedups in parallel. 
+- Projection of points back onto level set.
 
 ## [3.1.4] - 2020-11-15
 - Laplacian smoothing at termination for 2D meshing...significantly improves minimum cell quality.

SeismicMesh/generation/mesh_generator.py

@@ -170,6 +170,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
 
@@ -247,6 +248,7 @@ def print_msg2(msg):
                 + " iterations...no slivers detected!",
             )
             p, t, _ = geometry.fix_mesh(p, t, dim=dim, delete_unused=True)
+            p = _improve_level_set_newton(p, t, fd, deps, deps * 1000)
             return p, t
 
         p0, p1, p2, p3 = (
@@ -715,7 +717,8 @@ def _improve_level_set_newton(p, t, fd, deps, tol):
     """Reduce level set error by using Newton's minimization method"""
     dim = p.shape[1]
     bid = geometry.get_boundary_vertices(t, dim)
-    for _ in range(5):
+    alpha = 1
+    for iteration in range(5):
         d = fd(p[bid])
 
         def _deps_vec(i):
@@ -726,7 +729,8 @@ def _deps_vec(i):
         dgrads = [(fd(p[bid] + _deps_vec(i)) - d) / deps for i in range(dim)]
         dgrad2 = sum(dgrad ** 2 for dgrad in dgrads)
         dgrad2 = np.where(dgrad2 < deps, deps, dgrad2)
-        p[bid] -= (d * np.vstack(dgrads) / dgrad2).T  # Project
+        p[bid] -= alpha * (d * np.vstack(dgrads) / dgrad2).T  # Project
+        alpha /= iteration + 1
     return p
 
 

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = SeismicMesh
-version = 3.1.4
+version = 3.1.5
 url = https://github.com/krober10nd/SeismicMesh
 author = Keith Roberts
 email = keithrbt0@gmail.com