Feature: polyhedral cross-section current sources and a general VIM formulation of 3-D magnetostatics

bluemira/magnetostatics/baseclass.py

@@ -21,9 +21,10 @@
 
 from bluemira.geometry.bound_box import BoundingBox
 from bluemira.geometry.coordinates import rotation_matrix
+from bluemira.magnetostatics.error import MagnetostaticsError
 from bluemira.utilities.plot_tools import Plot3D
 
-__all__ = ["CurrentSource", "RectangularCrossSectionCurrentSource", "SourceGroup"]
+__all__ = ["CurrentSource", "CrossSectionCurrentSource", "SourceGroup"]
 
 
 class CurrentSource(ABC):
@@ -99,17 +100,16 @@
         return deepcopy(self)
 
 
-class RectangularCrossSectionCurrentSource(CurrentSource):
+class CrossSectionCurrentSource(CurrentSource):
     """
-    Abstract base class for a current source with a rectangular cross-section.
+    Abstract class for a current source with a cross-section
     """
 
-    origin: np.array
-    dcm: np.array
-    points: np.array
-    breadth: float
-    depth: float
-    length: float
+    _origin: np.array
+    _dcm: np.array
+    _points: np.array
+    _rho: float
+    _area: float
 
     def set_current(self, current: float):
         """
@@ -121,7 +121,7 @@
             The current of the source [A]
         """
         super().set_current(current)
-        self.rho = current / (4 * self.breadth * self.depth)
+        self._rho = current / self._area
 
     def rotate(self, angle: float, axis: Union[np.ndarray, str]):
         """
@@ -135,21 +135,21 @@
             The axis of rotation
         """
         r = rotation_matrix(np.deg2rad(angle), axis).T
-        self.origin = self.origin @ r
-        self.points = np.array([p @ r for p in self.points], dtype=object)
-        self.dcm = self.dcm @ r
+        self._origin = self._origin @ r
+        self._points = np.array([p @ r for p in self._points], dtype=object)
+        self._dcm = self._dcm @ r
 
     def _local_to_global(self, points: np.ndarray) -> np.ndarray:
         """
         Convert local x', y', z' point coordinates to global x, y, z point coordinates.
         """
-        return np.array([self.origin + self.dcm.T @ p for p in points])
+        return np.array([self._origin + self._dcm.T @ p for p in points])
 
     def _global_to_local(self, points: np.ndarray) -> np.ndarray:
         """
         Convert global x, y, z point coordinates to local x', y', z' point coordinates.
         """
