Add edge grading with Chop.preserve keyword

Currently all block gradings are specified as `edgeGrading` (TODO:
switch only when necessary).

- Add an example to demonstrate the 'preserve' keyword
- Remove 'take' keyword - always calculate cell count from average block
  size in specified axis (TODO: add an option to grade a single edge)
- Add keyword hints to Operation.chop() method

Author: FranzBangar

CHANGELOG.md

@@ -4,14 +4,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-# [1.6.4]
+## Unreleased
+### Added
+- **Edge grading**
+
+## [1.6.4]
 ### Added
 - MappedSketch.merge() method
 
 ### Changed
 - Improved blocking in cyclone example
 
-# [1.6.3]
+## [1.6.3]
 ### Changed
 - Bugfix: sorting of cells by sensitivity
 - Bugfix: sensitivity calculation moved clamps around
@@ -21,7 +25,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Removed
 - Quality caching (produced invalid results and did not speed up optimization)
 
-# [1.6.2]
+## [1.6.2]
 ### Added
 - QuarterSplineRound and QuarterSplineRoundRing; cylinders and rings with an arbitrary, parametrized spline cross-section
 - Lofts and Shapes, created with spline side-edges (when multiple sketches for mid-sections are provided)
@@ -30,15 +34,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Bugfixes
 
-# [1.6.1]
+## [1.6.1]
 ### Added
 - Assemblies and *Joints
   - Examples thereof
 
 ### Changed
 - Some bugfixes
 
-# [1.6.0]
+## [1.6.0]
 ### Added
 - Array element; handling multiple points at once (makes transforms faster)
 - Complete overhaul of Optimization:
@@ -55,14 +59,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Removed
 - QuadMap is no longer needed (handled by Grid classes)
 
-# [1.5.3]
+## [1.5.3]
 ### Added
 - RoundSolidShape.remove_inner_edges() can now remove edges from a specific face (start, end or both)
 
 ### Changed
 - Bugfix: invalid cells in round shapes (due to wrong spline point calculation)
 
-# [1.5.2]
+## [1.5.2]
 ### Added
 - Spline edges on QuarterDisk, HalfDisk and Disk (FourCoreDisk) to improve mesh quality
 - Face.remove_edges() to reset all edges to simple lines
@@ -71,7 +75,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Updated affected tutorials (diffusers, cyclone)
 
-# [1.5.1]
+## [1.5.1]
 ### Added
 - Optimization:
   - Choice of solver (different problems require - work best - with different solvers)
@@ -86,7 +90,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Improved optimization speed
 
-# [1.5.0]
+## [1.5.0]
 ### Added
 - Curves:
   - get_tangent()
@@ -114,15 +118,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Calling .transform() with a Mirror transformation will warn about creating an inverted block
 
 
-# [1.4.1]
+## [1.4.1]
 ### Changed
 - Improved optimization output
 
 ### Removed
 - Relaxation within optimization
 - Numerical integration of analytic curve lengths; use discretization instead
 
-# [1.4.0]
+## [1.4.0]
 ### Added
 - Channel example
 - Cyclone example
@@ -157,7 +161,7 @@ that parameter is of a similar magnitudes than points' coordinates.
 ### Removed
 - Junction.delta() is now handled by optimization automatically
 
-# [1.3.3]
+## [1.3.3]
 ### Added
 - Airfoil example
 - Translation and Rotational Link: move together with vertices being optimized (see the airfoil example)
@@ -172,19 +176,19 @@ that parameter is of a similar magnitudes than points' coordinates.
 ### Removed
 - Curve.get_closest_param() now finds initial_param automatically
 
-# [1.3.2] Curves
+## [1.3.2] Curves
 ### Added
 - *Curve objects for dealing with edge specification and optimization
 
-# [1.3.1] Optimization/Backport
+## [1.3.1] Optimization/Backport
 ### Added
 - mesh.clear() removes lists of all items that were populated during mesh.assemble()
 - mesh.backport() updates user supplied operations' points with results of optimization/modification
 
 ### Changed
 - Optimizer: under-relaxation for the first optimization iterations
 
-# [1.3.0] Blocking Optimization
+## [1.3.0] Blocking Optimization
 ### Added
 - **Blocking Optimization**
   - Finders for easier fetching vertices, generated by mesh.assemble()
@@ -208,7 +212,7 @@ that parameter is of a similar magnitudes than points' coordinates.
   - Calling .project_side() on an Operation will add the new geometry to edges instead of replacing them (but will replace existing label for the side)
 
 
-# [1.2.0] Shell Shape
+## [1.2.0] Shell Shape
 ### Added
 - A Shell Shape, created from arbitrary set of faces
 ### Changed

examples/advanced/edge_grading.py

