Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement alternating division schemes between adjacent hexahedra and wedges #172

Merged
merged 11 commits into from
Nov 11, 2024
Merged

Implement alternating division schemes between adjacent hexahedra and wedges #172

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
fb4960d
Implement alternation flag
connoramoreno Nov 1, 2024
a2ae206
Change alternate scheme flag to local variable and change name
connoramoreno Nov 2, 2024
d99f5fe
Introduce additional alternate scheme flags to ensure face alignment …
connoramoreno Nov 2, 2024
4d339b3
Fix which hexahedron division scheme is considered the alternate
connoramoreno Nov 2, 2024
da08ac4
Incorporate review suggestions
connoramoreno Nov 7, 2024
1857293
Fix num_theta and num_phi checks and update tests
connoramoreno Nov 7, 2024
ed787fa
Change naming convention for references to flux-coordinate axes
connoramoreno Nov 8, 2024
6b26f2c
Finalize changes to naming convention
connoramoreno Nov 8, 2024
939e946
Modify vertex indexing scheme
connoramoreno Nov 8, 2024
4270974
Update parameter references in source mesh test
connoramoreno Nov 8, 2024
12f55ff
Incorporate review suggestions
connoramoreno Nov 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 82 additions & 17 deletions parastell/source_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
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:
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
Expand Down Expand Up @@ -342,7 +378,9 @@ def _create_tets_from_hex(self, s_idx, theta_idx, phi_idx):
(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.
"""

# relative offsets of vertices in a 3-D index space
Expand Down Expand Up @@ -375,23 +413,38 @@ 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]],
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]],
],
[
[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]],
],
]
gonuke marked this conversation as resolved.
Show resolved Hide resolved

for vertex_ids in hex_canon_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
connoramoreno marked this conversation as resolved.
Show resolved Hide resolved

for vertex_ids in canonical_ordering_schemes[scheme_idx]:
self._create_tet(vertex_ids)

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:
idx_list (list of int): list of wedge vertex indices.
theta_idx (int): index defining location along poloidal angle axis.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand the stencil below. Maybe it's been this way for a while, but I don't see any points in this stencil at a s_idx=1 which I expect.
I'm also not sure why this goes from [1,num_phi] istead of [0,num_phi-1]

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I agree that it's slightly confusing, and our indexing scheme could use some changes to make it more intuitive. The way we've defined things is that the vertex at the magnetic axis and those on the first closed flux surface all have a s_idx of 0. This is also why theta_idx begins at 1, to differentiate the first vertex on the first CFS from that on the magnetic axis.

It would make much more sense to have s_idx = 0 to always refer to the magnetic axis vertices, and s_idx = 1 refer to the first CFS. As such, we could then have theta_idx start from 0. I think I'll also change the naming convention from using to theta_idx and phi_idx to using poloidal_idx and toroidal_idx, respectively, to avoid confusion.

phi_idx (int): index defining location along toroidal angle axis.
"""

# relative offsets of wedge vertices in a 3-D index space
Expand Down Expand Up @@ -420,13 +473,25 @@ 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]],
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]],
],
[
[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:
# 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

for vertex_ids in canonical_ordering_schemes[scheme_idx]:
self._create_tet(vertex_ids)

def create_mesh(self):
Expand Down
2 changes: 1 addition & 1 deletion tests/test_parastell.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
12 changes: 6 additions & 6 deletions tests/test_source_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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()
Expand All @@ -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)

Expand All @@ -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
Expand Down