Repository: fusion-energy/paramak (language: python)

Issue:
Add new component for dished vacuum vessel heads
I would be handy to have a few curved shapes for vacuum vessel heads.

The straight section can be added afterwards.

Then vacuum vessels could be made from a mixture of heads




Discussion / hints:

Fix patch:
diff --git a/docs/source/API-Reference.rst b/docs/source/API-Reference.rst
index dedbf4cfa..aa38f3b5d 100644
--- a/docs/source/API-Reference.rst
+++ b/docs/source/API-Reference.rst
@@ -1516,6 +1516,43 @@ CapsuleVacuumVessel()
    :show-inheritance:
 
 
+DishedVacuumVessel()
+^^^^^^^^^^^^^^^^^^^^
+.. cadquery::
+   :select: cadquery_object
+   :gridsize: 0
+
+   import paramak
+   my_component = paramak.DishedVacuumVessel(rotation_angle=180)
+
+   cadquery_object = my_component.solid
+
+.. image:: https://user-images.githubusercontent.com/8583900/160503281-3aabf145-22b0-4953-bc4f-3f75a6696b5e.png
+
+.. automodule:: paramak.parametric_components.dished_vacuum_vessel
+   :members:
+   :show-inheritance:
+
+ConstantThicknessDome()
+^^^^^^^^^^^^^^^^^^^^^^^
+.. cadquery::
+   :select: cadquery_object
+   :gridsize: 0
+
+   import paramak
+   my_component = paramak.ConstantThicknessDome(rotation_angle=180)
+
+   cadquery_object = my_component.solid
+
+.. image:: https://user-images.githubusercontent.com/8583900/160503286-71a6d771-1d47-476e-a0cf-85fb94c9389c.png
+
+.. automodule:: paramak.parametric_components.constant_thickness_dome
+   :members:
+   :show-inheritance:
+
+
+
+
 Other components
 ----------------
 
diff --git a/paramak/__init__.py b/paramak/__init__.py
index 7f8a70a5b..b631f8399 100644
--- a/paramak/__init__.py
+++ b/paramak/__init__.py
@@ -13,7 +13,15 @@
 
 from .shape import Shape
 from .reactor import Reactor
-from .utils import rotate, extend, distance_between_two_points, diff_between_angles
+from .utils import (
+    rotate,
+    extend,
+    distance_between_two_points,
+    diff_between_angles,
+    find_center_point_of_circle,
+    angle_between_two_points_on_circle,
+    find_radius_of_circle,
+)
 from .utils import EdgeLengthSelector, FaceAreaSelector
 
 from .parametric_shapes.extruded_mixed_shape import ExtrudeMixedShape
@@ -113,7 +121,10 @@
     ToroidalFieldCoilRectangleRoundCorners,
 )
 
+from .parametric_components.constant_thickness_dome import ConstantThicknessDome
 from .parametric_components.vacuum_vessel import VacuumVessel
+
+from .parametric_components.dished_vacuum_vessel import DishedVacuumVessel
 from .parametric_components.vacuum_vessel_inner_leg import VacuumVesselInnerLeg
 from .parametric_components.capsule_vacuum_vessel import CapsuleVacuumVessel
 from .parametric_components.hollow_cube import HollowCube