@@ -0,0 +1,35 @@
+import os
+
+import numpy as np
+
+import classy_blocks as cb
+
+mesh = cb.Mesh()
+
+start = cb.Box([0, 0, 0], [1, 1, 0.1])
+
+start.chop(0, start_size=0.1)
+start.chop(1, length_ratio=0.5, start_size=0.01, c2c_expansion=1.2, preserve="start_size")
+# start.chop(1, length_ratio=0.5, end_size=0.01, c2c_expansion=1 / 1.2, preserve="end_size")
+start.chop(2, count=1)
+mesh.add(start)
+
+expand_start = start.get_face("right")
+expand = cb.Loft(expand_start, expand_start.copy().translate([1, 0, 0]).scale(2))
+expand.chop(2, start_size=0.1)
+mesh.add(expand)
+
+contract_start = expand.get_face("top")
+contract = cb.Loft(contract_start, contract_start.copy().translate([1, 0, 0]).scale(0.25))
+contract.chop(2, start_size=0.1)
+mesh.add(contract)
+
+end = cb.Extrude(contract.get_face("top"), 1)
+end.rotate(np.pi, [0, 0, 1])
+end.chop(2, start_size=0.1)
+
+mesh.add(end)
+
+mesh.set_default_patch("walls", "wall")
+
+mesh.write(os.path.join("..", "case", "system", "blockMeshDict"))

examples/complex/airfoil/airfoil.py

@@ -135,11 +135,12 @@ def get_loft(indexes):
 # Keep in mind that not all blocks need exact specification as
 # chopping will propagate automatically through blocking
 lofts[0].chop(2, count=1)  # 1 cell in the 3rd dimension
-lofts[1].chop(1, start_size=BL_THICKNESS, c2c_expansion=C2C_EXPANSION, take="max")
-lofts[8].chop(1, start_size=CELL_SIZE, take="max")
+# keep consistent first cell thickness by using edge grading
+lofts[1].chop(1, start_size=BL_THICKNESS, c2c_expansion=C2C_EXPANSION, preserve="start_size")
+lofts[8].chop(1, start_size=CELL_SIZE)
 
 for i in (0, 1, 2, 3, 4, 5, 6):
-    lofts[i].chop(0, start_size=CELL_SIZE, take="max")
+    lofts[i].chop(0, start_size=CELL_SIZE)
 
 for loft in lofts:
     mesh.add(loft)
@@ -195,7 +196,9 @@ def optimize_along_line(point_index, line_index_1, line_index_2):
 optimize_along_line(5, 3, 6)
 
 if OPTIMIZE:
-    optimizer.optimize(tolerance=1e-7, method="SLSQP")
+    optimizer.optimize(tolerance=1e-3, method="SLSQP")
+    mesh.backport()
+    mesh.clear()
 
 ### Write the mesh
 mesh.modify_patch("topAndBottom", "empty")

examples/complex/heater/heater.py

@@ -43,7 +43,7 @@ def set_cell_zones(shape: cb.Shape):
 
 # The curved part of heater (and fluid around it); constructed from 4 revolves
 heater_arch = cb.RevolvedStack(straight_1.sketch_2, np.pi, [0, 0, 1], [0, 0, 0], 4)
-heater_arch.chop(start_size=p.solid_cell_size, take="min")
+heater_arch.chop(start_size=p.solid_cell_size)
 for shape in heater_arch.shapes:
     set_cell_zones(shape)
 mesh.add(heater_arch)

showcase/classy_classic_grading.png

@@ -18,6 +18,10 @@ class UndefinedGradingsError(Exception):
     define all blocks in the mesh"""
 
 
+class InconsistentGradingsError(Exception):
+    """Raised when cell counts for edges on the same axis is not consistent"""
+
+
 class ShapeCreationError(Exception):
     """Base class for shape creation errors (invalid parameters/types to
     shape constructors)"""

showcase/classy_edges_grading.png

@@ -1,7 +1,8 @@
 import warnings
-from typing import Dict, List, Optional, TypeVar, Union, get_args
+from typing import Dict, List, Optional, Union, get_args
 
 import numpy as np
+from typing_extensions import Unpack
 
 from classy_blocks.base.element import ElementBase
 from classy_blocks.base.exceptions import EdgeCreationError
@@ -10,15 +11,13 @@
 from classy_blocks.construct.flat.face import Face
 from classy_blocks.construct.point import Point
 from classy_blocks.grading.chop import Chop
-from classy_blocks.types import AxisType, NPPointType, OrientType, PointType, ProjectToType, VectorType
+from classy_blocks.types import AxisType, ChopArgs, NPPointType, OrientType, PointType, ProjectToType, VectorType
 from classy_blocks.util import constants
 from classy_blocks.util import functions as f
 from classy_blocks.util.constants import SIDES_MAP
 from classy_blocks.util.frame import Frame
 from classy_blocks.util.tools import edge_map
 
-OperationT = TypeVar("OperationT", bound="Operation")
-
 
 class Operation(ElementBase):
     """A base class for all single-block operations
@@ -59,7 +58,7 @@ def add_side_edge(self, corner_idx: int, edge_data: EdgeData) -> None:
 
         self.side_edges[corner_idx] = edge_data
 
-    def chop(self, axis: AxisType, **kwargs) -> None:
+    def chop(self, axis: AxisType, **kwargs: Unpack[ChopArgs]) -> None:
         """Chop the operation (count/grading) in given axis:
         0: along first edge of a face
         1: along second edge of a face

src/classy_blocks/base/exceptions.py

@@ -1,9 +1,9 @@
 import dataclasses
 from functools import lru_cache
-from typing import Callable, List, Optional, Set, Tuple, Union
+from typing import Callable, Dict, List, Optional, Set, Tuple, Union
 
 from classy_blocks.grading import relations as rel
-from classy_blocks.types import ChopTakeType
+from classy_blocks.types import ChopPreserveType
 
 
 @dataclasses.dataclass
@@ -58,9 +58,9 @@ class Chop:
     end_size: Optional[float] = None
     total_expansion: Optional[float] = None
     invert: bool = False
-    take: ChopTakeType = "avg"
+    preserve: ChopPreserveType = "c2c_expansion"
 
-    def __post_init__(self):
+    def __post_init__(self) -> None:
         # default: take c2c_expansion=1 if there's less than 2 parameters given
         grading_params = [self.start_size, self.end_size, self.count, self.total_expansion, self.c2c_expansion]
         if len(grading_params) - grading_params.count(None) < 2:
@@ -71,25 +71,31 @@ def __post_init__(self):
         if self.count is not None:
             self.count = max(int(self.count), 1)
 
+        # stored results from self.calculate() method;
+        # this will be used for creating new chops taking self.preserve
+        # into account. User input (chop arguments) should not be overridden because
+        # if mesh is modified (optimization etc.) all numbers will be wrong.
+        self.results: Dict[Union[ChopPreserveType, str], Union[float, int]] = dict()
+
     def calculate(self, length: float) -> Tuple[int, float]:
         """Calculates cell count and total expansion ratio for this chop
         by calling functions that take known variables and return new values"""
         data = dataclasses.asdict(self)
+        self.results = data
         calculated: Set[str] = set()
 
-        for key in data.keys():
+        for key in self.results.keys():
             if data[key] is not None:
                 calculated.add(key)
 
-        for _ in range(20):
-            if {"count", "total_expansion"}.issubset(calculated):
-                count = int(data["count"])
-                total_expansion = data["total_expansion"]
+        for _ in range(12):
+            if {"count", "total_expansion", "c2c_expansion", "start_size", "end_size"}.issubset(calculated):
+                self.results["count"] = int(self.results["count"])
 
                 if self.invert:
-                    return count, 1 / total_expansion
+                    data["total_expansion"] = 1 / data["total_expansion"]
 
-                return count, total_expansion
+                return data["count"], data["total_expansion"]
 
             for chop_rel in ChopRelation.get_possible_combinations():
                 output = chop_rel.output
@@ -106,3 +112,32 @@ def calculate(self, length: float) -> Tuple[int, float]:
                     calculated.add(output)
 
         raise ValueError(f"Could not calculate count and grading for given parameters: {data}")
+
+    def copy_preserving(self, inverted: bool = False) -> "Chop":
+        """Creates a copy of this Chop with equal count but
+        sets other parameters from current data so that
+        the correct start/end size or c2c is maintained"""
+        args = dataclasses.asdict(self)
+        args["count"] = self.results["count"]
+
+        for arg in ["total_expansion", "c2c_expansion", "start_size", "end_size"]:
+            args[arg] = None
+
+        # TODO: get rid of if-ing
+        if self.preserve == "start_size":
+            if inverted:
+                args["end_size"] = self.results["start_size"]
+            else:
+                args["start_size"] = self.results["start_size"]
+        elif self.preserve == "end_size":
+            if inverted:
+                args["start_size"] = self.results["end_size"]
+            else:
+                args["end_size"] = self.results["end_size"]
+        else:
+            if inverted:
+                args["c2c_expansion"] = 1 / self.results["c2c_expansion"]
+            else:
+                args["c2c_expansion"] = self.results["c2c_expansion"]
+
+        return Chop(**args)

src/classy_blocks/construct/operations/operation.py

@@ -49,10 +49,10 @@ class Grading:
     def __init__(self, length: float):
         # "multi-grading" specification according to:
         # https://cfd.direct/openfoam/user-guide/v9-blockMesh/#multi-grading
-        self.specification: List[List] = []  # a list of lists [length ratio, count ratio, total expansion]
-
         self.length = length
 
+        self.specification: List[List] = []  # a list of lists [length ratio, count ratio, total expansion]
+
     def add_chop(self, chop: Chop) -> None:
         """Add a grading division to block specification.
         Use length_ratio for multigrading (see documentation).
@@ -101,10 +101,15 @@ def inverted(self) -> "Grading":
 
         return g_inv
 
+    @property
+    def counts(self) -> List[int]:
+        """Counts per chop"""
+        return [d[1] for d in self.specification]
+
     @property
     def count(self) -> int:
         """Return number of cells, summed over all sub-divisions"""
-        return sum(d[1] for d in self.specification)
+        return sum(self.counts)
 
     @property
     def is_defined(self) -> bool: