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

Merge pull request #2362 from eepeterson/polygon_fix

Browse source 
try to constrain triangulation
This commit is contained in:
Paul Romano 2023-03-25 15:30:42 -05:00 • committed by GitHub
commit 2630dbf4bc
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 200 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

172
openmc/model/surface_composite.py
View file

@ -659,21 +659,10 @@ class Polygon(CompositeSurface):
def __init__(self, points, basis='rz'):
check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
self._basis = basis
points = np.asarray(points, dtype=float)
check_iterable_type('points', points, float, min_depth=2, max_depth=2)
check_length('points', points[0, :], 2, 2)
# If the last point is the same as the first, remove it and make sure
# there are still at least 3 points for a valid polygon.
if np.allclose(points[0, :], points[-1, :]):
points = points[:-1, :]
check_length('points', points, 3)
# Order the points counter-clockwise (necessary for offset method)
self._points = self._make_ccw(points)
# Create a triangulation of the points.
self._tri = Delaunay(self._points, qhull_options='QJ')
# Create a constrained triangulation of the validated points.
# The constrained triangulation is set to the _tri attribute
self._constrain_triangulation(self._validate_points(points))
# Decompose the polygon into groups of simplices forming convex subsets
# and get the sets of (surface, operator) pairs defining the polygon
@ -723,7 +712,7 @@ class Polygon(CompositeSurface):
@property
def points(self):
return self._points
return self._tri.points
@property
def basis(self):
@ -738,7 +727,7 @@ class Polygon(CompositeSurface):
# Get the unit vectors that point from one point in the polygon to the
# next given that they are ordered counterclockwise and that the final
# point is connected to the first point
tangents = np.diff(self._points, axis=0, append=[self._points[0, :]])
tangents = np.diff(self.points, axis=0, append=[self.points[0, :]])
tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
# Rotate the tangent vectors clockwise by 90 degrees, which for a
# counter-clockwise ordered polygon will produce the outward normal
@ -761,8 +750,9 @@ class Polygon(CompositeSurface):
def region(self):
return self._region
def _make_ccw(self, points):
"""Order a set of points counter-clockwise.
def _validate_points(self, points):
"""Ensure the closed path defined by points does not intersect and is
oriented counter-clockwise.
Parameters
----------
@ -773,15 +763,145 @@ class Polygon(CompositeSurface):
-------
ordered_points : the input points ordered counter-clockwise
"""
points = np.asarray(points, dtype=float)
check_iterable_type('points', points, float, min_depth=2, max_depth=2)
check_length('points', points[0, :], 2, 2)
# If the last point is the same as the first, remove it and make sure
# there are still at least 3 points for a valid polygon.
if np.allclose(points[0, :], points[-1, :]):
points = points[:-1, :]
check_length('points', points, 3)
if len(points) != len(np.unique(points, axis=0)):
raise ValueError('Duplicate points were detected in the Polygon input')
# Order the points counter-clockwise (necessary for offset method)
# Calculates twice the signed area of the polygon using the "Shoelace
# Formula" https://en.wikipedia.org/wiki/Shoelace_formula
# If signed area is positive the curve is oriented counter-clockwise.
# If the signed area is negative the curve is oriented clockwise.
xpts, ypts = points.T
if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) < 0:
points = points[::-1, :]
# If signed area is positive the curve is oriented counter-clockwise
if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) > 0:
return points
# Check if polygon is self-intersecting by comparing edges pairwise
n = len(points)
for i in range(n):
p0 = points[i, :]
p1 = points[(i + 1) % n, :]
for j in range(i + 1, n):
p2 = points[j, :]
p3 = points[(j + 1) % n, :]
# Compute orientation of p0 wrt p2->p3 line segment
cp0 = np.cross(p3-p0, p2-p0)
# Compute orientation of p1 wrt p2->p3 line segment
cp1 = np.cross(p3-p1, p2-p1)
# Compute orientation of p2 wrt p0->p1 line segment
cp2 = np.cross(p1-p2, p0-p2)
# Compute orientation of p3 wrt p0->p1 line segment
cp3 = np.cross(p1-p3, p0-p3)
return points[::-1, :]
# Group cross products in an array and find out how many are 0
cross_products = np.array([[cp0, cp1], [cp2, cp3]])
cps_near_zero = np.isclose(cross_products, 0).astype(int)
num_zeros = np.sum(cps_near_zero)
# Topologies of 2 finite line segments categorized by the number
# of zero-valued cross products:
#
# 0: No 3 points lie on the same line
# 1: 1 point lies on the same line defined by the other line
# segment, but is not coincident with either of the points
# 2: 2 points are coincident, but the line segments are not
# collinear which guarantees no intersection
# 3: not possible, except maybe floating point issues?
# 4: Both line segments are collinear, simply need to check if
# they overlap or not
# adapted from algorithm linked below and modified to only
# consider intersections on the interior of line segments as
# proper intersections: i.e. segments sharing end points do not
# count as intersections.
# https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
if num_zeros == 0:
# If the orientations of p0 and p1 have opposite signs
# and the orientations of p2 and p3 have opposite signs
# then there is an intersection.
if all(np.prod(cross_products, axis=-1) < 0):
raise ValueError('Polygon cannot be self-intersecting')
continue
elif num_zeros == 1:
# determine which line segment has 2 out of the 3 collinear
# points
idx = np.argwhere(np.sum(cps_near_zero, axis=-1) == 0)
if np.prod(cross_products[idx, :]) < 0:
raise ValueError('Polygon cannot be self-intersecting')
continue
elif num_zeros == 2:
continue
elif num_zeros == 3:
warnings.warn('Unclear if Polygon is self-intersecting')
continue
else:
# All 4 cross products are zero
# Determine number of unique points, x span and y span for
# both line segments
xmin1, xmax1 = min(p0[0], p1[0]), max(p0[0], p1[0])
ymin1, ymax1 = min(p0[1], p1[1]), max(p0[1], p1[1])
xmin2, xmax2 = min(p2[0], p3[0]), max(p2[0], p3[0])
ymin2, ymax2 = min(p2[1], p3[1]), max(p2[1], p3[1])
xlap = xmin1 < xmax2 and xmin2 < xmax1
ylap = ymin1 < ymax2 and ymin2 < ymax1
if xlap or ylap:
raise ValueError('Polygon cannot be self-intersecting')
continue
return points
def _constrain_triangulation(self, points, depth=0):
"""Generate a constrained triangulation by ensuring all edges of the
Polygon are contained within the simplices.
Parameters
----------
points : np.ndarray (Nx2)
An Nx2 array of coordinate pairs describing the vertices. These
points represent a planar straight line graph.
Returns
-------
None
"""
# Only attempt the triangulation up to 3 times.
if depth > 2:
raise RuntimeError('Could not create a valid triangulation after 3'
' attempts')
tri = Delaunay(points, qhull_options='QJ')
# Loop through the boundary edges of the polygon. If an edge is not
# included in the triangulation, break it into two line segments.
n = len(points)
new_pts = []
for i, j in zip(range(n), range(1, n +1)):
# If both vertices of any edge are not found in any simplex, insert
# a new point between them.
if not any([i in s and j % n in s for s in tri.simplices]):
newpt = (points[i, :] + points[j % n, :]) / 2
new_pts.append((j, newpt))
# If all the edges are included in the triangulation set it, otherwise
# try again with additional points inserted on offending edges.
if not new_pts:
self._tri = tri
else:
for i, pt in new_pts[::-1]:
points = np.insert(points, i, pt, axis=0)
self._constrain_triangulation(points, depth=depth + 1)
def _group_simplices(self, neighbor_map, group=None):
"""Generate a convex grouping of simplices.
@ -819,7 +939,7 @@ class Polygon(CompositeSurface):
continue
test_group = group + [n]
test_point_idx = np.unique(self._tri.simplices[test_group, :])
test_points = self._tri.points[test_point_idx]
test_points = self.points[test_point_idx]
# If test_points are convex keep adding to this group
if len(test_points) == len(ConvexHull(test_points).vertices):
group = self._group_simplices(neighbor_map, group=test_group)
@ -902,8 +1022,8 @@ class Polygon(CompositeSurface):
# Get centroids of all the simplices and determine if they are inside
# the polygon defined by input vertices or not.
centroids = np.mean(self._points[self._tri.simplices], axis=1)
in_polygon = Path(self._points).contains_points(centroids)
centroids = np.mean(self.points[self._tri.simplices], axis=1)
in_polygon = Path(self.points).contains_points(centroids)
self._in_polygon = in_polygon
# Build a map with keys of simplex indices inside the polygon whose
@ -931,7 +1051,7 @@ class Polygon(CompositeSurface):
# generate the convex hull and find the resulting surfaces and
# unary operators that represent this convex subset of the polygon.
idx = np.unique(self._tri.simplices[group, :])
qhull = ConvexHull(self._tri.points[idx, :])
qhull = ConvexHull(self.points[idx, :])
surf_ops = self._get_convex_hull_surfs(qhull)
surfsets.append(surf_ops)
return surfsets

