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 : #pragma once
11 :
12 : // MOOSE includes
13 : #include "StaticallyAllocatedSet.h"
14 : #include "Conversion.h"
15 : #include "MooseTypes.h"
16 : #include "libMeshReducedNamespace.h"
17 :
18 : // Local includes
19 : #include "DebugRay.h"
20 : #include "ElemExtrema.h"
21 : #include "RayTracingCommon.h"
22 : #include "NeighborInfo.h"
23 :
24 : // libMesh includes
25 : #include "libmesh/face.h"
26 : #include "libmesh/cell_hex.h"
27 : #include "libmesh/cell_prism.h"
28 : #include "libmesh/cell_pyramid.h"
29 : #include "libmesh/cell_tet.h"
30 : #include "libmesh/remote_elem.h"
31 :
32 : namespace libMesh
33 : {
34 : class Mesh;
35 : class Edge;
36 : class Cell;
37 : }
38 :
39 : namespace TraceRayTools
40 : {
41 :
42 : /// The element types that are traceable
43 : extern const std::set<int> TRACEABLE_ELEMTYPES;
44 : /// The element types that are traceable with adaptivity
45 : extern const std::set<int> ADAPTIVITY_TRACEABLE_ELEMTYPES;
46 :
47 : /// The standard tolerance to use in tracing
48 : const Real TRACE_TOLERANCE = 1e-8;
49 : /// Looser tolerance for use in error checking in difficult situations
50 : const Real LOOSE_TRACE_TOLERANCE = 1e-5;
51 :
52 : /// Value for an invalid integer
53 : const int INVALID_INT = std::numeric_limits<int>::max();
54 :
55 : /**
56 : * Enum for the possible vertices on a segment used in lineLineIntersect2D()
57 : */
58 : enum SegmentVertices
59 : {
60 : SEGMENT_VERTEX_0 = 0,
61 : SEGMENT_VERTEX_1 = 1,
62 : SEGMENT_VERTEX_NONE = 7
63 : };
64 :
65 : /**
66 : * Checks for the intersection of the line u0 -> u1 with the line v0 -> v1,
67 : * where u0 = start and u1 = start + direction * length.
68 : *
69 : * From: https://stackoverflow.com/a/565282
70 : *
71 : * @return Whether or not the lines intersect
72 : *
73 : * @param start The start point, u0
74 : * @param direction The direction used to define u1
75 : * @param length The length used to define u1
76 : * @param v0 The point v0
77 : * @param v1 The point v1
78 : * @param intersection_point To be filled with the point of intersection
79 : * @param intersection_distance To be filled with the distance of intersection
80 : * @param segment_vertex To be filled with the intersected vertex, if any
81 : */
82 : bool lineLineIntersect2D(const Point & start,
83 : const Point & direction,
84 : const Real length,
85 : const Point & v0,
86 : const Point & v1,
87 : Point & intersection_point,
88 : Real & intersection_distance,
89 : SegmentVertices & segment_vertex
90 : #ifdef DEBUG_RAY_INTERSECTIONS
91 : ,
92 : const bool debug
93 : #endif
94 : );
95 :
96 : /**
97 : * Checks whether or not a point is within a line segment
98 : * @param segment1 The first point on the segment
99 : * @param segment2 The second point on the segment
100 : * @param point The point
101 : * @param tolerance The tolerance to use
102 : * @return If point is within the segment defined by [segment1, segment2]
103 : */
104 : bool isWithinSegment(const Point & segment1,
105 : const Point & segment2,
106 : const Point & point,
107 : const Real tolerance = TRACE_TOLERANCE);
108 :
109 : /**
110 : * Checks whether or not a point is within a line segment
111 : * @param segment1 The first point on the segment
112 : * @param segment2 The second point on the segment
113 : * @param segment_length The segment length (for optimization if it's already computed)
114 : * @param point The point
115 : * @param tolerance The tolerance to use
116 : * @return If point is within the segment defined by [segment1, segment2]
117 : */
118 : bool isWithinSegment(const Point & segment1,
119 : const Point & segment2,
120 : const Real segment_length,
121 : const Point & point,
122 : const Real tolerance = TRACE_TOLERANCE);
123 :
124 : /**
125 : * Rewrite of the find_point_neighbors function in libMesh, instead using a statically allocated
126 : * set: returns the active point neighbors at p within elem.
127 : *
128 : * @param elem The element
129 : * @param p The point
130 : * @param neighbor_set The set to fill the neighbors into
131 : * @param untested_set Set for internal use
132 : * @param next_untested_set Set for internal use
133 : * @param active_neighbor_children Temporary vector for use in the search
134 : */
135 : void findPointNeighbors(
136 : const Elem * const elem,
137 : const Point & point,
138 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
139 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
140 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
141 : std::vector<const Elem *> active_neighbor_children,
142 : std::vector<NeighborInfo> & info);
143 :
144 : void findNodeNeighbors(
145 : const Elem * const elem,
146 : const Node * const node,
147 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
148 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
149 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
150 : std::vector<const Elem *> active_neighbor_children,
151 : std::vector<NeighborInfo> & info);
152 :
153 : void findEdgeNeighbors(
154 : const Elem * const elem,
155 : const Node * const node1,
156 : const Node * const node2,
157 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
158 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
159 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
160 : std::vector<const Elem *> active_neighbor_children,
161 : std::vector<NeighborInfo> & info);
162 :
163 : /**
164 : * More generalized form of the find_point_neighbors function in libMesh. Instead uses a statically
165 : * allocated set and accepts a functor for the rejection/acceptance of an element.
166 : *
167 : * Returns the active neighbors that fit the criteria of keep_functor.
168 : * Does not the return the current element (elem)
169 : *
170 : * @param elem The element for which we are searching for the neighbors
171 : * @param neighbor_set The set to fill the neighbors into
172 : * @param untested_set Set for internal use
173 : * @param next_untested_set Set for internal use
174 : * @param active_neighbor_children Vector to use in the search
175 : * @param keep_functor The functor that is used to accept or reject an element
176 : */
177 : template <typename KeepFunctor>
178 : void
179 334448 : findNeighbors(
180 : const Elem * const elem,
181 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
182 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
183 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
184 : std::vector<const Elem *> active_neighbor_children,
185 : KeepFunctor & keep_functor)
186 : {
187 : mooseAssert(elem->active(), "Inactive element");
188 :
189 : neighbor_set.clear();
190 : untested_set.clear();
191 : next_untested_set.clear();
192 :
193 334448 : untested_set.insert(elem);
194 :
195 1364027 : while (!untested_set.empty())
196 : {
197 : // Loop over all the elements in the patch that haven't already been tested
198 2662803 : for (const Elem * const test_elem : untested_set)
199 9426509 : for (const auto * const current_neighbor : test_elem->neighbor_ptr_range())
200 7793285 : if (current_neighbor && current_neighbor != remote_elem &&
201 6994339 : current_neighbor != elem) // we have a real neighbor on elem side
202 : {
203 : if (current_neighbor->active()) // ... if it is active
204 : {
205 6337056 : if (!neighbor_set.contains(current_neighbor) && keep_functor(current_neighbor))
206 : {
207 1296504 : next_untested_set.insert(current_neighbor);
208 1296504 : neighbor_set.insert(current_neighbor);
209 : }
210 : }
211 : #ifdef LIBMESH_ENABLE_AMR
212 : else // ... the neighbor is *not* active,
213 : { // ... so add *all* neighboring
214 : // active children that touch p
215 8576 : current_neighbor->active_family_tree_by_neighbor(active_neighbor_children, test_elem);
216 :
217 36155 : for (const Elem * current_child : active_neighbor_children)
218 27579 : if (!neighbor_set.contains(current_child) && keep_functor(current_child) &&
219 : current_child != elem)
220 : {
221 2272 : next_untested_set.insert(current_child);
222 2272 : neighbor_set.insert(current_child);
223 : }
224 : }
225 : #endif // #ifdef LIBMESH_ENABLE_AMR
226 : }
227 1029579 : untested_set.swap(next_untested_set);
228 : next_untested_set.clear();
229 : }
230 334448 : }
231 :
232 : template <typename T>
233 : bool
234 2897955 : findEdgeNeighborsWithinEdgeInternal(const Elem * const candidate,
235 : const Elem * const elem,
236 : const Node * const vertex1,
237 : const Node * const vertex2,
238 : const Real edge_length,
239 : std::vector<NeighborInfo> & info)
240 : {
241 : // If we have the type already: use it so we can avoid all of the virtual
242 : // calls that we would be making on Elem
243 : mooseAssert(dynamic_cast<const T *>(candidate), "Incorrect elem type");
244 2897955 : const T * const T_candidate = static_cast<const T *>(candidate);
245 :
246 : // Local index of our two vertices of interest (if any)
247 2897955 : auto v1 = T_candidate->get_node_index(vertex1);
248 2897955 : auto v2 = T_candidate->get_node_index(vertex2);
249 : // Whether or not we have said vertices
250 2897955 : const bool has_v1 = v1 != invalid_uint;
251 2897955 : const bool has_v2 = v2 != invalid_uint;
252 :
253 : // If the candidate has both vertices, it shares the complete edge
254 2897955 : if (has_v1 && has_v2)
255 : {
256 : // Add the sides that contain said edge
257 3903501 : for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
258 3903501 : if (T_candidate->is_node_on_edge(v1, e) && T_candidate->is_node_on_edge(v2, e))
259 : {
260 828829 : std::vector<unsigned short> sides(2);
261 828829 : sides[0] = T::edge_sides_map[e][0];
262 828829 : sides[1] = T::edge_sides_map[e][1];
263 828829 : info.emplace_back(candidate, std::move(sides), 0, 1);
264 : return true;
265 828829 : }
266 :
267 0 : mooseError("Failed to find a side that the vertices are on");
268 : }
269 :
270 : const auto n_vertices = T_candidate->n_vertices();
271 :
272 : // If we only have one of the vertices, we can still be contained within the edge if we have
273 : // another vertex that is within the edge
274 2069126 : if (has_v1 || has_v2)
275 : {
276 : // Local index of the vertex that the candidate and the target edge have in common
277 1636068 : const auto common_v = has_v1 ? v1 : v2;
278 :
279 : // See if another vertex that isn't the common node is contained
280 : MooseIndex(n_vertices) other_v;
281 10072356 : for (other_v = 0; other_v < n_vertices; ++other_v)
282 15239518 : if (other_v != common_v &&
283 6802026 : isWithinSegment(*vertex1, *vertex2, edge_length, T_candidate->point(other_v)))
284 : break;
285 :
286 : // If we have the common vertex and another vertex within the target edge, use the sides
287 : // that contain both of those vertices
288 1636068 : if (other_v != n_vertices)
289 : {
290 8140 : for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
291 8140 : if (T_candidate->is_node_on_edge(common_v, e) && T_candidate->is_node_on_edge(other_v, e))
292 : {
293 1204 : std::vector<unsigned short> sides(2);
294 1204 : sides[0] = T::edge_sides_map[e][0];
295 1204 : sides[1] = T::edge_sides_map[e][1];
296 :
297 1204 : info.emplace_back(
298 : candidate,
299 : std::move(sides),
300 1204 : has_v1 ? 0 : (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length,
301 1204 : has_v1 ? (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length : 1);
302 : return true;
303 1204 : }
304 0 : mooseError("Failed to find a side that the vertices are on");
305 : }
306 1634864 : else if (T_candidate->level() < elem->level())
307 : {
308 2961 : for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
309 3581 : if (T_candidate->is_node_on_edge(common_v, e) &&
310 770 : isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
311 770 : T_candidate->point(T::edge_nodes_map[e][1]),
312 770 : has_v1 ? *vertex2 : *vertex1))
313 : {
314 150 : other_v = T::edge_nodes_map[e][0] == common_v ? T::edge_nodes_map[e][1]
315 : : T::edge_nodes_map[e][0];
316 150 : std::vector<unsigned short> sides(2);
317 150 : sides[0] = T::edge_sides_map[e][0];
318 150 : sides[1] = T::edge_sides_map[e][1];
319 :
320 150 : info.emplace_back(
321 : candidate,
322 : std::move(sides),
323 150 : has_v1 ? 0 : (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length,
324 150 : has_v1 ? (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length : 1);
325 : return true;
326 150 : }
327 :
328 : return false;
329 : }
330 : }
331 : // When the candidate level is less, one of our vertices could be in a candidate's face and the
332 : // other could be within one of the candidate's edges
333 433058 : else if (T_candidate->level() < elem->level())
334 : {
335 : auto v1_edge = RayTracingCommon::invalid_edge;
336 : auto v2_edge = RayTracingCommon::invalid_edge;
337 :
338 : // See if any of the edges contain one of the vertices
339 57148 : for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
340 : {
341 105049 : if (v1_edge == RayTracingCommon::invalid_edge &&
342 52297 : isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
343 52297 : T_candidate->point(T::edge_nodes_map[e][1]),
344 : *vertex1))
345 82 : v1_edge = e;
346 102456 : if (v2_edge == RayTracingCommon::invalid_edge &&
347 49704 : isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
348 49704 : T_candidate->point(T::edge_nodes_map[e][1]),
349 : *vertex2))
350 590 : v2_edge = e;
351 : }
352 :
353 4396 : const auto v1_within = v1_edge != RayTracingCommon::invalid_edge;
354 4396 : const auto v2_within = v2_edge != RayTracingCommon::invalid_edge;
355 :
356 4396 : if (v1_within && v2_within)
357 : {
358 : mooseAssert(v1_edge == v2_edge, "Vertices should be contained in same edge");
359 :
360 0 : std::vector<unsigned short> sides(2);
361 0 : sides[0] = T::edge_sides_map[v1_edge][0];
362 0 : sides[1] = T::edge_sides_map[v1_edge][1];
363 0 : info.emplace_back(candidate, std::move(sides), 0, 1);
364 : return true;
365 0 : }
366 4396 : else if (v1_within || v2_within)
367 : {
368 672 : const auto in_edge = v1_within ? v1_edge : v2_edge;
369 672 : const Point & check_point = v1_within ? *vertex2 : *vertex1;
370 :
371 672 : std::vector<unsigned short> sides(1);
372 672 : Real lower_bound = 0;
373 672 : Real upper_bound = 1;
374 :
375 672 : if (candidate->side_ptr(T::edge_sides_map[in_edge][0])->contains_point(check_point))
376 9 : sides[0] = T::edge_sides_map[in_edge][0];
377 663 : else if (candidate->side_ptr(T::edge_sides_map[in_edge][1])->contains_point(check_point))
378 55 : sides[0] = T::edge_sides_map[in_edge][1];
379 : else
380 : {
381 608 : const Real point_bound = v1_within ? 0 : 1;
382 608 : lower_bound = point_bound;
383 608 : upper_bound = point_bound;
384 608 : sides[0] = T::edge_sides_map[in_edge][0];
385 608 : sides.push_back(T::edge_sides_map[in_edge][1]);
386 : }
387 :
388 672 : info.emplace_back(candidate, std::move(sides), lower_bound, upper_bound);
389 : return true;
390 672 : }
391 :
392 : return false;
393 : }
394 : // If we have neither of the vertices, and the candidate's edge can't fully contain [vertex1,
395 : // vertex2], check for other vertices that are contained within the edge
396 : else
397 : {
398 3247709 : for (v1 = 0; v1 < n_vertices; ++v1)
399 2819047 : if (isWithinSegment(*vertex1, *vertex2, edge_length, candidate->point(v1)))
400 : {
401 0 : for (v2 = v1 + 1; v2 < n_vertices; ++v2)
402 0 : if (isWithinSegment(*vertex1, *vertex2, edge_length, candidate->point(v2)))
403 : break;
404 : break;
405 : }
406 :
407 : // No vertices are contained within the edge
408 428662 : if (v1 == n_vertices)
409 : return false;
410 :
411 : // Only one vertex contained: add the sides associated with that vertex
412 0 : if (v2 >= n_vertices)
413 : {
414 : std::vector<unsigned short> sides;
415 0 : sides.reserve(2);
416 :
417 0 : for (MooseIndex(T::num_sides) s = 0; s < T::num_sides; ++s)
418 0 : if (T_candidate->is_node_on_side(v1, s))
419 0 : sides.push_back(s);
420 :
421 0 : if (sides.empty())
422 0 : mooseError("Failed to find sides that the vertex is on");
423 :
424 : // Bounds are the same because we only touch at a point
425 0 : const auto bound = (T_candidate->point(v1) - (Point)*vertex1).norm() / edge_length;
426 0 : info.emplace_back(candidate, std::move(sides), bound, bound);
427 : return true;
428 0 : }
429 : // Two vertices contained: find the sides that contain both vertices
430 : else
431 : {
432 0 : for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
433 0 : if (T_candidate->is_node_on_edge(v1, e) && T_candidate->is_node_on_edge(v2, e))
434 : {
435 0 : auto lower_bound = (T_candidate->point(v1) - (Point)*vertex1).norm() / edge_length;
436 0 : auto upper_bound = (T_candidate->point(v2) - (Point)*vertex1).norm() / edge_length;
437 0 : if (lower_bound > upper_bound)
438 : {
439 : const auto temp = lower_bound;
440 0 : lower_bound = upper_bound;
441 0 : upper_bound = temp;
442 : }
443 :
444 0 : std::vector<unsigned short> sides(2);
445 0 : sides[0] = T::edge_sides_map[e][0];
446 0 : sides[1] = T::edge_sides_map[e][1];
447 0 : info.emplace_back(candidate, std::move(sides), lower_bound, upper_bound);
448 : return true;
449 0 : }
450 :
451 0 : mooseError("Failed to find sides that the vertices are on");
452 : }
453 : }
454 :
455 : return false;
456 : }
457 :
458 : /**
459 : * Find the child of an elem that contains a point on a specified side of \p elem.
460 : * @param elem The parent element
461 : * @param point The point that is in elem
462 : * @param side The side that point is on in elem
463 : * @return The child that contains the point
464 : */
465 : const Elem *
466 : childContainingPointOnSide(const Elem * elem, const Point & point, const unsigned short side);
467 :
468 : /**
469 : * Get the active neighbor on side of elem that contains point.
470 : *
471 : * @param elem The element
472 : * @param side The side of elem that point is on
473 : * @param point The point on the side you want the neighbor for
474 : * @return The active neighbor that contains the point
475 : */
476 : const Elem * getActiveNeighbor(const Elem * elem, const unsigned short side, const Point & point);
477 :
478 : /**
479 : * Checks for the intersection of a ray and a triangular face
480 : *
481 : * Adapted from:
482 : * webserver2.tecgraf.puc-rio.br/~mgattass/cg/trbRR/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
483 : *
484 : * @param start Start point of the ray
485 : * @param direction Direction of the ray
486 : * @param elem The elem that contains the triangular face (used to get the vertex points)
487 : * @param v0 Vertex index on elem that is the zeroth vertex of the triangular face
488 : * @param v1 Vertex index on elem that is the first vertex of the triangular face
489 : * @param v2 Vertex index on elem that is the second vertex of the triangular face
490 : * @param intersection_distance If intersected, storage for the intersection distance
491 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
492 : * @return Whether or not the ray intersected the triangular face
493 : */
494 : bool intersectTriangle(const Point & start,
495 : const Point & direction,
496 : const Elem * const elem,
497 : const unsigned short v0,
498 : const unsigned short v1,
499 : const unsigned short v2,
500 : Real & intersection_distance,
501 : ElemExtrema & intersected_extrema,
502 : const Real hmax
503 : #ifdef DEBUG_RAY_INTERSECTIONS
504 : ,
505 : const bool debug
506 : #endif
507 : );
508 :
509 : /**
510 : * Checks for the intersection of a ray and a quadrilateral, numbered as such:
511 : *
512 : * v01 v11
513 : * o--------o
514 : * | |
515 : * | |
516 : * o--------o
517 : * v00 v10
518 : *
519 : * Uses intersectTriangle() to check the possible intersection with the two triangles that make
520 : * up the quad
521 : *
522 : * @param start Start point of the ray
523 : * @param direction Direction of the ray
524 : * @param elem The elem that contains the quad face (used to get the vertex points)
525 : * @param v00 Vertex index on elem that represents v00 in the method description
526 : * @param v10 Vertex index on elem that represents v10 in the method description
527 : * @param v11 Vertex index on elem that represents v11 in the method description
528 : * @param v01 Vertex index on elem that represents v01 in the method description
529 : * @param intersection_distance If intersected, storage for the intersection distance
530 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
531 : * @return Whether or not the ray intersected the quadrilateral
532 : */
533 : bool intersectQuad(const Point & start,
534 : const Point & direction,
535 : const Elem * const elem,
536 : const unsigned short v00,
537 : const unsigned short v10,
538 : const unsigned short v11,
539 : const unsigned short v01,
540 : Real & intersection_distance,
541 : ElemExtrema & intersected_extrema,
542 : const Real hmax
543 : #ifdef DEBUG_RAY_INTERSECTIONS
544 : ,
545 : const bool debug
546 : #endif
547 : );
548 :
549 : /**
550 : * SFINAE typed object for checking the intersection of a line segment with the face of a Face
551 : * typed element, i.e., the intersection with the edge of a 2D element.
552 : * @param elem The element to check
553 : * @param start_point Start point of the segment
554 : * @param direction Direction of the segment
555 : * @param side The side of elem to check
556 : * @param max_length A length that is longer than the length of the line segment
557 : * @param intersection_point If intersected, storage for the intersection point
558 : * @param intersection_distance If intersected, storage for the intersection distance
559 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
560 : * @return Whether or not the line segment intersected with the face
561 : */
562 : template <typename T>
563 : typename std::enable_if<std::is_base_of<libMesh::Face, T>::value, bool>::type
564 3084612 : sideIntersectedByLine(const Elem * elem,
565 : const Point & start_point,
566 : const Point & direction,
567 : const unsigned short side,
568 : const Real max_length,
569 : Point & intersection_point,
570 : Real & intersection_distance,
571 : ElemExtrema & intersected_extrema,
572 : const Real /* hmax */
573 : #ifdef DEBUG_RAY_INTERSECTIONS
574 : ,
575 : const bool debug
576 : #endif
577 : )
578 : {
579 : mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
580 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
581 :
582 3084612 : SegmentVertices segment_vertex = SEGMENT_VERTEX_NONE;
583 :
584 3084612 : const bool intersected = lineLineIntersect2D(start_point,
585 : direction,
586 : max_length,
587 3084612 : elem->point(T::side_nodes_map[side][0]),
588 3084612 : elem->point(T::side_nodes_map[side][1]),
589 : intersection_point,
590 : intersection_distance,
591 : segment_vertex
592 : #ifdef DEBUG_RAY_INTERSECTIONS
593 : ,
594 : debug
595 : #endif
596 : );
597 :
598 3084612 : if (segment_vertex != SEGMENT_VERTEX_NONE)
599 : {
600 657969 : intersected_extrema.setVertex(T::side_nodes_map[side][segment_vertex]);
601 : mooseAssert(
602 : intersected_extrema.vertexPoint(elem).absolute_fuzzy_equals(intersection_point,
603 : 3 * TRACE_TOLERANCE),
604 : "Doesn't intersect vertex at: " + Moose::stringify(intersected_extrema.vertexPoint(elem)) +
605 : " tentative intersection point: " + Moose::stringify(intersection_point));
606 : }
607 :
608 3084612 : return intersected;
609 : }
610 :
611 : /**
612 : * SFINAE typed object for checking the intersection of a line segment with the face of a Hex
613 : * typed element
614 : * @param elem The element to check
615 : * @param start_point Start point of the segment
616 : * @param direction Direction of the segment
617 : * @param side The side of elem to check
618 : * @param intersection_point If intersected, storage for the intersection point
619 : * @param intersection_distance If intersected, storage for the intersection distance
620 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
621 : * @return Whether or not the line segment intersected with the face
622 : */
623 : template <typename T>
624 : typename std::enable_if<std::is_base_of<Hex, T>::value, bool>::type
625 51905308 : sideIntersectedByLine(const Elem * elem,
626 : const Point & start_point,
627 : const Point & direction,
628 : const unsigned short side,
629 : const Real,
630 : Point & intersection_point,
631 : Real & intersection_distance,
632 : ElemExtrema & intersected_extrema,
633 : const Real hmax
634 : #ifdef DEBUG_RAY_INTERSECTIONS
635 : ,
636 : const bool debug
637 : #endif
638 : )
639 : {
640 : mooseAssert(elem->first_order_equivalent_type(elem->type()) == HEX8, "Not a Hex");
641 : mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
642 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
643 :
644 51905308 : const bool intersected = intersectQuad(start_point,
645 : direction,
646 : elem,
647 51905308 : T::side_nodes_map[side][3],
648 51905308 : T::side_nodes_map[side][2],
649 51905308 : T::side_nodes_map[side][1],
650 51905308 : T::side_nodes_map[side][0],
651 : intersection_distance,
652 : intersected_extrema,
653 : hmax
654 : #ifdef DEBUG_RAY_INTERSECTIONS
655 : ,
656 : debug
657 : #endif
658 : );
659 :
660 51905308 : if (intersected)
661 9316983 : intersection_point = start_point + intersection_distance * direction;
662 :
663 51905308 : return intersected;
664 : }
665 :
666 : /**
667 : * SFINAE typed object for checking the intersection of a line segment with the face of a Tet
668 : * typed element
669 : * @param elem The element to check
670 : * @param start_point Start point of the segment
671 : * @param direction Direction of the segment
672 : * @param side The side of elem to check
673 : * @param max_length A length that is longer than the length of the line segment
674 : * @param intersection_point If intersected, storage for the intersection point
675 : * @param intersection_distance If intersected, storage for the intersection distance
676 : * @param intersected_vertex If a vertex is intersected, storage for said intersection
677 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
678 : * @return Whether or not the line segment intersected with the face
679 : */
680 : template <typename T>
681 : typename std::enable_if<std::is_base_of<Tet, T>::value, bool>::type
682 71302236 : sideIntersectedByLine(const Elem * elem,
683 : const Point & start_point,
684 : const Point & direction,
685 : const unsigned short side,
686 : const Real /* max_length */,
687 : Point & intersection_point,
688 : Real & intersection_distance,
689 : ElemExtrema & intersected_extrema,
690 : const Real hmax
691 : #ifdef DEBUG_RAY_INTERSECTIONS
692 : ,
693 : const bool debug
694 : #endif
695 : )
696 : {
697 : mooseAssert(elem->first_order_equivalent_type(elem->type()) == TET4, "Not a Tet");
698 : mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
699 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
700 :
701 71302236 : const bool intersected = intersectTriangle(start_point,
702 : direction,
703 : elem,
704 71302236 : T::side_nodes_map[side][2],
705 71302236 : T::side_nodes_map[side][1],
706 71302236 : T::side_nodes_map[side][0],
707 : intersection_distance,
708 : intersected_extrema,
709 : hmax
710 : #ifdef DEBUG_RAY_INTERSECTIONS
711 : ,
712 : debug
713 : #endif
714 : );
715 :
716 71302236 : if (intersected)
717 20479564 : intersection_point = start_point + intersection_distance * direction;
718 :
719 71302236 : return intersected;
720 : }
721 :
722 : /**
723 : * SFINAE typed object for checking the intersection of a line segment with the face of a Pyramid
724 : * typed element
725 : * @param elem The element to check
726 : * @param start_point Start point of the segment
727 : * @param direction Direction of the segment
728 : * @param side The side of elem to check
729 : * @param max_length A length that is longer than the length of the line segment
730 : * @param intersection_point If intersected, storage for the intersection point
731 : * @param intersection_distance If intersected, storage for the intersection distance
732 : * @param intersected_vertex If a vertex is intersected, storage for said intersection
733 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
734 : * @return Whether or not the line segment intersected with the face
735 : */
736 : template <typename T>
737 : typename std::enable_if<std::is_base_of<libMesh::Pyramid, T>::value, bool>::type
738 30757692 : sideIntersectedByLine(const Elem * elem,
739 : const Point & start_point,
740 : const Point & direction,
741 : const unsigned short side,
742 : const Real /* max_length */,
743 : Point & intersection_point,
744 : Real & intersection_distance,
745 : ElemExtrema & intersected_extrema,
746 : const Real hmax
747 : #ifdef DEBUG_RAY_INTERSECTIONS
748 : ,
749 : const bool debug
750 : #endif
751 : )
752 : {
753 : mooseAssert(elem->first_order_equivalent_type(elem->type()) == PYRAMID5, "Not a Pyramid");
754 : mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
755 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
756 :
757 30757692 : const bool intersected = side < 4 ? intersectTriangle(start_point,
758 : direction,
759 : elem,
760 25552944 : T::side_nodes_map[side][2],
761 25552944 : T::side_nodes_map[side][1],
762 25552944 : T::side_nodes_map[side][0],
763 : intersection_distance,
764 : intersected_extrema,
765 : hmax
766 : #ifdef DEBUG_RAY_INTERSECTIONS
767 : ,
768 : debug
769 : #endif
770 : )
771 5204748 : : intersectQuad(start_point,
772 : direction,
773 : elem,
774 5204748 : T::side_nodes_map[side][3],
775 5204748 : T::side_nodes_map[side][2],
776 5204748 : T::side_nodes_map[side][1],
777 5204748 : T::side_nodes_map[side][0],
778 : intersection_distance,
779 : intersected_extrema,
780 : hmax
781 : #ifdef DEBUG_RAY_INTERSECTIONS
782 : ,
783 : debug
784 : #endif
785 : );
786 :
787 30757692 : if (intersected)
788 6452758 : intersection_point = start_point + intersection_distance * direction;
789 :
790 30757692 : return intersected;
791 : }
792 :
793 : /**
794 : * SFINAE typed object for checking the intersection of a line segment with the face of a Prism
795 : * typed element
796 : * @param elem The element to check
797 : * @param start_point Start point of the segment
798 : * @param direction Direction of the segment
799 : * @param side The side of elem to check
800 : * @param max_length A length that is longer than the length of the line segment
801 : * @param intersection_point If intersected, storage for the intersection point
802 : * @param intersection_distance If intersected, storage for the intersection distance
803 : * @param intersected_vertex If a vertex is intersected, storage for said intersection
804 : * @param intersected_extrema If a vertex/edge is intersected, will be filled with the intersection
805 : * @return Whether or not the line segment intersected with the face
806 : */
807 : template <typename T>
808 : typename std::enable_if<std::is_base_of<libMesh::Prism, T>::value, bool>::type
809 23412576 : sideIntersectedByLine(const Elem * elem,
810 : const Point & start_point,
811 : const Point & direction,
812 : const unsigned short side,
813 : const Real /* max_length */,
814 : Point & intersection_point,
815 : Real & intersection_distance,
816 : ElemExtrema & intersected_extrema,
817 : const Real hmax
818 : #ifdef DEBUG_RAY_INTERSECTIONS
819 : ,
820 : const bool debug
821 : #endif
822 : )
823 : {
824 : mooseAssert(elem->first_order_equivalent_type(elem->type()) == PRISM6, "Not a Prism");
825 : mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
826 : mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");
827 :
828 23412576 : const bool intersected = (side == 0 || side == 4) ? intersectTriangle(start_point,
829 : direction,
830 : elem,
831 10044637 : T::side_nodes_map[side][2],
832 10044637 : T::side_nodes_map[side][1],
833 10044637 : T::side_nodes_map[side][0],
834 : intersection_distance,
835 : intersected_extrema,
836 : hmax
837 : #ifdef DEBUG_RAY_INTERSECTIONS
838 : ,
839 : debug
840 : #endif
841 : )
842 13367939 : : intersectQuad(start_point,
843 : direction,
844 : elem,
845 13367939 : T::side_nodes_map[side][3],
846 13367939 : T::side_nodes_map[side][2],
847 13367939 : T::side_nodes_map[side][1],
848 13367939 : T::side_nodes_map[side][0],
849 : intersection_distance,
850 : intersected_extrema,
851 : hmax
852 : #ifdef DEBUG_RAY_INTERSECTIONS
853 : ,
854 : debug
855 : #endif
856 : );
857 :
858 23412576 : if (intersected)
859 4474773 : intersection_point = start_point + intersection_distance * direction;
860 :
861 23412576 : return intersected;
862 : }
863 :
864 : /**
865 : * @return Whether not the element is traceable
866 : */
867 : bool isTraceableElem(const Elem * elem);
868 : /**
869 : * @return Whether or not the element is traceable with adaptivity
870 : */
871 : bool isAdaptivityTraceableElem(const Elem * elem);
872 :
873 : /**
874 : * Determines if a point is at a vertex of an element.
875 : * @param elem The element
876 : * @param point The point
877 : * @return The local vertex ID the point is at, if any (invalid_point otherwise)
878 : */
879 : unsigned short atVertex(const Elem * elem, const Point & point);
880 :
881 : /**
882 : * Determines if a point is at a vertex on the side of en element.
883 : * @param elem The element
884 : * @param point The point
885 : * @param side The side
886 : * @return The local vertex ID the point is at, if any (invalid_point otherwise)
887 : */
888 : unsigned short atVertexOnSide(const Elem * elem, const Point & point, const unsigned short side);
889 :
890 : /**
891 : * Returns the number of nodes on a side for an Elem that is not a Pyramid or Prism.
892 : */
893 : template <typename T>
894 : inline typename std::enable_if<!std::is_base_of<libMesh::Pyramid, T>::value &&
895 : !std::is_base_of<libMesh::Prism, T>::value,
896 : unsigned short>::type
897 : nodesPerSide(const unsigned short)
898 : {
899 : return T::nodes_per_side;
900 : }
901 :
902 : /**
903 : * Returns the number of nodes on a side on a Pyramid elem.
904 : */
905 : template <typename T>
906 : inline typename std::enable_if<std::is_base_of<Pyramid, T>::value, unsigned short>::type
907 : nodesPerSide(const unsigned short side)
908 : {
909 1072083 : return T::nodes_per_side - (side != 4);
910 : }
911 :
912 : /**
913 : * Returns the number of nodes on a side on a Prism elem.
914 : */
915 : template <typename T>
916 : inline typename std::enable_if<std::is_base_of<libMesh::Prism, T>::value, unsigned short>::type
917 : nodesPerSide(const unsigned short side)
918 : {
919 2139448 : return T::nodes_per_side - (side == 0 || side == 4);
920 : }
921 :
922 : /**
923 : * Determines if a point is at a vertex on the side of an element.
924 : *
925 : * SFINAE here makes this method available for only 2D/3D elem types (not edges)
926 : *
927 : * The template argument should be the first order type of the elem - this is used to
928 : * access the vertex mappings for said element without calling any virtuals
929 : *
930 : * @param elem The element
931 : * @param point The point
932 : * @param side The side
933 : * @return The local vertex ID the point is at, if any (invalid_point otherwise)
934 : */
935 : template <typename T>
936 : typename std::enable_if<!std::is_base_of<Edge, T>::value, unsigned short>::type
937 : atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side);
938 :
939 : /**
940 : * Determines if a point is at a vertex on the side of an element.
941 : *
942 : * SFINAE here makes this method available for only 1D elem types (edges)
943 : *
944 : * The template argument should be the first order type of the elem - this is used to
945 : * access the vertex mappings for said element without calling any virtuals
946 : *
947 : * @param elem The element
948 : * @param point The point
949 : * @param side The side
950 : * @return The local vertex ID the point is at, if any (invalid_point otherwise)
951 : */
952 : template <typename T>
953 : typename std::enable_if<std::is_base_of<Edge, T>::value, unsigned short>::type
954 : atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side);
955 :
956 : /**
957 : * Determines if a point is within edge on an element.
958 : *
959 : * The template argument should be a derived element type that corresponds to the elem->type(),
960 : * and is used to obtain edge/node mapping without using virtuals.
961 : *
962 : * @param elem The element
963 : * @param point The point
964 : * @param extrema To be filled with the local vertex IDs that contain the edge the point is
965 : * within, if any
966 : * @param tolerance The tolerance to use
967 : * @return If the point is within an edge of the element
968 : */
969 : template <typename T>
970 : bool withinEdgeTempl(const Elem * elem,
971 : const Point & point,
972 : ElemExtrema & extrema,
973 : const Real tolerance = TRACE_TOLERANCE);
974 :
975 : /**
976 : * Determines if a point is within an edge on an element.
977 : * @param elem The element
978 : * @param point The point
979 : * @param extrema To be filled with the local vertex IDs that contain the edge the point is
980 : * within, if any
981 : * @param tolerance The tolerance to use
982 : * @return If the point is within an edge of the element
983 : */
984 : bool withinEdge(const Elem * elem,
985 : const Point & point,
986 : ElemExtrema & extrema,
987 : const Real tolerance = TRACE_TOLERANCE);
988 :
989 : /**
990 : * Determines if a point is within an edge on the side of an element.
991 : * @param elem The element
992 : * @param point The point
993 : * @param side The side
994 : * @param extrema To be filled with the local vertex IDs that contain the edge the point is
995 : * within, if any
996 : * @return If the point is within an edge on the side of the element
997 : */
998 : bool withinEdgeOnSide(const Elem * const elem,
999 : const Point & point,
1000 : const unsigned short side,
1001 : ElemExtrema & extrema);
1002 :
1003 : /**
1004 : * Determines if a point is within an Elem's extrema (at vertex/within edge) on a side
1005 : * @param elem The element
1006 : * @param point The point
1007 : * @param side The side
1008 : * @param dim The element dimension
1009 : * @param extrema To be filled with the extrema if any are found
1010 : * @return If the point is at a vertex/within an edge on the side
1011 : */
1012 : bool withinExtremaOnSide(const Elem * const elem,
1013 : const Point & point,
1014 : const unsigned short side,
1015 : const unsigned int dim,
1016 : ElemExtrema & extrema);
1017 :
1018 : /**
1019 : * Determines if a point is within an edge on the side of an element.
1020 : *
1021 : * SFINAE here makes this method available for only 3D elem types (Cell)
1022 : *
1023 : * The template argument should be the first order type of the elem - this is used to
1024 : * access the vertex mappings for said element without calling any virtuals
1025 : *
1026 : * @param elem The element
1027 : * @param point The point
1028 : * @param side The side
1029 : * @param extrema To be filled with the local vertex IDs that contain the edge the point is
1030 : * within, if any
1031 : * @return If the point is within an edge on the side of the element
1032 : */
1033 : template <typename T>
1034 : typename std::enable_if<std::is_base_of<libMesh::Cell, T>::value, bool>::type withinEdgeOnSideTempl(
1035 : const Elem * const elem, const Point & point, const unsigned short side, ElemExtrema & extrema);
1036 :
1037 : /**
1038 : * Determines if a point is within an edge on the side of an element.
1039 : *
1040 : * SFINAE here makes this method available for only 1D/2D elem types (not Cell), which do not have
1041 : * edges, therefore this function errors.
1042 : */
1043 : template <typename T>
1044 : typename std::enable_if<!std::is_base_of<libMesh::Cell, T>::value, bool>::type
1045 : withinEdgeOnSideTempl(const Elem * const, const Point &, const unsigned short, ElemExtrema &)
1046 : {
1047 : mooseError("Should not call withinEdgeOnSideTempl() with a non-Cell derived Elem");
1048 : }
1049 :
1050 : /**
1051 : * Whether or not \p point is on the boundary (min/max) of \p bbox.
1052 : *
1053 : * Checks \p dim dimensions.
1054 : */
1055 : bool onBoundingBoxBoundary(const BoundingBox & bbox,
1056 : const Point & point,
1057 : const unsigned int dim,
1058 : const Real tolerance);
1059 :
1060 : }
|
