Merge branch 'development' into feature/autochop

README.md

@@ -367,6 +367,13 @@ Airfoil core with blunt trailing edge (imported points from NACA generator) and
 (see `examples/complex/airfoil.py`). A simulation-ready mesh needs additional blocks to expand domain further away from the airfoil.
 ![Airfoil](showcase/airfoil.png "Airfoil core mesh")
 
+## Automatic Edge Grading
+When setting cell counts and expansion ratios, it is possible to specify which value to keep constant. Mostly this will be used for keeping thickness of the first cell at the wall consistent to maintain desired `y+` throughout the mesh. This is done by simple specifying a `preserve="..."` keyword.
+
+Example: a block chopped in a classic way where cell sizes will increase when block size increases:
+![Classic block grading](showcase/classy_classic_grading.png "Classic block grading")
+The same case but with a specified `preserve="start_size"` keyword for the bottom and `preserve="end_size"` for the top edge:
+![Grading for consistent cell size](showcase/classy_edges_grading.png "Classy block grading")
 
 ## Debugging
 

examples/optimization/duct.py

@@ -0,0 +1,39 @@
+# An example where a Shape is optimized *before* it is added to mesh, using ShapeOptimizer
+import os
+
+import classy_blocks as cb
+from classy_blocks.optimize.junction import ClampExistsError
+from classy_blocks.optimize.optimizer import ShapeOptimizer
+
+mesh = cb.Mesh()
+
+start_sketch = cb.SplineDisk([0, 0, 0], [2.5, 0, 0], [0, 1, 0], 0, 0)
+end_sketch = cb.SplineDisk([0, 0, 0], [1, 0, 0], [0, 2.5, 0], 0, 0).translate([0, 0, 1])
+
+shape = cb.LoftedShape(start_sketch, end_sketch)
+
+optimizer = ShapeOptimizer(shape.operations)
+
+for operation in shape.operations[:4]:
+    # remove edges because inner splines will ruin the day
+    # TODO: make edges move with points too
+    operation.top_face.remove_edges()
+    operation.bottom_face.remove_edges()
+
+    for point in operation.point_array:
+        try:
+            optimizer.add_clamp(cb.FreeClamp(point))
+        except ClampExistsError:
+            pass
+
+optimizer.optimize(tolerance=0.01)
+
+# Quick'n'dirty chopping, don't do this at home
+for operation in shape.operations:
+    for axis in range(3):
+        operation.chop(axis, 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/construct/curves/interpolated.py

@@ -25,9 +25,9 @@ class InterpolatedCurveBase(FunctionCurveBase, abc.ABC):
 
     _interpolator: Type[InterpolatorBase]
 
-    def __init__(self, points: PointListType):
+    def __init__(self, points: PointListType, extrapolate: bool = False, equalize: bool = True):
         self.array = Array(points)
-        self.function = self._interpolator(self.array, False)
+        self.function = self._interpolator(self.array, extrapolate, equalize)
         self.bounds = (0, 1)
 
     @property

src/classy_blocks/construct/curves/interpolators.py

@@ -2,10 +2,9 @@
 
 import numpy as np
 import scipy.interpolate
-from numpy.typing import NDArray
 
 from classy_blocks.construct.array import Array
-from classy_blocks.types import NPPointType, ParamCurveFuncType
+from classy_blocks.types import FloatListType, NPPointType, ParamCurveFuncType
 
 
 class InterpolatorBase(abc.ABC):
@@ -19,9 +18,10 @@ class InterpolatorBase(abc.ABC):
     def _get_function(self) -> ParamCurveFuncType:
         """Returns an interpolation function from stored points"""
 
-    def __init__(self, points: Array, extrapolate: bool):
+    def __init__(self, points: Array, extrapolate: bool, equalize: bool = True):
         self.points = points
         self.extrapolate = extrapolate
+        self.equalize = equalize
 
         self.function = self._get_function()
         self._valid = True
@@ -37,7 +37,16 @@ def invalidate(self) -> None:
         self._valid = False
 
     @property
-    def params(self) -> NDArray:
+    def params(self) -> FloatListType:
+        """A list of parameters for the interpolation curve.
+        If not equalized, it's just linearly spaced floats;
+        if equalized, scaled distances between provided points are taken so that
+        evenly spaced parameters will produce evenly spaced points even if
+        interpolation points are unequally spaced."""
+        if self.equalize:
+            lengths = np.cumsum(np.sqrt(np.sum((self.points[:-1] - self.points[1:]) ** 2, axis=1)))
+            return np.concatenate(([0], lengths / lengths[-1]))
+
         return np.linspace(0, 1, num=len(self.points))
 
 

src/classy_blocks/construct/flat/face.py

@@ -105,6 +105,8 @@ def project_edge(self, corner: int, label: ProjectToType) -> None:
     def invert(self) -> "Face":
         """Reverses the order of points in this face."""
         self.points.reverse()
+        self.edges.reverse()
+        self.edges = [self.edges[i] for i in (1, 2, 3, 0)]
 
         return self
 
@@ -168,8 +170,7 @@ def project(self, label: str, edges: bool = False, points: bool = False) -> None
                 self.points[i].project(label)
 
     def shift(self, count: int) -> "Face":
-        """Shifts points of this face by 'count', changing its
-        starting point"""
+        """Shifts points of this face by 'count', changing its starting point"""
         indexes = collections.deque(range(4))
         indexes.rotate(count)
 

src/classy_blocks/construct/flat/sketches/spline_round.py

@@ -336,7 +336,7 @@ def __init__(
     @property
     def grid(self) -> List[List[Face]]:
         if len(self.faces) > 6:
-            return [self.faces[:4], self.faces[4:]]
+            return [self.faces[::3], [face for i, face in enumerate(self.faces) if not i % 3 == 0]]
         else:
             return super().grid
 

src/classy_blocks/grading/autograding/probe.py

@@ -1,12 +1,15 @@
 import functools
-from typing import Dict, List, Optional, get_args
+from typing import Dict, List, Optional, Set, get_args
 
 from classy_blocks.base.exceptions import BlockNotFoundError, NoInstructionError
 from classy_blocks.items.block import Block
+from classy_blocks.items.vertex import Vertex
 from classy_blocks.items.wires.axis import Axis
 from classy_blocks.items.wires.wire import Wire
 from classy_blocks.mesh import Mesh
-from classy_blocks.types import ChopTakeType, DirectionType
+from classy_blocks.optimize.grid import HexGrid
+from classy_blocks.types import ChopTakeType, DirectionType, OrientType
+from classy_blocks.util.constants import FACE_MAP
 
 
 @functools.lru_cache(maxsize=3000)  # that's for 1000 blocks
@@ -18,6 +21,20 @@ def get_block_from_axis(mesh: Mesh, axis: Axis) -> Block:
     raise RuntimeError("Block for Axis not found!")
 
 
+@functools.lru_cache(maxsize=2)
+def get_defined_wall_vertices(mesh: Mesh) -> Set[Vertex]:
+    """Returns vertices that are on the 'wall' patches"""
+    wall_vertices: set[Vertex] = set()
+
+    # explicitly defined walls
+    for patch in mesh.patches:
+        if patch.kind == "wall":
+            for side in patch.sides:
+                wall_vertices.update(set(side.vertices))
+
+    return wall_vertices
+
+
 class Instruction:
     """A descriptor that tells in which direction the specific block can be chopped."""
 
@@ -151,10 +168,46 @@ class Probe:
 
     def __init__(self, mesh: Mesh):
         self.mesh = mesh
+
+        # maps blocks to rows
         self.catalogue = Catalogue(self.mesh)
 
+        # finds blocks' neighbours
+        self.grid = HexGrid.from_mesh(self.mesh)
+
     def get_row_blocks(self, block: Block, direction: DirectionType) -> List[Block]:
         return self.catalogue.get_row_blocks(block, direction)
 
     def get_rows(self, direction: DirectionType) -> List[Row]:
         return self.catalogue.rows[direction]
+
+    def get_explicit_wall_vertices(self, block: Block) -> Set[Vertex]:
+        """Returns vertices from a block that lie on explicitly defined wall patches"""
+        mesh_vertices = get_defined_wall_vertices(self.mesh)
+        block_vertices = set(block.vertices)
+
+        return block_vertices.intersection(mesh_vertices)
+
+    def get_default_wall_vertices(self, block: Block) -> Set[Vertex]:
+        """Returns vertices that lie on default 'wall' patch"""
+        wall_vertices: Set[Vertex] = set()
+
+        # other sides when mesh has a default wall patch
+        if self.mesh.patch_list.default["kind"] == "wall":
+            # find block boundaries
+            block_index = self.mesh.blocks.index(block)
+            cell = self.grid.cells[block_index]
+
+            # sides with no neighbours are on boundary
+            boundaries: List[OrientType] = [
+                orient for orient, neighbours in cell.neighbours.items() if neighbours is None
+            ]
+
+            for orient in boundaries:
+                wall_vertices.union({block.vertices[i] for i in FACE_MAP[orient]})
+
+        return wall_vertices
+
+    def get_wall_vertices(self, block: Block) -> Set[Vertex]:
+        """Returns vertices that are on the 'wall' patches"""
+        return self.get_explicit_wall_vertices(block).union(self.get_default_wall_vertices(block))

src/classy_blocks/mesh.py

@@ -8,6 +8,7 @@
 from classy_blocks.construct.shape import Shape
 from classy_blocks.construct.stack import Stack
 from classy_blocks.items.block import Block
+from classy_blocks.items.patch import Patch
 from classy_blocks.items.vertex import Vertex
 from classy_blocks.lists.block_list import BlockList
 from classy_blocks.lists.edge_list import EdgeList
@@ -243,6 +244,10 @@ def is_assembled(self) -> bool:
     def vertices(self) -> List[Vertex]:
         return self.vertex_list.vertices
 
+    @property
+    def patches(self) -> List[Patch]:
+        return list(self.patch_list.patches.values())
+
     @property
     def operations(self) -> List[Operation]:
         """Returns a list of operations from all entities in depot"""

src/classy_blocks/optimize/cell.py

@@ -164,12 +164,19 @@ def q_scale(base, exponent, factor, value):
             quality += np.sum(q_scale(3, 2.5, 3, aspect_factor))
 
         except RuntimeWarning:
-            raise ValueError("Degenerate Cell") from RuntimeWarning
+            raise ValueError(f"Degenerate Cell: {self}") from RuntimeWarning
         finally:
             warnings.resetwarnings()
 
         return quality
 
+    @property
+    def min_length(self) -> float:
+        return min(self.get_edge_lengths())
+
+    def __str__(self):
+        return "-".join([str(index) for index in self.indexes])
+
 
 class QuadCell(CellBase):
     # Like constants.FACE_MAP but for quadrangle sides as line segments

src/classy_blocks/optimize/grid.py

@@ -1,14 +1,20 @@
-from typing import List, Type
+from typing import List, Type, Union
 
 import numpy as np
 
 from classy_blocks.base.exceptions import InvalidLinkError, NoJunctionError
+from classy_blocks.construct.assemblies.assembly import Assembly
+from classy_blocks.construct.flat.sketch import Sketch
 from classy_blocks.construct.flat.sketches.mapped import MappedSketch
+from classy_blocks.construct.operations.operation import Operation
+from classy_blocks.construct.shape import Shape
+from classy_blocks.construct.stack import Stack
 from classy_blocks.mesh import Mesh
 from classy_blocks.optimize.cell import CellBase, HexCell, QuadCell
 from classy_blocks.optimize.clamps.clamp import ClampBase
 from classy_blocks.optimize.junction import Junction
 from classy_blocks.optimize.links import LinkBase
+from classy_blocks.optimize.mapper import Mapper
 from classy_blocks.types import IndexType, NPPointListType, NPPointType
 from classy_blocks.util import functions as f
 from classy_blocks.util.constants import TOL
@@ -124,16 +130,41 @@ class QuadGrid(GridBase):
     cell_class = QuadCell
 
     @classmethod
-    def from_sketch(cls, sketch: MappedSketch) -> "QuadGrid":
-        # TODO: make grids from ANY sketch
-        return QuadGrid(sketch.positions, sketch.indexes)
+    def from_sketch(cls, sketch: Sketch) -> "QuadGrid":
+        if isinstance(sketch, MappedSketch):
+            # Use the mapper's indexes (provided by the user!)
+            return cls(sketch.positions, sketch.indexes)
+
+        # automatically create a mapping for arbitrary sketches
+        mapper = Mapper()
+        for face in sketch.faces:
+            mapper.add(face)
+
+        return cls(np.array(mapper.points), mapper.indexes)
 
 
 class HexGrid(GridBase):
     cell_class = HexCell
 
+    @classmethod
+    def from_elements(cls, elements: List[Union[Operation, Shape, Stack, Assembly]]) -> "HexGrid":
+        """Creates a grid from a list of elements"""
+        mapper = Mapper()
+
+        for element in elements:
+            if isinstance(element, Operation):
+                operations = [element]
+            else:
+                operations = element.operations
+
+            for operation in operations:
+                mapper.add(operation)
+
+        return cls(np.array(mapper.points), mapper.indexes)
+
     @classmethod
     def from_mesh(cls, mesh: Mesh) -> "HexGrid":
+        """Creates a grid from an assembled Mesh object"""
         points = np.array([vertex.position for vertex in mesh.vertices])
         addresses = [block.indexes for block in mesh.blocks]
 

src/classy_blocks/optimize/mapper.py

@@ -0,0 +1,54 @@
+from typing import List, Union
+
+from classy_blocks.construct.flat.face import Face
+from classy_blocks.construct.operations.operation import Operation
+from classy_blocks.types import IndexType, NPPointType
+from classy_blocks.util import functions as f
+from classy_blocks.util.constants import TOL
+
+ElementType = Union[Face, Operation]
+
+
+class Mapper:
+    """A helper that constructs mapped sketches/shapes
+    from arbitrary collection of faces/operations"""
+
+    def __init__(self) -> None:
+        self.points: List[NPPointType] = []
+        self.indexes: List[IndexType] = []
+        self.elements: List[Union[Face, Operation]] = []
+
+    def _add_point(self, point: NPPointType) -> int:
+        # TODO: this code is repeated several times all over;
+        # consolidate, unify, agglomerate, amass
+        # (especially in case one would need to utilize an octree or something)
+        for i, position in enumerate(self.points):
+            if f.norm(point - position) < TOL:
+                # reuse an existing point
+                index = i
+                break
+        else:
+            # no point found, create a new one
+            index = len(self.points)
+            self.points.append(point)
+
+        return index
+
+    def add(self, element: ElementType) -> None:
+        """Add Face's or Operation's points to the map"""
+        indexes = [self._add_point(point) for point in element.point_array]
+        self.indexes.append(indexes)
+        self.elements.append(element)
+
+    @classmethod
+    def from_map(cls, points: List[NPPointType], indexes: List[IndexType], elements: List[ElementType]) -> "Mapper":
+        """Creates a ready-made mapper from a sketch/shape that already has points/indexes defined"""
+        if len(indexes) != len(elements):
+            raise ValueError("Number of indexes and elements don't match!")
+
+        mapper = cls()
+        mapper.points = points
+        mapper.indexes = indexes
+        mapper.elements = elements
+
+        return mapper