-        return np.array([(self.dcm @ (p - self.origin)) for p in points])
+        return np.array([(self._dcm @ (p - self._origin)) for p in points])
 
     def plot(self, ax: Optional[Axes] = None, show_coord_sys: bool = False):
         """
@@ -166,21 +166,87 @@
             ax = Plot3D()
             # If no ax provided, we assume that we want to plot only this source,
             # and thus set aspect ratio equality on this term only
-            edge_points = np.concatenate(self.points)
+            edge_points = np.concatenate(self._points)
 
             # Invisible bounding box to set equal aspect ratio plot
             xbox, ybox, zbox = BoundingBox.from_xyz(*edge_points.T).get_box_arrays()
             ax.plot(1.1 * xbox, 1.1 * ybox, 1.1 * zbox, "s", alpha=0)
 
-        for points in self.points:
+        for points in self._points:
             ax.plot(*points.T, color="b", linewidth=1)
 
         # Plot local coordinate system
         if show_coord_sys:
-            ax.scatter(*self.origin, color="k")
-            ax.quiver(*self.origin, *self.dcm[0], length=self.breadth, color="r")
-            ax.quiver(*self.origin, *self.dcm[1], length=self.length, color="r")
-            ax.quiver(*self.origin, *self.dcm[2], length=self.depth, color="r")
+            ax.scatter(*self._origin, color="k")
+            ax.quiver(*self._origin, *self._dcm[0], length=1, color="r")
+            ax.quiver(*self._origin, *self._dcm[1], length=1, color="r")
+            ax.quiver(*self._origin, *self._dcm[2], length=1, color="r")
+
+
+class PolyhedralCrossSectionCurrentSource(CrossSectionCurrentSource):
+    """
+    Abstract base class for a current source with a polyhedral cross-section.
+    """
+
+    _face_points: np.ndarray
+    _face_normals: np.ndarray
+    _mid_points: np.ndarray
+
+    def rotate(self, angle: float, axis: Union[np.ndarray, str]):
+        """
+        Rotate the CurrentSource about an axis.
+
+        Parameters
+        ----------
+        angle:
+            The rotation degree [degree]
+        axis:
+            The axis of rotation
+        """
+        super().rotate(angle, axis)
+        r = rotation_matrix(np.deg2rad(angle), axis).T
+        self._face_normals = np.array([n @ r for n in self._face_normals])
+        self._face_points = np.array([p @ r for p in self._face_points])
+        self._mid_points = np.array([p @ r for p in self._mid_points])
+
+
+class PrismEndCapMixin:
+    def _check_angle_values(self, alpha, beta):
+        """
+        Check that end-cap angles are acceptable.
+        """
+        sign_alpha = np.sign(alpha)
+        sign_beta = np.sign(beta)
+        one_zero = np.any(np.array([sign_alpha, sign_beta]) == 0)
+        if not one_zero and sign_alpha != sign_beta:
+            raise MagnetostaticsError(
+                f"{self.__class__.__name__} instantiation error: end-cap angles "
+                f"must have the same sign {alpha=:.3f}, {beta=:.3f}."
+            )
+        if not (0 <= abs(alpha) < 0.5 * np.pi):
+            raise MagnetostaticsError(
+                f"{self.__class__.__name__} instantiation error: {alpha=:.3f} is outside"
+                " bounds of [0, 180°)."
+            )
+        if not (0 <= abs(beta) < 0.5 * np.pi):
+            raise MagnetostaticsError(
+                f"{self.__class__.__name__} instantiation error: {beta=:.3f} is outside "
+                "bounds of [0, 180°)."
+            )
+
+    def _check_raise_self_intersection(
+        self, length: float, breadth: float, alpha: float, beta: float
+    ):
+        """
+        Check for bad combinations of source length and end-cap angles.
+        """
+        a = np.tan(alpha) * breadth
+        b = np.tan(beta) * breadth
+        if (a + b) > length:
+            raise MagnetostaticsError(
+                f"{self.__class__.__name__} instantiation error: source length and "
+                "angles imply a self-intersecting trapezoidal prism."
+            )
 
 
 class SourceGroup(ABC):
@@ -189,11 +255,11 @@
     """
 
     sources: List[CurrentSource]
-    points: np.array
+    _points: np.array
 
     def __init__(self, sources: List[CurrentSource]):
         self.sources = sources
-        self.points = np.vstack([np.vstack(s.points) for s in self.sources])
+        self._points = np.vstack([np.vstack(s._points) for s in self.sources])
 
     def set_current(self, current: float):
         """
@@ -244,7 +310,7 @@
         """
         for source in self.sources:
             source.rotate(angle, axis)
-        self.points = self.points @ rotation_matrix(angle, axis)
+        self._points = self._points @ rotation_matrix(angle, axis)
 
     def plot(self, ax: Optional[Axes] = None, show_coord_sys: bool = False):
         """
@@ -261,7 +327,7 @@
             ax = Plot3D()
 
         # Invisible bounding box to set equal aspect ratio plot
-        xbox, ybox, zbox = BoundingBox.from_xyz(*self.points.T).get_box_arrays()
+        xbox, ybox, zbox = BoundingBox.from_xyz(*self._points.T).get_box_arrays()
         ax.plot(1.1 * xbox, 1.1 * ybox, 1.1 * zbox, "s", alpha=0)
 
         for source in self.sources:

bluemira/magnetostatics/biot_savart.py

@@ -79,12 +79,12 @@
                 self.length_scale = np.min(lengths)
 
         # Assemble arrays and vector
-        self.d_l = np.vstack(d_ls)
-        self.d_l_hat = np.linalg.norm(self.d_l, axis=1)
-        self.mid_points = np.vstack(mids_points)
-        self.points = np.vstack(points)
+        self._d_l = np.vstack(d_ls)
+        self._d_l_hat = np.linalg.norm(self._d_l, axis=1)
+        self._mid_points = np.vstack(mids_points)
+        self._points = np.vstack(points)
         self._arrays = arrays
-        self.radius = radius
+        self._radius = radius
         self.current = current
 
     @staticmethod
@@ -124,16 +124,16 @@
         The vector potential at the point due to the arbitrarily shaped Coordinates
         """
         point = np.array([x, y, z])
