Speed up assembly with registry (and other goodies)

- Improve grading error messages
- Change the way deleted operations work (just comment out block output
  in blockMesh but still take all the data from deleted operations*)
- Cache wire length and store wires and edges in a dict for faster addressing
- Utilize scipy's KDTree to get rid of plain linear searches for block
  neighbours (and in many other parts in the code**)

* TODO: might result in an error if the block carries an edge that is
  not shared by any other block
** TODO: get rid of sketch mapper that still uses that old linear search

TODO! Face merging is broken! Put the functionality back and write tests
that confirm that.

Author: FranzBangar

examples/advanced/smooth_grader.py

@@ -28,7 +28,7 @@
     vertex.translate([0, 0.8, 0])
 
 mesh.set_default_patch("walls", "wall")
-mesh.block_list.update_lengths()
+
 grader = cb.SmoothGrader(mesh, 0.05)
 grader.grade()
 

examples/shape/cylinder.py

@@ -4,9 +4,9 @@
 
 mesh = cb.Mesh()
 
-axis_point_1 = [0.0, 0.0, 0.0]
-axis_point_2 = [5.0, 5.0, 0.0]
-radius_point_1 = [0.0, 0.0, 2.0]
+axis_point_1 = [0, 0, 0]
+axis_point_2 = [0, 0, 1]
+radius_point_1 = [1, 0, 0]
 
 cylinder = cb.Cylinder(axis_point_1, axis_point_2, radius_point_1)
 

src/classy_blocks/assemble/connection_registry.py

@@ -0,0 +1,65 @@
+from typing import ClassVar, Dict, List, Set, Tuple
+
+from classy_blocks.cbtyping import IndexType, NPPointListType
+from classy_blocks.util.constants import EDGE_PAIRS
+
+ConnectionType = Tuple[int, int]
+
+
+def get_key(index_1: int, index_2: int):
+    return (min(index_1, index_2), max(index_1, index_2))
+
+
+class Connection:
+    """Connects two points by their index in unique point list"""
+
+    def __init__(self, index_1: int, index_2: int):
+        self.indexes = get_key(index_1, index_2)
+
+    def get_other_index(self, index):
+        if self.indexes[0] == index:
+            return self.indexes[1]
+
+        return self.indexes[0]
+
+    def __repr__(self):
+        return f"Connection {self.indexes[0]}-{self.indexes[1]}"
+
+    def __hash__(self):
+        return hash(self.indexes)
+
+
+class ConnectionRegistryBase:
+    edge_pairs: ClassVar[List[ConnectionType]]
+
+    def __init__(self, points: NPPointListType, addressing: List[IndexType]):
+        self.points = points
+        self.addressing = addressing
+
+        self.connections: Dict[ConnectionType, Connection] = {}
+        self.nodes: Dict[int, Set[Connection]] = {i: set() for i in range(len(self.points))}
+
+        for i in range(len(self.addressing)):
+            cell_indexes = self.addressing[i]
+            for pair in self.edge_pairs:
+                index_1 = cell_indexes[pair[0]]
+                index_2 = cell_indexes[pair[1]]
+
+                edge_key = get_key(index_1, index_2)
+
+                if edge_key not in self.connections:
+                    connection = Connection(index_1, index_2)
+                    self.connections[edge_key] = connection
+                    self.nodes[index_1].add(connection)
+                    self.nodes[index_2].add(connection)
+
+    def get_connected_indexes(self, index: int) -> Set[int]:
+        return {c.get_other_index(index) for c in self.nodes[index]}
+
+
+class QuadConnectionRegistry(ConnectionRegistryBase):
+    edge_pairs: ClassVar = [(0, 1), (1, 2), (2, 3), (3, 0)]
+
+
+class HexConnectionRegistry(ConnectionRegistryBase):
+    edge_pairs: ClassVar = EDGE_PAIRS

src/classy_blocks/assemble/face_registry.py

@@ -0,0 +1,18 @@
+from typing import List, Tuple
+
+FaceType = Tuple[int, int, int, int]
+
+
+def get_key(indexes: List[int]):
+    return tuple(sorted(indexes))
+
+
+class Face:
+    def __init__(self, indexes: FaceType):
+        self.indexes = indexes
+
+    def __hash__(self):
+        return hash(self.indexes)
+
+    def __repr__(self):
+        return f"Face {self.indexes}"

src/classy_blocks/assemble/navigator.py renamed to src/classy_blocks/assemble/point_registry.py

@@ -1,23 +1,27 @@
-from typing import List, Set, Tuple, Type, TypeVar
+from typing import List, Set, Type, TypeVar
 
 import numpy as np
 import scipy.spatial
 
 from classy_blocks.base.exceptions import VertexNotFoundError
 from classy_blocks.cbtyping import IndexType, NPPointListType, PointType
