Partially optimize mesh assembly

ACHTUNG! Still not fixed: face merging & duplicate vertices

Author: FranzBangar

pyproject.toml

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

src/classy_blocks/assemble/cell_registry.py

@@ -0,0 +1,22 @@
+from typing import Dict, List, Set
+
+from classy_blocks.cbtyping import IndexType
+from classy_blocks.util import functions as f
+
+
+class CellRegistry:
+    def __init__(self, addressing: List[IndexType]):
+        self.addressing = addressing
+
+        # build a map of cells that belong to each point;
+        # first, find the maximum point index
+        max_index = max(f.flatten_2d_list(addressing))
+
+        self.near_cells: Dict[int, Set[int]] = {i: set() for i in range(max_index + 1)}
+
+        for i_cell, cell_indexes in enumerate(self.addressing):
+            for i_point in cell_indexes:
+                self.near_cells[i_point].add(i_cell)
+
+    def get_near_cells(self, point_index: int) -> Set[int]:
+        return self.near_cells[point_index]

src/classy_blocks/assemble/face_registry.py

@@ -1,18 +1,49 @@
-from typing import List, Tuple
+import abc
+from typing import ClassVar, Dict, List, Set, Tuple, Type
 
-FaceType = Tuple[int, int, int, int]
+from classy_blocks.cbtyping import IndexType, OrientType
+from classy_blocks.optimize.cell import CellBase, HexCell, QuadCell
 
 
-def get_key(indexes: List[int]):
-    return tuple(sorted(indexes))
+class FaceRegistryBase(abc.ABC):
+    cell_type: ClassVar[Type[CellBase]]
 
+    def get_key_from_side(self, cell: int, side: OrientType) -> Tuple:
+        cell_indexes = self.addressing[cell]
+        face_corners = self.orient_indexes[side]
+        face_indexes = [cell_indexes[c] for c in face_corners]
 
-class Face:
-    def __init__(self, indexes: FaceType):
-        self.indexes = indexes
+        return tuple(sorted(face_indexes))
 
-    def __hash__(self):
-        return hash(self.indexes)
+    def __init__(self, addressing: List[IndexType]):
+        self.addressing = addressing
 
-    def __repr__(self):
-        return f"Face {self.indexes}"
+        # this is almost equal to constants.FACE_MAP except
+        # faces are oriented so that they point towards cell centers
+        self.orient_indexes: Dict[OrientType, IndexType] = {
+            name: self.cell_type.side_indexes[i] for i, name in enumerate(self.cell_type.side_names)
+        }
+
+        # will hold faces, accessible in O(1) time, and their coincident cell(s);
+        # boundary faces will have a single coincident cell, internal two
+        self.faces: Dict[Tuple, Set[int]] = {}
+
+        for cell in range(len(self.addressing)):
+            for side in self.orient_indexes.keys():
+                face_key = self.get_key_from_side(cell, side)
+
+                if face_key not in self.faces:
+                    self.faces[face_key] = set()
+
+                self.faces[face_key].add(cell)
+
+    def get_cells(self, cell: int, side: OrientType) -> Set[int]:
+        return self.faces[self.get_key_from_side(cell, side)]
+
+
+class QuadFaceRegistry(FaceRegistryBase):
+    cell_type = QuadCell
+
+
+class HexFaceRegistry(FaceRegistryBase):
+    cell_type = HexCell

src/classy_blocks/assemble/point_registry.py

@@ -41,7 +41,7 @@ def _compile_unique(self) -> NPPointListType:
         handled_indexes: Set[int] = set()
 
         for point in self._repeated_points:
-            coincident_indexes = self._repeated_point_tree.query_ball_point(point, r=self.merge_tol, workers=-1)
+            coincident_indexes = self._repeated_point_tree.query_ball_point(point, r=self.merge_tol, workers=1)
 
             if set(coincident_indexes).isdisjoint(handled_indexes):
                 # this vertex hasn't been handled yet
@@ -52,15 +52,15 @@ def _compile_unique(self) -> NPPointListType:
 
     def _query_unique(self, positions) -> List[int]:
         """A shortcut to KDTree.query_ball_point()"""
-        result = self._unique_point_tree.query_ball_point(positions, r=self.merge_tol, workers=-1)
+        result = self._unique_point_tree.query_ball_point(positions, r=self.merge_tol, workers=1)
 
         if len(np.shape(positions)) > 1:
             return f.flatten_2d_list(result)
 
         return result
 
     def _query_repeated(self, positions) -> List[int]:
-        result = self._repeated_point_tree.query_ball_point(positions, r=self.merge_tol, workers=-1)
+        result = self._repeated_point_tree.query_ball_point(positions, r=self.merge_tol, workers=1)
 
         if len(np.shape(positions)) > 1:
             return f.flatten_2d_list(result)

src/classy_blocks/construct/assemblies/joints.py

