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 "TraceRayTools.h"
11 :
12 : // MOOSE includes
13 : #include "MooseTypes.h"
14 : #include "MooseUtils.h"
15 :
16 : // libMesh includes
17 : #include "libmesh/cell_hex8.h"
18 : #include "libmesh/cell_hex20.h"
19 : #include "libmesh/cell_hex27.h"
20 : #include "libmesh/cell_prism6.h"
21 : #include "libmesh/cell_prism15.h"
22 : #include "libmesh/cell_prism18.h"
23 : #include "libmesh/cell_pyramid5.h"
24 : #include "libmesh/cell_pyramid13.h"
25 : #include "libmesh/cell_pyramid14.h"
26 : #include "libmesh/cell_tet4.h"
27 : #include "libmesh/cell_tet10.h"
28 : #include "libmesh/cell_tet14.h"
29 : #include "libmesh/edge_edge2.h"
30 : #include "libmesh/enum_to_string.h"
31 : #include "libmesh/face_quad4.h"
32 : #include "libmesh/face_tri3.h"
33 : #include "libmesh/mesh.h"
34 : #include "libmesh/remote_elem.h"
35 : #include "libmesh/tensor_value.h"
36 :
37 : using namespace libMesh;
38 :
39 : namespace TraceRayTools
40 : {
41 :
42 : const std::set<int> TRACEABLE_ELEMTYPES = {
43 : HEX8, HEX20, HEX27, QUAD4, QUAD8, QUAD9, TET4, TET10, TET14, TRI3, TRI6,
44 : TRI7, EDGE2, EDGE3, EDGE4, PYRAMID5, PYRAMID13, PYRAMID14, PRISM6, PRISM15, PRISM18};
45 : const std::set<int> ADAPTIVITY_TRACEABLE_ELEMTYPES = {QUAD4, HEX8, TRI3, TET4, EDGE2};
46 :
47 : bool
48 3084612 : lineLineIntersect2D(const Point & start,
49 : const Point & direction,
50 : const Real length,
51 : const Point & v0,
52 : const Point & v1,
53 : Point & intersection_point,
54 : Real & intersection_distance,
55 : SegmentVertices & segment_vertex
56 : #ifdef DEBUG_RAY_INTERSECTIONS
57 : ,
58 : const bool debug
59 : #endif
60 : )
61 :
62 : {
63 : // TODO: consider using hmax scaling here
64 : mooseAssert(segment_vertex == SEGMENT_VERTEX_NONE, "Vertex should be none");
65 : debugRaySimple("Called lineLineIntersect2D()");
66 : debugRaySimple(" start = ", start);
67 : debugRaySimple(" direction = ", direction);
68 : debugRaySimple(" length = ", length);
69 : debugRaySimple(" v0 = ", v0);
70 : debugRaySimple(" v1 = ", v1);
71 :
72 : const auto r = direction * length;
73 : const auto s = v1 - v0;
74 :
75 3084612 : const auto rxs = r(0) * s(1) - r(1) * s(0);
76 : debugRaySimple(" rxs = ", rxs);
77 :
78 : // Lines are parallel or colinear
79 3084612 : if (std::abs(rxs) < TRACE_TOLERANCE)
80 : return false;
81 :
82 : const auto v0mu0 = v0 - start;
83 :
84 2695725 : const auto t = (v0mu0(0) * s(1) - v0mu0(1) * s(0)) / rxs;
85 : debugRaySimple(" t = ", t);
86 2695725 : if (0 >= t + TRACE_TOLERANCE || t - TRACE_TOLERANCE > 1.0)
87 : {
88 : debugRaySimple("lineLineIntersect2D did not intersect: t out of range");
89 : return false;
90 : }
91 :
92 1995229 : const auto u = (v0mu0(0) * r(1) - v0mu0(1) * r(0)) / rxs;
93 : debugRaySimple(" u = ", u);
94 1995229 : if (0 < u + TRACE_TOLERANCE && u - TRACE_TOLERANCE <= 1.0)
95 : {
96 1382509 : intersection_point = start + r * t;
97 1382509 : intersection_distance = t * length;
98 :
99 1382509 : if (u < TRACE_TOLERANCE)
100 359438 : segment_vertex = SEGMENT_VERTEX_0;
101 1023071 : else if (u > 1.0 - TRACE_TOLERANCE)
102 298531 : segment_vertex = SEGMENT_VERTEX_1;
103 :
104 : debugRaySimple("lineLineIntersect2D intersected with:");
105 : debugRaySimple(" intersection_distance = ", intersection_point);
106 : debugRaySimple(" intersection_distance = ", intersection_distance);
107 : debugRaySimple(" segment_vertex = ", Utility::enum_to_string(segment_vertex));
108 :
109 1382509 : return true;
110 : }
111 :
112 : // Not parallel, but don't intersect
113 : debugRaySimple("lineLineIntersect2d() did not intersect: u out of range");
114 : return false;
115 : }
116 :
117 : void
118 17104 : findPointNeighbors(
119 : const Elem * const elem,
120 : const Point & point,
121 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
122 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
123 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
124 : std::vector<const Elem *> active_neighbor_children,
125 : std::vector<NeighborInfo> & info)
126 : {
127 : mooseAssert(elem->contains_point(point), "Doesn't contain point");
128 :
129 17104 : info.clear();
130 :
131 : // Helper for avoiding extraneous allocation when building side elements
132 17104 : std::unique_ptr<const Elem> side_helper;
133 :
134 232856 : auto contains_point = [&point, &info, &side_helper, &elem](const Elem * const candidate)
135 : {
136 232856 : if (candidate->contains_point(point))
137 : {
138 : std::vector<unsigned short> sides;
139 492386 : for (const auto s : candidate->side_index_range())
140 : {
141 398470 : candidate->build_side_ptr(side_helper, s);
142 398470 : if (side_helper->contains_point(point))
143 229704 : sides.push_back(s);
144 : }
145 :
146 : // Dont add the local element
147 93916 : if (!sides.empty() && candidate != elem)
148 : {
149 76812 : info.emplace_back(candidate, std::move(sides));
150 : return true;
151 : }
152 93916 : }
153 :
154 : return false;
155 17104 : };
156 :
157 : // Fill info for the element that was passed in
158 17104 : contains_point(elem);
159 :
160 17104 : findNeighbors(elem,
161 : neighbor_set,
162 : untested_set,
163 : next_untested_set,
164 : active_neighbor_children,
165 : contains_point);
166 :
167 : #ifndef NDEBUG
168 : // In non-opt modes, verify that we found all of the correct neighbors
169 : // using the more expensive libMesh routine
170 : std::set<const Elem *> point_neighbors;
171 : elem->find_point_neighbors(point, point_neighbors);
172 : for (const auto & point_neighbor : point_neighbors)
173 : if (!neighbor_set.contains(point_neighbor) && point_neighbor != elem)
174 : mooseError("Missed a point neighbor");
175 : #endif
176 17104 : }
177 :
178 : void
179 105739 : findNodeNeighbors(
180 : const Elem * const elem,
181 : const Node * const node,
182 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
183 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
184 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
185 : std::vector<const Elem *> active_neighbor_children,
186 : std::vector<NeighborInfo> & info)
187 : {
188 : mooseAssert(elem->get_node_index(node) != libMesh::invalid_uint, "Doesn't contain node");
189 :
190 105739 : info.clear();
191 :
192 : // Helper for avoiding extraneous allocations when building side elements
193 105739 : std::unique_ptr<const Elem> side_helper;
194 :
195 1793872 : auto contains_node = [&node, &elem, &info, &side_helper](const Elem * const candidate)
196 : {
197 : // Candidate has this node and it is a vertex - add sides that contain said node
198 1793872 : const auto n = candidate->get_node_index(node);
199 1793872 : if (n != invalid_uint && candidate->is_vertex(n))
200 : {
201 : std::vector<unsigned short> sides;
202 4119134 : for (const auto s : candidate->side_index_range())
203 3410002 : if (candidate->is_node_on_side(n, s))
204 2135910 : sides.push_back(s);
205 :
206 709132 : if (sides.empty())
207 0 : mooseError("Failed to find a side containing node");
208 :
209 709132 : info.emplace_back(candidate, std::move(sides));
210 : return true;
211 709132 : }
212 : // In the case of a less refined candidate, the node can be a hanging node. The candidate
213 : // will only ever have the hanging node if it is less refined.
214 1084740 : if (candidate->level() < elem->level() && candidate->contains_point(*node))
215 : {
216 : std::vector<unsigned short> sides;
217 6368 : for (const auto s : candidate->side_index_range())
218 : {
219 5302 : candidate->build_side_ptr(side_helper, s);
220 5302 : if (side_helper->contains_point(*node))
221 1356 : sides.push_back(s);
222 : }
223 :
224 1066 : if (!sides.empty())
225 : {
226 1066 : info.emplace_back(candidate, std::move(sides));
227 : return true;
228 : }
229 1066 : }
230 :
231 : return false;
232 105739 : };
233 :
234 : // Fill info for the element that was passed in
235 105739 : contains_node(elem);
236 :
237 105739 : findNeighbors(
238 : elem, neighbor_set, untested_set, next_untested_set, active_neighbor_children, contains_node);
239 :
240 : #ifndef NDEBUG
241 : // In non-opt modes, verify that we found all of the correct neighbors
242 : // using the more expensive libMesh routine
243 : std::set<const Elem *> point_neighbors;
244 : elem->find_point_neighbors(*node, point_neighbors);
245 : for (const auto & point_neighbor : point_neighbors)
246 : for (const auto & neighbor_node : point_neighbor->node_ref_range())
247 : if (node == &neighbor_node && !neighbor_set.contains(point_neighbor) &&
248 : point_neighbor != elem)
249 : mooseError("Missed a node neighbor");
250 : #endif
251 105739 : }
252 :
253 : void
254 211605 : findEdgeNeighbors(
255 : const Elem * const elem,
256 : const Node * const node1,
257 : const Node * const node2,
258 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
259 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
260 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
261 : std::vector<const Elem *> active_neighbor_children,
262 : std::vector<NeighborInfo> & info)
263 : {
264 : mooseAssert(elem->get_node_index(node1) != libMesh::invalid_uint, "Doesn't contain node");
265 : mooseAssert(elem->get_node_index(node2) != libMesh::invalid_uint, "Doesn't contain node");
266 :
267 211605 : info.clear();
268 :
269 : // The length for this edge used for checking if a point is contained within said edge
270 211605 : const Real edge_length = ((Point)*node1 - (Point)*node2).norm();
271 :
272 : // Lambda that returns whether or not a candidate element contains an edge that is within the
273 : // target edge defined by node1 and node2. Also fills "info" if a match is found
274 2897955 : auto within_edge = [&elem, &node1, &node2, &edge_length, &info](const Elem * const candidate)
275 : {
276 2897955 : switch (candidate->type())
277 : {
278 283564 : case HEX8:
279 283564 : return findEdgeNeighborsWithinEdgeInternal<Hex8>(
280 283564 : candidate, elem, node1, node2, edge_length, info);
281 275726 : case TET4:
282 275726 : return findEdgeNeighborsWithinEdgeInternal<Tet4>(
283 275726 : candidate, elem, node1, node2, edge_length, info);
284 396728 : case PYRAMID5:
285 396728 : return findEdgeNeighborsWithinEdgeInternal<Pyramid5>(
286 396728 : candidate, elem, node1, node2, edge_length, info);
287 120409 : case PRISM6:
288 120409 : return findEdgeNeighborsWithinEdgeInternal<Prism6>(
289 120409 : candidate, elem, node1, node2, edge_length, info);
290 124108 : case HEX20:
291 124108 : return findEdgeNeighborsWithinEdgeInternal<Hex20>(
292 124108 : candidate, elem, node1, node2, edge_length, info);
293 124091 : case HEX27:
294 124091 : return findEdgeNeighborsWithinEdgeInternal<Hex27>(
295 124091 : candidate, elem, node1, node2, edge_length, info);
296 271006 : case TET10:
297 271006 : return findEdgeNeighborsWithinEdgeInternal<Tet10>(
298 271006 : candidate, elem, node1, node2, edge_length, info);
299 270876 : case TET14:
300 270876 : return findEdgeNeighborsWithinEdgeInternal<Tet14>(
301 270876 : candidate, elem, node1, node2, edge_length, info);
302 395462 : case PYRAMID13:
303 395462 : return findEdgeNeighborsWithinEdgeInternal<Pyramid13>(
304 395462 : candidate, elem, node1, node2, edge_length, info);
305 395242 : case PYRAMID14:
306 395242 : return findEdgeNeighborsWithinEdgeInternal<Pyramid14>(
307 395242 : candidate, elem, node1, node2, edge_length, info);
308 120413 : case PRISM15:
309 120413 : return findEdgeNeighborsWithinEdgeInternal<Prism15>(
310 120413 : candidate, elem, node1, node2, edge_length, info);
311 120330 : case PRISM18:
312 120330 : return findEdgeNeighborsWithinEdgeInternal<Prism18>(
313 120330 : candidate, elem, node1, node2, edge_length, info);
314 0 : default:
315 0 : mooseError("Element type ",
316 0 : Utility::enum_to_string(candidate->type()),
317 : " not supported in TraceRayTools::findEdgeNeighbors()");
318 : }
319 211605 : };
320 :
321 : // Fill info for the element that was passed in
322 211605 : within_edge(elem);
323 :
324 211605 : findNeighbors(
325 : elem, neighbor_set, untested_set, next_untested_set, active_neighbor_children, within_edge);
326 211605 : }
327 :
328 : const Elem *
329 9273 : childContainingPointOnSide(const Elem * elem, const Point & point, const unsigned short side)
330 : {
331 : mooseAssert(!elem->active(), "Should be inactive");
332 : mooseAssert(elem->side_ptr(side)->contains_point(point), "Side should contain point");
333 :
334 21958 : for (unsigned int c = 0; c < elem->n_children(); ++c)
335 : {
336 21958 : if (!elem->is_child_on_side(c, side))
337 6615 : continue;
338 :
339 : const auto child = elem->child_ptr(c);
340 : // Experience shows that we need to loosen this tolerance just a little
341 : // The default is libMesh::TOLERANCE = 1e-6
342 15343 : if (child->close_to_point(point, 5.e-5))
343 : {
344 : if (child->active())
345 : return child;
346 : else
347 0 : return childContainingPointOnSide(child, point, side);
348 : }
349 : }
350 :
351 0 : mooseError("Failed to find child containing point on side");
352 : }
353 :
354 : const Elem *
355 189283 : getActiveNeighbor(const Elem * elem, const unsigned short side, const Point & point)
356 : {
357 189283 : const auto neighbor = elem->neighbor_ptr(side);
358 189283 : if (!neighbor || neighbor->active())
359 : return neighbor;
360 :
361 : // There is adaptivity... need to find the active child that contains the point
362 9273 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
363 9273 : return childContainingPointOnSide(neighbor, point, neighbor_side);
364 : }
365 :
366 : bool
367 239549538 : intersectTriangle(const Point & start,
368 : const Point & direction,
369 : const Elem * const elem,
370 : const unsigned short v0,
371 : const unsigned short v1,
372 : const unsigned short v2,
373 : Real & intersection_distance,
374 : ElemExtrema & intersected_extrema,
375 : const Real hmax
376 : #ifdef DEBUG_RAY_INTERSECTIONS
377 : ,
378 : const bool debug
379 : #endif
380 : )
381 : {
382 : debugRaySimple("intersectTriangle() called:");
383 : debugRaySimple(" start = ", start);
384 : debugRaySimple(" direction = ", direction);
385 : debugRaySimple(" elem->id() = ", elem->id());
386 : debugRaySimple(" v0 = ", v0, " at ", elem->point(v0));
387 : debugRaySimple(" v1 = ", v1, " at ", elem->point(v1));
388 : debugRaySimple(" v2 = ", v2, " at ", elem->point(v2));
389 : debugRaySimple(" hmax = ", hmax);
390 : mooseAssert(elem->is_vertex(v0), "Not a vertex");
391 : mooseAssert(elem->is_vertex(v1), "Not a vertex");
392 : mooseAssert(elem->is_vertex(v2), "Not a vertex");
393 :
394 : // We are scaling the whole element (start, v0, v1, v2) by 1 / hmax as an alternative to scaling
395 : // the tolerance by hmax. If an intersection is found, the resulting intersection distance is
396 : // then scaled by hmax to reverse the original scaling.
397 239549538 : const auto inv_hmax = 1.0 / hmax;
398 :
399 : const auto & v0_point = elem->point(v0);
400 :
401 : const auto edge1 = (elem->point(v1) - v0_point) * inv_hmax;
402 : const auto edge2 = (elem->point(v2) - v0_point) * inv_hmax;
403 :
404 239549538 : const auto pvec = direction.cross(edge2);
405 :
406 : auto det = edge1 * pvec;
407 : debugRaySimple(" det = ", det);
408 239549538 : if (det < TRACE_TOLERANCE)
409 : {
410 : debugRaySimple("intersectTriangle() did not intersect: det < tol");
411 : return false;
412 : }
413 :
414 : const auto tvec = (start - v0_point) * inv_hmax;
415 : const auto u = tvec * pvec;
416 : debugRaySimple(" u = ", u);
417 : debugRaySimple(" u / det = ", u / det);
418 124012804 : if (u < -TRACE_TOLERANCE || u > det + TRACE_TOLERANCE)
419 : {
420 : debugRaySimple("intersectTriangle() did not intersect: u out of range");
421 : return false;
422 : }
423 :
424 70369057 : const auto qvec = tvec.cross(edge1);
425 : const auto v = direction * qvec;
426 : debugRaySimple(" v = ", v);
427 : debugRaySimple(" v / det = ", v / det);
428 : debugRaySimple(" (u + v) / det = ", (u + v) / det);
429 70369057 : if (v < -TRACE_TOLERANCE || u + v > det + TRACE_TOLERANCE)
430 : {
431 : debugRaySimple("intersectTriangle() did not intersect: v out of range");
432 : return false;
433 : }
434 :
435 52293630 : const auto possible_distance = (edge2 * qvec) / det;
436 : debugRaySimple(" possible_distance = ", possible_distance);
437 52293630 : if (possible_distance <= TRACE_TOLERANCE)
438 : {
439 : debugRaySimple("intersectTriangle() did not intersect: distance too small");
440 : return false;
441 : }
442 :
443 : // Recall that the element was scaled by (1 / hmax), reverse this scaling by
444 : // scaling the intersection distance by hmax
445 40724078 : intersection_distance = possible_distance * hmax;
446 :
447 : // Here, u and v aren't truly u and v. The actual u and v are obtained with:
448 : // u = u / det and v = v / det -> move det to the RHS to avoid division
449 40724078 : if (u < TRACE_TOLERANCE * det)
450 : {
451 3861991 : if (v < TRACE_TOLERANCE * det) // u = 0, v = 0
452 : intersected_extrema.setVertex(v0);
453 2668717 : else if (v > (1.0 - TRACE_TOLERANCE) * det) // u = 0, v = 1
454 : intersected_extrema.setVertex(v2);
455 : else // u = 0
456 : intersected_extrema.setEdge(v0, v2);
457 : }
458 36862087 : else if (v < TRACE_TOLERANCE * det)
459 : {
460 4156255 : if (u > (1.0 - TRACE_TOLERANCE) * det) // u = 1, v = 0
461 : intersected_extrema.setVertex(v1);
462 : else // v = 0
463 : intersected_extrema.setEdge(v0, v1);
464 : }
465 32705832 : else if ((u + v > (1.0 - TRACE_TOLERANCE) * det)) // u + v = 1
466 : intersected_extrema.setEdge(v1, v2);
467 :
468 : debugRaySimple("intersectTriangle() intersected with:");
469 : debugRaySimple(" intersection_distance = ", intersection_distance);
470 : debugRaySimple(" intersected_extrema = ", intersected_extrema);
471 :
472 : return true;
473 : }
474 :
475 : bool
476 70477995 : intersectQuad(const Point & start,
477 : const Point & direction,
478 : const Elem * const elem,
479 : const unsigned short v00,
480 : const unsigned short v10,
481 : const unsigned short v11,
482 : const unsigned short v01,
483 : Real & intersection_distance,
484 : ElemExtrema & intersected_extrema,
485 : const Real hmax
486 : #ifdef DEBUG_RAY_INTERSECTIONS
487 : ,
488 : const bool debug
489 : #endif
490 : )
491 : {
492 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
493 : debugRaySimple("intersectQuad() called:");
494 : debugRaySimple(" start = ", start);
495 : debugRaySimple(" direction = ", direction);
496 : debugRaySimple(" elem->id() = ", elem->id());
497 : debugRaySimple(" v00 = ", v00, " at ", elem->point(v00));
498 : debugRaySimple(" v10 = ", v10, " at ", elem->point(v10));
499 : debugRaySimple(" v11 = ", v11, " at ", elem->point(v11));
500 : debugRaySimple(" v01 = ", v01, " at ", elem->point(v01));
501 :
502 : // NOTE discovered by @GiudGiud: In the case that you have
503 : // a glancing intersection (the direction is within the plane
504 : // of the face), you could possibly miss a further intersection
505 : // on the second triangle. In reality, this should not be
506 : // a problem because we check all sides of the element, so
507 : // we would find a further intersection in the future.
508 :
509 : // First check the triangle contained by v00, v10, v11
510 70477995 : bool intersects = intersectTriangle(start,
511 : direction,
512 : elem,
513 : v00,
514 : v10,
515 : v11,
516 : intersection_distance,
517 : intersected_extrema,
518 : hmax
519 : #ifdef DEBUG_RAY_INTERSECTIONS
520 : ,
521 : debug
522 : #endif
523 : );
524 : // If no intersection, check the triangle contained by v11, v01, v00
525 70477995 : if (!intersects)
526 62171726 : intersects = intersectTriangle(start,
527 : direction,
528 : elem,
529 : v11,
530 : v01,
531 : v00,
532 : intersection_distance,
533 : intersected_extrema,
534 : hmax
535 : #ifdef DEBUG_RAY_INTERSECTIONS
536 : ,
537 : debug
538 : #endif
539 : );
540 :
541 : // Because we split the quad into two triangles, we could intersect the edge v00 - v11. However,
542 : // that really isn't an edge - it's a diagonal across the quad. If we intersect this edge, be sure
543 : // to invalidate it
544 62171726 : if (intersects && intersected_extrema.atEdge(v00, v11))
545 : intersected_extrema.invalidate();
546 :
547 70477995 : return intersects;
548 : }
549 :
550 : bool
551 456868 : isTraceableElem(const Elem * elem)
552 : {
553 456868 : return TRACEABLE_ELEMTYPES.count(elem->type());
554 : }
555 :
556 : bool
557 5670 : isAdaptivityTraceableElem(const Elem * elem)
558 : {
559 5670 : return ADAPTIVITY_TRACEABLE_ELEMTYPES.count(elem->type());
560 : }
561 :
562 : unsigned short
563 5987 : atVertex(const Elem * elem, const Point & point)
564 : {
565 19025 : for (unsigned int v = 0; v < elem->n_vertices(); ++v)
566 18161 : if (elem->point(v).absolute_fuzzy_equals(point, TRACE_TOLERANCE))
567 5123 : return v;
568 :
569 : return RayTracingCommon::invalid_vertex;
570 : }
571 :
572 : template <typename T>
573 : bool
574 8833 : withinEdgeTempl(const Elem * elem,
575 : const Point & point,
576 : ElemExtrema & extrema,
577 : const Real tolerance /* = TRACE_TOLERANCE */)
578 : {
579 : mooseAssert(extrema.isInvalid(), "Should be invalid");
580 :
581 29353 : for (int e = 0; e < T::num_edges; ++e)
582 29353 : if (isWithinSegment(elem->point(T::edge_nodes_map[e][0]),
583 29353 : elem->point(T::edge_nodes_map[e][1]),
584 : point,
585 : tolerance))
586 : {
587 8833 : extrema.setEdge(T::edge_nodes_map[e][0], T::edge_nodes_map[e][1]);
588 8833 : return true;
589 : }
590 :
591 : return false;
592 : }
593 :
594 : bool
595 8833 : withinEdge(const Elem * elem,
596 : const Point & point,
597 : ElemExtrema & extrema,
598 : const Real tolerance /* = TRACE_TOLERANCE */)
599 : {
600 8833 : switch (elem->type())
601 : {
602 481 : case HEX8:
603 : case HEX20:
604 : case HEX27:
605 481 : return withinEdgeTempl<Hex8>(elem, point, extrema, tolerance);
606 5760 : case TET4:
607 : case TET10:
608 : case TET14:
609 5760 : return withinEdgeTempl<Tet4>(elem, point, extrema, tolerance);
610 1872 : case PYRAMID5:
611 : case PYRAMID13:
612 : case PYRAMID14:
613 1872 : return withinEdgeTempl<Pyramid5>(elem, point, extrema, tolerance);
614 720 : case PRISM6:
615 : case PRISM15:
616 : case PRISM18:
617 720 : return withinEdgeTempl<Prism6>(elem, point, extrema, tolerance);
618 0 : default:
619 0 : mooseError("Element type ",
620 0 : Utility::enum_to_string(elem->type()),
621 : " not supported in TraceRayTools::withinEdge()");
622 : }
623 :
624 : return false;
625 : }
626 :
627 : unsigned short
628 2177867 : atVertexOnSide(const Elem * elem, const Point & point, const unsigned short side)
629 : {
630 2177867 : switch (elem->type())
631 : {
632 1030639 : case HEX8:
633 : case HEX20:
634 : case HEX27:
635 1030639 : return atVertexOnSideTempl<Hex8>(elem, point, side);
636 76649 : case QUAD4:
637 : case QUAD8:
638 : case QUAD9:
639 76649 : return atVertexOnSideTempl<Quad4>(elem, point, side);
640 5739 : case TRI3:
641 : case TRI6:
642 : case TRI7:
643 5739 : return atVertexOnSideTempl<Tri3>(elem, point, side);
644 317316 : case TET4:
645 : case TET10:
646 : case TET14:
647 317316 : return atVertexOnSideTempl<Tet4>(elem, point, side);
648 220777 : case PYRAMID5:
649 : case PYRAMID13:
650 : case PYRAMID14:
651 220777 : return atVertexOnSideTempl<Pyramid5>(elem, point, side);
652 526735 : case PRISM6:
653 : case PRISM15:
654 : case PRISM18:
655 526735 : return atVertexOnSideTempl<Prism6>(elem, point, side);
656 12 : case EDGE2:
657 : case EDGE3:
658 : case EDGE4:
659 12 : return atVertexOnSideTempl<Edge2>(elem, point, side);
660 0 : default:
661 0 : mooseError("Element type ",
662 0 : Utility::enum_to_string(elem->type()),
663 : " not supported in TraceRayTools::atVertexOnSide()");
664 : }
665 :
666 : return RayTracingCommon::invalid_vertex;
667 : }
668 :
669 : template <typename T>
670 : typename std::enable_if<!std::is_base_of<Edge, T>::value, unsigned short>::type
671 2177855 : atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side)
672 : {
673 : mooseAssert(side < elem->n_sides(), "Invalid side");
674 : mooseAssert(elem->side_ptr(side)->close_to_point(point, LOOSE_TRACE_TOLERANCE),
675 : "Side does not contain point");
676 :
677 8344092 : for (int i = 0; i < nodesPerSide<T>(side); ++i)
678 6465987 : if (elem->point(T::side_nodes_map[side][i]).absolute_fuzzy_equals(point, TRACE_TOLERANCE))
679 1221073 : return T::side_nodes_map[side][i];
680 :
681 : return RayTracingCommon::invalid_vertex;
682 : }
683 :
684 : template <typename T>
685 : typename std::enable_if<std::is_base_of<Edge, T>::value, unsigned short>::type
686 12 : atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side)
687 : {
688 : mooseAssert(side < elem->n_sides(), "Invalid side");
689 : mooseAssert(elem->side_ptr(side)->close_to_point(point, LOOSE_TRACE_TOLERANCE),
690 : "Side does not contain point");
691 :
692 12 : if (elem->point(side).absolute_fuzzy_equals(point, TRACE_TOLERANCE))
693 12 : return side;
694 :
695 : return RayTracingCommon::invalid_vertex;
696 : }
697 :
698 : bool
699 977974 : withinEdgeOnSide(const Elem * const elem,
700 : const Point & point,
701 : const unsigned short side,
702 : ElemExtrema & extrema)
703 : {
704 977974 : switch (elem->type())
705 : {
706 398934 : case HEX8:
707 : case HEX20:
708 : case HEX27:
709 398934 : return withinEdgeOnSideTempl<Hex8>(elem, point, side, extrema);
710 154580 : case TET4:
711 : case TET10:
712 : case TET14:
713 154580 : return withinEdgeOnSideTempl<Tet4>(elem, point, side, extrema);
714 150760 : case PYRAMID5:
715 : case PYRAMID13:
716 : case PYRAMID14:
717 150760 : return withinEdgeOnSideTempl<Pyramid5>(elem, point, side, extrema);
718 273700 : case PRISM6:
719 : case PRISM15:
720 : case PRISM18:
721 273700 : return withinEdgeOnSideTempl<Prism6>(elem, point, side, extrema);
722 0 : default:
723 0 : mooseError("Element type ",
724 0 : Utility::enum_to_string(elem->type()),
725 : " not supported in TraceRayTools::withinEdgeOnSide()");
726 : }
727 :
728 : return false;
729 : }
730 :
731 : template <typename T>
732 : typename std::enable_if<std::is_base_of<Cell, T>::value, bool>::type
733 977974 : withinEdgeOnSideTempl(const Elem * const elem,
734 : const Point & point,
735 : const unsigned short side,
736 : ElemExtrema & extrema)
737 : {
738 : mooseAssert(side < elem->n_sides(), "Invalid side");
739 : mooseAssert(elem->side_ptr(side)->close_to_point(point, LOOSE_TRACE_TOLERANCE),
740 : "Side does not contain point");
741 : mooseAssert(extrema.isInvalid(), "Should be invalid");
742 :
743 977974 : int last_n = T::side_nodes_map[side][nodesPerSide<T>(side) - 1];
744 :
745 2300102 : for (int side_v = 0; side_v < nodesPerSide<T>(side); ++side_v)
746 2296694 : if (isWithinSegment(elem->point(last_n), elem->point(T::side_nodes_map[side][side_v]), point))
747 : {
748 974566 : extrema.setEdge(last_n, T::side_nodes_map[side][side_v]);
749 : mooseAssert(extrema.buildEdge(elem)->close_to_point(point, LOOSE_TRACE_TOLERANCE),
750 : "Edge doesn't contain point");
751 974566 : return true;
752 : }
753 : else
754 1322128 : last_n = T::side_nodes_map[side][side_v];
755 :
756 : return false;
757 : }
758 :
759 : bool
760 2148372 : withinExtremaOnSide(const Elem * const elem,
761 : const Point & point,
762 : const unsigned short side,
763 : const unsigned int dim,
764 : ElemExtrema & extrema)
765 : {
766 : mooseAssert(extrema.isInvalid(), "Extrema should be invalid");
767 : mooseAssert(dim == elem->dim(), "Incorrect dim");
768 :
769 2148372 : extrema.first = atVertexOnSide(elem, point, side);
770 : if (extrema.atVertex())
771 : return true;
772 953254 : if (dim == 3 && withinEdgeOnSide(elem, point, side, extrema))
773 : return true;
774 :
775 : return false;
776 : }
777 :
778 : bool
779 2429038 : isWithinSegment(const Point & segment1,
780 : const Point & segment2,
781 : const Point & point,
782 : const Real tolerance /* = TRACE_TOLERANCE */)
783 : {
784 : mooseAssert(!segment1.absolute_fuzzy_equals(segment2, TRACE_TOLERANCE), "Same endpoints");
785 :
786 2429038 : const auto segment_length = (segment1 - segment2).norm();
787 2429038 : return isWithinSegment(segment1, segment2, segment_length, point, tolerance);
788 : }
789 :
790 : bool
791 12050111 : isWithinSegment(const Point & segment1,
792 : const Point & segment2,
793 : const Real segment_length,
794 : const Point & point,
795 : const Real tolerance /* = TRACE_TOLERANCE */)
796 : {
797 : mooseAssert(!segment1.absolute_fuzzy_equals(segment2, TRACE_TOLERANCE), "Same endpoints");
798 : mooseAssert(MooseUtils::absoluteFuzzyEqual((segment1 - segment2).norm(), segment_length),
799 : "Invalid segment length");
800 :
801 : const auto diff1 = point - segment1;
802 : const auto diff2 = point - segment2;
803 :
804 12050111 : if (diff1 * diff2 > tolerance * segment_length)
805 : return false;
806 :
807 1443029 : return std::abs(diff1.norm() + diff2.norm() - segment_length) < tolerance * segment_length;
808 : }
809 :
810 : bool
811 4116843 : onBoundingBoxBoundary(const BoundingBox & bbox,
812 : const Point & point,
813 : const unsigned int dim,
814 : const Real tolerance)
815 : {
816 14039555 : for (unsigned int d = 0; d < dim; ++d)
817 11039902 : if (MooseUtils::absoluteFuzzyEqual(point(d), bbox.min()(d), tolerance) ||
818 : MooseUtils::absoluteFuzzyEqual(point(d), bbox.max()(d), tolerance))
819 : return true;
820 :
821 : return false;
822 : }
823 : }
|