+from classy_blocks.construct.flat.sketch import Sketch
 from classy_blocks.construct.operations.operation import Operation
 from classy_blocks.util import functions as f
 from classy_blocks.util.constants import TOL, vector_format
 
-NavigatorType = TypeVar("NavigatorType", bound="NavigatorBase")
+PointRegistryType = TypeVar("PointRegistryType", bound="PointRegistryBase")
 
 
-class NavigatorBase:
-    """Searches for, connects and creates unique indexes points, taken from lists of operations"""
+class PointRegistryBase:
+    """Searches for, connects and creates unique points, taken from lists of elements (quads, hexas, ...)"""
 
     cell_size: int  # 4 for quads, 8 for hexas
 
     def __init__(self, flattened_points: NPPointListType, merge_tol: float) -> None:
+        if len(flattened_points) % self.cell_size != 0:
+            raise ValueError(f"Number of points not divisible by cell_size: {len(flattened_points)} % {self.cell_size}")
+
         self.merge_tol = merge_tol
 
         # a flattened list of all points, possibly multiple at the same spot;
@@ -29,6 +33,8 @@ def __init__(self, flattened_points: NPPointListType, merge_tol: float) -> None:
         self.unique_points = self._compile_unique()
         self._unique_point_tree = scipy.spatial.KDTree(self.unique_points)
 
+        self.cell_addressing = [self._query_unique(self.find_cell_points(i)) for i in range(self.cell_count)]
+
     def _compile_unique(self) -> NPPointListType:
         # create a list of unique vertices, taken from the list of operations
         unique_points = []
@@ -83,10 +89,11 @@ def find_cell_points(self, cell: int) -> NPPointListType:
 
     def find_cell_indexes(self, cell: int) -> List[int]:
         """Returns indexes of points that define this cell"""
-        return self._query_unique(self.find_cell_points(cell))
+        return self.cell_addressing[cell]
 
     def find_cell_neighbours(self, cell: int) -> List[int]:
         """Returns indexes of this and every touching cell"""
+        # TODO: remove
         cell_points = self.find_cell_points(cell)
 
         indexes = []
@@ -96,38 +103,50 @@ def find_cell_neighbours(self, cell: int) -> List[int]:
 
         return list(set(indexes))
 
-    def get_map(self) -> Tuple[NPPointListType, List[IndexType]]:
-        cell_indexes: List[IndexType] = []
-
-        for i in range(len(self._repeated_points) // self.cell_size):
-            cell_indexes.append(self.find_cell_indexes(i))
-
-        return self.unique_points, cell_indexes
-
     @staticmethod
     def flatten(points, length) -> NPPointListType:
         return np.reshape(points, (length, 3))
 
     @classmethod
     def from_addresses(
-        cls: Type[NavigatorType], points: NPPointListType, addressing: List[IndexType], merge_tol: float = TOL
-    ) -> NavigatorType:
+        cls: Type[PointRegistryType], points: NPPointListType, addressing: List[IndexType], merge_tol: float = TOL
+    ) -> PointRegistryType:
         all_points = cls.flatten(
             [np.take(points, addr, axis=0) for addr in addressing], len(addressing) * cls.cell_size
         )
 
         return cls(all_points, merge_tol)
 
+    @property
+    def cell_count(self) -> int:
+        return len(self._repeated_points) // self.cell_size
+
+    @property
+    def point_count(self) -> int:
+        return len(self.unique_points)
+
+
+class QuadPointRegistry(PointRegistryBase):
+    """A registry of points, taken from a list of quads"""
 
-class QuadNavigator(NavigatorBase):
     cell_size = 4
 
+    @classmethod
+    def from_sketch(cls: Type[PointRegistryType], sketch: Sketch, merge_tol: float = TOL) -> PointRegistryType:
+        return cls(
+            cls.flatten([face.point_array for face in sketch.faces], len(sketch.faces) * cls.cell_size), merge_tol
+        )
+
+
+class HexPointRegistry(PointRegistryBase):
+    """A registry of points, taken from a list of hexas"""
 
-class HexNavigator(NavigatorBase):
     cell_size = 8
 
     @classmethod
-    def from_operations(cls: Type[NavigatorType], operations: List[Operation], merge_tol: float = TOL) -> NavigatorType:
+    def from_operations(
+        cls: Type[PointRegistryType], operations: List[Operation], merge_tol: float = TOL
+    ) -> PointRegistryType:
         all_points = cls.flatten([op.point_array for op in operations], len(operations) * cls.cell_size)
 
         return cls(all_points, merge_tol)

