Merge pull request #151 from svalinn/filament-data-refactor

Examples/radial_build_distance_example.py

@@ -1,29 +1,34 @@
 from parastell.radial_distance_utils import *
 
-coils_file = "coils_wistelld.txt"
-circ_cross_section = ["circle", 25]
+coils_file = "coils.example"
+width = 40.0
+thickness = 50.0
 toroidal_extent = 90.0
-vmec_file = "plasma_wistelld.nc"
+vmec_file = "wout_vmec.nc"
 vmec = read_vmec.VMECData(vmec_file)
 
 magnet_set = magnet_coils.MagnetSet(
-    coils_file, circ_cross_section, toroidal_extent
+    coils_file, width, thickness, toroidal_extent
 )
 
 reordered_filaments = get_reordered_filaments(magnet_set)
 
 build_magnet_surface(reordered_filaments)
 
-toroidal_angles = np.linspace(0, 90, 72)
-poloidal_angles = np.linspace(0, 360, 96)
+toroidal_angles = np.linspace(0, 90, num=4)
+poloidal_angles = np.linspace(0, 360, num=4)
 wall_s = 1.08
 
 radial_build_dict = {}
 
 radial_build = ivb.RadialBuild(
-    toroidal_angles, poloidal_angles, wall_s, radial_build_dict
+    toroidal_angles,
+    poloidal_angles,
+    wall_s,
+    radial_build_dict,
+    split_chamber=False,
 )
-build = ivb.InVesselBuild(vmec, radial_build, split_chamber=False)
+build = ivb.InVesselBuild(vmec, radial_build, num_ribs=72, num_rib_pts=96)
 build.populate_surfaces()
 build.calculate_loci()
 ribs = build.Surfaces["chamber"].Ribs

parastell/magnet_coils.py

@@ -8,7 +8,7 @@
 
 from . import log
 from . import cubit_io as cubit_io
-from .utils import normalize, read_yaml_config, filter_kwargs, m2cm
+from .utils import read_yaml_config, filter_kwargs, m2cm
 
 export_allowed_kwargs = ["step_filename", "export_mesh", "mesh_filename"]
 
@@ -112,15 +112,16 @@ def logger(self):
     def logger(self, logger_object):
         self._logger = log.check_init(logger_object)
 
-    def _extract_filament_data(self):
-        """Extracts filament coordinate numerical data from input data file.
+    def _instantiate_coils(self):
+        """Extracts filament coordinate data from input data file and
+        instantiates MagnetCoil class objects.
         (Internal function not intended to be called externally)
         """
         with open(self.coils_file, "r") as file:
             data = file.readlines()[self.start_line :]
 
         coords = []
-        filament_coords = []
+        self.magnet_coils = []
 
         for line in data:
             columns = line.strip().split()
@@ -139,11 +140,17 @@ def _extract_filament_data(self):
 
             else:
                 coords.append(coords[0])
-                filament_coords.append(np.array(coords))
+                self.magnet_coils.append(
+                    MagnetCoil(
+                        np.array(coords),
+                        np.average(coords[:-1], axis=0),
+                        self._width,
+                        self._thickness,
+                        self.sample_mod,
+                    )
+                )
                 coords.clear()
 
-        self.filament_coords = filament_coords
-
     def _compute_radial_distance_data(self):
         """Computes average and maximum radial distance of filament points.
         (Internal function not intended to be called externally)
@@ -152,8 +159,8 @@ def _compute_radial_distance_data(self):
         self.average_radial_distance = 0
         self.max_radial_distance = -1
 
