From fb4960d0fc0ab66318ac9263f2b39771100250c8 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 1 Nov 2024 14:26:15 -0500 Subject: [PATCH 01/11] Implement alternation flag --- parastell/source_mesh.py | 49 ++++++++++++++++++++++++++++++---------- 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index ae10f06..b5b376f 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -375,13 +375,25 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): # first, ordered clockwise relative to the thumb, followed by the # remaining vertex at the end of the thumb. # See Moreno, Bader, Wilson 2024 for hexahedron splitting - hex_canon_ids = [ - [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], - [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], - [idx_list[0], idx_list[7], idx_list[5], idx_list[4]], - [idx_list[7], idx_list[2], idx_list[5], idx_list[6]], - [idx_list[0], idx_list[2], idx_list[5], idx_list[7]], - ] + # Conditionally alternate ordering of vertices defining hexahedron + # splitting to avoid gaps and overlaps between non-planar hexahedron + # faces + if self.alt_flag: + hex_canon_ids = [ + [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], + [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], + [idx_list[0], idx_list[7], idx_list[5], idx_list[4]], + [idx_list[7], idx_list[2], idx_list[5], idx_list[6]], + [idx_list[0], idx_list[2], idx_list[5], idx_list[7]], + ] + else: + hex_canon_ids = [ + [idx_list[0], idx_list[3], idx_list[1], idx_list[4]], + [idx_list[1], idx_list[3], idx_list[2], idx_list[6]], + [idx_list[1], idx_list[4], idx_list[6], idx_list[5]], + [idx_list[3], idx_list[6], idx_list[4], idx_list[7]], + [idx_list[1], idx_list[3], idx_list[6], idx_list[4]], + ] for vertex_ids in hex_canon_ids: self._create_tet(vertex_ids) @@ -420,11 +432,20 @@ def _create_tets_from_wedge(self, theta_idx, phi_idx): # first, ordered clockwise relative to the thumb, followed by the # remaining vertex at the end of the thumb. # See Moreno, Bader, Wilson 2024 for wedge splitting - wedge_canon_ids = [ - [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], - [idx_list[3], idx_list[2], idx_list[4], idx_list[5]], - [idx_list[3], idx_list[2], idx_list[1], idx_list[4]], - ] + # Conditionally alternate ordering of vertices defining wedge splitting + # to avoid gaps and overlaps between non-planar wedge faces + if self.alt_flag: + wedge_canon_ids = [ + [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], + [idx_list[1], idx_list[3], idx_list[5], idx_list[4]], + [idx_list[1], idx_list[3], idx_list[2], idx_list[5]], + ] + else: + wedge_canon_ids = [ + [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], + [idx_list[3], idx_list[2], idx_list[4], idx_list[5]], + [idx_list[3], idx_list[2], idx_list[1], idx_list[4]], + ] for vertex_ids in wedge_canon_ids: self._create_tet(vertex_ids) @@ -437,14 +458,18 @@ def create_mesh(self): self.mbc.add_entity(self.mesh_set, self.verts) for phi_idx in range(self.num_phi - 1): + # Set alternation flag to true at beginning of each toroidal block + self.alt_flag = True # Create tetrahedra for wedges at center of plasma for theta_idx in range(1, self.num_theta): self._create_tets_from_wedge(theta_idx, phi_idx) + self.alt_flag = not self.alt_flag # Create tetrahedra for hexahedra beyond center of plasma for s_idx in range(self.num_s - 2): for theta_idx in range(1, self.num_theta): self._create_tets_from_hex(s_idx, theta_idx, phi_idx) + self.alt_flag = not self.alt_flag def export_mesh(self, filename="source_mesh", export_dir=""): """Use PyMOAB interface to write source mesh with source strengths From a2ae20643fc7f9a15cc717066132a246fa09ca4e Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Sat, 2 Nov 2024 12:31:46 -0500 Subject: [PATCH 02/11] Change alternate scheme flag to local variable and change name --- parastell/source_mesh.py | 54 +++++++++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index b5b376f..9a1994a 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -337,12 +337,21 @@ def _get_vertex_id(self, vertex_idx): return id - def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): + def _create_tets_from_hex( + self, s_idx, theta_idx, phi_idx, alternate_scheme + ): """Creates five tetrahedra from defined hexahedron. (Internal function not intended to be called externally) Arguments: - idx_list (list of int): list of hexahedron vertex indices. + s_idx (int): index defining location along CFS axis. + theta_idx (int): index defining location along poloidal angle axis. + phi_idx (int): index defining location along toroidal angle axis. + alternate_scheme (bool): flag indicating whether alternate + hexahedron division scheme should be used. + + Returns: + alternate_scheme (bool): switched flag. """ # relative offsets of vertices in a 3-D index space @@ -375,10 +384,9 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): # first, ordered clockwise relative to the thumb, followed by the # remaining vertex at the end of the thumb. # See Moreno, Bader, Wilson 2024 for hexahedron splitting - # Conditionally alternate ordering of vertices defining hexahedron - # splitting to avoid gaps and overlaps between non-planar hexahedron - # faces - if self.alt_flag: + # Conditionally use alternate scheme defining hexahedron splitting to + # avoid gaps and overlaps between non-planar hexahedron faces + if alternate_scheme: hex_canon_ids = [ [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], @@ -398,12 +406,20 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): for vertex_ids in hex_canon_ids: self._create_tet(vertex_ids) - def _create_tets_from_wedge(self, theta_idx, phi_idx): + return not alternate_scheme + + def _create_tets_from_wedge(self, theta_idx, phi_idx, alternate_scheme): """Creates three tetrahedra from defined wedge. (Internal function not intended to be called externally) Arguments: - idx_list (list of int): list of wedge vertex indices. + theta_idx (int): index defining location along poloidal angle axis. + phi_idx (int): index defining location along toroidal angle axis. + alternate_scheme (bool): flag indicating whether alternate + hexahedron division scheme should be used. + + Returns: + alternate_scheme (bool): switched flag. """ # relative offsets of wedge vertices in a 3-D index space @@ -434,7 +450,7 @@ def _create_tets_from_wedge(self, theta_idx, phi_idx): # See Moreno, Bader, Wilson 2024 for wedge splitting # Conditionally alternate ordering of vertices defining wedge splitting # to avoid gaps and overlaps between non-planar wedge faces - if self.alt_flag: + if alternate_scheme: wedge_canon_ids = [ [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], [idx_list[1], idx_list[3], idx_list[5], idx_list[4]], @@ -450,6 +466,8 @@ def _create_tets_from_wedge(self, theta_idx, phi_idx): for vertex_ids in wedge_canon_ids: self._create_tet(vertex_ids) + return not alternate_scheme + def create_mesh(self): """Creates volumetric source mesh in real space.""" self._logger.info("Constructing source mesh...") @@ -458,18 +476,24 @@ def create_mesh(self): self.mbc.add_entity(self.mesh_set, self.verts) for phi_idx in range(self.num_phi - 1): - # Set alternation flag to true at beginning of each toroidal block - self.alt_flag = True + # Initialize alternate scheme flag at beginning of each toroidal + # block + alternate_scheme = False # Create tetrahedra for wedges at center of plasma for theta_idx in range(1, self.num_theta): - self._create_tets_from_wedge(theta_idx, phi_idx) - self.alt_flag = not self.alt_flag + alternate_scheme = self._create_tets_from_wedge( + theta_idx, phi_idx, alternate_scheme + ) # Create tetrahedra for hexahedra beyond center of plasma for s_idx in range(self.num_s - 2): + # Initialize alternate scheme flag at beginning of each + # CFS block + alternate_scheme = False for theta_idx in range(1, self.num_theta): - self._create_tets_from_hex(s_idx, theta_idx, phi_idx) - self.alt_flag = not self.alt_flag + alternate_scheme = self._create_tets_from_hex( + s_idx, theta_idx, phi_idx, alternate_scheme + ) def export_mesh(self, filename="source_mesh", export_dir=""): """Use PyMOAB interface to write source mesh with source strengths From d99f5fe839c4b5f6ec47230f1711663f4ecb9f5c Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Sat, 2 Nov 2024 12:42:55 -0500 Subject: [PATCH 03/11] Introduce additional alternate scheme flags to ensure face alignment between adjacent toroidal and CFS blocks --- parastell/source_mesh.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index 9a1994a..c30d784 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -475,26 +475,38 @@ def create_mesh(self): self.mesh_set = self.mbc.create_meshset() self.mbc.add_entity(self.mesh_set, self.verts) + # Initialize alternate scheme flag for adjacent toroidal blocks + toroidal_alt_scheme = False + for phi_idx in range(self.num_phi - 1): # Initialize alternate scheme flag at beginning of each toroidal # block - alternate_scheme = False + alternate_scheme = toroidal_alt_scheme + # Create tetrahedra for wedges at center of plasma for theta_idx in range(1, self.num_theta): alternate_scheme = self._create_tets_from_wedge( theta_idx, phi_idx, alternate_scheme ) + # Initialize alternate scheme flag for adjacent CFS blocks + cfs_alt_scheme = not toroidal_alt_scheme + # Create tetrahedra for hexahedra beyond center of plasma for s_idx in range(self.num_s - 2): # Initialize alternate scheme flag at beginning of each # CFS block - alternate_scheme = False + alternate_scheme = cfs_alt_scheme + for theta_idx in range(1, self.num_theta): alternate_scheme = self._create_tets_from_hex( s_idx, theta_idx, phi_idx, alternate_scheme ) + cfs_alt_scheme = not cfs_alt_scheme + + toroidal_alt_scheme = not toroidal_alt_scheme + def export_mesh(self, filename="source_mesh", export_dir=""): """Use PyMOAB interface to write source mesh with source strengths tagged. From 4d339b3e0143c1dc88bf75a0b7cf07ebe1c4c730 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Sat, 2 Nov 2024 13:06:13 -0500 Subject: [PATCH 04/11] Fix which hexahedron division scheme is considered the alternate --- parastell/source_mesh.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index c30d784..8490004 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -387,14 +387,6 @@ def _create_tets_from_hex( # Conditionally use alternate scheme defining hexahedron splitting to # avoid gaps and overlaps between non-planar hexahedron faces if alternate_scheme: - hex_canon_ids = [ - [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], - [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], - [idx_list[0], idx_list[7], idx_list[5], idx_list[4]], - [idx_list[7], idx_list[2], idx_list[5], idx_list[6]], - [idx_list[0], idx_list[2], idx_list[5], idx_list[7]], - ] - else: hex_canon_ids = [ [idx_list[0], idx_list[3], idx_list[1], idx_list[4]], [idx_list[1], idx_list[3], idx_list[2], idx_list[6]], @@ -402,6 +394,14 @@ def _create_tets_from_hex( [idx_list[3], idx_list[6], idx_list[4], idx_list[7]], [idx_list[1], idx_list[3], idx_list[6], idx_list[4]], ] + else: + hex_canon_ids = [ + [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], + [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], + [idx_list[0], idx_list[7], idx_list[5], idx_list[4]], + [idx_list[7], idx_list[2], idx_list[5], idx_list[6]], + [idx_list[0], idx_list[2], idx_list[5], idx_list[7]], + ] for vertex_ids in hex_canon_ids: self._create_tet(vertex_ids) From da08ac4a3f388200e7a71dc143bcbb87f57632de Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 7 Nov 2024 15:02:14 -0600 Subject: [PATCH 05/11] Incorporate review suggestions --- parastell/source_mesh.py | 130 ++++++++++++++++++++------------------- 1 file changed, 67 insertions(+), 63 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index 8490004..9caefe2 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -133,18 +133,54 @@ def __init__( self._create_mbc() + @property + def num_theta(self): + return self._num_theta + + @num_theta.setter + def num_theta(self, value): + if value % 2 != 1.0: + e = AttributeError( + "To ensure that tetrahedral faces are coincident at the end of " + "the closed poloidal loop, the number of poloidal intervals " + "must be even. To ensure this, the number of poloidal grid " + "points must be odd." + ) + self._logger.error(e.args[0]) + raise e + self._num_theta = value + + @property + def num_phi(self): + return self._num_phi + + @num_phi.setter + def num_phi(self, value): + self._num_phi = value + @property def toroidal_extent(self): return self._toroidal_extent @toroidal_extent.setter def toroidal_extent(self, angle): - self._toroidal_extent = np.deg2rad(angle) - if self._toroidal_extent > 360.0: - e = ValueError("Toroidal extent cannot exceed 360.0 degrees.") + if angle > 360.0: + e = AttributeError("Toroidal extent cannot exceed 360.0 degrees.") + self._logger.error(e.args[0]) + raise e + + if angle == 360.0 and self._num_phi % 2 != 1.0: + e = AttributeError( + "To ensure that tetrahedral faces are coincident at the end of " + "the closed toroidal loop, the number of toroidal intervals " + "must be even. To ensure this, the number of toroidal grid " + "points must be odd." + ) self._logger.error(e.args[0]) raise e + self._toroidal_extent = np.deg2rad(angle) + @property def logger(self): return self._logger @@ -337,9 +373,7 @@ def _get_vertex_id(self, vertex_idx): return id - def _create_tets_from_hex( - self, s_idx, theta_idx, phi_idx, alternate_scheme - ): + def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): """Creates five tetrahedra from defined hexahedron. (Internal function not intended to be called externally) @@ -347,11 +381,6 @@ def _create_tets_from_hex( s_idx (int): index defining location along CFS axis. theta_idx (int): index defining location along poloidal angle axis. phi_idx (int): index defining location along toroidal angle axis. - alternate_scheme (bool): flag indicating whether alternate - hexahedron division scheme should be used. - - Returns: - alternate_scheme (bool): switched flag. """ # relative offsets of vertices in a 3-D index space @@ -384,42 +413,38 @@ def _create_tets_from_hex( # first, ordered clockwise relative to the thumb, followed by the # remaining vertex at the end of the thumb. # See Moreno, Bader, Wilson 2024 for hexahedron splitting - # Conditionally use alternate scheme defining hexahedron splitting to - # avoid gaps and overlaps between non-planar hexahedron faces - if alternate_scheme: - hex_canon_ids = [ + canonical_ordering_schemes = [ + [ [idx_list[0], idx_list[3], idx_list[1], idx_list[4]], [idx_list[1], idx_list[3], idx_list[2], idx_list[6]], [idx_list[1], idx_list[4], idx_list[6], idx_list[5]], [idx_list[3], idx_list[6], idx_list[4], idx_list[7]], [idx_list[1], idx_list[3], idx_list[6], idx_list[4]], - ] - else: - hex_canon_ids = [ + ], + [ [idx_list[0], idx_list[2], idx_list[1], idx_list[5]], [idx_list[0], idx_list[3], idx_list[2], idx_list[7]], [idx_list[0], idx_list[7], idx_list[5], idx_list[4]], [idx_list[7], idx_list[2], idx_list[5], idx_list[6]], [idx_list[0], idx_list[2], idx_list[5], idx_list[7]], - ] + ], + ] - for vertex_ids in hex_canon_ids: - self._create_tet(vertex_ids) + # Alternate canonical ordering schemes defining hexahedron splitting to + # avoid gaps and overlaps between non-planar hexahedron faces + # s_idx begins at 0 and theta_idx begins at 1 + scheme_idx = ((s_idx + 1) + (theta_idx - 1) + phi_idx) % 2 - return not alternate_scheme + for vertex_ids in canonical_ordering_schemes[scheme_idx]: + self._create_tet(vertex_ids) - def _create_tets_from_wedge(self, theta_idx, phi_idx, alternate_scheme): + def _create_tets_from_wedge(self, theta_idx, phi_idx): """Creates three tetrahedra from defined wedge. (Internal function not intended to be called externally) Arguments: theta_idx (int): index defining location along poloidal angle axis. phi_idx (int): index defining location along toroidal angle axis. - alternate_scheme (bool): flag indicating whether alternate - hexahedron division scheme should be used. - - Returns: - alternate_scheme (bool): switched flag. """ # relative offsets of wedge vertices in a 3-D index space @@ -448,25 +473,26 @@ def _create_tets_from_wedge(self, theta_idx, phi_idx, alternate_scheme): # first, ordered clockwise relative to the thumb, followed by the # remaining vertex at the end of the thumb. # See Moreno, Bader, Wilson 2024 for wedge splitting - # Conditionally alternate ordering of vertices defining wedge splitting - # to avoid gaps and overlaps between non-planar wedge faces - if alternate_scheme: - wedge_canon_ids = [ + canonical_ordering_schemes = [ + [ [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], [idx_list[1], idx_list[3], idx_list[5], idx_list[4]], [idx_list[1], idx_list[3], idx_list[2], idx_list[5]], - ] - else: - wedge_canon_ids = [ + ], + [ [idx_list[0], idx_list[2], idx_list[1], idx_list[3]], [idx_list[3], idx_list[2], idx_list[4], idx_list[5]], [idx_list[3], idx_list[2], idx_list[1], idx_list[4]], - ] + ], + ] - for vertex_ids in wedge_canon_ids: - self._create_tet(vertex_ids) + # Alternate canonical ordering schemes defining wedge splitting to + # avoid gaps and overlaps between non-planar wedge faces + # s_idx begins at 0 and theta_idx begins at 1 + scheme_idx = ((theta_idx - 1) + phi_idx) % 2 - return not alternate_scheme + for vertex_ids in canonical_ordering_schemes[scheme_idx]: + self._create_tet(vertex_ids) def create_mesh(self): """Creates volumetric source mesh in real space.""" @@ -475,37 +501,15 @@ def create_mesh(self): self.mesh_set = self.mbc.create_meshset() self.mbc.add_entity(self.mesh_set, self.verts) - # Initialize alternate scheme flag for adjacent toroidal blocks - toroidal_alt_scheme = False - for phi_idx in range(self.num_phi - 1): - # Initialize alternate scheme flag at beginning of each toroidal - # block - alternate_scheme = toroidal_alt_scheme - # Create tetrahedra for wedges at center of plasma for theta_idx in range(1, self.num_theta): - alternate_scheme = self._create_tets_from_wedge( - theta_idx, phi_idx, alternate_scheme - ) - - # Initialize alternate scheme flag for adjacent CFS blocks - cfs_alt_scheme = not toroidal_alt_scheme + self._create_tets_from_wedge(theta_idx, phi_idx) # Create tetrahedra for hexahedra beyond center of plasma for s_idx in range(self.num_s - 2): - # Initialize alternate scheme flag at beginning of each - # CFS block - alternate_scheme = cfs_alt_scheme - for theta_idx in range(1, self.num_theta): - alternate_scheme = self._create_tets_from_hex( - s_idx, theta_idx, phi_idx, alternate_scheme - ) - - cfs_alt_scheme = not cfs_alt_scheme - - toroidal_alt_scheme = not toroidal_alt_scheme + self._create_tets_from_hex(s_idx, theta_idx, phi_idx) def export_mesh(self, filename="source_mesh", export_dir=""): """Use PyMOAB interface to write source mesh with source strengths From 18572936fa8efae48216e7b523542079e5afceb3 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 7 Nov 2024 15:16:29 -0600 Subject: [PATCH 06/11] Fix num_theta and num_phi checks and update tests --- parastell/source_mesh.py | 4 ++-- tests/test_parastell.py | 2 +- tests/test_source_mesh.py | 12 ++++++------ 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index 9caefe2..a2484bf 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -139,7 +139,7 @@ def num_theta(self): @num_theta.setter def num_theta(self, value): - if value % 2 != 1.0: + if value % 2 != 1: e = AttributeError( "To ensure that tetrahedral faces are coincident at the end of " "the closed poloidal loop, the number of poloidal intervals " @@ -169,7 +169,7 @@ def toroidal_extent(self, angle): self._logger.error(e.args[0]) raise e - if angle == 360.0 and self._num_phi % 2 != 1.0: + if angle == 360.0 and self._num_phi % 2 != 1: e = AttributeError( "To ensure that tetrahedral faces are coincident at the end of " "the closed toroidal loop, the number of toroidal intervals " diff --git a/tests/test_parastell.py b/tests/test_parastell.py index 70a6f40..b05095d 100644 --- a/tests/test_parastell.py +++ b/tests/test_parastell.py @@ -103,7 +103,7 @@ def test_parastell(stellarator): assert Path(step_filename_exp).with_suffix(".step").exists() assert Path(mesh_filename_exp).with_suffix(".h5m").exists() - mesh_size = (4, 8, 4) + mesh_size = (6, 41, 9) toroidal_extent = 15.0 stellarator.construct_source_mesh(mesh_size, toroidal_extent) diff --git a/tests/test_source_mesh.py b/tests/test_source_mesh.py index f0a2fbb..dab9372 100644 --- a/tests/test_source_mesh.py +++ b/tests/test_source_mesh.py @@ -24,8 +24,8 @@ def source_mesh(): # Set mesh size to minimum that maintains element aspect ratios that do not # result in negative volumes - mesh_size = (6, 41, 51) - toroidal_extent = 90.0 + mesh_size = (6, 41, 9) + toroidal_extent = 15.0 source_mesh_obj = sm.SourceMesh(vmec_obj, mesh_size, toroidal_extent) @@ -36,8 +36,8 @@ def test_mesh_basics(source_mesh): num_s_exp = 6 num_theta_exp = 41 - num_phi_exp = 51 - tor_ext_exp = 90.0 + num_phi_exp = 9 + tor_ext_exp = 15.0 scale_exp = 100 remove_files() @@ -55,7 +55,7 @@ def test_vertices(source_mesh): num_s = 6 num_theta = 41 - num_phi = 51 + num_phi = 9 num_verts_exp = num_phi * ((num_s - 1) * (num_theta - 1) + 1) @@ -74,7 +74,7 @@ def test_mesh_generation(source_mesh): num_s = 6 num_theta = 41 - num_phi = 51 + num_phi = 9 tets_per_wedge = 3 tets_per_hex = 5 From ed787fa17a8f3b72e4fee02260340c43aee70bf2 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 8 Nov 2024 10:09:20 -0600 Subject: [PATCH 07/11] Change naming convention for references to flux-coordinate axes --- parastell/source_mesh.py | 133 +++++++++++++++++++++------------------ 1 file changed, 73 insertions(+), 60 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index a2484bf..6b1fbcd 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -85,10 +85,10 @@ class SourceMesh(object): any closed flux surface label, s, poloidal angle, theta, and toroidal angle, phi. mesh_size (tuple of int): number of grid points along each axis of - flux-coordinate space, in the order (num_s, num_theta, num_phi). - 'num_s' is the number of closed flux surfaces for vertex locations - in each toroidal plane. 'num_theta' is the number of poloidal - angles for vertex locations in each toroidal plane. 'num_phi' is + flux-coordinate space, in the order (num_cfs_pts, num_poloidal_pts, num_toroidal_pts). + 'num_cfs_pts' is the number of closed flux surfaces for vertex locations + in each toroidal plane. 'num_poloidal_pts' is the number of poloidal + angles for vertex locations in each toroidal plane. 'num_toroidal_pts' is the number of toroidal angles for planes of vertices. toroidal_extent (float): extent of source mesh in toroidal direction [deg]. @@ -112,9 +112,9 @@ def __init__( self.logger = logger self.vmec_obj = vmec_obj - self.num_s = mesh_size[0] - self.num_theta = mesh_size[1] - self.num_phi = mesh_size[2] + self.num_cfs_pts = mesh_size[0] + self.num_poloidal_pts = mesh_size[1] + self.num_toroidal_pts = mesh_size[2] self.toroidal_extent = toroidal_extent self.scale = m2cm @@ -134,11 +134,11 @@ def __init__( self._create_mbc() @property - def num_theta(self): - return self._num_theta + def num_poloidal_pts(self): + return self._num_poloidal_pts - @num_theta.setter - def num_theta(self, value): + @num_poloidal_pts.setter + def num_poloidal_pts(self, value): if value % 2 != 1: e = AttributeError( "To ensure that tetrahedral faces are coincident at the end of " @@ -148,15 +148,15 @@ def num_theta(self, value): ) self._logger.error(e.args[0]) raise e - self._num_theta = value + self._num_poloidal_pts = value @property - def num_phi(self): - return self._num_phi + def num_toroidal_pts(self): + return self._num_toroidal_pts - @num_phi.setter - def num_phi(self, value): - self._num_phi = value + @num_toroidal_pts.setter + def num_toroidal_pts(self, value): + self._num_toroidal_pts = value @property def toroidal_extent(self): @@ -169,7 +169,7 @@ def toroidal_extent(self, angle): self._logger.error(e.args[0]) raise e - if angle == 360.0 and self._num_phi % 2 != 1: + if angle == 360.0 and self._num_toroidal_pts % 2 != 1: e = AttributeError( "To ensure that tetrahedral faces are coincident at the end of " "the closed toroidal loop, the number of toroidal intervals " @@ -228,44 +228,51 @@ def create_vertices(self): """ self._logger.info("Computing source mesh point cloud...") - phi_list = np.linspace(0, self._toroidal_extent, num=self.num_phi) + toroidal_grid_pts = np.linspace( + 0, self._toroidal_extent, num=self._num_toroidal_pts + ) # don't include magnetic axis in list of s values - s_list = np.linspace(0.0, 1.0, num=self.num_s)[1:] + cfs_grid_pts = np.linspace(0.0, 1.0, num=self.num_cfs_pts)[1:] # don't include repeated entry at 0 == 2*pi - theta_list = np.linspace(0, 2 * np.pi, num=self.num_theta)[:-1] + poloidal_grid_pts = np.linspace( + 0, 2 * np.pi, num=self._num_poloidal_pts + )[:-1] # don't include repeated entry at 0 == 2*pi if self._toroidal_extent == 2 * np.pi: - phi_list = phi_list[:-1] + toroidal_grid_pts = toroidal_grid_pts[:-1] - self.verts_per_ring = theta_list.shape[0] + self.verts_per_ring = poloidal_grid_pts.shape[0] # add one vertex per plane for magenetic axis - self.verts_per_plane = s_list.shape[0] * self.verts_per_ring + 1 + self.verts_per_plane = cfs_grid_pts.shape[0] * self.verts_per_ring + 1 - num_verts = phi_list.shape[0] * self.verts_per_plane + num_verts = toroidal_grid_pts.shape[0] * self.verts_per_plane self.coords = np.zeros((num_verts, 3)) self.coords_s = np.zeros(num_verts) # Initialize vertex index vert_idx = 0 - for phi in phi_list: + for toroidal_ang in toroidal_grid_pts: # vertex coordinates on magnetic axis self.coords[vert_idx, :] = ( - np.array(self.vmec_obj.vmec2xyz(0, 0, phi)) * self.scale + np.array(self.vmec_obj.vmec2xyz(0, 0, toroidal_ang)) + * self.scale ) self.coords_s[vert_idx] = 0 vert_idx += 1 # vertex coordinate away from magnetic axis - for s in s_list: - for theta in theta_list: + for cfs in cfs_grid_pts: + for poloidal_ang in poloidal_grid_pts: self.coords[vert_idx, :] = ( - np.array(self.vmec_obj.vmec2xyz(s, theta, phi)) + np.array( + self.vmec_obj.vmec2xyz(cfs, poloidal_ang, phi) + ) * self.scale ) - self.coords_s[vert_idx] = s + self.coords_s[vert_idx] = cfs vert_idx += 1 @@ -352,35 +359,38 @@ def _get_vertex_id(self, vertex_idx): id (int): vertex index in row-major order as stored by MOAB """ - s_idx, theta_idx, phi_idx = vertex_idx + cfs_idx, poloidal_idx, toroidal_idx = vertex_idx - ma_offset = phi_idx * self.verts_per_plane + ma_offset = toroidal_idx * self.verts_per_plane # Wrap around if final plane and it is 2*pi - if self._toroidal_extent == 2 * np.pi and phi_idx == self.num_phi - 1: + if ( + self._toroidal_extent == 2 * np.pi + and toroidal_idx == self._num_toroidal_pts - 1 + ): ma_offset = 0 # Compute index offset from closed flux surface - s_offset = s_idx * self.verts_per_ring + s_offset = cfs_idx * self.verts_per_ring - theta_offset = theta_idx + theta_offset = poloidal_idx # Wrap around if theta is 2*pi - if theta_idx == self.num_theta: + if poloidal_idx == self._num_poloidal_pts: theta_offset = 1 id = ma_offset + s_offset + theta_offset return id - def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): + def _create_tets_from_hex(self, cfs_idx, poloidal_idx, toroidal_idx): """Creates five tetrahedra from defined hexahedron. (Internal function not intended to be called externally) Arguments: - s_idx (int): index defining location along CFS axis. - theta_idx (int): index defining location along poloidal angle axis. - phi_idx (int): index defining location along toroidal angle axis. + cfs_idx (int): index defining location along CFS axis. + poloidal_idx (int): index defining location along poloidal angle axis. + toroidal_idx (int): index defining location along toroidal angle axis. """ # relative offsets of vertices in a 3-D index space @@ -399,7 +409,8 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): # Ids of hex vertices applying offset stencil to current point hex_idx_data = ( - np.array([s_idx, theta_idx, phi_idx]) + hex_vertex_stencil + np.array([cfs_idx, poloidal_idx, toroidal_idx]) + + hex_vertex_stencil ) idx_list = [ @@ -432,35 +443,35 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx): # Alternate canonical ordering schemes defining hexahedron splitting to # avoid gaps and overlaps between non-planar hexahedron faces - # s_idx begins at 0 and theta_idx begins at 1 - scheme_idx = ((s_idx + 1) + (theta_idx - 1) + phi_idx) % 2 + # cfs_idx begins at 0 and poloidal_idx begins at 1 + scheme_idx = ((cfs_idx + 1) + (poloidal_idx - 1) + toroidal_idx) % 2 for vertex_ids in canonical_ordering_schemes[scheme_idx]: self._create_tet(vertex_ids) - def _create_tets_from_wedge(self, theta_idx, phi_idx): + def _create_tets_from_wedge(self, poloidal_idx, toroidal_idx): """Creates three tetrahedra from defined wedge. (Internal function not intended to be called externally) Arguments: - theta_idx (int): index defining location along poloidal angle axis. - phi_idx (int): index defining location along toroidal angle axis. + poloidal_idx (int): index defining location along poloidal angle axis. + toroidal_idx (int): index defining location along toroidal angle axis. """ # relative offsets of wedge vertices in a 3-D index space wedge_vertex_stencil = np.array( [ [0, 0, 0], - [0, theta_idx, 0], - [0, theta_idx + 1, 0], + [0, poloidal_idx, 0], + [0, poloidal_idx + 1, 0], [0, 0, 1], - [0, theta_idx, 1], - [0, theta_idx + 1, 1], + [0, poloidal_idx, 1], + [0, poloidal_idx + 1, 1], ] ) # Ids of wedge vertices applying offset stencil to current point - wedge_idx_data = np.array([0, 0, phi_idx]) + wedge_vertex_stencil + wedge_idx_data = np.array([0, 0, toroidal_idx]) + wedge_vertex_stencil idx_list = [ self._get_vertex_id(vertex_idx) for vertex_idx in wedge_idx_data @@ -488,8 +499,8 @@ def _create_tets_from_wedge(self, theta_idx, phi_idx): # Alternate canonical ordering schemes defining wedge splitting to # avoid gaps and overlaps between non-planar wedge faces - # s_idx begins at 0 and theta_idx begins at 1 - scheme_idx = ((theta_idx - 1) + phi_idx) % 2 + # cfs_idx begins at 0 and poloidal_idx begins at 1 + scheme_idx = ((poloidal_idx - 1) + toroidal_idx) % 2 for vertex_ids in canonical_ordering_schemes[scheme_idx]: self._create_tet(vertex_ids) @@ -501,15 +512,17 @@ def create_mesh(self): self.mesh_set = self.mbc.create_meshset() self.mbc.add_entity(self.mesh_set, self.verts) - for phi_idx in range(self.num_phi - 1): + for toroidal_idx in range(self._num_toroidal_pts - 1): # Create tetrahedra for wedges at center of plasma - for theta_idx in range(1, self.num_theta): - self._create_tets_from_wedge(theta_idx, phi_idx) + for poloidal_idx in range(1, self._num_poloidal_pts): + self._create_tets_from_wedge(poloidal_idx, toroidal_idx) # Create tetrahedra for hexahedra beyond center of plasma - for s_idx in range(self.num_s - 2): - for theta_idx in range(1, self.num_theta): - self._create_tets_from_hex(s_idx, theta_idx, phi_idx) + for cfs_idx in range(self.num_cfs_pts - 2): + for poloidal_idx in range(1, self._num_poloidal_pts): + self._create_tets_from_hex( + cfs_idx, poloidal_idx, toroidal_idx + ) def export_mesh(self, filename="source_mesh", export_dir=""): """Use PyMOAB interface to write source mesh with source strengths From 6b26f2ccc147b3af715433d4410a2c5cc932cec3 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 8 Nov 2024 10:17:28 -0600 Subject: [PATCH 08/11] Finalize changes to naming convention --- parastell/source_mesh.py | 61 +++++++++++++++++++++------------------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index 6b1fbcd..f618597 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -66,30 +66,31 @@ class SourceMesh(object): neutrons in a magnetically confined plasma described by a VMEC plasma equilibrium. - The mesh will be defined on a regular grid in the plasma coordinates of s, - theta, phi. Mesh vertices will be defined on circular grid at each toroidal - plane, and connected between toroidal planes. This results in wedge elements - along the magnetic axis and hexagonal elements throughout the remainder of - the mesh. Each of these elements will be subdivided into tetrahedra (4 for - the wedges and 5 for the hexahedra) to result in a mesh that is simpler to - use. - - Each tetrahedron will be tagged with the volumetric neutron source intensity - in n/cm3/s, using on a finite-element based quadrature of the source - intensity evaluated at each vertex. + The mesh will be defined on a regular grid in the flux coordinates of + (CFS value, poloidal angle, toroidal angle). Mesh vertices will be defined + on circular grid at each toroidal plane, and connected between toroidal + planes. This results in wedge elements along the magnetic axis and + hexagonal elements throughout the remainder of the mesh. Each of these + elements will be subdivided into tetrahedra (3 for the wedges and 5 for the + hexahedra) to result in a mesh that is simpler to use. + + Each tetrahedron will be tagged with the volumetric neutron source + intensity in n/cm3/s, using on a finite-element based quadrature of the + source intensity evaluated at each vertex. Arguments: vmec_obj (object): plasma equilibrium VMEC object as defined by the PyStell-UW VMEC reader. Must have a method - 'vmec2xyz(s, theta, phi)' that returns an (x,y,z) coordinate for - any closed flux surface label, s, poloidal angle, theta, and - toroidal angle, phi. + 'vmec2xyz(cfs, poloidal_ang, toroidal_ang)' that returns + a (x,y,z) coordinate for any closed flux surface value, cfs, + poloidal angle, poloidal_ang, and toroidal angle, toroidal_ang. mesh_size (tuple of int): number of grid points along each axis of - flux-coordinate space, in the order (num_cfs_pts, num_poloidal_pts, num_toroidal_pts). - 'num_cfs_pts' is the number of closed flux surfaces for vertex locations - in each toroidal plane. 'num_poloidal_pts' is the number of poloidal - angles for vertex locations in each toroidal plane. 'num_toroidal_pts' is - the number of toroidal angles for planes of vertices. + flux-coordinate space, in the order (num_cfs_pts, num_poloidal_pts, + num_toroidal_pts). 'num_cfs_pts' is the number of closed flux + surfaces for vertex locations in each toroidal plane. + 'num_poloidal_pts' is the number of poloidal angles for vertex + locations in each toroidal plane. 'num_toroidal_pts' is the number + of toroidal angles for planes of vertices. toroidal_extent (float): extent of source mesh in toroidal direction [deg]. logger (object): logger object (optional, defaults to None). If no @@ -248,7 +249,7 @@ def create_vertices(self): num_verts = toroidal_grid_pts.shape[0] * self.verts_per_plane self.coords = np.zeros((num_verts, 3)) - self.coords_s = np.zeros(num_verts) + self.coords_cfs = np.zeros(num_verts) # Initialize vertex index vert_idx = 0 @@ -259,7 +260,7 @@ def create_vertices(self): np.array(self.vmec_obj.vmec2xyz(0, 0, toroidal_ang)) * self.scale ) - self.coords_s[vert_idx] = 0 + self.coords_cfs[vert_idx] = 0 vert_idx += 1 @@ -268,11 +269,13 @@ def create_vertices(self): for poloidal_ang in poloidal_grid_pts: self.coords[vert_idx, :] = ( np.array( - self.vmec_obj.vmec2xyz(cfs, poloidal_ang, phi) + self.vmec_obj.vmec2xyz( + cfs, poloidal_ang, toroidal_ang + ) ) * self.scale ) - self.coords_s[vert_idx] = cfs + self.coords_cfs[vert_idx] = cfs vert_idx += 1 @@ -295,7 +298,7 @@ def _source_strength(self, tet_ids): # Initialize list of source strengths for each tetrahedron vertex vertex_strengths = [ - self.reaction_rate(*self.plasma_conditions(self.coords_s[id])) + self.reaction_rate(*self.plasma_conditions(self.coords_cfs[id])) for id in tet_ids ] @@ -371,15 +374,15 @@ def _get_vertex_id(self, vertex_idx): ma_offset = 0 # Compute index offset from closed flux surface - s_offset = cfs_idx * self.verts_per_ring + cfs_offset = cfs_idx * self.verts_per_ring - theta_offset = poloidal_idx + poloidal_offset = poloidal_idx - # Wrap around if theta is 2*pi + # Wrap around if poloidal angle is 2*pi if poloidal_idx == self._num_poloidal_pts: - theta_offset = 1 + poloidal_offset = 1 - id = ma_offset + s_offset + theta_offset + id = ma_offset + cfs_offset + poloidal_offset return id From 939e9468b319ad515d050d9ee41bae24d6f19ca3 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 8 Nov 2024 11:54:28 -0600 Subject: [PATCH 09/11] Modify vertex indexing scheme --- parastell/source_mesh.py | 51 +++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 27 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index f618597..b9ff8d9 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -229,17 +229,18 @@ def create_vertices(self): """ self._logger.info("Computing source mesh point cloud...") - toroidal_grid_pts = np.linspace( - 0, self._toroidal_extent, num=self._num_toroidal_pts - ) - # don't include magnetic axis in list of s values + # Exclude entry at magnetic axis cfs_grid_pts = np.linspace(0.0, 1.0, num=self.num_cfs_pts)[1:] - # don't include repeated entry at 0 == 2*pi + + # Exclude repeated entry at 0 == 2*pi poloidal_grid_pts = np.linspace( 0, 2 * np.pi, num=self._num_poloidal_pts )[:-1] - # don't include repeated entry at 0 == 2*pi + toroidal_grid_pts = np.linspace( + 0, self._toroidal_extent, num=self._num_toroidal_pts + ) + # Conditionally exclude repeated entry at 0 == 2*pi if self._toroidal_extent == 2 * np.pi: toroidal_grid_pts = toroidal_grid_pts[:-1] @@ -292,7 +293,6 @@ def _source_strength(self, tet_ids): Returns: ss (float): integrated source strength for tetrahedron. """ - # Initialize list of vertex coordinates for each tetrahedron vertex tet_coords = [self.coords[id] for id in tet_ids] @@ -335,7 +335,6 @@ def _create_tet(self, tet_ids): Arguments: tet_ids (list of int): tetrahedron vertex indices. """ - tet_verts = [self.verts[int(id)] for id in tet_ids] tet = self.mbc.create_element(types.MBTET, tet_verts) self.mbc.add_entity(self.mesh_set, tet) @@ -361,7 +360,6 @@ def _get_vertex_id(self, vertex_idx): Returns: id (int): vertex index in row-major order as stored by MOAB """ - cfs_idx, poloidal_idx, toroidal_idx = vertex_idx ma_offset = toroidal_idx * self.verts_per_plane @@ -373,14 +371,17 @@ def _get_vertex_id(self, vertex_idx): ): ma_offset = 0 - # Compute index offset from closed flux surface - cfs_offset = cfs_idx * self.verts_per_ring + # Compute index offset from closed flux surface, taking single vertex + # at magnetic axis into account + if cfs_idx == 0: + cfs_offset = cfs_idx + else: + cfs_offset = (cfs_idx - 1) * self.verts_per_ring + 1 poloidal_offset = poloidal_idx - # Wrap around if poloidal angle is 2*pi - if poloidal_idx == self._num_poloidal_pts: - poloidal_offset = 1 + if poloidal_idx == self._num_poloidal_pts - 1: + poloidal_offset = 0 id = ma_offset + cfs_offset + poloidal_offset @@ -395,7 +396,6 @@ def _create_tets_from_hex(self, cfs_idx, poloidal_idx, toroidal_idx): poloidal_idx (int): index defining location along poloidal angle axis. toroidal_idx (int): index defining location along toroidal angle axis. """ - # relative offsets of vertices in a 3-D index space hex_vertex_stencil = np.array( [ @@ -446,8 +446,7 @@ def _create_tets_from_hex(self, cfs_idx, poloidal_idx, toroidal_idx): # Alternate canonical ordering schemes defining hexahedron splitting to # avoid gaps and overlaps between non-planar hexahedron faces - # cfs_idx begins at 0 and poloidal_idx begins at 1 - scheme_idx = ((cfs_idx + 1) + (poloidal_idx - 1) + toroidal_idx) % 2 + scheme_idx = (cfs_idx + poloidal_idx + toroidal_idx) % 2 for vertex_ids in canonical_ordering_schemes[scheme_idx]: self._create_tet(vertex_ids) @@ -460,16 +459,15 @@ def _create_tets_from_wedge(self, poloidal_idx, toroidal_idx): poloidal_idx (int): index defining location along poloidal angle axis. toroidal_idx (int): index defining location along toroidal angle axis. """ - # relative offsets of wedge vertices in a 3-D index space wedge_vertex_stencil = np.array( [ [0, 0, 0], - [0, poloidal_idx, 0], - [0, poloidal_idx + 1, 0], + [1, poloidal_idx, 0], + [1, poloidal_idx + 1, 0], [0, 0, 1], - [0, poloidal_idx, 1], - [0, poloidal_idx + 1, 1], + [1, poloidal_idx, 1], + [1, poloidal_idx + 1, 1], ] ) @@ -502,8 +500,7 @@ def _create_tets_from_wedge(self, poloidal_idx, toroidal_idx): # Alternate canonical ordering schemes defining wedge splitting to # avoid gaps and overlaps between non-planar wedge faces - # cfs_idx begins at 0 and poloidal_idx begins at 1 - scheme_idx = ((poloidal_idx - 1) + toroidal_idx) % 2 + scheme_idx = (poloidal_idx + toroidal_idx) % 2 for vertex_ids in canonical_ordering_schemes[scheme_idx]: self._create_tet(vertex_ids) @@ -517,12 +514,12 @@ def create_mesh(self): for toroidal_idx in range(self._num_toroidal_pts - 1): # Create tetrahedra for wedges at center of plasma - for poloidal_idx in range(1, self._num_poloidal_pts): + for poloidal_idx in range(self._num_poloidal_pts - 1): self._create_tets_from_wedge(poloidal_idx, toroidal_idx) # Create tetrahedra for hexahedra beyond center of plasma - for cfs_idx in range(self.num_cfs_pts - 2): - for poloidal_idx in range(1, self._num_poloidal_pts): + for cfs_idx in range(1, self.num_cfs_pts - 1): + for poloidal_idx in range(self._num_poloidal_pts - 1): self._create_tets_from_hex( cfs_idx, poloidal_idx, toroidal_idx ) From 427097463cfc69eb93406c1ff736868126da283b Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 8 Nov 2024 12:05:31 -0600 Subject: [PATCH 10/11] Update parameter references in source mesh test --- tests/test_source_mesh.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/tests/test_source_mesh.py b/tests/test_source_mesh.py index dab9372..2a64b7f 100644 --- a/tests/test_source_mesh.py +++ b/tests/test_source_mesh.py @@ -34,17 +34,17 @@ def source_mesh(): def test_mesh_basics(source_mesh): - num_s_exp = 6 - num_theta_exp = 41 - num_phi_exp = 9 + num_cfs_exp = 6 + num_poloidal_pts_exp = 41 + num_toroidal_pts_exp = 9 tor_ext_exp = 15.0 scale_exp = 100 remove_files() - assert source_mesh.num_s == num_s_exp - assert source_mesh.num_theta == num_theta_exp - assert source_mesh.num_phi == num_phi_exp + assert source_mesh.num_cfs_pts == num_cfs_exp + assert source_mesh.num_poloidal_pts == num_poloidal_pts_exp + assert source_mesh.num_toroidal_pts == num_toroidal_pts_exp assert source_mesh.toroidal_extent == np.deg2rad(tor_ext_exp) assert source_mesh.scale == scale_exp @@ -53,18 +53,20 @@ def test_mesh_basics(source_mesh): def test_vertices(source_mesh): - num_s = 6 - num_theta = 41 - num_phi = 9 + num_cfs_exp = 6 + num_poloidal_pts_exp = 41 + num_toroidal_pts_exp = 9 - num_verts_exp = num_phi * ((num_s - 1) * (num_theta - 1) + 1) + num_verts_exp = num_toroidal_pts_exp * ( + (num_cfs_exp - 1) * (num_poloidal_pts_exp - 1) + 1 + ) remove_files() source_mesh.create_vertices() assert source_mesh.coords.shape == (num_verts_exp, 3) - assert source_mesh.coords_s.shape == (num_verts_exp,) + assert source_mesh.coords_cfs.shape == (num_verts_exp,) assert len(source_mesh.verts) == num_verts_exp remove_files() From 12f55ff488f9717786e83aefe661727a58bb6455 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Mon, 11 Nov 2024 12:46:38 -0600 Subject: [PATCH 11/11] Incorporate review suggestions --- parastell/source_mesh.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/parastell/source_mesh.py b/parastell/source_mesh.py index b9ff8d9..0d7d932 100644 --- a/parastell/source_mesh.py +++ b/parastell/source_mesh.py @@ -165,12 +165,14 @@ def toroidal_extent(self): @toroidal_extent.setter def toroidal_extent(self, angle): - if angle > 360.0: + angle = np.deg2rad(angle) + + if angle > 2 * np.pi: e = AttributeError("Toroidal extent cannot exceed 360.0 degrees.") self._logger.error(e.args[0]) raise e - if angle == 360.0 and self._num_toroidal_pts % 2 != 1: + if angle == 2 * np.pi and self._num_toroidal_pts % 2 != 1: e = AttributeError( "To ensure that tetrahedral faces are coincident at the end of " "the closed toroidal loop, the number of toroidal intervals " @@ -180,7 +182,7 @@ def toroidal_extent(self, angle): self._logger.error(e.args[0]) raise e - self._toroidal_extent = np.deg2rad(angle) + self._toroidal_extent = angle @property def logger(self):