Improve optimization speed and add choices for solver

Quality is now being calculated on a static numpy array and is cached
inside Cell objects.
Optimization algorithm can now be selected for optimization and
reporting is improved (timing and overall improvement). Different
geometric situations work differently with each solver so best is to try
and see.

Timing and quality improvement is output for easier comparison between
choices.

Examples that demonstrate optimization have been corrected and slightly
improved.

examples/complex/airfoil/airfoil.py

@@ -166,7 +166,7 @@ def make_link(leader):
 
 
 # Points that slide along airfoil curve
-for index in (10, 11, 13, 14):
+for index in (10, 11, 12, 13, 14):
     opt_vertex = find_vertex(index)
     clamp = cb.CurveClamp(opt_vertex, foil_curve)
     clamp.add_link(make_link(opt_vertex))
@@ -195,7 +195,7 @@ def optimize_along_line(point_index, line_index_1, line_index_2):
 optimize_along_line(5, 3, 6)
 
 if OPTIMIZE:
-    optimizer.optimize(tolerance=1e-5, method="SLSQP")
+    optimizer.optimize(tolerance=1e-7, method="SLSQP")
 
 ### Write the mesh
 mesh.modify_patch("topAndBottom", "empty")

examples/shape/custom.py

@@ -63,5 +63,5 @@ def add_edges(self) -> None:
         op.chop(i, count=10)
 
 mesh.add(shape)
-
+mesh.set_default_patch("walls", "wall")
 mesh.write(os.path.join("..", "case", "system", "blockMeshDict"), debug_path="debug.vtk")

src/classy_blocks/mesh.py

@@ -95,7 +95,7 @@ def delete(self, operation: Operation) -> None:
         the data remains but it will not contribute to the mesh"""
         self.deleted.add(operation)
 
-    def assemble(self, skip_edges: bool = True) -> None:
+    def assemble(self, skip_edges: bool = False) -> None:
         """Converts classy_blocks entities (operations and shapes) to
         actual vertices, edges, blocks and other stuff to be inserted into
         blockMeshDict. After this has been done, the above objects
@@ -162,7 +162,7 @@ def backport(self) -> None:
             op.top_face.update(vertices[4:])
 
         self.clear()
-        self.assemble(True)
+        self.assemble()
 
     def format_settings(self) -> str:
         """Put self.settings in a proper, blockMesh-readable format"""
@@ -181,7 +181,7 @@ def write(self, output_path: str, debug_path: Optional[str] = None) -> None:
         a VTK file is created first where each block is a single cell, to see simplified
         blocking in case blockMesh fails with an unfriendly error message."""
         if not self.is_assembled:
-            self.assemble(False)
+            self.assemble()
 
         if debug_path is not None:
             write_vtk(debug_path, self.vertex_list.vertices, self.block_list.blocks)

src/classy_blocks/modify/grid.py

@@ -1,10 +1,17 @@
+from typing import List
+
 import numpy as np
 
 from classy_blocks.mesh import Mesh
 from classy_blocks.modify.cell import Cell
+from classy_blocks.modify.clamps.clamp import ClampBase
 from classy_blocks.modify.junction import Junction
 
 
+class NoJunctionError(Exception):
+    """Raised when there's a clamp defined for a vertex that doesn't exist"""
+
+
 class Grid:
     """A list of cells and junctions"""
 
@@ -34,13 +41,45 @@ def _bind_neighbours(self) -> None:
             for cell_2 in self.cells:
                 cell_1.add_neighbour(cell_2)
 
+    def get_junction_from_clamp(self, clamp: ClampBase) -> Junction:
+        for junction in self.junctions:
+            if junction.clamp == clamp:
+                return junction
+
+        raise NoJunctionError
+
     def update(self, junction: Junction) -> None:
         self.points[junction.vertex.index] = junction.vertex.position
-
         for cell in junction.cells:
             cell.invalidate()
 
+        # also update linked stuff
+        if junction.clamp is not None:
+            if junction.clamp.is_linked:
+                linked_junction = self.get_junction_from_clamp(junction.clamp)
+                for cell in linked_junction.cells:
+                    cell.invalidate()
+
+    def clear_cache(self):
+        for cell in self.cells:
+            cell._quality = None
+
+    def add_clamp(self, clamp: ClampBase) -> None:
+        for junction in self.junctions:
+            if junction.vertex == clamp.vertex:
+                junction.add_clamp(clamp)
+
+    @property
+    def clamps(self) -> List[ClampBase]:
+        clamps: List[ClampBase] = []
+
+        for junction in self.junctions:
+            if junction.clamp is not None:
+                clamps.append(junction.clamp)
+
+        return clamps
+
     @property
     def quality(self) -> float:
-        """Returns summed qualities of all cells in this grid"""
-        return sum([cell.quality for cell in self.cells])
+        """Returns summed qualities of all junctions"""
+        return sum([junction.quality for junction in self.junctions])

src/classy_blocks/modify/junction.py

@@ -1,7 +1,12 @@
-from typing import Set
+from typing import Optional, Set
 
 from classy_blocks.items.vertex import Vertex
 from classy_blocks.modify.cell import Cell
+from classy_blocks.modify.clamps.clamp import ClampBase
+
+
+class ClampExistsError(Exception):
+    """Raised when attempting to add a vertex that's already among existing clamps"""
 
 
 class Junction:
@@ -12,6 +17,8 @@ def __init__(self, vertex: Vertex):
         self.vertex = vertex
         self.cells: Set[Cell] = set()
 
+        self.clamp: Optional[ClampBase] = None
+
     def add_cell(self, cell: Cell) -> None:
         """Adds the given cell to the list if it is
         a part of this junction (one common vertex)"""
@@ -20,6 +27,12 @@ def add_cell(self, cell: Cell) -> None:
                 self.cells.add(cell)
                 return
 
+    def add_clamp(self, clamp: ClampBase) -> None:
+        if self.clamp is not None:
+            raise ClampExistsError(f"A clamp has already been defined for vertex {self.vertex}")
+
+        self.clamp = clamp
+
     @property
     def quality(self) -> float:
         """Returns average quality of all cells at this junction;

src/classy_blocks/modify/optimizer.py

@@ -1,30 +1,23 @@
 import copy
 import time
-from typing import List, Literal
+from typing import Literal
 
+import numpy as np
 import scipy.optimize
 
 from classy_blocks.mesh import Mesh
 from classy_blocks.modify.clamps.clamp import ClampBase
 from classy_blocks.modify.grid import Grid
 from classy_blocks.modify.iteration import ClampOptimizationData, IterationDriver
-from classy_blocks.modify.junction import Junction
+from classy_blocks.util.constants import TOL
 
 MinimizationMethodType = Literal["SLSQP", "L-BFGS-B", "Nelder-Mead", "Powell"]
 
 
-class NoJunctionError(Exception):
-    """Raised when there's a clamp defined for a vertex that doesn't exist"""
-
-
 class NoClampError(Exception):
     """Raised when there's no junction defined for a given Clamp"""
 
 
-class ClampExistsError(Exception):
-    """Raised when attempting to add a vertex that's already among existing clamps"""
-
-
 class Optimizer:
     """Provides tools for blocking optimization"""
 
@@ -33,28 +26,16 @@ def __init__(self, mesh: Mesh, report: bool = True):
         self.report = report
 
         self.grid = Grid(mesh)
-        self.clamps: List[ClampBase] = []
 
     def release_vertex(self, clamp: ClampBase) -> None:
         """Adds a clamp to optimization. Raises an exception if it already exists"""
-        for existing in self.clamps:
-            if existing.vertex == clamp.vertex:
-                raise ClampExistsError(f"A clamp has already been defined for vertex {existing}")
-        self.clamps.append(clamp)
+        self.grid.add_clamp(clamp)
 
-    def _get_junction(self, clamp: ClampBase) -> Junction:
-        """Returns a Junction that corresponds to clamp"""
-        for junction in self.grid.junctions:
-            if junction.vertex == clamp.vertex:
-                return junction
-
-        raise NoJunctionError
-
-    def optimize_clamp(self, clamp: ClampBase, iteration: int, method: MinimizationMethodType) -> None:
+    def optimize_clamp(self, clamp: ClampBase, method: MinimizationMethodType) -> 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)
-        junction = self._get_junction(clamp)
+        junction = self.grid.get_junction_from_clamp(clamp)
 
         reporter = ClampOptimizationData(clamp.vertex.index, self.grid.quality, junction.quality)
         reporter.report_start()
@@ -80,23 +61,26 @@ def fquality(params):
 
     def _get_sensitivity(self, clamp):
         """Returns maximum partial derivative at current params"""
-        junction = self._get_junction(clamp)
+        junction = self.grid.get_junction_from_clamp(clamp)
 
         def fquality(clamp, junction, params):
             clamp.update_params(params)
+            self.grid.update(junction)
             return junction.quality
 
-        # sensitivities = np.asarray(
-        #    scipy.optimize.approx_fprime(clamp.params, lambda p: fquality(clamp, junction, p), epsilon=10 * TOL)
-        # )
+        sensitivities = np.asarray(
+            scipy.optimize.approx_fprime(clamp.params, lambda p: fquality(clamp, junction, p), epsilon=10 * TOL)
+        )
+        return np.linalg.norm(sensitivities)
         # return np.max(np.abs(sensitivities.flatten()))
-        return junction.quality
 
-    def optimize_iteration(self, iteration: int, method: MinimizationMethodType) -> None:
-        self.clamps.sort(key=lambda c: self._get_sensitivity(c), reverse=True)
+    def optimize_iteration(self, method: MinimizationMethodType) -> None:
+        self.grid.clear_cache()
+
+        clamps = sorted(self.grid.clamps, key=lambda c: self._get_sensitivity(c), reverse=True)
 
