Merge with feature/shape_optimization

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/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

src/classy_blocks/optimize/optimizer.py

@@ -1,18 +1,20 @@
 import abc
 import copy
 import time
-from typing import Literal
+from typing import List, Literal
 
 import numpy as np
 import scipy.optimize
 
 from classy_blocks.construct.flat.sketches.mapped import MappedSketch
+from classy_blocks.construct.operations.operation import Operation
 from classy_blocks.mesh import Mesh
 from classy_blocks.optimize.clamps.clamp import ClampBase
 from classy_blocks.optimize.clamps.surface import PlaneClamp
 from classy_blocks.optimize.grid import GridBase, HexGrid, QuadGrid
 from classy_blocks.optimize.iteration import ClampOptimizationData, IterationDriver
 from classy_blocks.optimize.links import LinkBase
+from classy_blocks.optimize.mapper import Mapper
 from classy_blocks.util.constants import TOL
 
 MinimizationMethodType = Literal["SLSQP", "L-BFGS-B", "Nelder-Mead", "Powell"]
@@ -144,6 +146,26 @@ def backport(self):
             self.mesh.vertices[i].move_to(point)
 
 
+class ShapeOptimizer(OptimizerBase):
+    def __init__(self, operations: List[Operation], report: bool = True):
+        self.mapper = Mapper()
+
+        for operation in operations:
+            self.mapper.add(operation)
+
+        grid = HexGrid(np.array(self.mapper.points), self.mapper.indexes)
+
+        super().__init__(grid, report)
+
+    def backport(self) -> None:
+        # Move every point of every operation to wherever it is now
+        for iop, indexes in enumerate(self.mapper.indexes):
+            operation = self.mapper.elements[iop]
+
+            for ipnt, i in enumerate(indexes):
+                operation.points[ipnt].move_to(self.grid.points[i])
+
+
 class SketchOptimizer(OptimizerBase):
     def __init__(self, sketch: MappedSketch, report: bool = True):
         self.sketch = sketch

tests/test_optimize/test_mapper.py

@@ -0,0 +1,38 @@
+import unittest
+
+from classy_blocks.construct.flat.face import Face
+from classy_blocks.construct.operations.box import Box
+from classy_blocks.optimize.mapper import Mapper
+
+
+class MapperTests(unittest.TestCase):
+    def test_single_face(self):
+        mapper = Mapper()
+        mapper.add(Face([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]))
+
+        self.assertListEqual(mapper.indexes, [[0, 1, 2, 3]])
+
+    def test_two_faces(self):
+        mapper = Mapper()
+        face = Face([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]])
+        mapper.add(face)
+        mapper.add(face.copy().translate([1, 0, 0]))
+
+        self.assertListEqual(mapper.indexes, [[0, 1, 2, 3], [1, 4, 5, 2]])
+
+    def test_single_operation(self):
+        mapper = Mapper()
+        box = Box([0, 0, 0], [1, 1, 1])
+
+        mapper.add(box)
+
+        self.assertListEqual(mapper.indexes, [[0, 1, 2, 3, 4, 5, 6, 7]])
+
+    def test_two_operations(self):
+        mapper = Mapper()
+        box = Box([0, 0, 0], [1, 1, 1])
+
+        mapper.add(box)
+        mapper.add(box.copy().translate([1, 0, 0]))
+
+        self.assertListEqual(mapper.indexes, [[0, 1, 2, 3, 4, 5, 6, 7], [1, 8, 9, 2, 5, 10, 11, 6]])