Savegame slot 2

examples/advanced/inflation_grader.py

@@ -32,10 +32,6 @@
     vertex = list(finder.find_in_sphere(point))[0]
     vertex.translate([0, 0.8, 0])
 
-
-# TODO: Hack! mesh.assemble() won't work here but wires et. al. must be updated
-mesh.block_list.update()
-
 grader = InflationGrader(mesh, 1e-2, 0.1)
 grader.grade(take="max")
 

src/classy_blocks/grading/autograding/params/distributor.py

@@ -0,0 +1,171 @@
+import abc
+import dataclasses
+from typing import List
+
+import numpy as np
+import scipy.interpolate
+import scipy.optimize
+
+from classy_blocks.grading.autograding.params.base import sum_length
+from classy_blocks.grading.autograding.params.layers import InflationLayer
+from classy_blocks.grading.chop import Chop
+from classy_blocks.types import FloatListType
+from classy_blocks.util.constants import TOL
+
+
+@dataclasses.dataclass
+class DistributorBase(abc.ABC):
+    """Algorithm that creates chops from given count and sizes;
+    first distributes cells using a predefined 'ideal' sizes and ratios,
+    then optimizes their individual size to get as close to that ideal as possible.
+
+    Then, when cells are placed, Chops are created so that they produce as similar
+    cells to the calculated as possible. Since blockMesh only supports equal
+    cell-to-cell expansion for each chop, currently creating a limited number of Chops
+    is the only way to go. In a possible future (a.k.a. own meshing script),
+    calculated cell sizes could be used directly for meshing."""
+
+    count: int
+    size_before: float
+    length: float
+    size_after: float
+
+    @staticmethod
+    def get_actual_ratios(coords: FloatListType) -> FloatListType:
+        lengths = np.diff(coords)
+        return (lengths[:-1] / np.roll(lengths, -1)[:-1]) ** -1
+
+    @abc.abstractmethod
+    def get_ideal_ratios(self) -> FloatListType:
+        """Returns desired cell-to-cell ratios"""
+        return np.ones(self.count + 1)
+
+    @abc.abstractmethod
+    def get_ratio_weights(self) -> FloatListType:
+        """Returns weights of cell ratios"""
+
+    def get_raw_coords(self) -> FloatListType:
+        # 'count' denotes number of 'intervals' so there must be another point
+        return np.concatenate(
+            ([-self.size_before], np.linspace(0, self.length, num=self.count + 1), [self.length + self.size_after])
+        )
+
+    def get_smooth_coords(self) -> FloatListType:
+        coords = self.get_raw_coords()
+
+        # prepare a function for least squares minimization:
+        # f(coords) -> actual_ratios / ideal_ratios
+        # first two and last two points are fixed
+
+        def ratios(inner_coords):
+            coords[2:-2] = inner_coords
+
+            difference = self.get_actual_ratios(coords) - self.get_ideal_ratios()
+            return difference * self.get_ratio_weights()
+
+        scale = min(self.size_before, self.size_after, self.length / self.count) / 10
+        _ = scipy.optimize.least_squares(ratios, coords[2:-2], method="lm", ftol=scale / 100, x_scale=scale)
+
+        return coords
+
+    @property
+    def is_simple(self) -> bool:
+        # Don't overdo basic, simple-graded blocks
+        base_size = self.length / self.count
+
+        # TODO: use a more relaxed criterion?
+        return base_size - self.size_before < TOL and base_size - self.size_after < TOL
+
+    def get_chops(self, pieces: int) -> List[Chop]:
+        if self.is_simple:
+            return [Chop(count=self.count)]
+
+        coords = self.get_smooth_coords()
+        sizes = np.diff(coords)
+        ratios = self.get_actual_ratios(coords)
+
+        count = len(ratios) - 1
+
+        # create a piecewise linear function from a number of chosen indexes and their respective c2c ratios
+        # then optimize indexes to obtain best fit
+        # fratios = scipy.interpolate.interp1d(range(len(ratios)), ratios)
+
+        # def get_piecewise(indexes: List[int]) -> Callable:
+        #     values = np.take(ratios, indexes)
+        #     return scipy.interpolate.interp1d(indexes, values)
+
+        # def get_fitness(indexes: List[int]) -> float:
+        #     fitted = get_piecewise(indexes)(range(len(ratios)))
+
+        #     ss_tot = np.sum((ratios - np.mean(ratios)) ** 2)
+        #     ss_res = np.sum((ratios - fitted) ** 2)
+
+        #     return 1 - (ss_res / ss_tot)
+
+        # print(get_fitness([0, 9, 10, count]))
+        # print(get_fitness(np.linspace(0, count, num=pieces + 1, dtype=int)))
+
+        chops: List[Chop] = []
+        indexes = np.linspace(0, count, num=pieces + 1, dtype=int)
+
+        for i, index in enumerate(indexes[:-1]):
+            chop_ratios = ratios[index : indexes[i + 1]]
+            ratio = np.prod(chop_ratios)
+            chop_count = indexes[i + 1] - indexes[i]
+            avg_ratio = ratio ** (1 / chop_count)
+            length = sum_length(sizes[index], chop_count, avg_ratio)
+            chop = Chop(length_ratio=length, total_expansion=ratio, count=chop_count)
+            chops.append(chop)
+
+        # normalize length ratios
+        sum_ratio = sum([chop.length_ratio for chop in chops])
+        for chop in chops:
+            chop.length_ratio = chop.length_ratio / sum_ratio
+
+        return chops
+
+
+@dataclasses.dataclass
+class SmoothDistributor(DistributorBase):
+    def get_ideal_ratios(self):
+        # In a smooth grader, we want all cells to be as similar to their neighbours as possible
+        return super().get_ideal_ratios()
+
+    def get_ratio_weights(self):
+        # Enforce stricter policy on size_before and size_after
+        weights = np.ones(self.count + 1)
+        for i in (0, 1, 2, 3):
+            w = 2 ** (3 - i)
+            weights[i] = w
+            weights[-i - 1] = w
+
+        return weights
+
+
+@dataclasses.dataclass
+class InflationDistributor(SmoothDistributor):
+    c2c_expansion: float
+    bl_thickness_factor: int
+
+    @property
+    def is_simple(self) -> bool:
+        return False
+
+    def get_ideal_ratios(self):
+        # Ideal growth ratio in boundary layer is user-specified c2c_expansion;
+        inflation_layer = InflationLayer(self.size_before, self.c2c_expansion, self.bl_thickness_factor, 1e12)
+        inflation_count = inflation_layer.count
+        print(f"Inflation count: {inflation_count}")
+
+        ratios = super().get_ideal_ratios()
+
+        ratios[:inflation_count] = self.c2c_expansion
+        print(ratios)
+
+        return ratios
+
+    def get_ratio_weights(self):
+        return super().get_ratio_weights()
+
+    def _get_ratio_weights(self):
+        return np.ones(self.count + 1)