@@ -36,7 +36,7 @@ def __init__(
         for loft in self.shell:
             point_1 = loft.top_face.points[1].position
             point_2 = loft.top_face.points[2].position
-            edge_points = f.divide_arc(axis, axis_point_2, point_1, point_2, 5)
+            edge_points = f.divide_arc(axis_point_2, point_1, point_2, 5)
             loft.top_face.add_edge(1, Spline(edge_points))
 
         for loft in self.operations:

src/classy_blocks/construct/curves/curve.py

@@ -60,19 +60,7 @@ def length(self) -> float:
 
     @abc.abstractmethod
     def get_closest_param(self, point: PointType) -> float:
-        """Finds the parameter on curve where point is the closest to given point;
-        To improve search speed and reliability, an optional starting
-        estimation can be supplied."""
-        # because curves can have all sorts of shapes, find
-        # initial guess by checking distance to discretized points
-        point = np.array(point)
-        all_points = self.discretize()
-
-        distances = np.array([f.norm(p - point) for p in all_points])
-        params = np.linspace(self.bounds[0], self.bounds[1], num=len(distances))
-
-        i_distance = np.argmin(distances)
-        return params[i_distance]
+        """Finds the parameter on curve where point is the closest to given point."""
 
     def get_param_at_length(self, length: float) -> float:
         """Returns parameter at specified length along the curve"""
@@ -135,25 +123,28 @@ def discretize(
         return np.array([self.function(t) for t in params])
 
     def get_closest_param(self, point: PointType) -> float:
-        """Finds the param on curve where point is the closest to given point;
-        To improve search speed and reliability, an optional starting
-        estimation can be supplied."""
-        closest_param = super().get_closest_param(point)
+        """Finds the param on curve where point is the closest to given point"""
+        # because curves can have all sorts of shapes, find
+        # initial guess by checking distance to discretized points
         point = np.array(point)
+        all_points = self.discretize()
 
-        result = scipy.optimize.minimize(
-            lambda t: f.norm(self.get_point(t[0]) - point), (closest_param,), bounds=(self.bounds,)
-        )
+        distances = np.linalg.norm(all_points.T - point[:, None], axis=0)
+        params = np.linspace(self.bounds[0], self.bounds[1], num=len(distances))
 
-        return result.x[0]
+        i_smallest = int(np.argmin(distances))
+        i_prev = max(i_smallest - 1, 0)
+        i_next = min(i_smallest + 1, len(distances) - 1)
+
+        param_start = params[i_prev]
+        param_end = params[i_next]
+
+        result = scipy.optimize.minimize_scalar(
+            lambda t: f.norm(self.get_point(t) - point),
+            bounds=(param_start, param_end),
+        )
 
-        # TODO: optimize by using a simpler function (?)
-        # param_start = super().get_closest_param(point)
-        # param_end = (1 + closest_param) / 2
-        # result = scipy.optimize.minimize_scalar(
-        #     lambda t: f.norm(self.get_point(t) - point), bounds=(param_start, param_end), method="bounded"
-        # )
-        # return result.x
+        return result.x
 
     def get_point(self, param: float) -> NPPointType:
         self._check_param(param)

src/classy_blocks/construct/curves/discrete.py

@@ -46,10 +46,15 @@ def get_length(self, param_from: Optional[float] = None, param_to: Optional[floa
 
     def get_closest_param(self, point: PointType) -> float:
         """Returns the index of point on this curve where distance to supplied
-        point is the smallest.
+        point is the smallest."""
+        point = np.array(point)
+        all_points = self.discretize()
 
-        For DiscreteCurve, argument param_start is ignored since all points are checked."""
-        return super().get_closest_param(point)
+        distances = np.linalg.norm(all_points.T - point[:, None], axis=0)
+        params = np.linspace(self.bounds[0], self.bounds[1], num=len(distances))
+
+        i_distance = np.argmin(distances)
+        return params[i_distance]
 
     @property
     def center(self):

src/classy_blocks/grading/autograding/grader.py

@@ -15,6 +15,7 @@ def __init__(self, mesh: Mesh, rules: ChopRules):
         self.rules = rules
 
         self.mesh.assemble()
+        self.mesh.block_list.update_lengths()  # TODO: move to a better place
         self.probe = Probe(self.mesh)
 
     def _chop_wire(self, wire: Wire, chops: List[Chop]) -> None:

src/classy_blocks/items/edges/arcs/angle.py

@@ -37,7 +37,7 @@ def arc_from_theta(edge_point_1: PointType, edge_point_2: PointType, angle: floa
 
     center = pm - length * axis / 2 - rm * mag_chord / 2 / np.tan(angle / 2)
 
-    return f.arc_mid(axis, center, edge_point_1, edge_point_2)
+    return f.arc_mid(center, edge_point_1, edge_point_2)
 
 
 @dataclasses.dataclass

src/classy_blocks/items/edges/arcs/origin.py

@@ -73,7 +73,7 @@ def arc_from_origin(
         return arc_from_origin(p1, p3, new_center, False)
 
     # done, return the calculated point
-    return f.arc_mid(axis, center, edge_point_1, edge_point_2)
+    return f.arc_mid(center, edge_point_1, edge_point_2)
 
 
 @dataclasses.dataclass
@@ -82,6 +82,7 @@ class OriginEdge(ArcEdgeBase):
 
     data: edges.Origin
 
+    # TODO: make this accessible to the user
     adjust_center: ClassVar[bool] = True
 
     @property