Refactor SemiCylinders/Joints and fix bugs from other parts

Author: FranzBangar

examples/assembly/n_joint.py

@@ -0,0 +1,22 @@
+import os
+
+import classy_blocks as cb
+from classy_blocks.construct.assemblies.joints import NJoint
+
+mesh = cb.Mesh()
+
+axis_point_1 = [0.0, 0.0, 0.0]
+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)
+
+joint = NJoint([0, -1, 0], [0, 0, 0], [0.2, -1, 0], 3)
+
+for operation in joint.operations:
+    for i in range(3):
+        operation.chop(i, count=10)
+
+mesh.add(joint)
+
+mesh.write(os.path.join("..", "case", "system", "blockMeshDict"), debug_path="debug.vtk")

src/classy_blocks/construct/array.py

@@ -74,7 +74,7 @@ def shear(self, normal: VectorType, origin: PointType, direction: VectorType, an
             if f.point_to_plane_distance(origin, normal, point) > TOL:
                 distance = f.point_to_plane_distance(origin, normal, point)
                 direction = f.unit_vector(direction)
-                amount = distance * np.tan(angle)
+                amount = distance / np.tan(angle)
 
                 self.points[i] += direction * amount
         return self

src/classy_blocks/construct/assemblies/joints.py

@@ -1,11 +1,70 @@
+from typing import List
+
 import numpy as np
 
 from classy_blocks.construct.assemblies.assembly import Assembly
-from classy_blocks.construct.shapes.cylinder import SemiCylinder, SlashedCylinder
+from classy_blocks.construct.edges import Spline
+from classy_blocks.construct.flat.sketches.disk import HalfDisk
+from classy_blocks.construct.shape import Shape
+from classy_blocks.construct.shapes.cylinder import SemiCylinder
 from classy_blocks.types import PointType
 from classy_blocks.util import functions as f
 
 
+class CuspSemiCylinder(SemiCylinder):
+    sketch_class = HalfDisk
+
+    def __init__(
+        self, axis_point_1: PointType, axis_point_2: PointType, radius_point_1: PointType, end_angle: float = np.pi / 4
+    ):
+        axis_point_1 = np.asarray(axis_point_1)
+        axis_point_2 = np.asarray(axis_point_2)
+        axis = axis_point_2 - axis_point_1
+        radius_point_1 = np.asarray(radius_point_1)
+        radius_vector = radius_point_1 - axis_point_1
+
+        shear_normal = np.cross(-axis, radius_vector)
+
+        super().__init__(axis_point_1, axis_point_2, radius_point_1)
+
+        for loft in self.operations:
+            loft.top_face.remove_edges()
+
+        for loft in self.shell:
+            point_1 = loft.top_face.points[1].position
+            point_2 = loft.top_face.points[2].position
+            edge_points = f.divide_arc(axis, axis_point_2, point_1, point_2, 5)
+            loft.top_face.add_edge(1, Spline(edge_points))
+
+        for loft in self.operations:
+            loft.top_face.shear(shear_normal, axis_point_2, -axis, end_angle)
+
+        # TODO: include those inner edges (Disk.spline_ratios > Curve)
+        # self.remove_inner_edges(end=True)
+
+
+class CuspCylinder(Assembly):
+    def __init__(
+        self,
+        axis_point_1: PointType,
+        axis_point_2: PointType,
+        radius_point: PointType,
+        end_angle_left: float,
+        end_angle_right: float,
+    ):
+        axis_point_1 = np.asarray(axis_point_1)
+        radius_point_right = np.asarray(radius_point)
+        radius_vector_right = np.asarray(radius_point) - axis_point_1
+        radius_point_left = axis_point_1 - radius_vector_right
+
+        self.cusp_right = CuspSemiCylinder(axis_point_1, axis_point_2, radius_point_right, end_angle_right)
+        self.cusp_left = CuspSemiCylinder(axis_point_1, axis_point_2, radius_point_left, end_angle_left)
+
+    @property
+    def shapes(self):
+        return [self.cusp_right, self.cusp_left]
+
+
 class TJoint(Assembly):
     def __init__(self, start_point: PointType, center_point: PointType, right_point: PointType, radius: float):
         start_point = np.array(start_point)
@@ -19,21 +78,35 @@ def __init__(self, start_point: PointType, center_point: PointType, right_point:
 
         radius_vector = radius * f.unit_vector(np.cross(right_vector, start_vector))
 
-        # the right part
-        start_right_sc = SlashedCylinder(start_point, center_point, start_point + radius_vector)
-        right_sc = SlashedCylinder(right_point, center_point, right_point - radius_vector)
+        # middle
+        self.middle = CuspCylinder(start_point, center_point, start_point + radius_vector, np.pi / 4, np.pi / 4)
+
+        # the right 'arm'
+        self.right = CuspCylinder(right_point, center_point, right_point + radius_vector, np.pi / 4, np.pi / 2)
+
+        # the left 'arm'
+        self.left = CuspCylinder(left_point, center_point, left_point - radius_vector, np.pi / 4, np.pi / 2)
+
+        super().__init__([*self.middle.shapes, *self.right.shapes, *self.left.shapes])
+
 
-        # the left part
-        start_left_sc = SlashedCylinder(start_point, center_point, start_point - radius_vector)
-        left_sc = SlashedCylinder(left_point, center_point, left_point + radius_vector)
+class NJoint(Assembly):
+    def __init__(self, start_point: PointType, center_point: PointType, radius_point: PointType, branches: int = 4):
+        start_point = np.asarray(start_point)
+        center_point = np.asarray(center_point)
+        radius_point = np.asarray(radius_point)
 
-        # the top part
-        top_right = SemiCylinder(center_point, right_point, center_point - radius_vector)
-        top_left = SemiCylinder(center_point, left_point, center_point + radius_vector)
+        self.assemblies: List[Assembly] = []
+        shapes: List[Shape] = []
 
-        shapes = [start_right_sc, right_sc, start_left_sc, left_sc, top_right, top_left]
+        cusp_angle = np.pi / branches
+        base_asm = CuspCylinder(start_point, center_point, radius_point, cusp_angle, cusp_angle)
+        rotate_angles = np.linspace(0, 2 * np.pi, num=branches, endpoint=False)
+        rotation_axis = radius_point - start_point
 
-        for shape in shapes:
-            shape.remove_inner_edges(end=True)
+        for angle in rotate_angles:
+            asm = base_asm.copy().rotate(angle, rotation_axis, center_point)
+            self.assemblies.append(asm)
+            shapes += asm.shapes
 
         super().__init__(shapes)

src/classy_blocks/construct/shapes/cylinder.py

@@ -39,34 +39,6 @@ def __init__(self, axis_point_1: PointType, axis_point_2: PointType, radius_poin
         super().__init__(self.sketch_class(axis_point_1, radius_point_1, axis), transform_2, None)
 
 
-class SlashedCylinder(RoundSolidShape):
-    sketch_class = HalfDisk
-
-    def __init__(
-        self, axis_point_1: PointType, axis_point_2: PointType, radius_point_1: PointType, end_angle: float = np.pi / 4
-    ):
-        axis_point_1 = np.asarray(axis_point_1)
-        axis_point_2 = np.asarray(axis_point_2)
-        axis = axis_point_2 - axis_point_1
-        radius_point_1 = np.asarray(radius_point_1)
-        radius_vector = radius_point_1 - axis_point_1
-
-        shear_normal = np.cross(-axis, radius_vector)
-
-        diff = np.dot(axis, radius_point_1 - axis_point_1)
-        if diff > TOL:
-            raise CylinderCreationError(
-                "Axis and radius vectors are not perpendicular", f"Difference: {diff}, tolerance: {TOL}"
-            )
-
-        transform_2: List[tr.Transformation] = [
-            tr.Translation(axis),
-            tr.Shear(shear_normal, axis_point_2, -axis, end_angle),
-        ]
-
-        super().__init__(self.sketch_class(axis_point_1, radius_point_1, axis), transform_2, None)
-
-
 class Cylinder(SemiCylinder):
     sketch_class = Disk
 

src/classy_blocks/construct/shapes/sphere.py

@@ -101,16 +101,11 @@ class EighthSphere(Shape):
     def __init__(
         self, center_point: PointType, radius_point: PointType, normal: VectorType, diagonal_angle: float = np.pi / 5
     ):
-        # TODO: move these to properties
-        # BUG: move these to properties
-        # (will not change with transforms!)
-        self.center_point = np.asarray(center_point)
-        self.radius_point = np.asarray(radius_point)
-        self.normal = f.unit_vector(np.asarray(normal))
+        center_point = np.asarray(center_point)
+        radius_point = np.asarray(radius_point)
+        normal = f.unit_vector(np.asarray(normal))
 
-        self.lofts = eighth_sphere_lofts(
-            self.center_point, self.radius_point, self.normal, self.geometry_label, diagonal_angle
-        )
+        self.lofts = eighth_sphere_lofts(center_point, radius_point, normal, self.geometry_label, diagonal_angle)
 
     ### Chopping
     def chop_axial(self, **kwargs):
@@ -156,6 +151,18 @@ def shell(self):
     def grid(self):
         return [self.core, self.shell]
 
+    @property
+    def radius_point(self) -> NPPointType:
+        return self.lofts[1].bottom_face.points[1].position
+
+    @property
+    def center_point(self) -> NPPointType:
+        return self.lofts[0].bottom_face.points[0].position
+
+    @property
+    def normal(self) -> NPVectorType:
+        return self.lofts[0].bottom_face.normal
+
     @property
     def radius(self) -> float:
         """Radius of this sphere"""
@@ -166,6 +173,10 @@ def geometry_label(self) -> str:
         """Name of a unique geometry this will project to"""
         return f"sphere_{id(self)}"
 
+    @property
+    def center(self):
+        return self.center_point
+
     @property
     def geometry(self):
         return {

src/classy_blocks/items/edges/arcs/angle.py

@@ -36,9 +36,8 @@ def arc_from_theta(edge_point_1: PointType, edge_point_2: PointType, angle: floa
     mag_chord = f.norm(chord)
 
     center = pm - length * axis / 2 - rm * mag_chord / 2 / np.tan(angle / 2)
-    radius = f.norm(edge_point_1 - center)
 
-    return f.arc_mid(axis, center, radius, edge_point_1, edge_point_2)
+    return f.arc_mid(axis, center, edge_point_1, edge_point_2)
 
 
 @dataclasses.dataclass

src/classy_blocks/items/edges/arcs/origin.py

@@ -73,7 +73,7 @@ def arc_from_origin(
         return arc_from_origin(p1, p3, new_center, False)
 
     # done, return the calculated point
-    return f.arc_mid(axis, center, radius, edge_point_1, edge_point_2)
+    return f.arc_mid(axis, center, edge_point_1, edge_point_2)
 
 
 @dataclasses.dataclass

src/classy_blocks/util/functions.py

@@ -201,19 +201,27 @@ def arc_length_3point(p_start: NPPointType, p_btw: NPPointType, p_end: NPPointTy
     return angle * norm(radius)
 
 
-def arc_mid(axis: VectorType, center: PointType, radius: float, point_1: PointType, point_2: PointType) -> PointType:
-    """Returns the midpoint of the specified arc in 3D space"""
+def divide_arc(
+    axis: VectorType, center: PointType, point_1: PointType, point_2: PointType, count: int
+) -> NPPointListType:
     # Kudos to this guy for his shrewd solution
     # https://math.stackexchange.com/questions/3717427
+    # (extended here to create more than 1 point)
     axis = np.asarray(axis, dtype=constants.DTYPE)
     center = np.asarray(center, dtype=constants.DTYPE)
     point_1 = np.asarray(point_1, dtype=constants.DTYPE)
     point_2 = np.asarray(point_2, dtype=constants.DTYPE)
+    radius = norm(center - point_1)
+
+    secant_points = np.linspace(point_1, point_2, num=count + 2)[1:-1]
+    secant_vectors = [unit_vector(point - center) for point in secant_points]
 
-    sec = point_2 - point_1
-    sec_ort = np.cross(sec, axis)
+    return np.array([center + vector * radius for vector in secant_vectors])
 
-    return center + unit_vector(sec_ort) * radius
+
+def arc_mid(axis: VectorType, center: PointType, point_1: PointType, point_2: PointType) -> PointType:
+    """Returns the midpoint of the specified arc in 3D space"""
+    return divide_arc(axis, center, point_1, point_2, 1)[0]
 
 
 def mirror_matrix(normal: VectorType):

tests/test_construct/test_curves/test_interpolated.py

@@ -68,6 +68,21 @@ def test_transform(self):
     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.shear([0, 1, 0], [0, 0, 0], [1, 0, 0], np.pi / 4)
+
+        expected = [
+            [0, 0, 0],
+            [1, 1, 0],
+            [2, 1, 0],
+            [1, 0, 0],
+            [2, 0, 0],
+        ]
+
+        np.testing.assert_almost_equal(curve.discretize(count=5), expected)
+
 
 class SplineInterpolatedCurveTests(unittest.TestCase):
     def setUp(self):

tests/test_items/test_edge.py

@@ -269,13 +269,10 @@ class AlternativeArcTests(unittest.TestCase):
     def test_arc_mid(self):
         axis = f.vector(0, 0, 1)
         center = f.vector(0, 0, 0)
-        radius = 1
         edge_point_1 = f.vector(1, 0, 0)
         edge_point_2 = f.vector(0, 1, 0)
 
-        np.testing.assert_array_almost_equal(
-            f.arc_mid(axis, center, radius, edge_point_1, edge_point_2), self.unit_sq_corner
-        )
+        np.testing.assert_array_almost_equal(f.arc_mid(axis, center, edge_point_1, edge_point_2), self.unit_sq_corner)
 
     def test_arc_from_theta(self):
         edge_point_1 = f.vector(0, 1, 0)