src/classy_blocks/construct/curves/curve.py

@@ -138,15 +138,23 @@ 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."""
-        param_start = super().get_closest_param(point)
+        closest_param = super().get_closest_param(point)
         point = np.array(point)
 
         result = scipy.optimize.minimize(
-            lambda t: f.norm(self.get_point(t[0]) - point), (param_start,), bounds=(self.bounds,)
+            lambda t: f.norm(self.get_point(t[0]) - point), (closest_param,), bounds=(self.bounds,)
         )
 
         return result.x[0]
 
+        # 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
+
     def get_point(self, param: float) -> NPPointType:
         self._check_param(param)
         return self.function(param)

src/classy_blocks/items/block.py

@@ -38,6 +38,10 @@ def __init__(self, index: int, vertices: List[Vertex]):
         # (visible in blockMeshDict, useful for debugging)
         self.comment: str = ""
 
+        # 'hidden' blocks carry all of the data from the mesh
+        # but is not inserted into blockMeshDict;
+        self.visible = True
+
     def add_edge(self, corner_1: int, corner_2: int, edge: Edge):
         """Adds an edge between vertices at specified indexes."""
         if not (0 <= corner_1 < 8 and 0 <= corner_2 < 8):
@@ -146,13 +150,14 @@ def format_grading(self) -> str:
     @property
     def description(self) -> str:
         """hex definition for blockMesh"""
+        fmt_hidden = "" if self.visible else "// "
         fmt_vertices = "( " + " ".join(str(v.index) for v in self.vertices) + " )"
         fmt_count = "( " + " ".join([str(axis.count) for axis in self.axes]) + " )"
 
         fmt_grading = self.format_grading()
         fmt_comments = f"// {self.index} {self.comment}\n"
 
-        return f"\thex {fmt_vertices} {self.cell_zone} {fmt_count} {fmt_grading} {fmt_comments}"
+        return f"\t{fmt_hidden}hex {fmt_vertices} {self.cell_zone} {fmt_count} {fmt_grading} {fmt_comments}"
 
     def __hash__(self) -> int:
         return self.index

src/classy_blocks/lists/block_list.py

@@ -1,6 +1,6 @@
 from typing import List, Set
 
-from classy_blocks.assemble.navigator import HexNavigator
+from classy_blocks.assemble.point_registry import HexPointRegistry
 from classy_blocks.base.exceptions import UndefinedGradingsError
 from classy_blocks.items.block import Block
 
@@ -16,7 +16,7 @@ def add(self, block: Block) -> None:
         """Add blocks"""
         self.blocks.append(block)
 
-    def update_neighbours(self, navigator: HexNavigator) -> None:
+    def update_neighbours(self, navigator: HexPointRegistry) -> None:
         """Find and assign neighbours of a given block entry"""
         for block in self.blocks:
             neighbour_indexes = navigator.find_cell_neighbours(block.index)

src/classy_blocks/mesh.py

@@ -2,7 +2,7 @@
 
 from typing import List, Optional, Set, Union, get_args
 
-from classy_blocks.assemble.navigator import HexNavigator
+from classy_blocks.assemble.point_registry import HexPointRegistry
 from classy_blocks.base.exceptions import EdgeNotFoundError
 from classy_blocks.cbtyping import DirectionType
 from classy_blocks.construct.assemblies.assembly import Assembly
@@ -99,22 +99,6 @@ def delete(self, operation: Operation) -> None:
         the data remains but it will not contribute to the mesh"""
         self.deleted.add(operation)
 
-    @property
-    def _operations(self) -> List[Operation]:
-        operations: List[Operation] = []
-
-        for entity in self.depot:
-            if isinstance(entity, Operation):
-                ops_to_add = [entity]
-            else:
-                ops_to_add = entity.operations
-
-            for op in ops_to_add:
-                if op not in self.deleted:
-                    operations.append(op)
-
-        return operations
-
     def _add_geometry(self) -> None:
         for entity in self.depot:
             if entity.geometry is not None:
@@ -130,8 +114,8 @@ def assemble(self) -> None:
             self.block_list.update_lengths()
             return
 
-        operations = self._operations
-        navigator = HexNavigator.from_operations(operations)
+        operations = self.operations
+        navigator = HexPointRegistry.from_operations(operations)
 
         self.vertex_list.vertices = [Vertex(pos, i) for i, pos in enumerate(navigator.unique_points)]
 
@@ -164,6 +148,7 @@ def assemble(self) -> None:
                 block.add_chops(direction, operation.chops[direction])
 
             block.cell_zone = operation.cell_zone
+            block.visible = operation not in self.deleted
 
             self.block_list.add(block)