crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

fixed proper boundary_type setting?

Browse source
This commit is contained in:
Ethan Peterson 2022-10-20 12:46:29 -04:00
parent 299525984f
commit fdc3cfcfbb
1 changed files with 20 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

26
openmc/model/surface_composite.py
View file

@ -676,16 +676,18 @@ class Polygon(CompositeSurface):
surfnames = []
i = 0
for surfset in self._surfsets:
for surf, op in surfset:
setattr(self, f'surface_{i}', surf)
surfnames.append(f'surface_{i}')
i += 1
print(surfset)
for surf, op, on_boundary in surfset:
if on_boundary:
setattr(self, f'surface_{i}', surf)
surfnames.append(f'surface_{i}')
i += 1
self._surfnames = tuple(surfnames)
# Generate a list of regions whose union represents the polygon.
regions = []
for surfs_ops in self._surfsets:
regions.append([getattr(surf, op)() for surf, op in surfs_ops])
regions.append([getattr(surf, op)() for surf, op, _ in surfs_ops])
self._regions = [openmc.Intersection(regs) for regs in regions]
# Create the union of all the convex subsets
@ -717,6 +719,14 @@ class Polygon(CompositeSurface):
tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
return rotation.dot(tangents.T).T
@property
def _equations(self):
normals = self._normals
equations = np.empty((normals.shape[0], 3))
equations[:, 0:2] = normals
equations[:, 2] = -np.sum(normals*self.points, axis=-1)
return equations
@property
def regions(self):
return self._regions
@ -804,11 +814,15 @@ class Polygon(CompositeSurface):
"""
basis = self.basis
boundary_eqns = self._equations
# Collect surface/operator pairs such that the intersection of the
# regions defined by these pairs is the inside of the polygon.
surfs_ops = []
# hull facet equation: dx*x + dy*y + c = 0
for dx, dy, c in qhull.equations:
# check if this facet is on the boundary of the polygon
facet_eq = np.array([dx, dy, c])
on_boundary = any([np.allclose(facet_eq, eq) for eq in boundary_eqns])
# Check if the facet is horizontal
if isclose(dx, 0):
if basis in ('xz', 'yz', 'rz'):
@ -850,7 +864,7 @@ class Polygon(CompositeSurface):
else:
op = '__neg__'
surfs_ops.append((surf, op))
surfs_ops.append((surf, op, on_boundary))
return surfs_ops