diff --git a/paramak/parametric_components/capsule_vacuum_vessel.py b/paramak/parametric_components/capsule_vacuum_vessel.py
index bc77bfa2a..26e61432c 100644
--- a/paramak/parametric_components/capsule_vacuum_vessel.py
+++ b/paramak/parametric_components/capsule_vacuum_vessel.py
@@ -6,7 +6,7 @@
 
 class CapsuleVacuumVessel(RotateMixedShape):
     """A cylindrical vessel volume with constant thickness that has addition
-    spherical edges.
+    hemispherical head.
 
     Arguments:
         outer_start_point: the x,z coordinates of the outer bottom of the
diff --git a/paramak/parametric_components/constant_thickness_dome.py b/paramak/parametric_components/constant_thickness_dome.py
new file mode 100644
index 000000000..bff2acd3e
--- /dev/null
+++ b/paramak/parametric_components/constant_thickness_dome.py
@@ -0,0 +1,214 @@
+import math
+from paramak import RotateMixedShape, RotateStraightShape, Shape, CuttingWedge
+import cadquery as cq
+
+
+class ConstantThicknessDome(RotateMixedShape):
+    """A cylindrical vessel volume with constant thickness with a simple dished
+    head. This style of tank head has no knuckle radius or straight flange. The
+    dished shape is made from a chord of a circle.
+
+    Arguments:
+        thickness: the radial thickness of the dome.
+        chord_center_height: the vertical position of the chord center
+        chord_width: the width of the chord base
+        chord_height: the height of the chord which is also distance between
+            the chord_center_height and the inner surface of the dome
+        upper_or_lower: Curves the dish with a positive or negative direction
+            to allow the upper section or lower section of vacuum vessel
+            domes to be made.
+        name: the name of the shape, used in the graph legend and as a
+            filename prefix when exporting.
+    """
+
+    def __init__(
+        self,
+        thickness: float = 10,
+        chord_center_height: float = 0,
+        chord_width: float = 100,
+        chord_height: float = 20,
+        upper_or_lower: str = "upper",
+        name: str = "constant_thickness_dome",
+        **kwargs,
+    ):
+
+        self.thickness = thickness
+        self.chord_center_height = chord_center_height
+        self.chord_width = chord_width
+        self.chord_height = chord_height
+        self.upper_or_lower = upper_or_lower
+        self.name = name
+
+        super().__init__(name=name, **kwargs)
+
+    @property
+    def chord_width(self):
+        return self._chord_width
+
+    @chord_width.setter
+    def chord_width(self, value):
+        if not isinstance(value, (float, int)):
+            raise ValueError("ConstantThicknessDome.chord_width must be a number. Not", value)
+        if value <= 0:
+            msg = f"ConstantThicknessDome.chord_width must be a positive number above 0. Not {value}"
+            raise ValueError(msg)
+        self._chord_width = value
+
+    @property
+    def chord_height(self):
+        return self._chord_height
+
+    @chord_height.setter
+    def chord_height(self, value):
+        if not isinstance(value, (float, int)):
+            raise ValueError("ConstantThicknessDome.chord_height must be a number. Not", value)
+        if value <= 0:
+            msg = f"ConstantThicknessDome.chord_height must be a positive number above 0. Not {value}"
+            raise ValueError(msg)
+        self._chord_height = value
+
+    @property
+    def thickness(self):
+        return self._thickness
+
+    @thickness.setter
+    def thickness(self, value):
+        if not isinstance(value, (float, int)):
+            msg = f"VacuumVessel.thickness must be a number. Not {value}"
+            raise ValueError(msg)
+        if value <= 0:
+            msg = f"VacuumVessel.thickness must be a positive number above 0. Not {value}"
+            raise ValueError(msg)
+        self._thickness = value
+
+    def find_points(self):
+        """
+        Finds the XZ points joined by straight and circle connections that
+        describe the 2D profile of the vessel shape.
+        """
+
+        # Note these points are not used in the normal way when constructing
+        # the solid
+        #
+        #          6   -
+        #          |       -
+        #          7  -       -
+        #                -       -
+        #                  -       3
+        #                    -     |
+        #         cc          1 -- 2
+        #     chord center
+        #
+        #
+        #          cp
+        #     center point
+        #
+        #
+        #
+        #
+        #         cc           1 -- 2
+        #                    -      |
+        #                  -        3
+        #                -       -
+        #          7  -       -
+        #          |       -
+        #          6   -
+        #       far side
+
+        if self.chord_height * 2 >= self.chord_width:
+            msg = "ConstantThicknessDome requires that the self.chord_width is at least 2 times as large as the chord height"
+            raise ValueError(msg)
+
+        radius_of_sphere = ((math.pow(self.chord_width, 2)) + (4.0 * math.pow(self.chord_height, 2))) / (
+            8 * self.chord_height
+        )
+
+        # TODO set to 0 for now, add ability to shift the center of the chord left and right
+        chord_center = (0, self.chord_center_height)
+
+        point_1 = (chord_center[0] + (self.chord_width / 2), chord_center[1], "straight")
+
+        if self.upper_or_lower == "upper":
+            center_point = (chord_center[0], chord_center[1] + self.chord_height - radius_of_sphere)
+            inner_tri_angle = math.atan((center_point[1] - chord_center[1]) / (self.chord_width / 2))
+            outer_tri_adj = math.cos(inner_tri_angle) * self.thickness
+            # original ending type
+            # point_2 = (point_1[0] + outer_tri_adj, point_1[1], "straight")
+            point_2 = (point_1[0] + self.thickness, point_1[1], "straight")
+            outer_tri_opp = math.sqrt(math.pow(self.thickness, 2) - math.pow(outer_tri_adj, 2))
+            point_7 = (chord_center[0], chord_center[1] + radius_of_sphere, "straight")
+            point_6 = (chord_center[0], chord_center[1] + radius_of_sphere + self.thickness, "straight")
+            far_side = (center_point[0], center_point[1] - (radius_of_sphere + self.thickness))
+            point_3 = (point_2[0], point_2[1] + outer_tri_opp, "straight")
+        elif self.upper_or_lower == "lower":
+            center_point = (chord_center[0], chord_center[1] - self.chord_height + radius_of_sphere)
+            inner_tri_angle = math.atan((center_point[1] - chord_center[1]) / (self.chord_width / 2))
+            outer_tri_adj = math.cos(inner_tri_angle) * self.thickness
+            # original ending type
+            # point_2 = (point_1[0] + outer_tri_adj, point_1[1], "straight")
+            point_2 = (point_1[0] + self.thickness, point_1[1], "straight")
+            outer_tri_opp = math.sqrt(math.pow(self.thickness, 2) - math.pow(outer_tri_adj, 2))
+            point_7 = (chord_center[0], chord_center[1] - radius_of_sphere, "straight")
+            point_6 = (chord_center[0], chord_center[1] - (radius_of_sphere + self.thickness), "straight")
+            far_side = (center_point[0], center_point[1] + radius_of_sphere + self.thickness)
+            point_3 = (point_2[0], point_2[1] - outer_tri_opp, "straight")
+        else:
+            msg = f'upper_or_lower should be either "upper"  or "lower". Not {self.upper_or_lower}'
+            raise ValueError(msg)
+
+        self.points = [point_1, point_2, point_3, point_6, point_7]
+
+    def create_solid(self):
+        """Creates a rotated 3d solid using points with circular edges.
+
+        Returns:
+           A CadQuery solid: A 3D solid volume
+        """
+
+        radius_of_sphere = ((math.pow(self.chord_width, 2)) + (4.0 * math.pow(self.chord_height, 2))) / (
+            8 * self.chord_height
+        )
+
+        # TODO set to 0 for now, add ability to shift the center of the chord left and right
+        chord_center = (0, self.chord_center_height)
+
+        if self.upper_or_lower == "upper":
+            center_point = (chord_center[0], chord_center[1] + self.chord_height - radius_of_sphere)
+            far_side = (center_point[0], center_point[1] - (radius_of_sphere + self.thickness))
+        elif self.upper_or_lower == "lower":
+            center_point = (chord_center[0], chord_center[1] - self.chord_height + radius_of_sphere)
+            far_side = (center_point[0], center_point[1] + radius_of_sphere + self.thickness)
+        else:
+            raise ValueError("self.upper_or_lower")
+
+        big_sphere = (
+            cq.Workplane(self.workplane)
+            .moveTo(center_point[0], center_point[1])
+            .sphere(radius_of_sphere + self.thickness)
+        )
+        small_sphere = cq.Workplane(self.workplane).moveTo(center_point[0], center_point[1]).sphere(radius_of_sphere)
+
+        outer_cylinder_cutter = RotateStraightShape(
+            workplane=self.workplane,
+            points=(
+                (chord_center[0], chord_center[1]),  # cc
+                (self.points[1][0], self.points[1][1]),  # point 2
+                (self.points[2][0], self.points[2][1]),  # point 3
+                (self.points[2][0] + radius_of_sphere, self.points[2][1]),  # point 3 wider
+                (self.points[2][0] + radius_of_sphere, far_side[1]),
+                far_side,
+            ),
+            rotation_angle=360,
+        )
+
+        cap = Shape()
+        cap.solid = big_sphere.cut(small_sphere)
+
+        height = 2 * (radius_of_sphere + abs(center_point[1]) + self.thickness)
+        radius = 2 * (radius_of_sphere + abs(center_point[0]) + self.thickness)
+        cutter = CuttingWedge(height=height, radius=radius, rotation_angle=self.rotation_angle)
+
+        cap.solid = cap.solid.intersect(cutter.solid)
+        cap.solid = cap.solid.cut(outer_cylinder_cutter.solid)
+
+        self.solid = cap.solid
diff --git a/paramak/parametric_components/cutting_wedge_fs.py b/paramak/parametric_components/cutting_wedge_fs.py
index f732e8850..18f592e45 100644
--- a/paramak/parametric_components/cutting_wedge_fs.py
+++ b/paramak/parametric_components/cutting_wedge_fs.py
@@ -94,7 +94,10 @@ def azimuth_placement_angle(self, value):
     def find_radius_height(self):
         shape = self.shape
         if shape.rotation_angle == 360:
-            msg = "cutting_wedge cannot be created, rotation_angle must be " "less than 360 degrees"
+            msg = (
+                "cutting_wedge cannot be created, the shapes rotation_angle "
+                f"must be less than 360 degrees. shape.rotation_angle is {shape.rotation_angle}"
+            )
             raise ValueError(msg)
         shape_points = shape.points
         if hasattr(shape, "radius") and len(shape_points) == 1:
diff --git a/paramak/parametric_components/dished_vacuum_vessel.py b/paramak/parametric_components/dished_vacuum_vessel.py
new file mode 100644
index 000000000..068a74537
--- /dev/null
+++ b/paramak/parametric_components/dished_vacuum_vessel.py
@@ -0,0 +1,123 @@
+
+from paramak import RotateMixedShape, CenterColumnShieldCylinder, ConstantThicknessDome
+
+
+class DishedVacuumVessel(RotateMixedShape):
+    """A cylindrical vessel volume with constant thickness with a simple dished
+    head. This style of tank head has no knuckle radius or straight flange.
+
+    Arguments:
+        radius: the radius from which the centres of the vessel meets the outer
+            circumference.
+        center_point: the x,z coordinates of the center of the vessel
+        dish_height: the height of the dish section. This is also the chord
+            heigh of the circle used to make the dish.
+        cylinder_height: the height of the cylindrical section of the vacuum
+            vessel.
+        thickness: the radial thickness of the vessel in cm.
+    """
+
+    def __init__(
+        self,
+        radius: float=300,
+        center_point: float=0,
+        dish_height: float=50,
+        cylinder_height: float=400,
+        thickness: float=15,
+        name: str='dished_vacuum_vessel',
+        **kwargs,
+    ):
+        self.radius = radius
+        self.center_point = center_point
+        self.dish_height = dish_height
+        self.cylinder_height = cylinder_height
+        self.thickness = thickness
+        self.name = name
+
+        super().__init__(name=name, **kwargs)
+
+    @property
+    def radius(self):
+        return self._radius
+
+    @radius.setter
+    def radius(self, value):
+        if not isinstance(value, (float, int)):
+            raise ValueError("VacuumVessel.radius must be a number. Not", value)
+        if value <= 0:
+            msg = "VacuumVessel.radius must be a positive number above 0. " f"Not {value}"
+            raise ValueError(msg)
+        self._radius = value
+
+    @property
+    def thickness(self):
+        return self._thickness
+
+    @thickness.setter
+    def thickness(self, value):
+        if not isinstance(value, (float, int)):
+            msg = f"VacuumVessel.thickness must be a number. Not {value}"
+            raise ValueError(msg)
+        if value <= 0:
+            msg = f"VacuumVessel.thickness must be a positive number above 0. Not {value}"
+            raise ValueError(msg)
+        self._thickness = value
+
+
+    def create_solid(self):
+        """Creates a rotated 3d solid using points with circular edges.
+
+        Returns:
+           A CadQuery solid: A 3D solid volume
+        """
+
+        #
+        #          -   -
+        #                  -
+        #          -  -       -
+        #                -       -
+        #                  -       -
+        #                    -     |
+        #                     |    |
+        #                     |    |
+        #                     |    |
+        #          c,p        |    |
+        #                     |    |
+        #                     |    |
+        #                     |    |
+        #                    -     |
+        #                  -      -
+        #                -       -
+        #          -  -       -
+        #                  -
+        #          -   -
+        #
+
+        cylinder_section =CenterColumnShieldCylinder(
+            height=self.cylinder_height,
+            inner_radius=self.radius - self.thickness,
+            outer_radius=self.radius,
+            center_height=self.center_point,
+            rotation_angle=self.rotation_angle,
+        )
+
+        upper_dome_section = ConstantThicknessDome(
+            thickness=self.thickness,
+            chord_center_height=self.center_point + 0.5 * self.cylinder_height,
+            chord_width=(self.radius - self.thickness) * 2,
+            chord_height=self.dish_height,
+            upper_or_lower="upper",
+            rotation_angle=self.rotation_angle,
+        )
+
+        lower_dome_section = ConstantThicknessDome(
+            thickness=self.thickness,
+            chord_center_height=self.center_point - 0.5 * self.cylinder_height,
+            chord_width=(self.radius - self.thickness) * 2,
+            chord_height=self.dish_height,
+            upper_or_lower="lower",
+            rotation_angle=self.rotation_angle,
+        )
+
+        upper_dome_section.solid = upper_dome_section.solid.union(cylinder_section.solid)
+        self.solid = lower_dome_section.solid.union(upper_dome_section.solid)
diff --git a/paramak/shape.py b/paramak/shape.py
index ec9647b82..9233e3e44 100644
--- a/paramak/shape.py
+++ b/paramak/shape.py
@@ -1061,16 +1061,9 @@ def export_html(
             facet_splines=facet_splines,
             facet_circles=facet_circles,
             tolerance=tolerance,
-            title=(f"coordinates of {self.__class__.__name__} shape, viewed " "from the {view_plane} plane"),
+            title=f"coordinates of {self.__class__.__name__} shape, viewed from the {view_plane} plane",
         )
 
-        if self.points is not None:
-            fig.add_trace(plotly_trace(points=self.points, mode="markers", name="Shape.points"))
-
-        # sweep shapes have .path_points but not .points attribute
-        if self.path_points:
-            fig.add_trace(plotly_trace(points=self.path_points, mode="markers", name="Shape.path_points"))
-
         if filename is not None:
 
             Path(filename).parents[0].mkdir(parents=True, exist_ok=True)
diff --git a/paramak/utils.py b/paramak/utils.py
index a79d52de8..cf7422c4b 100644
--- a/paramak/utils.py
+++ b/paramak/utils.py
@@ -146,6 +146,16 @@ def diff_between_angles(angle_a: float, angle_b: float) -> float:
     return delta_mod
 
 
+def angle_between_two_points_on_circle(
+    point_1: Tuple[float, float], point_2: Tuple[float, float], radius_of_circle: float
+):
+
+    separation = distance_between_two_points(point_1, point_2)
+    isos_tri_term = (2 * math.pow(radius_of_circle, 2) - math.pow(separation, 2)) / (2 * math.pow(radius_of_circle, 2))
+    angle = math.acos(isos_tri_term)
+    return angle
+
+
 def distance_between_two_points(point_a: Tuple[float, float], point_b: Tuple[float, float]) -> float:
     """Computes the distance between two points.
 