54
tests/unit_tests/test_surface_composite.py
View file

@ -339,3 +339,57 @@ def test_polygon():
offset_star = star_poly.offset(.6)
assert (0, 0, 0) in -offset_star
assert any([(0, 0, 0) in reg for reg in offset_star.regions])
# check invalid Polygon input points
# duplicate points not just at start and end
rz_points = np.array([[6.88, 3.02],
[6.88, 2.72],
[6.88, 3.02],
[7.63, 0.0],
[5.75, 0.0],
[5.75, 1.22],
[7.63, 0.0],
[6.30, 1.22],
[6.30, 3.02],
[6.88, 3.02]])
with pytest.raises(ValueError):
openmc.model.Polygon(rz_points)
# segment traces back on previous segment
rz_points = np.array([[6.88, 3.02],
[6.88, 2.72],
[6.88, 2.32],
[6.88, 2.52],
[7.63, 0.0],
[5.75, 0.0],
[6.75, 0.0],
[5.75, 1.22],
[6.30, 1.22],
[6.30, 3.02],
[6.88, 3.02]])
with pytest.raises(ValueError):
openmc.model.Polygon(rz_points)
# segments intersect (line-line)
rz_points = np.array([[6.88, 3.02],
[5.88, 2.32],
[7.63, 0.0],
[5.75, 0.0],
[5.75, 1.22],
[6.30, 1.22],
[6.30, 3.02],
[6.88, 3.02]])
with pytest.raises(ValueError):
openmc.model.Polygon(rz_points)
# segments intersect (line-point)
rz_points = np.array([[6.88, 3.02],
[6.3, 2.32],
[7.63, 0.0],
[5.75, 0.0],
[5.75, 1.22],
[6.30, 1.22],
[6.30, 3.02],
[6.88, 3.02]])
with pytest.raises(ValueError):
openmc.model.Polygon(rz_points)