Improve optimization

- handle concave quads (TODO: refactor the function)
- use under-relaxation (TODO: rewrite to use actual number instead of
  bool)
- use RMS value when summing up quality of junctions
- handle crashes when calculating sensitivity outside bounds (TODO:
  handle bounds within the calculation)
- use a different transform function for angles (does not go to infinity
  and is non-linear below some threshold)

Author: FranzBangar

docs/optimization/quality-criteria.ods

@@ -5,6 +5,8 @@
 import classy_blocks as cb
 from classy_blocks.util import functions as f
 
+# TODO! Rewrite this to use a MappedSketch and get rid of 3D optimization of a 2D domain
+
 # This rather lengthy tutorial does the following:
 # - reads 2D airfoil points*
 # - rotates the airfoil to a desired angle of attack

examples/complex/airfoil/airfoil.py

@@ -80,7 +80,7 @@ def add_regions(regions: List[Region]) -> None:
     for clamp in region.get_clamps(mesh):
         optimizer.add_clamp(clamp)
 
-optimizer.optimize(tolerance=1e-3, method="SLSQP")
+optimizer.optimize(tolerance=1e-3, relax=True, max_iterations=10)
 
 # Now Block objects contain optimization results but those are not reflected in
 # user-created Operations. Mesh.backport() will copy the data back.

examples/complex/cyclone/cyclone.py

@@ -127,7 +127,7 @@ def quality(self) -> float:
         """Returns summed qualities of all junctions"""
         # It is only called when optimizing linked clamps
         # or at the end of an iteration.
-        return sum(junction.quality for junction in self.junctions)
+        return sum(junction.quality**2 for junction in self.junctions) ** 0.5
 
     def update(self, index: int, position: NPPointType) -> float:
         self.points[index] = position

src/classy_blocks/optimize/grid.py

@@ -62,4 +62,7 @@ def quality(self) -> float:
         else:
             quality_function = get_quad_quality
 
-        return sum(quality_function(self.points, np.array(cell.indexes, dtype=np.int32)) for cell in self.cells)
+        return (
+            sum(quality_function(self.points, np.array(cell.indexes, dtype=np.int32)) ** 2 for cell in self.cells)
+            ** 0.5
+        )

src/classy_blocks/optimize/junction.py

@@ -34,7 +34,7 @@ def add_clamp(self, clamp: ClampBase) -> None:
     def add_link(self, link: LinkBase) -> None:
         self.grid.add_link(link)
 
-    def optimize_clamp(self, clamp: ClampBase, method: MinimizationMethodType) -> None:
+    def optimize_clamp(self, clamp: ClampBase, method: MinimizationMethodType, relax: bool, iteration_no: int) -> None:
         """Move clamp.vertex so that quality at junction is improved;
         rollback changes if grid quality decreased after optimization"""
         initial_params = copy.copy(clamp.params)
@@ -46,11 +46,23 @@ def optimize_clamp(self, clamp: ClampBase, method: MinimizationMethodType) -> No
         def fquality(params):
             # move all vertices according to X
             clamp.update_params(params)
-            return self.grid.update(junction.index, clamp.position)
+            quality = self.grid.update(junction.index, clamp.position)
+            return quality
 
         try:
             scipy.optimize.minimize(fquality, clamp.params, bounds=clamp.bounds, method=method, tol=1e-2)
 
+            if relax:
+                relaxation_factor = 1 - 2 ** -(iteration_no + 1)
+                if relaxation_factor > 0.9:
+                    relaxation_factor = 1
+
+                for i, param in enumerate(clamp.params):
+                    clamp.params[i] = initial_params[i] + relaxation_factor * (param - initial_params[i])
+
+                # update the position after relaxation
+                fquality(clamp.params)
+
             reporter.junction_final = junction.quality
             reporter.grid_final = self.grid.quality
 
@@ -78,36 +90,45 @@ def fquality(clamp, junction, params):
             self.grid.update(junction.index, clamp.position)
             return junction.quality
 
-        sensitivities = np.asarray(
-            scipy.optimize.approx_fprime(clamp.params, lambda p: fquality(clamp, junction, p), epsilon=10 * TOL)
-        )
+        try:
+            sensitivities = np.asarray(
+                scipy.optimize.approx_fprime(clamp.params, lambda p: fquality(clamp, junction, p), epsilon=10 * TOL)
+            )
+        except ValueError:
+            clamp.update_params(initial_params)
+            self.grid.update(junction.index, clamp.position)
+            return np.ones(len(clamp.params))
 
         clamp.update_params(initial_params)
         self.grid.update(junction.index, clamp.position)
 
         return np.linalg.norm(sensitivities)
 
-    def optimize_iteration(self, method: MinimizationMethodType) -> None:
+    def optimize_iteration(self, method: MinimizationMethodType, relax: bool, iteration_no: int) -> None:
         clamps = sorted(self.grid.clamps, key=lambda c: self._get_sensitivity(c), reverse=True)
 
         for clamp in clamps:
-            self.optimize_clamp(clamp, method)
+            self.optimize_clamp(clamp, method, relax, iteration_no)
 
     def optimize(
-        self, max_iterations: int = 20, tolerance: float = 0.1, method: MinimizationMethodType = "SLSQP"
+        self,
+        max_iterations: int = 20,
+        tolerance: float = 0.1,
+        method: MinimizationMethodType = "SLSQP",
+        relax: bool = False,
     ) -> IterationDriver:
         """Move vertices, defined and restrained with Clamps
         so that better mesh quality is obtained.
 
         Within each iteration, all vertices will be moved, starting with the one with the most influence on quality.
-        Lower tolerance values"""
+        Lower tolerance values."""
         driver = IterationDriver(max_iterations, tolerance)
 
         start_time = time.time()
 
         while not driver.converged:
             driver.begin_iteration(self.grid.quality)
-            self.optimize_iteration(method)
+            self.optimize_iteration(method, relax, len(driver.iterations))
             driver.end_iteration(self.grid.quality)
 
         end_time = time.time()

src/classy_blocks/optimize/optimizer.py

@@ -17,17 +17,17 @@ def scale_quality(base: float, exponent: float, factor: float, value: float) ->
 
 @numba.jit(nopython=True, cache=True)
 def scale_non_ortho(angle: float) -> float:
-    return scale_quality(1.5, 0.25, 0.5, angle)
+    return scale_quality(1.4, 0.25, 0.5, angle)
 
 
 @numba.jit(nopython=True, cache=True)
 def scale_inner_angle(angle: float) -> float:
-    return scale_quality(1.5, 0.25, 0.5, np.abs(angle))
+    return scale_quality(1.4, 0.25, 0.5, np.abs(angle))
 
 
 @numba.jit(nopython=True, cache=True)
 def scale_aspect(ratio: float) -> float:
-    return scale_quality(3, 2.5, 3, np.log10(ratio))
+    return scale_quality(4, 3, 3, np.log10(ratio))
 
 
 @numba.jit(nopython=True, cache=True)
@@ -73,6 +73,41 @@ def get_quad_normal(points: NPPointListType) -> Tuple[NPVectorType, NPVectorType
     return center, normal / np.linalg.norm(normal), max_length / (min_length + VSMALL)
 
 
+@numba.jit(nopython=True, cache=True)
+def is_quad_convex(points: NPPointListType) -> bool:
+    # Compute normal using the first three points
+    normal = np.cross(points[1] - points[0], points[2] - points[1])
+    normal /= np.linalg.norm(normal)
+
+    sign = 0
+    for i in range(4):
+        prev_leg = points[i] - points[(i - 1) % 4]
+        next_leg = points[(i + 1) % 4] - points[i]
+        cross = np.cross(prev_leg, next_leg)
+        dot = np.dot(cross, normal)
+        if i == 0:
+            sign = np.sign(dot)
+            if sign == 0:
+                return False  # Degenerate
+        else:
+            if np.sign(dot) != sign:
+                return False
+    return True
+
+
+@numba.jit(nopython=True, cache=True)
+def scale_angle(angle: float) -> float:
+    n = 4
+    m = 10
+    threshold = 65
+    a = m / (n * threshold ** (n - 1))
+
+    if angle <= threshold:
+        return a * angle**n
+
+    return a * threshold**n + m * (angle - threshold)
+
+
 @numba.jit(nopython=True, cache=True)
 def get_quad_non_ortho_quality(
     quad_points: NPPointListType, quad_center: NPPointType, quad_normal: NPPointType
@@ -92,7 +127,9 @@ def get_quad_non_ortho_quality(
         center_vector /= np.linalg.norm(center_vector)
 
         angle = 180 * np.arccos(np.dot(side_normal, center_vector)) / np.pi
-        quality += scale_non_ortho(angle)
+        if not is_quad_convex(quad_points):
+            angle += 180
+        quality += scale_angle(angle)
 
     return quality
 
@@ -113,7 +150,10 @@ def get_quad_angle_quality(quad_points: NPPointListType) -> float:
         side_2 /= np.linalg.norm(side_2) + VSMALL
 
         angle = 180 * np.arccos(np.dot(side_1, side_2)) / np.pi - 90
-        quality += scale_inner_angle(angle)
+        if not is_quad_convex(quad_points):
+            angle += 180
+
+        quality += scale_angle(angle)
 
     return quality
 
@@ -129,6 +169,7 @@ def get_quad_quality(grid_points: NPPointListType, cell_indexes: NPIndexType) ->
     # inner angles
     quality += get_quad_angle_quality(cell_points)
 
+    # aspect ratio
     quality += scale_aspect(cell_aspect)
 
     return quality
@@ -141,8 +182,6 @@ def get_hex_quality(grid_points: NPPointListType, cell_indexes: NPIndexType) ->
 
     side_indexes = np.array([[0, 1, 2, 3], [7, 6, 5, 4], [4, 0, 3, 7], [6, 2, 1, 5], [0, 4, 5, 1], [7, 3, 2, 6]])
 
-    max_aspect = 1
-
     quality = 0
 
     for side in side_indexes:
@@ -158,14 +197,14 @@ def get_hex_quality(grid_points: NPPointListType, cell_indexes: NPIndexType) ->
 
         center_vector /= np.linalg.norm(center_vector)
 
-        angle = 180 * np.arccos(np.dot(side_normal, center_vector)) / np.pi
-        quality += scale_non_ortho(angle)
-
-        max_aspect = max(max_aspect, side_aspect)
+        angle = 180 * np.arccos(min(1 - VSMALL, np.dot(side_normal, center_vector))) / np.pi
+        if not is_quad_convex(side_points):
+            angle = 180 - angle
+        quality += scale_angle(angle)
 
-        # take inner angles and aspect from quad calculation
+        # take inner angles and aspect from quad calculation;
         quality += get_quad_angle_quality(side_points)
 
-    quality += scale_aspect(max_aspect)
+        quality += scale_aspect(side_aspect)
 
     return quality