src/classy_blocks/grading/autograding/params/inflation.py

@@ -1,146 +1,10 @@
-import abc
-from typing import List, Optional, Tuple
+from typing import List, Optional
 
-from classy_blocks.grading import relations as gr
+from classy_blocks.grading.autograding.params.distributor import InflationDistributor, SmoothDistributor
+from classy_blocks.grading.autograding.params.layers import BufferLayer, BulkLayer, InflationLayer, LayerStack
 from classy_blocks.grading.autograding.params.smooth import SmoothGraderParams
 from classy_blocks.grading.autograding.probe import WireInfo
 from classy_blocks.grading.chop import Chop
-from classy_blocks.util.constants import VBIG
-
-
-class Layer(abc.ABC):
-    """A common interface to all layers of a grading (inflation, buffer, bulk)"""
-
-    # Defines one chop and tools to handle it
-    start_size: float
-    c2c_expansion: float
-    length: float
-    count: int
-    end_size: float
-
-    def _construct(
-        self, length_limit: float = VBIG, count_limit: int = 10**12, size_limit: float = VBIG
-    ) -> Tuple[float, float, int]:
-        # stop construction of the layer when it hits any of the above limits
-        count = 1
-        size = self.start_size
-        length = self.start_size
-
-        for i in range(count_limit):
-            if length >= length_limit:
-                break
-
-            if count >= count_limit:
-                break
-
-            if size >= size_limit:
-                break
-
-            length += size
-            size *= self.c2c_expansion
-            count = i
-
-        return length, size, count
-
-    def __init__(self, length_limit: float = VBIG, count_limit: int = 10**12, size_limit: float = VBIG):
-        # stop construction of the layer when it hits any of the above limits
-        self.length, self.end_size, self.count = self._construct(length_limit, count_limit, size_limit)
-
-    def get_chop_count(self, total_count: int) -> int:
-        return min(self.count, total_count)
-
-    def get_chop(self, actual_length: float, remaining_count: int, invert: bool) -> Tuple[Chop, int]:
-        """Returns a Chop, adapter to a given length; count will not exceed remaining"""
-        # length ratios will be normalized later
-        length, end_size, count = self._construct(actual_length, remaining_count)
-        chop = Chop(length_ratio=length, end_size=end_size, count=count)
-
-        if invert:
-            chop.start_size = None
-            chop.end_size = self.start_size
-
-        return chop, count
-
-    def __repr__(self):
-        return f"{self.length}-{self.count}"
-
-
-class InflationLayer(Layer):
-    def __init__(self, wall_size: float, c2c_expansion: float, thickness_factor: int, _max_length: float):
-        self.start_size = wall_size
-        self.c2c_expansion = c2c_expansion
-
-        super().__init__(length_limit=thickness_factor * wall_size)
-
-
-class BufferLayer(Layer):
-    def __init__(self, start_size: float, c2c_expansion: float, bulk_size: float, _max_length: float):
-        self.start_size = start_size
-        self.c2c_expansion = c2c_expansion
-
-        super().__init__(size_limit=bulk_size)
-
-    def get_chop_count(self, total_count: int) -> int:
-        # use all remaining cells
-        return max(self.count, total_count)
-
-
-class BulkLayer(Layer):
-    # Must be able to adapt end_size to match size_after
-    def __init__(self, start_size: float, end_size: float, remainder: float):
-        self.start_size = start_size
-        self.end_size = end_size
-        self.length = remainder
-
-        total_expansion = gr.get_total_expansion__start_size__end_size(self.length, self.start_size, self.end_size)
-        self.count = gr.get_count__total_expansion__start_size(self.length, total_expansion, self.start_size)
-        self.c2c_expansion = gr.get_c2c_expansion__count__end_size(self.length, self.count, self.end_size)
-
-
-class LayerStack:
-    """A collection of one, two or three layers (chops) for InflationGrader"""
-
-    def __init__(self, length: float):
-        self.length = length
-        self.layers: List[Layer] = []
-
-    @property
-    def count(self) -> int:
-        return sum(layer.count for layer in self.layers)
-
-    @property
-    def remaining_length(self) -> float:
-        return self.length - sum(layer.length for layer in self.layers)
-
-    def add(self, layer: Layer) -> bool:
-        self.layers.append(layer)
-        return self.is_done
-
-    @property
-    def is_done(self) -> bool:
-        """Returns True if no more layers need to be added"""
-        if len(self.layers) == 0:
-            return False
-
-        if len(self.layers) == 3:
-            # nothing more to be added?
-            return True
-
-        return self.remaining_length <= self.layers[0].start_size
-
-    @property
-    def last_size(self) -> float:
-        return self.layers[-1].end_size
-
-    def get_chops(self, _total_count: int, _invert: bool) -> List[Chop]:
-        # normalize length_ratios
-        # ratios = [chop.length_ratio for chop in chops]
-
-        # for chop in chops:
-        #     chop.length_ratio = chop.length_ratio / sum(ratios)
-
-        # return chops
-        raise NotImplementedError
 
 
 class InflationGraderParams(SmoothGraderParams):
