Add chop smoothing via grading for HighRe grader

Need a better objective function because this doesn't help much.

Author: FranzBangar

examples/advanced/autograding_highre.py

@@ -0,0 +1,35 @@
+import os
+
+import classy_blocks as cb
+from classy_blocks.grading.autograding.grader import HighReGrader, SimpleGrader
+from classy_blocks.grading.autograding.params import HighReChopParams, SimpleChopParams
+
+mesh = cb.Mesh()
+
+base = cb.Grid([0, 0, 0], [3, 2, 0], 3, 2)
+
+shape = cb.ExtrudedShape(base, 1)
+mesh.add(shape)
+mesh.assemble()
+finder = cb.GeometricFinder(mesh)
+
+move_points = [[0, 1, 1], [2, 1, 1], [3, 1, 1]]
+
+for point in move_points:
+    vertex = list(finder.find_in_sphere(point))[0]
+    vertex.translate([0, 0.8, 0])
+
+# TODO! Un-hack
+mesh.block_list.update()
+
+mesh.set_default_patch("walls", "wall")
+
+params = HighReChopParams(0.1)
+grader = HighReGrader(mesh, params)
+grader.grade(take="max")
+
+params = SimpleChopParams(0.1)
+grader = SimpleGrader(mesh, params)
+# grader.grade()
+
+mesh.write(os.path.join("..", "case", "system", "blockMeshDict"), debug_path="debug.vtk")

src/classy_blocks/grading/autograding/grader.py

@@ -1,7 +1,8 @@
-from typing import get_args
+from typing import Set, get_args
 
 from classy_blocks.grading.autograding.params import ChopParams, FixedCountParams, HighReChopParams, SimpleChopParams
 from classy_blocks.grading.autograding.probe import Probe
+from classy_blocks.items.wires.wire import Wire
 from classy_blocks.mesh import Mesh
 from classy_blocks.types import ChopTakeType, DirectionType
 
@@ -14,13 +15,45 @@ def __init__(self, mesh: Mesh, params: ChopParams):
         self.mesh.assemble()
         self.probe = Probe(self.mesh)
 
+    def _get_end_size(self, wires: Set[Wire]):
+        """Returns average size of wires' last cell"""
+        if len(wires) == 0:
+            return 0
+
+        return sum(wire.grading.end_size for wire in wires) / len(wires)
+
+    def _get_start_size(self, wires: Set[Wire]):
+        """Returns average size of wires' first cell"""
+        if len(wires) == 0:
+            return 0
+
+        return sum(wire.grading.start_size for wire in wires) / len(wires)
+
     def grade_axis(self, axis: DirectionType, take: ChopTakeType) -> None:
         for row in self.probe.get_rows(axis):
-            length = row.get_length(take)
-            count = self.params.get_count(length)
+            # determine count
+            wires = row.get_wires()
+
+            for wire in wires:
+                if wire.is_defined:
+                    # there's a wire with a defined count already, use that
+                    count = wire.grading.count
+                    break
+            else:
+                # take length from a row, as requested
+                length = row.get_length(take)
+                count = self.params.get_count(length)
 
             for wire in row.get_wires():
-                chops = self.params.get_chops(count, wire.length)
+                # don't touch defined wires
+                if wire.is_defined:
+                    # TODO: test
+                    continue
+
+                size_before = self._get_end_size(wire.before)
+                size_after = self._get_start_size(wire.after)
+                chops = self.params.get_chops(count, wire.length, size_before, size_after)
+
                 for chop in chops:
                     wire.grading.add_chop(chop)
 

src/classy_blocks/grading/autograding/params.py

@@ -3,7 +3,9 @@
 import warnings
 from typing import List, Tuple
 
-from classy_blocks.grading import relations as gr
+import scipy.optimize
+
+import classy_blocks.grading.relations as gr
 from classy_blocks.grading.chop import Chop
 from classy_blocks.types import ChopTakeType
 
@@ -20,23 +22,13 @@ def sum_length(start_size: float, count: int, c2c_expansion: float) -> float:
     return length
 
 
-def sum_count(lengths: List[float], chops: List[Chop]):
-    count = 0
-
-    for i, chop in enumerate(chops):
-        length = lengths[i]
-        count += chop.calculate(length).count
-
-    return count
-
-
 class ChopParams(abc.ABC):
     @abc.abstractmethod
     def get_count(self, length: float) -> int:
         """Calculates count based on given length - used once only"""
 
     @abc.abstractmethod