-        r = point - self.points
+        r = point - self._points
         r_mag = tools.norm(r, axis=1)
         r_mag[r_mag < EPS] = EPS
-        core = r_mag / self.radius
-        core[r_mag > self.radius] = 1
+        core = r_mag / self._radius
+        core[r_mag > self._radius] = 1
 
         # The below einsum operation is equivalent to:
         # self.current * np.sum(core * self.d_l.T / r_mag, axis=0) / (4 * np.pi)
         return np.einsum(
-            "i, ji, ... -> j", core, self.d_l / r_mag[None], ONE_4PI * self.current
+            "i, ji, ... -> j", core, self._d_l / r_mag[None], ONE_4PI * self.current
         )
 
     @process_xyz_array
@@ -167,18 +167,18 @@
         smoothing. Do not use for values near the coil current centreline.
         """  # noqa: W505, E501
         point = np.array([x, y, z])
-        r = point - self.mid_points
+        r = point - self._mid_points
         r3 = np.linalg.norm(r, axis=1) ** 3
 
-        ds = np.cross(self.d_l, r)
+        ds = np.cross(self._d_l, r)
 
         # Coil core correction
-        d_l_hat = self.d_l_hat[:, None]
+        d_l_hat = self._d_l_hat[:, None]
         ds_mag = np.linalg.norm(ds / d_l_hat, axis=1)
         ds_mag = np.tile(ds_mag, (3, 1)).T
         ds_mag[ds_mag < EPS] = EPS
-        core = ds_mag**2 / self.radius**2
-        core[ds_mag > self.radius] = 1
+        core = ds_mag**2 / self._radius**2
+        core[ds_mag > self._radius] = 1
         return MU_0_4PI * self.current * np.sum(core * ds / r3[:, np.newaxis], axis=0)
 
     def inductance(self) -> float:
@@ -203,16 +203,16 @@
         inductance = 0
         for _i, (x1, dx1) in enumerate(zip(self.ref_mid_points, self.ref_d_l)):
             # We create a mask to drop the point where x1 == x2
-            r = x1 - self.mid_points
+            r = x1 - self._mid_points
             mask = np.sum(r**2, axis=1) > 0.5 * self.length_scale
             inductance += np.sum(
-                np.dot(dx1, self.d_l[mask].T) / np.linalg.norm(r[mask], axis=1)
+                np.dot(dx1, self._d_l[mask].T) / np.linalg.norm(r[mask], axis=1)
             )
 
         # Self-inductance correction (Y = 0.5 for homogenous current distribution)
         # Equation 6 of https://arxiv.org/pdf/1204.1486.pdf
         error_tail = 0
-        a, b = self.radius, 0.5 * self.length_scale
+        a, b = self._radius, 0.5 * self.length_scale
         if b > 10 * a:
             # Equation A.4 of https://arxiv.org/pdf/1204.1486.pdf
             error_tail = a**2 / b**2 - 3 / (8 * b**4) * (a**4 - 2 * a**2)
@@ -232,9 +232,9 @@
             The axis of rotation
         """
         r = rotation_matrix(np.deg2rad(angle), axis).T
-        self.points = self.points @ r
-        self.d_l = self.d_l @ r
-        self.mid_points = self.mid_points @ r
+        self._points = self._points @ r
+        self._d_l = self._d_l @ r
+        self._mid_points = self._mid_points @ r
         self.ref_d_l = self.ref_d_l @ r
         self.ref_mid_points = self.ref_mid_points @ r
         self._arrays = [array @ r for array in self._arrays]
@@ -255,7 +255,7 @@
             # If no ax provided, we assume that we want to plot only this source,
             # and thus set aspect ratio equality on this term only
             # Invisible bounding box to set equal aspect ratio plot
-            xbox, ybox, zbox = BoundingBox.from_xyz(*self.points.T).get_box_arrays()
+            xbox, ybox, zbox = BoundingBox.from_xyz(*self._points.T).get_box_arrays()
             ax.plot(1.1 * xbox, 1.1 * ybox, 1.1 * zbox, "s", alpha=0)
 
         for array in self._arrays: