diff --git a/README.md b/README.md index f8871c4..cb62e9d 100644 --- a/README.md +++ b/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 diff --git a/examples/optimization/duct.py b/examples/optimization/duct.py new file mode 100644 index 0000000..120b010 --- /dev/null +++ b/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") diff --git a/src/classy_blocks/construct/curves/interpolated.py b/src/classy_blocks/construct/curves/interpolated.py index e7eba40..51dfd61 100644 --- a/src/classy_blocks/construct/curves/interpolated.py +++ b/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 diff --git a/src/classy_blocks/construct/curves/interpolators.py b/src/classy_blocks/construct/curves/interpolators.py index c25e266..81ce759 100644 --- a/src/classy_blocks/construct/curves/interpolators.py +++ b/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)) diff --git a/src/classy_blocks/construct/flat/face.py b/src/classy_blocks/construct/flat/face.py index b35a262..9581c4b 100644 --- a/src/classy_blocks/construct/flat/face.py +++ b/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) diff --git a/src/classy_blocks/construct/flat/sketches/spline_round.py b/src/classy_blocks/construct/flat/sketches/spline_round.py index 61b91b5..e02967f 100644 --- a/src/classy_blocks/construct/flat/sketches/spline_round.py +++ b/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 diff --git a/src/classy_blocks/grading/autograding/probe.py b/src/classy_blocks/grading/autograding/probe.py index a115672..b1979f9 100644 --- a/src/classy_blocks/grading/autograding/probe.py +++ b/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)) diff --git a/src/classy_blocks/mesh.py b/src/classy_blocks/mesh.py index f967094..489a2ab 100644 --- a/src/classy_blocks/mesh.py +++ b/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""" diff --git a/src/classy_blocks/optimize/cell.py b/src/classy_blocks/optimize/cell.py index 988fbc9..4e07c58 100644 --- a/src/classy_blocks/optimize/cell.py +++ b/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 diff --git a/src/classy_blocks/optimize/grid.py b/src/classy_blocks/optimize/grid.py index e5022a1..70fa84d 100644 --- a/src/classy_blocks/optimize/grid.py +++ b/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] diff --git a/src/classy_blocks/optimize/mapper.py b/src/classy_blocks/optimize/mapper.py new file mode 100644 index 0000000..76db6db --- /dev/null +++ b/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 diff --git a/src/classy_blocks/optimize/optimizer.py b/src/classy_blocks/optimize/optimizer.py index a1e81e0..1dd4a76 100644 --- a/src/classy_blocks/optimize/optimizer.py +++ b/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 diff --git a/tests/test_construct/test_curves/test_interpolated.py b/tests/test_construct/test_curves/test_interpolated.py index 4514830..b6c9ece 100644 --- a/tests/test_construct/test_curves/test_interpolated.py +++ b/tests/test_construct/test_curves/test_interpolated.py @@ -69,7 +69,7 @@ def test_param_at_length(self): self.assertAlmostEqual(self.curve.get_param_at_length(1), 0.25) def test_shear(self): - curve = self.curve + curve = LinearInterpolatedCurve(self.points, equalize=False) curve.shear([0, 1, 0], [0, 0, 0], [1, 0, 0], np.pi / 4) diff --git a/tests/test_construct/test_shell.py b/tests/test_construct/test_shell.py index d6ec679..83ed161 100644 --- a/tests/test_construct/test_shell.py +++ b/tests/test_construct/test_shell.py @@ -263,7 +263,7 @@ def test_chop(self): self.assertEqual(len(shell.operations[0].chops[2]), 1) - def test_set_outer_patch(self): + def test_set_outer_patch(self) -> None: orients: List[OrientType] = ["front", "right"] shell = self.get_shell(orients) shell.set_outer_patch("roof") diff --git a/tests/test_grading/test_probe.py b/tests/test_grading/test_probe.py index 24a6ab9..55fdc6b 100644 --- a/tests/test_grading/test_probe.py +++ b/tests/test_grading/test_probe.py @@ -1,9 +1,11 @@ -from typing import get_args +from typing import Set, get_args from parameterized import parameterized from classy_blocks.grading.autograding.probe import Probe, get_block_from_axis +from classy_blocks.items.vertex import Vertex from classy_blocks.mesh import Mesh +from classy_blocks.modify.find.shape import RoundSolidFinder from classy_blocks.types import DirectionType from tests.test_grading.test_autograde import AutogradeTestsBase @@ -98,3 +100,64 @@ def test_get_blocks_cylinder(self, axis, row, blocks): indexes.add(block.index) self.assertSetEqual(indexes, blocks) + + def test_wall_vertices_defined(self) -> None: + """Catch wall vertices from explicitly defined wall patches""" + cylinder = self.get_cylinder() + cylinder.set_outer_patch("outer") + + self.mesh.add(cylinder) + self.mesh.modify_patch("outer", "wall") + self.mesh.assemble() + + probe = Probe(self.mesh) + + finder = RoundSolidFinder(self.mesh, cylinder) + shell_vertices = finder.find_shell(True).union(finder.find_shell(False)) + wall_vertices: Set[Vertex] = set() + + for block in self.mesh.blocks: + wall_vertices.update(probe.get_explicit_wall_vertices(block)) + + self.assertSetEqual(shell_vertices, wall_vertices) + + def test_wall_vertices_default(self) -> None: + """Catch wall vertices from default patch""" + cylinder = self.get_cylinder() + cylinder.set_start_patch("inlet") + cylinder.set_end_patch("outlet") + + self.mesh.set_default_patch("outer", "wall") + self.mesh.assemble() + + probe = Probe(self.mesh) + + finder = RoundSolidFinder(self.mesh, cylinder) + shell_vertices = finder.find_shell(True).union(finder.find_shell(False)) + wall_vertices: Set[Vertex] = set() + + for block in self.mesh.blocks: + wall_vertices.update(probe.get_default_wall_vertices(block)) + + self.assertSetEqual(shell_vertices, wall_vertices) + + def test_wall_vertices_combined(self) -> None: + cylinder = self.get_cylinder() + cylinder.set_end_patch("outlet") + + cylinder.set_start_patch("bottom") + self.mesh.modify_patch("bottom", "wall") + + self.mesh.set_default_patch("outer", "wall") + self.mesh.assemble() + + probe = Probe(self.mesh) + + finder = RoundSolidFinder(self.mesh, cylinder) + shell_vertices = finder.find_shell(True).union(finder.find_shell(False)).union(finder.find_core(False)) + wall_vertices: Set[Vertex] = set() + + for block in self.mesh.blocks: + wall_vertices.update(probe.get_default_wall_vertices(block)) + + self.assertSetEqual(shell_vertices, wall_vertices) diff --git a/tests/test_items/test_wire.py b/tests/test_items/test_wire.py index aea1945..cb6947b 100644 --- a/tests/test_items/test_wire.py +++ b/tests/test_items/test_wire.py @@ -13,7 +13,7 @@ class WireTests(DataTestCase): """Tests of Pair object""" - def setUp(self): + def setUp(self) -> None: block = self.get_single_data(0) self.vertices = [Vertex(block.points[i], i) for i in range(8)] diff --git a/tests/test_optimize/test_mapper.py b/tests/test_optimize/test_mapper.py new file mode 100644 index 0000000..a83cd9c --- /dev/null +++ b/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]])