-    def get_chops(self, count: int, length: float) -> List[Chop]:
+    def get_chops(self, count: int, length: float, size_before: float = 0, size_after: float = 0) -> List[Chop]:
         """Fixes cell count but modifies chops so that proper cell sizing will be obeyed"""
         # That depends on inherited classes' philosophy
 
@@ -48,7 +40,7 @@ class FixedCountParams(ChopParams):
     def get_count(self, _length):
         return self.count
 
-    def get_chops(self, count, _length) -> List[Chop]:
+    def get_chops(self, count, _length, _size_before=0, _size_after=0) -> List[Chop]:
         return [Chop(count=count)]
 
 
@@ -60,7 +52,7 @@ class SimpleChopParams(ChopParams):
     def get_count(self, length: float):
         return int(length / self.cell_size)
 
-    def get_chops(self, count, _length):
+    def get_chops(self, count, _length, _size_before=0, _size_after=0):
         return [Chop(count=count)]
 
 
@@ -69,17 +61,61 @@ class HighReChopParams(ChopParams):
     cell_size: float
 
     def get_count(self, length: float):
-        # the first chop defines the count; it's a very simple one
-        return int(length / self.cell_size)
-
-    def get_chops(self, count: int, _):
-        # TODO: adjust length ratio for smoothest transition in the middle of block
-        return [
-            Chop(length_ratio=0.5, count=count // 2, start_size=self.cell_size),
-            Chop(length_ratio=0.5, count=count // 2, end_size=self.cell_size),
+        # the first chop defines the count;
+        count = int(length / self.cell_size)
+        # it must be divisible by 2
+        if count % 2 != 0:
+            count += 1
+
+        return count
+
+    def get_chops(self, count, length, size_before=0, size_after=0):
+        # length of the wire that was used to set count
+        if size_before == 0:
+            size_before = self.cell_size
+        if size_after == 0:
+            size_after = self.cell_size
+
+        chops = [
+            Chop(count=count // 2),
+            Chop(count=count // 2),
         ]
 
+        def objfun(params):
+            chops[0].length_ratio = params[0]
+            chops[1].length_ratio = 1 - params[0]
+
+            chops[0].total_expansion = params[1]
+            chops[1].total_expansion = params[2]
+
+            data_1 = chops[0].calculate(length)
+            data_2 = chops[1].calculate(length)
 
+            ofstart = (size_before - data_1.start_size) ** 2
+            ofmid1 = (data_1.end_size - self.cell_size) ** 2
+            ofmid2 = (data_2.start_size - self.cell_size) ** 2
+            ofend = (data_2.end_size - size_after) ** 2
+
+            return max([ofstart, ofmid1, ofmid2, ofend])
+
+        initial = [0.5, 1, 1]
+        bounds = (
+            (0.1, 0.9),
+            (0.1, 10),
+            (0.1, 10),
+        )
+        result = scipy.optimize.minimize(objfun, initial, bounds=bounds).x
+
+        chops[0].length_ratio = result[0]
+        chops[1].length_ratio = 1 - result[0]
+
+        chops[0].total_expansion = result[1]
+        chops[1].total_expansion = result[2]
+
+        return chops
+
+
+# INVALID! Next on list
 @dataclasses.dataclass
 class LowReChopParams(ChopParams):
     """Parameters for mesh grading for Low-Re cases.
@@ -164,16 +200,17 @@ def get_count(self, length: float):
 
         if remaining_length <= 0:
             warnings.warn("Stopping chops at boundary layer (not enough space)!", stacklevel=1)
-            return sum_count([length], chops)
+            # return chops
+            return 0
 
         # buffer
         buffer, buffer_size = self._get_buffer_chop(last_bl_size)
         buffer.length_ratio = buffer_size / length
         chops.append(buffer)
         if buffer_size >= remaining_length:
             warnings.warn("Stopping chops at buffer layer (not enough space)!", stacklevel=1)
-
-            return sum_count([self.boundary_layer_thickness, buffer_size], chops)
+            # return chops
+            return 1
 
         # bulk
         remaining_length = remaining_length - buffer_size
@@ -182,7 +219,7 @@ def get_count(self, length: float):
         chops.append(bulk)
 
         # return chops
-        return sum_count([self.boundary_layer_thickness, buffer_size, remaining_length], chops)
+        return 1
 
-    def get_chops(self, count: int, length: float) -> List[Chop]:
+    def get_chops(self, count, length, size_before=0, size_after=0) -> List[Chop]:
         raise NotImplementedError("TODO!")

src/classy_blocks/mesh.py

@@ -258,4 +258,3 @@ def operations(self) -> List[Operation]:
     @property
     def blocks(self) -> List[Block]:
         return self.block_list.blocks
-        return self.block_list.blocks