@@ -157,8 +167,8 @@ def distance_between_two_points(point_a: Tuple[float, float], point_b: Tuple[flo
         float: distance between A and B
     """
 
-    xa, ya = point_a
-    xb, yb = point_b
+    xa, ya = point_a[0], point_a[1]
+    xb, yb = point_b[0], point_b[1]
     u_vec = [xb - xa, yb - ya]
     return np.linalg.norm(u_vec)
 
@@ -174,8 +184,8 @@ def extend(point_a: Tuple[float, float], point_b: Tuple[float, float], L: float)
         float, float: point C coordinates
     """
 
-    xa, ya = point_a
-    xb, yb = point_b
+    xa, ya = point_a[0], point_a[1]
+    xb, yb = point_b[0], point_b[1]
     u_vec = [xb - xa, yb - ya]
     u_vec /= np.linalg.norm(u_vec)
 
@@ -270,7 +280,7 @@ def rotate(origin: Tuple[float, float], point: Tuple[float, float], angle: float
     """
 
     ox, oy = origin
-    px, py = point
+    px, py = point[0], point[1]
 
     qx = ox + math.cos(angle) * (px - ox) - math.sin(angle) * (py - oy)
     qy = oy + math.sin(angle) * (px - ox) + math.cos(angle) * (py - oy)
@@ -448,9 +458,9 @@ def plotly_trace(
     text_values = []
 
     for i, point in enumerate(points):
-        text = "point number= {i} <br> x={point[0]} <br> y= {point[1]}"
+        text = f"point number= {i} <br> x={point[0]} <br> y= {point[1]}"
         if len(point) == 3:
-            text = text + "<br> z= {point[2]} <br>"
+            text = text + f"<br> z= {point[2]} <br>"
 
         text_values.append(text)
 


Test patch:
diff --git a/tests/test_parametric_components/test_constant_thickness_dome.py b/tests/test_parametric_components/test_constant_thickness_dome.py
new file mode 100644
index 000000000..d62e6021b
--- /dev/null
+++ b/tests/test_parametric_components/test_constant_thickness_dome.py
@@ -0,0 +1,49 @@
+
+import unittest
+import pytest
+
+import paramak
+
+class TestConstantThicknessDome(unittest.TestCase):
+    """tests for the ConstantThicknessDome class"""
+    def test_volume_increases_with_rotation_angle(self):
+        """Tests that the volume doubles when rotation angle doubles"""
+
+        test_shape_1 = paramak.ConstantThicknessDome(rotation_angle=180)
+        test_shape_2 = paramak.ConstantThicknessDome(rotation_angle=360)
+        
+        assert test_shape_1.volume() *2 == pytest.approx(test_shape_2.volume())
+
+    def test_upper_lower_flips_points(self):
+        """Checks that the coords of the flips version are the same for p1 and p2
+        and negative for part of p3"""
+        test_shape_1 = paramak.ConstantThicknessDome(upper_or_lower="upper")
+        test_shape_2 = paramak.ConstantThicknessDome(upper_or_lower="lower")
+        assert test_shape_1.points[0] == test_shape_2.points[0]
+        assert test_shape_1.points[1] == test_shape_2.points[1]
+        assert test_shape_1.points[2][0] == test_shape_2.points[2][0]
+        assert test_shape_1.points[2][1] == -test_shape_2.points[2][1]
+
+        assert test_shape_1.volume() == pytest.approx(test_shape_2.volume())
+
+    def test_invalid_parameters_errors(self):
+        """Checks that the correct errors are raised when invalid arguments are
+        input as shape parameters."""
+
+        def incorrect_shape_height_width_ratio():
+            my_shape = paramak.ConstantThicknessDome(chord_width=10, chord_height=40)
+            my_shape.solid
+
+        def incorrect_thickness():
+            paramak.ConstantThicknessDome(thickness=-1)
+
+        def incorrect_chord_height():
+            paramak.ConstantThicknessDome(chord_height=-1)
+
+        def incorrect_chord_width():
+            paramak.ConstantThicknessDome(chord_width=-1)
+
+        self.assertRaises(ValueError, incorrect_shape_height_width_ratio)
+        self.assertRaises(ValueError, incorrect_thickness)
+        self.assertRaises(ValueError, incorrect_chord_height)
+        self.assertRaises(ValueError, incorrect_chord_width)
diff --git a/tests/test_parametric_components/test_dished_vacuum_vessel.py b/tests/test_parametric_components/test_dished_vacuum_vessel.py
new file mode 100644
index 000000000..d30bc76ca
--- /dev/null
+++ b/tests/test_parametric_components/test_dished_vacuum_vessel.py
@@ -0,0 +1,10 @@
+
+import pytest
+
+import paramak
+
+
+def test_volume_increases_with_rotation_angle():
+    test_shape_1 = paramak.DishedVacuumVessel(rotation_angle=180)
+    test_shape_2 = paramak.DishedVacuumVessel(rotation_angle=360)
+    assert test_shape_1.volume() *2 == pytest.approx(test_shape_2.volume())