From d48ca7c3ecbe6946d0ba3b3e0d02a361f5660d1a Mon Sep 17 00:00:00 2001
From: Keith Roberts <krober@usp.br>
Date: Tue, 24 Nov 2020 16:29:33 -0300
Subject: [PATCH] dampening in Newton optmi. (#157)

* dampening in Newton optmi.
* boundary projection occurs after `sliver_removal` too
---
 README.md                                | 6 +++++-
 SeismicMesh/generation/mesh_generator.py | 8 ++++++--
 setup.cfg                                | 2 +-
 3 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index 5e90017b..01478ed9 100644
--- a/README.md
+++ b/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.
diff --git a/SeismicMesh/generation/mesh_generator.py b/SeismicMesh/generation/mesh_generator.py
index f1c58e65..84ce8cbb 100644
--- a/SeismicMesh/generation/mesh_generator.py
+++ b/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
 
 
diff --git a/setup.cfg b/setup.cfg
index 2e9d7b93..0c0a51ed 100644
--- a/setup.cfg
+++ b/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