Improve speed and quality of optimization

Author: FranzBangar

.vscode/settings.json

@@ -7,8 +7,6 @@
             "source.organizeImports.ruff": "always"
         },
     },
-    "python.analysis.typeCheckingMode": "basic",
-    "python.analysis.diagnosticMode": "openFilesOnly",
     "python.testing.pytestEnabled": false,
     "files.exclude": {
         "**/.git": true,

classy_blocks.code-workspace

@@ -13,16 +13,14 @@
             "scipy"
         ],
         "editor.autoClosingBrackets": "never",
-        "isort.importStrategy": "fromEnvironment",
-        "editor.defaultFormatter": null,
+        "editor.defaultFormatter": "charliermarsh.ruff",
         "ruff.organizeImports": true
     },
     "extensions": {
         "recommendations": [
             "ms-python.vscode-pylance",
             "ms-python.python",
-            "charliermarsh.ruff",
-            "ms-python.black-formatter"
+            "charliermarsh.ruff"
         ]
     },
 }

docs/optimization/quality-criteria.ods

@@ -7,7 +7,7 @@ license = "MIT"
 keywords = ["classy_blocks", "OpenFOAM", "blockMesh"]
 authors = [{ name = "Nejc Jurkovic", email = "kandelabr@gmail.com" }]
 requires-python = ">=3.8"
-dependencies = ["numpy", "scipy", "nptyping", "numba"]
+dependencies = ["numpy", "scipy", "nptyping", "numba", "parameterized"]
 
 [project.urls]
 "Homepage" = "https://github.com/damogranlabs/classy_blocks"

pyproject.toml

@@ -49,7 +49,7 @@ def fquality(params):
             return self.grid.update(junction.index, clamp.position)
 
         try:
-            scipy.optimize.minimize(fquality, clamp.params, bounds=clamp.bounds, method=method)
+            scipy.optimize.minimize(fquality, clamp.params, bounds=clamp.bounds, method=method, tol=1e-2)
 
             reporter.junction_final = junction.quality
             reporter.grid_final = self.grid.quality

src/classy_blocks/optimize/optimizer.py

@@ -17,12 +17,12 @@ 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.25, 0.35, 0.8, angle)
+    return scale_quality(1.5, 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.15, np.abs(angle))
+    return scale_quality(1.5, 0.25, 0.5, np.abs(angle))
 
 
 @numba.jit(nopython=True, cache=True)
@@ -146,12 +146,12 @@ def get_hex_quality(grid_points: NPPointListType, cell_indexes: NPIndexType) ->
     quality = 0
 
     for side in side_indexes:
-        # Non-ortho angle is measured between two lines;
-        # neighbour_center-to-this_center and face_center-to-this_center,
-        # but for cells on the boundary simply face center is taken.
-        # For this kind of optimization it is quite sufficient to take
-        # only the latter as it's not much different and we'll optimize
-        # other cells too.
+        # Non-ortho angle in a hexahedron is measured between two vectors:
+        # <neighbour_center-this_center> and <face_normal>
+        # but since cells on the boundary don't have a neighbour
+        # simply <face_center> is taken.
+        # For this kind of optimization it is quite sufficient to
+        # take <face_center> for all cells.
         side_points = take(cell_points, side)
         side_center, side_normal, side_aspect = get_quad_normal(side_points)
         center_vector = cell_center - side_center

src/classy_blocks/optimize/quality.py

@@ -61,27 +61,40 @@ def angle_between(vect_1: VectorType, vect_2: VectorType) -> float:
     return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
 
 