@@ -212,16 +76,39 @@ def is_squeezed(self, count: int, info: WireInfo) -> bool:
         # TODO: replace 0.9 with something less arbitrary (a better rule)
         return self.get_stack(info.length).last_size < 0.9 * self.bulk_cell_size
 
-    def get_chops(self, count, info: WireInfo) -> List[Chop]:
-        if not (info.starts_at_wall or info.ends_at_wall):
-            return super().get_chops(count, info)
-
-        stack = self.get_stack(info.length, info.size_after)
+    def get_squeezed_chops(self, count: int, info: WireInfo) -> List[Chop]:
+        return self.get_chops(count, info)
 
+    def get_chops(self, count, info: WireInfo) -> List[Chop]:
         if info.starts_at_wall and info.ends_at_wall:
             raise NotImplementedError
 
+        size_before = info.size_before
+        if size_before is None:
+            if info.starts_at_wall:
+                size_before = self.first_cell_size
+            else:
+                size_before = self.cell_size
+
+        size_after = info.size_after
+        if size_after is None:
+            if info.ends_at_wall:
+                size_after = self.first_cell_size
+            else:
+                size_after = self.cell_size
+
+        if not (info.starts_at_wall or info.ends_at_wall):
+            print("NOT AT WALL", info)
+            distributor = SmoothDistributor(count, size_before, info.length, size_after)
+        else:
+            print("AT WALL", info)
+            distributor = InflationDistributor(
+                count, size_before, info.length, size_after, self.c2c_expansion, self.bl_thickness_factor
+            )
+
+        chops = distributor.get_chops(3)
+
         if info.ends_at_wall:
-            return list(reversed(stack.get_chops(count, True)))
+            return list(reversed(chops))
 
-        return stack.get_chops(count, False)
+        return chops