-        for f in self.filament_coords:
-            radii = np.linalg.norm(f[:, :2], axis=1)
+        for coil in self.magnet_coils:
+            radii = np.linalg.norm(coil.coords[:-1, :2], axis=1)
             radii_count += len(radii)
             self.average_radial_distance += np.sum(radii)
             self.max_radial_distance = max(
@@ -162,18 +169,12 @@ def _compute_radial_distance_data(self):
 
         self.average_radial_distance /= radii_count
 
-    def _filter_filaments(self):
-        """Cleans filament data such that only filaments within the toroidal
-        extent of the model are included and filaments are sorted by toroidal
-        angle.
+    def _filter_coils(self):
+        """Filters list of MagnetCoil objects such that only those within the
+        toroidal extent of the model are included and coils are sorted by
+        center-of-mass toroidal angle.
         (Internal function not intended to be called externally)
         """
-        # Initialize coordinate data for filaments within toroidal extent of
-        # model
-        filtered_coords = []
-        # Initialize list of filament centers of mass
-        com_list = []
-
         # Define tolerance of toroidal extent to account for dimensionality of
         # coil cross-section
         # Multiply by factor of 2 to be conservative
@@ -183,39 +184,25 @@ def _filter_filaments(self):
         lower_bound = 2 * np.pi - tol
         upper_bound = self._toroidal_extent + tol
 
-        for coords in self.filament_coords:
-            # Compute filament center of mass
-            com = np.average(coords, axis=0)
-            # Compute toroidal angle of each point in filament
-            toroidal_angles = np.arctan2(coords[:, 1], coords[:, 0])
-            # Ensure angles are positive
-            toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi)
-            # Compute bounds of toroidal extent of filament
-            min_tor_ang = np.min(toroidal_angles)
-            max_tor_ang = np.max(toroidal_angles)
-
-            # Determine if filament toroidal extent overlaps with that of model
-            if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or (
-                max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
-            ):
-                filtered_coords.append(coords)
-                com_list.append(com)
-
-        filtered_coords = np.array(filtered_coords)
-        com_list = np.array(com_list)
+        # Create filter determining whether each coil lies within model's
+        # toroidal extent
+        filtered_coils = [
+            coil
+            for coil in self.magnet_coils
+            if coil.in_toroidal_extent(lower_bound, upper_bound)
+        ]
 
         # Compute toroidal angles of filament centers of mass
+        com_list = np.array([coil.center_of_mass for coil in filtered_coils])
         com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0])
+        # Ensure angles are positive
         com_toroidal_angles = (com_toroidal_angles + 2 * np.pi) % (2 * np.pi)
 
-        # Sort filaments by toroidal angle and overwrite coordinate data
-        self.filament_coords = np.array(
-            [x for _, x in sorted(zip(com_toroidal_angles, filtered_coords))]
-        )
-        # Store filtered centers of mass since they have already been computed
-        self.filament_com = np.array(
-            [x for _, x in sorted(zip(com_toroidal_angles, com_list))]
-        )
+        # Sort coils by center-of-mass toroidal angle and overwrite stored list
+        self.magnet_coils = [
+            coil
+            for _, coil in sorted(zip(com_toroidal_angles, filtered_coils))
+        ]
 
     def _cut_magnets(self):
         """Cuts the magnets at the planes defining the toriodal extent.
@@ -239,29 +226,22 @@ def _cut_magnets(self):
             cut_coil = coil.solid.intersect(toroidal_region)
             coil.solid = cut_coil
 
+    def populate_magnet_coils(self):
+        """Populates MagnetCoil class objects representing each of the magnetic
+        coils that lie within the specified toroidal extent.
+        """
+        self._logger.info("Populating magnet coils...")
+
+        self._instantiate_coils()
+        self._compute_radial_distance_data()
+        self._filter_coils()
+
     def build_magnet_coils(self):
         """Builds each filament in self.filtered_filaments in cubit, then cuts
         to the toroidal extent using self._cut_magnets().
         """
         self._logger.info("Constructing magnet coils...")
 
-        self._extract_filament_data()
-        self._compute_radial_distance_data()
-        self._filter_filaments()
-
-        self.magnet_coils = [
-            MagnetCoil(
-                coords,
-                center_of_mass,
-                self._width,
-                self._thickness,
-                self.sample_mod,
-            )
-            for coords, center_of_mass in zip(
-                self.filament_coords, self.filament_com
-            )
-        ]
-
         [magnet_coil.create_magnet() for magnet_coil in self.magnet_coils]
 
         self._cut_magnets()
@@ -349,57 +329,89 @@ def coords(self):
 
     @coords.setter
     def coords(self, data):
+        self._coords = data
+
         # Compute tangents
         tangents = np.subtract(
             np.append(data[1:], [data[1]], axis=0),
             np.append([data[-2]], data[0:-1], axis=0),
         )
-        tangents = tangents / np.linalg.norm(tangents, axis=1)[:, np.newaxis]
+        self.tangents = (
+            tangents / np.linalg.norm(tangents, axis=1)[:, np.newaxis]
+        )
 
-        # Sample filament coordinates and tangents by modifier
-        self._coords = data[0 : -1 : self.sample_mod]
-        self._coords = np.append(self._coords, [self._coords[0]], axis=0)
-        self.tangents = tangents[0 : -1 : self.sample_mod]
-        self.tangents = np.append(self.tangents, [self.tangents[0]], axis=0)
+    def in_toroidal_extent(self, lower_bound, upper_bound):
+        """Determines if the coil lies within a given toroidal angular extent,
+        based on filament coordinates.
+
+        Arguments:
+            lower_bound (float): lower bound of toroidal extent [rad].
+            upper_bound (float): upper bound of toroidal extent [rad].
+
+        Returns:
+            in_toroidal_extent (bool): flag to indicate whether coil lies
+                within toroidal bounds.
+        """
+        # Compute toroidal angle of each point in filament
+        toroidal_angles = np.arctan2(self._coords[:, 1], self._coords[:, 0])
+        # Ensure angles are positive
+        toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi)
+        # Compute bounds of toroidal extent of filament
+        min_tor_ang = np.min(toroidal_angles)
+        max_tor_ang = np.max(toroidal_angles)
+
+        # Determine if filament toroidal extent overlaps with that of model
+        if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or (
+            max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
+        ):
+            in_toroidal_extent = True
+        else:
+            in_toroidal_extent = False
+
+        return in_toroidal_extent
 
     def create_magnet(self):
         """Creates a single magnet coil CAD solid in CadQuery.
 
         Returns:
             coil (object): cq.Solid object representing a single magnet coil.
         """
-        tangent_vectors = [
-            cq.Vector(tuple(tangent)) for tangent in self.tangents
-        ]
+        # Sample filament coordinates and tangents by modifier
+        coords = self._coords[0 : -1 : self.sample_mod]
+        coords = np.append(coords, [self._coords[0]], axis=0)
+        tangents = self.tangents[0 : -1 : self.sample_mod]
+        tangents = np.append(tangents, [self.tangents[0]], axis=0)
+
+        tangent_vectors = [cq.Vector(tuple(tangent)) for tangent in tangents]
 
         # Define coil filament path normals such that they face the filament
         # center of mass
         # Compute "outward" direction as difference between filament positions
         # and filament center of mass
-        outward_dirs = self._coords - self.center_of_mass
+        outward_dirs = coords - self.center_of_mass
         outward_dirs = (
             outward_dirs / np.linalg.norm(outward_dirs, axis=1)[:, np.newaxis]
         )
 
         # Project outward directions onto desired coil cross-section (CS) plane
         # at each filament position to define filament path normals
         parallel_parts = np.diagonal(
-            np.matmul(outward_dirs, self.tangents.transpose())
+            np.matmul(outward_dirs, tangents.transpose())
         )
 
-        normals = outward_dirs - parallel_parts[:, np.newaxis] * self.tangents
+        normals = outward_dirs - parallel_parts[:, np.newaxis] * tangents
         normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis]
 
         # Compute binormals projected onto CS plane at each position
-        binormals = np.cross(self.tangents, normals)
+        binormals = np.cross(tangents, normals)
 
         # Compute coordinates of edges of rectangular coils
         edge_offsets = np.array([[-1, -1], [-1, 1], [1, 1], [1, -1]])
 
         coil_edge_coords = []
         for edge_offset in edge_offsets:
             coil_edge = (
-                self._coords
+                coords
                 + edge_offset[0] * binormals * (self.width / 2)
                 + edge_offset[1] * normals * (self.thickness / 2)
             )

parastell/parastell.py

@@ -224,6 +224,7 @@ def construct_magnets(
             **kwargs,
         )
 
+        self.magnet_set.populate_magnet_coils()
         self.magnet_set.build_magnet_coils()
 
     def export_magnets(

parastell/radial_distance_utils.py

@@ -90,9 +90,10 @@ def reorder_filaments(filaments):
             reordered_filament = np.flip(reordered_filament, axis=0)
 
         # ensure filament is a closed loop
-        reordered_filament = np.concatenate(
-            [reordered_filament, [reordered_filament[0]]]
-        )
+        if max_z_index != 0:
+            reordered_filament = np.concatenate(
+                [reordered_filament, [reordered_filament[0]]]
+            )
 
         filaments[filament_index] = reordered_filament
 
@@ -105,11 +106,11 @@ def get_reordered_filaments(magnet_set):
     """
     Convenience function to get the reordered filament data from a magnet
     """
-    magnet_set._extract_filaments()
-    magnet_set._set_average_radial_distance()
-    magnet_set._set_filtered_filaments()
+    magnet_set.populate_magnet_coils()
 
-    filtered_filaments = magnet_set.filtered_filaments
+    filtered_filaments = np.array(
+        [coil.coords for coil in magnet_set.magnet_coils]
+    )
     filaments = reorder_filaments(filtered_filaments)
 
     return filaments

tests/test_magnet_coils.py

@@ -42,16 +42,16 @@ def test_magnet_construction(coil_set):
     thickness_exp = 50.0
     toroidal_extent_exp = np.deg2rad(90.0)
     max_cs_len_exp = 50.0
-    average_radial_distance_exp = 1028.5044183872055
+    average_radial_distance_exp = 1023.7170384211436
     max_radial_distance_exp = 1646.3258131460148
-    len_filaments_exp = 1
-    len_coords_exp = 14
+    len_coils_exp = 1
+    len_coords_exp = 129
 
     remove_files()
 
-    coil_set.build_magnet_coils()
+    coil_set.populate_magnet_coils()
 
-    assert len(coil_set.filament_coords) == len_filaments_exp
+    assert len(coil_set.magnet_coils) == len_coils_exp
     assert coil_set.width == width_exp
     assert coil_set.thickness == thickness_exp
     assert coil_set.toroidal_extent == toroidal_extent_exp
@@ -69,6 +69,7 @@ def test_magnet_exports(coil_set):
 
     remove_files()
 
+    coil_set.populate_magnet_coils()
     coil_set.build_magnet_coils()
     coil_set.export_step()
     assert Path("magnet_set.step").exists()