Fix typing/format/test issues

examples/shape/cylinder.py

@@ -12,11 +12,12 @@
 axis_point_2 = [5.0, 5.0, 0.0]
 radius_point_1 = [0.0, 0.0, 2.0]
 
-cylinder = cb.Cylinder(axis_point_1, axis_point_2, radius_point_1)
+cylinder = cb.QuarterCylinder(axis_point_1, axis_point_2, radius_point_1)
 
 cylinder.set_start_patch("inlet")
 cylinder.set_end_patch("outlet")
 cylinder.set_outer_patch("walls")
+cylinder.set_symmetry_patch('sym')
 
 # if curved core edges get in the way (when moving vertices, optimization, ...),
 # remove them with this method:

examples/shape/frustum.py

@@ -1,11 +1,12 @@
 import os
+import numpy as np
 
 import classy_blocks as cb
 
 mesh = cb.Mesh()
 
 axis_point_1 = [0.0, 0.0, 0.0]
-axis_point_2 = [2.0, 2.0, 0.0]
+axis_point_2 = [2.0, 0, 0.0]
 radius_point_1 = [0.0, 0.0, 2.0]
 radius_2 = 0.5
 
@@ -20,14 +21,32 @@
 # because of sharp edges; in those cases it is better to use RevolvedRing combined with
 # Cylinder/Frustum with non-flat start/end faces.
 frustum = cb.Frustum(axis_point_1, axis_point_2, radius_point_1, radius_2, radius_mid=1.1)
+cylinder = cb.Cylinder.chain(frustum, 6, start_face=True)
 
-frustum.set_start_patch("inlet")
+pos = cylinder.sketch_1.positions
+pos[:9] += np.array([-0.3,0,0])
+cylinder.sketch_1.update(pos)
+cylinder.sketch_1.add_edges()
+
+pos = frustum.sketch_1.positions
+pos[:9] += np.array([-0.3,0,0])
+frustum.sketch_1.update(pos)
+frustum.sketch_1.add_edges()
+
+
+
+
+cylinder.set_end_patch("inlet")
 frustum.set_outer_patch("walls")
+cylinder.set_outer_patch("walls")
 frustum.set_end_patch("outlet")
 
+cylinder.chop_axial(count=90)
 frustum.chop_axial(count=30)
 frustum.chop_radial(start_size=core_size, end_size=bl_thickness)
 frustum.chop_tangential(start_size=core_size)
 
+mesh.add(cylinder)
 mesh.add(frustum)
 mesh.write(os.path.join("..", "case", "system", "blockMeshDict"), debug_path="debug.vtk")
+

examples/shape/splined/combined_example.py

@@ -33,16 +33,7 @@
                               corner_2_point=np.asarray([0, 0, 2.2]) + direction * 3.4,
                               side_1=0.3, side_2=0.3)
 
-# It is possible to define the oval rom the radius instead of the side lengths.
-# It is sometimes necessary to increase the number of spline points
-# in the transition from curved to straight for good mesh quality.
-# This is done using the n_straight_spline_points key word.
-"""
-oval_sketch_2 = cb.SplineDisk.init_from_radius(center_point=center_point + direction * 4,
-                                               corner_1_point=np.asarray([0, 1.5, 0]) + direction * 4,
-                                               corner_2_point=np.asarray([0, 0, 2.5]) + direction * 4,
-                                               r_1=0.7, r_2=0.7, n_straight_spline_points=100)
-"""
+
 
 oval_sketch_2 = cb.SplineDisk(center_point=center_point + direction * 4,
                               corner_1_point=np.asarray([0, 1.5, 0]) + direction * 4,
@@ -52,16 +43,16 @@
 # The ring is defined as the disk,
 # where a ring with same center, corner_1... as a disk will fit on the outside of the disk.
 oval_ring_1 = cb.SplineRing(center_point=oval_sketch_2.center,
-                            corner_1_point=oval_sketch_2.corner_1,
-                            corner_2_point=oval_sketch_2.corner_2,
+                            corner_1_point=oval_sketch_2.radius_1_point,
+                            corner_2_point=oval_sketch_2.radius_2_point,
                             side_1=oval_sketch_2.side_1,
                             side_2=oval_sketch_2.side_2,
                             width_1=0.2, width_2=0.2, n_outer_spline_points=100)
 
 # Note it is possible to access the center and corners of a defined sketch. These are stable on transformation.
 oval_ring_2 = cb.SplineRing(center_point=oval_sketch_2.center + direction,
-                            corner_1_point=oval_sketch_2.corner_1 + direction,
-                            corner_2_point=oval_sketch_2.corner_2 + direction,
+                            corner_1_point=oval_sketch_2.radius_1_point + direction,
+                            corner_2_point=oval_sketch_2.radius_2_point + direction,
                             side_1=0,
                             side_2=0,
                             width_1=0.2, width_2=0.2, n_outer_spline_points=100)

src/classy_blocks/__init__.py

@@ -25,7 +25,7 @@
 from .construct.operations.revolve import Revolve
 from .construct.operations.wedge import Wedge
 from .construct.shape import ExtrudedShape, LoftedShape, RevolvedShape, Shape
-from .construct.shapes.cylinder import Cylinder, SemiCylinder
+from .construct.shapes.cylinder import Cylinder, QuarterCylinder, SemiCylinder
 from .construct.shapes.elbow import Elbow
 from .construct.shapes.frustum import Frustum
 from .construct.shapes.rings import ExtrudedRing, RevolvedRing
@@ -100,6 +100,7 @@
     "Frustum",
     "Cylinder",
     "SemiCylinder",
+    "QuarterCylinder",
     "ExtrudedRing",
     "RevolvedRing",
     "Hemisphere",

src/classy_blocks/construct/flat/sketches/annulus.py

@@ -24,12 +24,13 @@ def __init__(
         normal: VectorType,
         inner_radius: float,
         n_segments: int = 8,
+        angle: float = 2 * np.pi,
     ):
         center_point = np.asarray(center_point)
         normal = f.unit_vector(np.asarray(normal))
         outer_radius_point = np.asarray(outer_radius_point)
         inner_radius_point = center_point + f.unit_vector(outer_radius_point - center_point) * inner_radius
-        segment_angle = 2 * np.pi / n_segments
+        segment_angle = angle / n_segments
 
         face = Face(
             [  # points
@@ -66,7 +67,17 @@ def grid(self):
 
     @property
     def center(self) -> NPPointType:
-        return np.average([f.center for f in self.faces], axis=0)
+        """Return center of sketch by assuming radial sides of faces intersect in the center"""
+        return np.average(
+            [
+                p[0]
+                - f.norm(np.cross(p[2] - p[3], p[3] - p[0]))
+                / f.norm(np.cross(p[2] - p[3], p[1] - p[0]))
+                * (p[1] - p[0])
+                for p in (face.point_array for face in self.faces)
+            ],
+            axis=0,
+        )
 
     @property
     def n_segments(self):

src/classy_blocks/construct/flat/sketches/disk.py

@@ -1,12 +1,21 @@
 import abc
-from typing import ClassVar, List
+from typing import ClassVar, List, Optional
 
 import numpy as np
 
 from classy_blocks.construct.edges import Origin, Spline
 from classy_blocks.construct.flat.face import Face
 from classy_blocks.construct.flat.sketches.mapped import MappedSketch
-from classy_blocks.types import NPPointListType, NPPointType, NPVectorType, PointType, VectorType
+from classy_blocks.construct.point import Point
+from classy_blocks.types import (
+    IndexType,
+    NPPointListType,
+    NPPointType,
+    NPVectorType,
+    PointListType,
+    PointType,
+    VectorType,
+)
 from classy_blocks.util import functions as f
 
 
@@ -51,8 +60,8 @@ class DiskBase(MappedSketch, abc.ABC):
         0.421177,
         0.502076,
         0.566526,
-        0.610535,
         0.610532,
+        0.610535,
         0.625913,
         0.652300,
         0.686516,
@@ -61,52 +70,84 @@ class DiskBase(MappedSketch, abc.ABC):
         0.792025,
     )
 
+    def __init__(self, positions: PointListType, quads: List[IndexType]):
+        # Center point as a constant.
+        self.origo_point = Point(positions[0])
+
+        super().__init__(positions, quads)
+
+    @property
+    def origo(self):
+        return self.origo_point.position
+
     # Relative size of the inner square (O-1-2-3) in a single core cylinder:
     # - too small will cause unnecessary high number of small cells in the square;
     # - too large will prevent creating large numbers of boundary layers
     @property
     def diagonal_ratio(self) -> float:
-        return 2**0.5 * self.spline_ratios[8] / 0.8 * self.core_ratio
-
-    def add_spline_edges(self) -> None:
+        return 2**0.5 * self.spline_ratios[7] / 0.8 * self.core_ratio
+
+    def circular_core_spline(
+        self,
+        p_core_ratio: PointType,
+        p_diagonal_ratio: PointType,
+        reverse: bool = False,
+        center: Optional[PointType] = None,
+    ) -> NPPointListType:
+        """Creates the spline points for the core."""
+        p_0 = np.asarray(p_core_ratio)
+        p_1 = np.asarray(p_diagonal_ratio)
+        if center is None:
+            center = self.center
+
+        # Spline points in unitary coordinates
+        spline_points_u = np.array([self.spline_ratios[-1:6:-1]]).T * np.array([0, 1, 0]) + np.array(
+            [self.spline_ratios[:7]]
+        ).T * np.array([0, 0, 1])
+
+        # p_1 and p_2 in unitary coordinates
+        p_0_u = np.array([0, 0.8, 0])
+        p_1_u = np.array([0, 6.10535e-01, 6.10535e-01])
+
+        # orthogonal vectors based on p_0_u and p_1_u
+        u_0_org = p_0_u
+        u_1_org = p_1_u - np.dot(p_1_u, f.unit_vector(u_0_org)) * f.unit_vector(u_0_org)
+
+        # Spline points in u_0_org and u_1_org
+        spline_d_0_org = np.dot(spline_points_u, f.unit_vector(u_0_org)).reshape((-1, 1)) / f.norm(u_0_org)
+        spline_d_1_org = np.dot(spline_points_u, f.unit_vector(u_1_org)).reshape((-1, 1)) / f.norm(u_1_org)
+
+        # New plane defined by new points
+        u_0 = p_0 - center
+        u_1 = p_1 - center - np.dot(p_1 - center, f.unit_vector(u_0)) * f.unit_vector(u_0)
+
+        spline_points_new = center + spline_d_0_org * u_0 + spline_d_1_org * u_1
+        if reverse:
+            return spline_points_new[::-1]
+        else:
+            return spline_points_new
+
+    def add_core_spline_edges(self) -> None:
         """Add a spline to the core blocks for an optimized mesh."""
-        spline_ratios = np.array(self.spline_ratios) / 0.8 * self.core_ratio
-        spl_len = len(spline_ratios)
-
-        spline_points = [
-            self.center
-            + self.radius_vector * spline_ratios[i]
-            + self.perp_radius_vector * spline_ratios[spl_len - i - 1]
-            for i in range(spl_len)
-        ]
-        spline_points.reverse()
-
-        for i in range(len(self.core)):
-            angle = i * np.pi / 2
-            points = [f.rotate(p, angle, self.normal, self.center) for p in spline_points]
+        for i, face in enumerate(self.core):
+            p_0 = face.point_array[(i + 1) % 4]  # Core point on radius vector
+            p_1 = face.point_array[(i + 2) % 4]  # Core point on diagonal
+            p_2 = face.point_array[(i + 3) % 4]  # Core point on perpendicular radius vector
 
-            points_1 = points[: spl_len // 2 - 1]
-            points_2 = points[1 + spl_len // 2 :]
-
-            if i == 2:
-                points_1.reverse()
-            if i == 1:
-                points_2.reverse()
-
-            curve_1 = Spline(points_1)
-            curve_2 = Spline(points_2)
+            spline_curve_0_1 = Spline(self.circular_core_spline(p_0, p_1, reverse=i == 2))
+            spline_curve_1_2 = Spline(self.circular_core_spline(p_2, p_1, reverse=i != 1))
 
+            # Add curves to edges
             edge_1 = (i + 1) % 4
             edge_2 = (i + 2) % 4
-
-            self.grid[0][i].add_edge(edge_1, curve_1)
-            self.grid[0][i].add_edge(edge_2, curve_2)
+            face.add_edge(edge_1, spline_curve_0_1)
+            face.add_edge(edge_2, spline_curve_1_2)
 
     def add_edges(self):
-        for face in self.grid[-1]:
-            face.add_edge(1, Origin(self.center))
+        for face in self.shell:
+            face.add_edge(1, Origin(self.origo))
 
-        self.add_spline_edges()
+        self.add_core_spline_edges()
 
     @property
     def center(self) -> NPPointType:
@@ -116,19 +157,15 @@ def center(self) -> NPPointType:
     @property
     def radius_point(self) -> NPPointType:
         """Point at outer radius"""
-        return self.grid[1][0].points[1].position
+        return self.shell[0].points[1].position
 
     @property
     def radius_vector(self) -> NPVectorType:
         """Vector that points from center of this
-        *Circle to its (first) radius point"""
-        return self.radius_point - self.center
-
-    @property
-    def perp_radius_vector(self) -> NPVectorType:
-        """Vector that points from center of this
-        *Circle to its (second) radius point"""
-        return f.unit_vector(np.cross(self.normal, self.radius_vector)) * f.norm(self.radius_vector)
+        *Circle to its (first) radius point
+        Origo is used instead of center to ensure outside is constant, when moving the core,
+        does not change the outer shape."""
+        return self.radius_point - self.origo
 
     @property
     def radius(self) -> float:
@@ -147,6 +184,10 @@ def core(self) -> List[Face]:
     def shell(self) -> List[Face]:
         return self.grid[-1]
 
+    @property
+    def parts(self):
+        return [self.origo_point, *super().parts]
+
 
 class OneCoreDisk(DiskBase):
     """A disk with a single block in  the center and four blocks around;
@@ -182,15 +223,15 @@ def center(self):
     def grid(self):
         return [self.faces[:1], self.faces[1:]]
 
-    def add_spline_edges(self):
+    def add_core_spline_edges(self):
         pass
 
 
 class QuarterDisk(DiskBase):
     """A quarter of a four-core disk; see docs/blocking for point numbers and faces/grid indexing"""
 
     chops: ClassVar = [
-        [0],  # axis 0
+        [1],  # axis 0
         [1, 2],  # axis 1
     ]
 
@@ -333,7 +374,7 @@ def grid(self):
 
     def add_edges(self):
         for face in self.grid[1]:
-            face.add_edge(1, Origin(self.center))
+            face.add_edge(1, Origin(self.origo))
 
     @property
     def center(self):