-def rotation_matrix(axis: VectorType, theta: float):
-    """
-    Return the rotation matrix associated with counterclockwise rotation about
-    the given axis by theta radians.
-    # Kudos to
+@jit(nopython=True, cache=True)
+def _rotation_matrix(axis: NPVectorType, angle: float):
+    axis = axis / np.sqrt(np.dot(axis, axis))
+    a = np.cos(angle / 2.0)
+    b, c, d = -axis * np.sin(angle / 2.0)
+    aa, bb, cc, dd = a * a, b * b, c * c, d * d
+    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
+
+    return np.array(
+        [
+            [aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
+            [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
+            [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc],
+        ]
+    )
+
+
+def rotation_matrix(axis: VectorType, angle: float):
+    """Returns the rotation matrix associated with counterclockwise rotation about
+    the given axis by theta radians."""
+    # Kudos to https://stackoverflow.com/questions/6802577/rotation-of-3d-vector
+    axis = np.asarray(axis, dtype=constants.DTYPE)
+
+    return _rotation_matrix(axis, angle)
+
+
+@jit(nopython=True, cache=True)
+def _rotate(point: NPPointType, angle: float, axis: NPVectorType, origin: NPPointType) -> NPPointType:
+    # Rotation matrix associated with counterclockwise rotation about
+    # the given axis by theta radians. Kudos to
     # https://stackoverflow.com/questions/6802577/rotation-of-3d-vector
-    #import math
-    #
-    # axis = np.asarray(axis)
-    # axis = axis / math.sqrt(np.dot(axis, axis))
-    # a = math.cos(theta / 2.0)
-    # b, c, d = -axis * math.sin(theta / 2.0)
-    # aa, bb, cc, dd = a * a, b * b, c * c, d * d
-    # bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
-    # return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
-    #                  [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
-    #                  [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
-
-    # Also Kudos to the guy with another answer for the same question (used here):"""
-    axis = np.asarray(axis)
-    return scipy.linalg.expm(np.cross(np.eye(3), axis / norm(axis) * theta))
+
+    rotated_point = np.dot(_rotation_matrix(axis, angle), point - origin)
+    return rotated_point + origin
 
 
 def rotate(point: PointType, angle: float, axis: VectorType, origin: PointType) -> NPPointType:
@@ -90,8 +103,7 @@ def rotate(point: PointType, angle: float, axis: VectorType, origin: PointType)
     axis = np.asarray(axis, dtype=constants.DTYPE)
     origin = np.asarray(origin, dtype=constants.DTYPE)
 
-    rotated_point = np.dot(rotation_matrix(axis, angle), point - origin)
-    return rotated_point + origin
+    return _rotate(point, angle, axis, origin)
 
 
 def scale(point: PointType, ratio: float, origin: Optional[PointType]) -> NPPointType:
@@ -201,14 +213,6 @@ def arc_length_3point(p_start: NPPointType, p_btw: NPPointType, p_end: NPPointTy
     return angle * np.linalg.norm(radius)
 
 
-# def arc_length_3point(p_start: NPPointType, p_btw: NPPointType, p_end: NPPointType) -> float:
-#     p_start = np.asarray(p_start, dtype=np.float64)
-#     p_btw = np.asarray(p_btw, dtype=np.float64)
-#     p_end = np.asarray(p_end, dtype=np.float64)
-
-#     return _arc_length_3point(p_start, p_btw, p_end)
-
-
 @jit(nopython=True, cache=True)
 def divide_arc(center: NPPointType, point_1: NPPointType, point_2: NPPointType, count: int) -> NPPointListType:
     radius = np.linalg.norm(center - point_1)
@@ -225,17 +229,6 @@ def divide_arc(center: NPPointType, point_1: NPPointType, point_2: NPPointType,
     return result
 
 
-# def divide_arc(center: PointType, point_1: PointType, point_2: PointType, count: int) -> NPPointListType:
-#     # Kudos to this guy for his shrewd solution
-#     # https://math.stackexchange.com/questions/3717427
-#     # (extended here to create more than 1 point)
-#     center = np.asarray(center, dtype=np.float64)
-#     point_1 = np.asarray(point_1, dtype=np.float64)
-#     point_2 = np.asarray(point_2, dtype=np.float64)
-
-#     return _divide_arc(center, point_1, point_2, count)
-
-
 def arc_mid(center: PointType, point_1: PointType, point_2: PointType) -> PointType:
     """Returns the midpoint of the specified arc in 3D space"""
     return divide_arc(center, point_1, point_2, 1)[0]

src/classy_blocks/util/functions.py

@@ -45,7 +45,7 @@ def check_success(self, chained_shape, end_center):
         self.assertEqual(len(self.mesh.blocks), 24)
         self.assertEqual(len(self.mesh.vertices), 3 * 17)
 
-        np.testing.assert_allclose(chained_shape.sketch_2.center, end_center)
+        np.testing.assert_allclose(chained_shape.sketch_2.center, end_center, atol=1e-7)
 
     def test_to_elbow_end(self):
         """Chain an elbow to an elbow on an end sketch"""

tests/test_construct/test_chaining.py

@@ -28,14 +28,14 @@ def test_arc_edge_translate(self):
         arc_edge = ArcEdge(self.vertex_1, self.vertex_2, edges.Arc([0.5, 0, 0]))
 
         arc_edge.translate([1, 1, 1])
-        np.testing.assert_array_equal(arc_edge.data.point.position, [1.5, 1, 1])
+        np.testing.assert_almost_equal(arc_edge.data.point.position, [1.5, 1, 1])
 
     def test_angle_edge_rotate_1(self):
         angle_edge = AngleEdge(self.vertex_1, self.vertex_2, edges.Angle(np.pi / 2, [0, 0, 1]))
 
         angle_edge.rotate(np.pi / 2, [0, 0, 1], [0, 0, 0])
 
-        np.testing.assert_array_equal(angle_edge.data.axis.components, [0, 0, 1])
+        np.testing.assert_almost_equal(angle_edge.data.axis.components, [0, 0, 1])
 
         self.assertEqual(angle_edge.data.angle, np.pi / 2)
 
@@ -44,7 +44,7 @@ def test_angle_edge_rotate_2(self):
 
         angle_edge.rotate(np.pi / 2, [1, 0, 0], [0, 0, 0])
 
-        np.testing.assert_array_equal(angle_edge.data.axis.components, [0, -1, 0])
+        np.testing.assert_almost_equal(angle_edge.data.axis.components, [0, -1, 0])
 
         self.assertEqual(angle_edge.data.angle, np.pi / 2)
 
@@ -63,7 +63,7 @@ def test_spline_edge_translate(self):
 
         spline_edge.translate([1, 1, 1])
 
-        np.testing.assert_array_equal(
+        np.testing.assert_almost_equal(
             spline_edge.point_array,
             [
                 [1.25, 1.1, 1],
@@ -119,21 +119,21 @@ def test_angle(self):
 
         self.assertIsInstance(edge, AngleEdge)
         self.assertEqual(angle, edge.data.angle)
-        np.testing.assert_array_equal(axis, edge.data.axis.components)
+        np.testing.assert_almost_equal(axis, edge.data.axis.components)
 
     def test_spline(self):
         points = [[0.3, 0.25, 0], [0.6, 0.1, 0], [0.3, 0.25, 0]]
         edg = factory.create(self.vertex_1, self.vertex_2, edges.Spline(points))
 
         self.assertIsInstance(edg, SplineEdge)
-        np.testing.assert_array_equal(points, edg.point_array)
+        np.testing.assert_almost_equal(points, edg.point_array)
 
     def test_polyline(self):
         points = [[0.3, 0.25, 0], [0.6, 0.1, 0], [0.3, 0.25, 0]]
         edg = factory.create(self.vertex_1, self.vertex_2, edges.PolyLine(points))
 
         self.assertIsInstance(edg, SplineEdge)
-        np.testing.assert_array_equal(points, edg.point_array)
+        np.testing.assert_almost_equal(points, edg.point_array)
 
     def test_project_edge_single(self):
         label = "terrain"

tests/test_items/test_edge.py

@@ -14,7 +14,7 @@ def test_translate(self):
         link.leader = np.array([3, 3, 3])
         link.update()
 
-        np.testing.assert_equal(link.follower, [4, 4, 4])
+        np.testing.assert_almost_equal(link.follower, [4, 4, 4])
 
 
 class RotationLinkTests(unittest.TestCase):
@@ -36,14 +36,14 @@ def test_radius(self, origin):
 
         link = RotationLink(leader, follower, [0, 0, 1], origin)
 
-        np.testing.assert_equal(link._get_radius(link.leader), [1, 0, 0])
+        np.testing.assert_almost_equal(link._get_radius(link.leader), [1, 0, 0])
 
     def test_rotate(self):
         link = RotationLink(self.leader, self.follower, [0, 0, 1], [0, 0, 0])
         link.leader = np.array([0, 1, 0])
         link.update()
 
-        np.testing.assert_equal(link.follower, [-1, 0, 0])
+        np.testing.assert_almost_equal(link.follower, [-1, 0, 0])
 
     def test_rotate_negative(self):
         """Rotate in negative direction"""
@@ -52,7 +52,7 @@ def test_rotate_negative(self):
         link.leader = np.array([0, -1, 0])
         link.update()
 
-        np.testing.assert_equal(link.follower, [1, 0, 0])
+        np.testing.assert_almost_equal(link.follower, [1, 0, 0])
 
     @parameterized.expand(
         (
@@ -108,4 +108,4 @@ def test_move(self):
         link.leader = np.array([0, 0, 2])
         link.update()
 
-        np.testing.assert_equal(link.follower, [0, 0, -2])
+        np.testing.assert_almost_equal(link.follower, [0, 0, -2])