-        for clamp in self.clamps:
-            self.optimize_clamp(clamp, iteration, method)
+        for clamp in clamps:
+            self.optimize_clamp(clamp, method)
 
     def optimize(
         self, max_iterations: int = 20, tolerance: float = 0.1, method: MinimizationMethodType = "SLSQP"
@@ -112,7 +96,7 @@ def optimize(
 
         while not driver.converged:
             driver.begin_iteration(self.grid.quality)
-            self.optimize_iteration(len(driver.iterations), method)
+            self.optimize_iteration(method)
             driver.end_iteration(self.grid.quality)
 
         end_time = time.time()
@@ -124,6 +108,7 @@ def optimize(
             rel_improvement = abs_improvement / start_quality
 
             print(
-                f"Overall improvement: {start_quality:.3e} > {end_quality:.3e} ({abs_improvement:.3e}, {rel_improvement*100:.0f}%)"
+                f"Overall improvement: {start_quality:.3e} > {end_quality:.3e}"
+                f"({abs_improvement:.3e}, {rel_improvement*100:.0f}%)"
             )
             print(f"Elapsed time: {end_time - start_time:.0f}s")

tests/test_modify/test_cell.py

@@ -1,10 +1,18 @@
+import numpy as np
 from parameterized import parameterized
 
 from classy_blocks.modify.cell import Cell, NoCommonSidesError
-from tests.fixtures.block import BlockTestCase
+from tests.fixtures.mesh import MeshTestCase
 
 
-class CellTests(BlockTestCase):
+class CellTests(MeshTestCase):
+    def setUp(self):
+        super().setUp()
+
+    @property
+    def mesh_points(self):
+        return np.array([vertex.position for vertex in self.mesh.vertices])
+
     @parameterized.expand(
         (
             (0, 1, 4),
@@ -14,43 +22,43 @@ class CellTests(BlockTestCase):
         )
     )
     def test_common_vertices(self, index_1, index_2, count):
-        block_1 = self.make_block(index_1)
-        block_2 = self.make_block(index_2)
+        block_1 = self.mesh.blocks[index_1]
+        block_2 = self.mesh.blocks[index_2]
 
-        cell_1 = Cell(block_1)
-        cell_2 = Cell(block_2)
+        cell_1 = Cell(block_1, self.mesh_points)
+        cell_2 = Cell(block_2, self.mesh_points)
 
         self.assertEqual(len(cell_1.get_common_vertices(cell_2)), count)
 
     @parameterized.expand(((0, 0, 0), (0, 1, 1), (1, 1, 0), (1, 8, 1)))
     def test_get_corner(self, block, vertex, corner):
-        cell = Cell(self.make_block(block))
+        cell = Cell(self.mesh.blocks[block], self.mesh_points)
 
         self.assertEqual(cell.get_corner(vertex), corner)
 
     @parameterized.expand(((0, 1, "right"), (1, 0, "left"), (1, 2, "back")))
     def test_get_common_side(self, index_1, index_2, orient):
-        cell_1 = Cell(self.make_block(index_1))
-        cell_2 = Cell(self.make_block(index_2))
+        cell_1 = Cell(self.mesh.blocks[index_1], self.mesh_points)
+        cell_2 = Cell(self.mesh.blocks[index_2], self.mesh_points)
 
         self.assertEqual(cell_1.get_common_side(cell_2), orient)
 
     def test_no_common_sides(self):
         with self.assertRaises(NoCommonSidesError):
-            cell_1 = Cell(self.make_block(0))
-            cell_2 = Cell(self.make_block(2))
+            cell_1 = Cell(self.mesh.blocks[0], self.mesh_points)
+            cell_2 = Cell(self.mesh.blocks[2], self.mesh_points)
 
             cell_1.get_common_side(cell_2)
 
     def test_quality_good(self):
-        cell = Cell(self.make_block(0))
+        cell = Cell(self.mesh.blocks[0], self.mesh_points)
 
         self.assertLess(cell.quality, 1)
 
     def test_quality_bad(self):
-        block = self.make_block(0)
+        block = self.mesh.blocks[0]
         block.vertices[0].move_to([-10, -10, -10])
 
-        cell = Cell(block)
+        cell = Cell(block, self.mesh_points)
 
         self.assertGreater(cell.quality, 100)

tests/test_optimize/test_optimizer.py

@@ -9,7 +9,8 @@
 from classy_blocks.modify.clamps.surface import PlaneClamp
 from classy_blocks.modify.find.geometric import GeometricFinder
 from classy_blocks.modify.iteration import IterationDriver
-from classy_blocks.modify.optimizer import ClampExistsError, Optimizer
+from classy_blocks.modify.junction import ClampExistsError
+from classy_blocks.modify.optimizer import Optimizer
 from classy_blocks.util import functions as f
 from classy_blocks.util.constants import VBIG