Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "GeometryUtils.h"
11 : #include "MooseUtils.h"
12 :
13 : namespace geom_utils
14 : {
15 :
16 : bool
17 2 : isPointZero(const Point & pt)
18 : {
19 2 : const Point zero(0.0, 0.0, 0.0);
20 4 : return pt.absolute_fuzzy_equals(zero);
21 : }
22 :
23 : Point
24 2 : unitVector(const Point & pt, const std::string & name)
25 : {
26 2 : if (isPointZero(pt))
27 3 : mooseError("'" + name + "' cannot have zero norm!");
28 :
29 1 : return pt.unit();
30 : }
31 :
32 : Real
33 0 : minDistanceToPoints(const Point & pt,
34 : const std::vector<Point> & candidates,
35 : const unsigned int axis)
36 : {
37 0 : const auto idx = projectedIndices(axis);
38 :
39 0 : Real distance = std::numeric_limits<Real>::max();
40 0 : for (const auto & c : candidates)
41 : {
42 0 : const Real dx = c(idx.first) - pt(idx.first);
43 0 : const Real dy = c(idx.second) - pt(idx.second);
44 0 : const Real d = std::sqrt(dx * dx + dy * dy);
45 0 : distance = std::min(d, distance);
46 : }
47 :
48 0 : return distance;
49 : }
50 :
51 : Point
52 12 : projectPoint(const Real x0, const Real x1, const unsigned int axis)
53 : {
54 12 : const auto i = projectedIndices(axis);
55 12 : Point point;
56 12 : point(i.first) = x0;
57 12 : point(i.second) = x1;
58 12 : point(axis) = 0.0;
59 :
60 24 : return point;
61 : }
62 :
63 : Real
64 51 : projectedLineHalfSpace(Point pt1, Point pt2, Point pt3, const unsigned int axis)
65 : {
66 : // project points onto plane perpendicular to axis
67 51 : pt1(axis) = 0.0;
68 51 : pt2(axis) = 0.0;
69 51 : pt3(axis) = 0.0;
70 :
71 51 : const auto i = projectedIndices(axis);
72 :
73 51 : return (pt1(i.first) - pt3(i.first)) * (pt2(i.second) - pt3(i.second)) -
74 51 : (pt2(i.first) - pt3(i.first)) * (pt1(i.second) - pt3(i.second));
75 : }
76 :
77 : bool
78 13 : pointInPolygon(const Point & point, const std::vector<Point> & corners, const unsigned int axis)
79 : {
80 13 : const auto n_pts = corners.size();
81 :
82 13 : std::vector<bool> negative_half_space;
83 13 : std::vector<bool> positive_half_space;
84 62 : for (const auto i : index_range(corners))
85 : {
86 49 : const int next = (i == n_pts - 1) ? 0 : i + 1;
87 49 : const auto half = projectedLineHalfSpace(point, corners[i], corners[next], axis);
88 49 : negative_half_space.push_back(half < 0);
89 49 : positive_half_space.push_back(half > 0);
90 : }
91 :
92 13 : const bool negative = std::find(negative_half_space.begin(), negative_half_space.end(), true) !=
93 13 : negative_half_space.end();
94 13 : const bool positive = std::find(positive_half_space.begin(), positive_half_space.end(), true) !=
95 13 : positive_half_space.end();
96 :
97 13 : const bool in_polygon = !(negative && positive);
98 13 : if (in_polygon)
99 8 : return true;
100 :
101 5 : if (pointOnEdge(point, corners, axis))
102 0 : return true;
103 :
104 5 : return false;
105 13 : }
106 :
107 : bool
108 5 : pointOnEdge(const Point & point, const std::vector<Point> & corners, const unsigned int axis)
109 : {
110 5 : const auto n_pts = corners.size();
111 5 : const auto idx = projectedIndices(axis);
112 :
113 5 : constexpr Real tol = 1e-8;
114 24 : for (const auto i : index_range(corners))
115 : {
116 19 : const int next = (i == n_pts - 1) ? 0 : i + 1;
117 19 : const auto & pt1 = corners[i];
118 19 : const auto & pt2 = corners[next];
119 19 : const bool close_to_line = projectedDistanceFromLine(point, pt1, pt2, axis) < tol;
120 :
121 : // we can stop early if we know we're not close to the line
122 19 : if (!close_to_line)
123 17 : continue;
124 :
125 : // check that the point is "between" the two points; TODO: first pass
126 : // we can just compare x and y coordinates
127 2 : const bool between_points = (point(idx.first) >= std::min(pt1(idx.first), pt2(idx.first))) &&
128 1 : (point(idx.first) <= std::max(pt1(idx.first), pt2(idx.first))) &&
129 3 : (point(idx.second) >= std::min(pt1(idx.second), pt2(idx.second))) &&
130 0 : (point(idx.second) <= std::max(pt1(idx.second), pt2(idx.second)));
131 :
132 : // point needs to be close to the line AND "between" the two points
133 2 : if (close_to_line && between_points)
134 0 : return true;
135 : }
136 :
137 5 : return false;
138 : }
139 :
140 : std::pair<unsigned int, unsigned int>
141 69 : projectedIndices(const unsigned int axis)
142 : {
143 69 : std::pair<unsigned int, unsigned int> indices;
144 :
145 69 : if (axis == 0)
146 : {
147 0 : indices.first = 1;
148 0 : indices.second = 2;
149 : }
150 69 : else if (axis == 1)
151 : {
152 0 : indices.first = 0;
153 0 : indices.second = 2;
154 : }
155 : else
156 : {
157 69 : indices.first = 0;
158 69 : indices.second = 1;
159 : }
160 :
161 69 : return indices;
162 : }
163 :
164 : Point
165 1 : projectedUnitNormal(Point pt1, Point pt2, const unsigned int axis)
166 : {
167 : // project the points to the plane perpendicular to the axis
168 1 : pt1(axis) = 0.0;
169 1 : pt2(axis) = 0.0;
170 :
171 1 : const auto i = projectedIndices(axis);
172 :
173 1 : const Real dx = pt2(i.first) - pt1(i.first);
174 1 : const Real dy = pt2(i.second) - pt1(i.second);
175 :
176 1 : const Point normal = projectPoint(dy, -dx, axis);
177 1 : const Point gap_line = pt2 - pt1;
178 :
179 1 : const auto cross_product = gap_line.cross(normal);
180 :
181 1 : if (cross_product(axis) > 0)
182 0 : return normal.unit();
183 : else
184 1 : return projectPoint(-dy, dx, axis).unit();
185 : }
186 :
187 : Real
188 28 : distanceFromLine(const Point & pt, const Point & line0, const Point & line1)
189 : {
190 28 : const Point a = pt - line0;
191 28 : const Point b = pt - line1;
192 28 : const Point c = line1 - line0;
193 :
194 28 : return (a.cross(b).norm()) / c.norm();
195 : }
196 :
197 : Real
198 22 : projectedDistanceFromLine(Point pt, Point line0, Point line1, const unsigned int axis)
199 : {
200 : // project all the points to the plane perpendicular to the axis
201 22 : pt(axis) = 0.0;
202 22 : line0(axis) = 0.0;
203 22 : line1(axis) = 0.0;
204 :
205 22 : return distanceFromLine(pt, line0, line1);
206 : }
207 :
208 : std::vector<Point>
209 2 : polygonCorners(const unsigned int num_sides, const Real radius, const unsigned int axis)
210 : {
211 2 : std::vector<Point> corners;
212 2 : const Real theta = 2.0 * M_PI / num_sides;
213 2 : const Real first_angle = M_PI / 2.0 - theta / 2.0;
214 :
215 12 : for (const auto i : make_range(num_sides))
216 : {
217 10 : const Real angle = first_angle + i * theta;
218 10 : const Real x = radius * cos(angle);
219 10 : const Real y = radius * sin(angle);
220 :
221 10 : corners.push_back(projectPoint(x, y, axis));
222 : }
223 :
224 2 : return corners;
225 0 : }
226 :
227 : Point
228 309 : rotatePointAboutAxis(const Point & p, const Real angle, const Point & axis)
229 : {
230 309 : const Real cos_theta = cos(angle);
231 309 : const Real sin_theta = sin(angle);
232 :
233 309 : Point pt;
234 309 : const Real xy = axis(0) * axis(1);
235 309 : const Real xz = axis(0) * axis(2);
236 309 : const Real yz = axis(1) * axis(2);
237 :
238 309 : const Point x_op(cos_theta + axis(0) * axis(0) * (1.0 - cos_theta),
239 618 : xy * (1.0 - cos_theta) - axis(2) * sin_theta,
240 618 : xz * (1.0 - cos_theta) + axis(1) * sin_theta);
241 :
242 618 : const Point y_op(xy * (1.0 - cos_theta) + axis(2) * sin_theta,
243 309 : cos_theta + axis(1) * axis(1) * (1.0 - cos_theta),
244 618 : yz * (1.0 - cos_theta) - axis(0) * sin_theta);
245 :
246 618 : const Point z_op(xz * (1.0 - cos_theta) - axis(1) * sin_theta,
247 618 : yz * (1.0 - cos_theta) + axis(0) * sin_theta,
248 309 : cos_theta + axis(2) * axis(2) * (1.0 - cos_theta));
249 :
250 309 : pt(0) = x_op * p;
251 309 : pt(1) = y_op * p;
252 309 : pt(2) = z_op * p;
253 618 : return pt;
254 : }
255 :
256 : std::vector<Point>
257 1 : boxCorners(const libMesh::BoundingBox & box, const Real factor)
258 : {
259 1 : Point diff = (box.max() - box.min()) / 2.0;
260 1 : const Point origin = box.min() + diff;
261 :
262 : // Rescale side length of box by specified factor
263 1 : diff *= factor;
264 :
265 : // Vectors for sides of box
266 1 : const Point dx(2.0 * diff(0), 0.0, 0.0);
267 1 : const Point dy(0.0, 2.0 * diff(1), 0.0);
268 1 : const Point dz(0.0, 0.0, 2.0 * diff(2));
269 :
270 1 : std::vector<Point> verts(8, origin - diff);
271 1 : const unsigned int pts_per_dim = 2;
272 3 : for (const auto z : make_range(pts_per_dim))
273 6 : for (const auto y : make_range(pts_per_dim))
274 12 : for (const auto x : make_range(pts_per_dim))
275 8 : verts[pts_per_dim * pts_per_dim * z + pts_per_dim * y + x] += x * dx + y * dy + z * dz;
276 :
277 2 : return verts;
278 0 : }
279 :
280 : } // end namespace geom_utils
|
