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 "TraceRay.h"
11 :
12 : // Local includes
13 : #include "Ray.h"
14 : #include "RayBoundaryConditionBase.h"
15 : #include "RayKernelBase.h"
16 : #include "RayTracingStudy.h"
17 : #include "TraceRayTools.h"
18 :
19 : // libMesh includes
20 : #include "libmesh/cell_tet4.h"
21 : #include "libmesh/cell_tet10.h"
22 : #include "libmesh/cell_tet14.h"
23 : #include "libmesh/cell_hex8.h"
24 : #include "libmesh/cell_hex20.h"
25 : #include "libmesh/cell_hex27.h"
26 : #include "libmesh/cell_prism6.h"
27 : #include "libmesh/cell_prism15.h"
28 : #include "libmesh/cell_prism18.h"
29 : #include "libmesh/cell_pyramid5.h"
30 : #include "libmesh/cell_pyramid13.h"
31 : #include "libmesh/cell_pyramid14.h"
32 : #include "libmesh/edge_edge2.h"
33 : #include "libmesh/edge_edge3.h"
34 : #include "libmesh/edge_edge4.h"
35 : #include "libmesh/face_quad4.h"
36 : #include "libmesh/face_quad8.h"
37 : #include "libmesh/face_quad9.h"
38 : #include "libmesh/face_tri3.h"
39 : #include "libmesh/face_tri6.h"
40 : #include "libmesh/face_tri7.h"
41 : #include "libmesh/enum_to_string.h"
42 : #include "libmesh/mesh.h"
43 :
44 : using namespace libMesh;
45 : using namespace TraceRayTools;
46 :
47 2366 : TraceRay::TraceRay(RayTracingStudy & study, const THREAD_ID tid)
48 2366 : : _study(study),
49 2366 : _mesh(study.getSubProblem().mesh()),
50 2366 : _dim(_mesh.dimension()),
51 2366 : _boundary_info(_mesh.getMesh().get_boundary_info()),
52 2366 : _pid(_study.comm().rank()),
53 2366 : _tid(tid),
54 2366 : _backface_culling(false),
55 2366 : _current_normals(nullptr),
56 4732 : _results(ENDED_STATIONARY + 1)
57 : {
58 2366 : }
59 :
60 : void
61 6201 : TraceRay::preExecute()
62 : {
63 6201 : _current_subdomain_id = Elem::invalid_subdomain_id;
64 6201 : _current_elem_type = INVALID_ELEM;
65 :
66 : // Zero out all results
67 99216 : for (auto & val : _results)
68 93015 : val = 0;
69 :
70 6201 : _has_ray_kernels = _study.hasRayKernels(_tid);
71 6201 : _is_rectangular_domain = _study.isRectangularDomain();
72 6201 : }
73 :
74 : void
75 150 : TraceRay::meshChanged()
76 : {
77 : // Invalidate the vertex and edge neighbor caches
78 : _vertex_neighbors.clear();
79 : _edge_neighbors.clear();
80 150 : }
81 :
82 : TraceRay::ExitsElemResult
83 31429648 : TraceRay::exitsElem(const Elem * elem,
84 : const ElemType elem_type,
85 : const unsigned short incoming_side,
86 : Point & intersection_point,
87 : unsigned short & intersected_side,
88 : ElemExtrema & intersected_extrema,
89 : Real & intersection_distance,
90 : const Point * normals)
91 : {
92 : debugRay("Called exitsElem()");
93 :
94 : traceAssert(elem_type == elem->type(), "elem_type incorrect");
95 : traceAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
96 : traceAssert(intersected_side == RayTracingCommon::invalid_side, "Side should be invalid");
97 : traceAssert(intersected_extrema.isInvalid(), "Extrema should be invalid");
98 : traceAssert(intersection_distance == RayTracingCommon::invalid_distance,
99 : "Distance should be invalid");
100 : if (_study.verifyRays() && incoming_side != RayTracingCommon::invalid_side &&
101 : !_study.sideIsNonPlanar(elem, incoming_side))
102 : traceAssert(_study.sideIsIncoming(elem, incoming_side, (*_current_ray)->direction(), _tid),
103 : "Incoming side is non-entrant");
104 :
105 : bool intersected;
106 31429648 : switch (elem_type)
107 : {
108 3894210 : case HEX8:
109 3894210 : intersected = exitsElem<Hex8, Hex8>(elem,
110 : incoming_side,
111 : intersection_point,
112 : intersected_side,
113 : intersected_extrema,
114 : intersection_distance,
115 : normals);
116 3894210 : break;
117 5182938 : case TET4:
118 5182938 : intersected = exitsElem<Tet4, Tet4>(elem,
119 : incoming_side,
120 : intersection_point,
121 : intersected_side,
122 : intersected_extrema,
123 : intersection_distance,
124 : normals);
125 5182938 : break;
126 1685591 : case PYRAMID5:
127 1685591 : intersected = exitsElem<Pyramid5, Pyramid5>(elem,
128 : incoming_side,
129 : intersection_point,
130 : intersected_side,
131 : intersected_extrema,
132 : intersection_distance,
133 : normals);
134 1685591 : break;
135 1276453 : case PRISM6:
136 1276453 : intersected = exitsElem<Prism6, Prism6>(elem,
137 : incoming_side,
138 : intersection_point,
139 : intersected_side,
140 : intersected_extrema,
141 : intersection_distance,
142 : normals);
143 1276453 : break;
144 588371 : case QUAD4:
145 588371 : intersected = exitsElem<Quad4, Quad4>(elem,
146 : incoming_side,
147 : intersection_point,
148 : intersected_side,
149 : intersected_extrema,
150 : intersection_distance,
151 : normals);
152 588371 : break;
153 36011 : case TRI3:
154 36011 : intersected = exitsElem<Tri3, Tri3>(elem,
155 : incoming_side,
156 : intersection_point,
157 : intersected_side,
158 : intersected_extrema,
159 : intersection_distance,
160 : normals);
161 36011 : break;
162 1191279 : case HEX20:
163 1191279 : intersected = exitsElem<Hex20, Hex8>(elem,
164 : incoming_side,
165 : intersection_point,
166 : intersected_side,
167 : intersected_extrema,
168 : intersection_distance,
169 : normals);
170 1191279 : break;
171 1191279 : case HEX27:
172 1191279 : intersected = exitsElem<Hex27, Hex8>(elem,
173 : incoming_side,
174 : intersection_point,
175 : intersected_side,
176 : intersected_extrema,
177 : intersection_distance,
178 : normals);
179 1191279 : break;
180 19460 : case QUAD8:
181 19460 : intersected = exitsElem<Quad8, Quad4>(elem,
182 : incoming_side,
183 : intersection_point,
184 : intersected_side,
185 : intersected_extrema,
186 : intersection_distance,
187 : normals);
188 19460 : break;
189 19460 : case QUAD9:
190 19460 : intersected = exitsElem<Quad9, Quad4>(elem,
191 : incoming_side,
192 : intersection_point,
193 : intersected_side,
194 : intersected_extrema,
195 : intersection_distance,
196 : normals);
197 19460 : break;
198 29506 : case TRI6:
199 29506 : intersected = exitsElem<Tri6, Tri3>(elem,
200 : incoming_side,
201 : intersection_point,
202 : intersected_side,
203 : intersected_extrema,
204 : intersection_distance,
205 : normals);
206 29506 : break;
207 29506 : case TRI7:
208 29506 : intersected = exitsElem<Tri7, Tri3>(elem,
209 : incoming_side,
210 : intersection_point,
211 : intersected_side,
212 : intersected_extrema,
213 : intersection_distance,
214 : normals);
215 29506 : break;
216 5177636 : case TET10:
217 5177636 : intersected = exitsElem<Tet10, Tet4>(elem,
218 : incoming_side,
219 : intersection_point,
220 : intersected_side,
221 : intersected_extrema,
222 : intersection_distance,
223 : normals);
224 5177636 : break;
225 5177636 : case TET14:
226 5177636 : intersected = exitsElem<Tet14, Tet4>(elem,
227 : incoming_side,
228 : intersection_point,
229 : intersected_side,
230 : intersected_extrema,
231 : intersection_distance,
232 : normals);
233 5177636 : break;
234 1685516 : case PYRAMID13:
235 1685516 : intersected = exitsElem<Pyramid13, Pyramid5>(elem,
236 : incoming_side,
237 : intersection_point,
238 : intersected_side,
239 : intersected_extrema,
240 : intersection_distance,
241 : normals);
242 1685516 : break;
243 1685516 : case PYRAMID14:
244 1685516 : intersected = exitsElem<Pyramid14, Pyramid5>(elem,
245 : incoming_side,
246 : intersection_point,
247 : intersected_side,
248 : intersected_extrema,
249 : intersection_distance,
250 : normals);
251 1685516 : break;
252 1276453 : case PRISM15:
253 1276453 : intersected = exitsElem<Prism15, Prism6>(elem,
254 : incoming_side,
255 : intersection_point,
256 : intersected_side,
257 : intersected_extrema,
258 : intersection_distance,
259 : normals);
260 1276453 : break;
261 1276453 : case PRISM18:
262 1276453 : intersected = exitsElem<Prism18, Prism6>(elem,
263 : incoming_side,
264 : intersection_point,
265 : intersected_side,
266 : intersected_extrema,
267 : intersection_distance,
268 : normals);
269 1276453 : break;
270 5966 : case EDGE2:
271 5966 : intersected = exitsElem<Edge2, Edge2>(elem,
272 : incoming_side,
273 : intersection_point,
274 : intersected_side,
275 : intersected_extrema,
276 : intersection_distance,
277 : normals);
278 5966 : break;
279 204 : case EDGE3:
280 204 : intersected = exitsElem<Edge3, Edge2>(elem,
281 : incoming_side,
282 : intersection_point,
283 : intersected_side,
284 : intersected_extrema,
285 : intersection_distance,
286 : normals);
287 204 : break;
288 204 : case EDGE4:
289 204 : intersected = exitsElem<Edge4, Edge2>(elem,
290 : incoming_side,
291 : intersection_point,
292 : intersected_side,
293 : intersected_extrema,
294 : intersection_distance,
295 : normals);
296 204 : break;
297 0 : default:
298 0 : mooseError(
299 0 : "Element type ", Utility::enum_to_string(elem->type()), " not supported by TraceRay");
300 : }
301 :
302 31429648 : if (intersected)
303 : {
304 24902341 : if (intersected_extrema.atExtrema())
305 : return intersected_extrema.atVertex() ? HIT_VERTEX : HIT_EDGE;
306 : return HIT_FACE;
307 : }
308 :
309 : return NO_EXIT;
310 : }
311 :
312 : template <typename T, typename FirstOrderT>
313 : typename std::enable_if<!std::is_base_of<Edge, T>::value, bool>::type
314 31423274 : TraceRay::exitsElem(const Elem * elem,
315 : const unsigned short incoming_side,
316 : Point & intersection_point,
317 : unsigned short & intersected_side,
318 : ElemExtrema & intersected_extrema,
319 : Real & intersection_distance,
320 : const Point * normals)
321 : {
322 31423274 : ++_results[INTERSECTION_CALLS];
323 :
324 : debugRay("Called exitsElem() in 2D or 3D");
325 :
326 31423274 : const auto hmax = subdomainHmax(elem);
327 :
328 : // Current point and distance (maybe not the best!)
329 : Point current_intersection_point;
330 : Real current_intersection_distance;
331 : ElemExtrema current_intersected_extrema;
332 : // Scale the minimum intersection distance required by hmax
333 31423274 : Real best_intersection_distance = TRACE_TOLERANCE * hmax;
334 : // Whether or not we are going to try backface culling the first time around
335 : // This depends on if the user set it to use culling
336 31423274 : bool use_backface_culling = _backface_culling;
337 : // If all sides have failed to find an intersection, whether or not we
338 : // are going to see if the nonplanar skip side (only if it IS nonplanar)
339 : // is also an exiting point
340 : bool try_nonplanar_incoming_side = false;
341 : // Easy access into the current direction
342 31423274 : const auto & direction = (*_current_ray)->direction();
343 : // The current side that we're checking
344 : unsigned short s = 0;
345 :
346 : // Loop over the side loop. We need this in the cases that:
347 : // - Backface culling is enabled and we were not able to find an intersection
348 : // so we will go through all culled sides and see if we can find one
349 : // - The incoming_side is non-planar and the ray also exits said side,
350 : // which will be checked after all other sides fail
351 : while (true)
352 : {
353 : debugRay(" use_backface_culling = ", use_backface_culling);
354 : debugRay(" try_nonplanar_incoming_side = ", try_nonplanar_incoming_side);
355 : // Loop over all of the sides
356 : while (true)
357 : {
358 : debugRay(" Side ", s, " with centroid ", _elem_side_builder(*elem, s).vertex_average());
359 :
360 : // All of the checks that follow are done while we're still searching through
361 : // all of the sides. try_nonplanar_incoming_side is our last possible check
362 : // and does not involve looping through all of the sides, so skip these
363 : // checks if we're at that point.
364 149031228 : if (!try_nonplanar_incoming_side)
365 : {
366 : // Don't search backwards. If we have a non-planar incoming side, we will
367 : // check it if all other sides have failed
368 148960087 : if (s == incoming_side)
369 : {
370 : debugRay(" Skipping due to incoming side");
371 29786003 : if (++s == T::num_sides)
372 : break;
373 : else
374 21993910 : continue;
375 : }
376 :
377 : // Using backface culling on this run through the sides
378 : // If the direction is non-entrant, skip this side
379 119174084 : if (use_backface_culling)
380 : {
381 : // Side is non-entrant per the culling, so skip
382 6110997 : if (normals[s] * direction < -LOOSE_TRACE_TOLERANCE)
383 : {
384 : debugRay(" Skipping due to backface culling dot = ", normals[s] * direction);
385 :
386 1972071 : if (++s == T::num_sides)
387 : break;
388 : else
389 1606975 : continue;
390 :
391 : ++_results[BACKFACE_CULLING_SUCCESSES];
392 : }
393 : }
394 : // We're not using backface culling but it is enabled, which means
395 : // we're on our second run through the sides because we could not
396 : // find any intersections while culling the sides. Try again and
397 : // check the sides that we previously skipped for intersections
398 113063087 : else if (_backface_culling)
399 : {
400 : // Side was non-entrant per the culling, so try again
401 1601280 : if (normals[s] * direction >= -LOOSE_TRACE_TOLERANCE)
402 : {
403 : debugRay(" Skipping because we already checked this side with culling enabled");
404 1045344 : if (++s == T::num_sides)
405 : break;
406 : else
407 855514 : continue;
408 : }
409 : else
410 : debugRay(" Side that was skipped due to culling");
411 :
412 555936 : ++_results[BACKFACE_CULLING_FAILURES];
413 : }
414 : }
415 :
416 : // Look for an intersection!
417 116227810 : current_intersection_point = RayTracingCommon::invalid_point;
418 : current_intersected_extrema.invalidate();
419 116227810 : const bool intersected = sideIntersectedByLine<FirstOrderT>(elem,
420 116227810 : _incoming_point,
421 : direction,
422 : s,
423 116227810 : _study.domainMaxLength(),
424 : current_intersection_point,
425 : current_intersection_distance,
426 : current_intersected_extrema,
427 : hmax
428 : #ifdef DEBUG_RAY_INTERSECTIONS
429 : ,
430 : DEBUG_RAY_IF
431 : #endif
432 : );
433 :
434 : // Do they intersect and is it further down the path than any other intersection?
435 : // If so, keep track of the intersection with the furthest distance
436 116227810 : if (intersected)
437 : {
438 : debugRay(" Intersected at point ",
439 : current_intersection_point,
440 : " with distance ",
441 : current_intersection_distance);
442 : debugRay(" Best intersection distance = ", best_intersection_distance);
443 :
444 : #ifndef NDEBUG
445 : // Only validate intersections if the side is planar
446 : if (_study.verifyTraceIntersections() && !_study.sideIsNonPlanar(elem, s) &&
447 : !_elem_side_builder(*elem, s).contains_point(current_intersection_point))
448 : failTrace("Intersected side does not contain intersection point",
449 : _study.tolerateFailure(),
450 : __LINE__);
451 : #endif
452 :
453 : // The intersection we just computed is further than any other intersection
454 : // that was found so far - mark it as the best
455 27316037 : if (current_intersection_distance > best_intersection_distance)
456 : {
457 : debugRay(" Best intersection so far");
458 :
459 25352837 : intersected_side = s;
460 25352837 : intersection_distance = current_intersection_distance;
461 25352837 : intersection_point = current_intersection_point;
462 : intersected_extrema = current_intersected_extrema;
463 : best_intersection_distance = current_intersection_distance;
464 : }
465 : }
466 :
467 116227810 : if (++s == T::num_sides || try_nonplanar_incoming_side)
468 : break;
469 : else
470 92759198 : continue;
471 : } // while(true)
472 :
473 : // Found an intersection
474 31815631 : if (intersected_side != RayTracingCommon::invalid_side)
475 : {
476 : debugRay(" Exiting with intersection");
477 : debugRay(" intersected_side = ", intersected_side);
478 : debugRay(" intersection_distance = ", intersection_distance);
479 : debugRay(" intersection_point = ", intersection_point);
480 : debugRay(" intersected_extrema = ", intersected_extrema);
481 :
482 : return true;
483 : }
484 :
485 : // This was our last possible check - no dice!
486 6919652 : if (try_nonplanar_incoming_side)
487 : return false;
488 :
489 : // We didn't find an intersection but we used backface culling. Try again without
490 : // it, checking only the sides that we skipped due to culling
491 6848511 : if (use_backface_culling)
492 : {
493 : debugRay(" Didn't find an intersection, retrying without backface culling");
494 : use_backface_culling = false;
495 : s = 0;
496 321216 : continue;
497 : }
498 :
499 : // Have tried all sides (potentially with and without culling). If the incoming
500 : // side is valid is non planar, see if we also exit out of it
501 6527295 : if (_incoming_side != RayTracingCommon::invalid_side &&
502 6524414 : _study.sideIsNonPlanar(elem, incoming_side))
503 : {
504 : debugRay(" Didn't find an intersection, trying non-planar incoming_side");
505 : try_nonplanar_incoming_side = true;
506 : s = incoming_side;
507 71141 : continue;
508 : }
509 :
510 : // No intersection found!
511 : return false;
512 : } // while(true)
513 : }
514 :
515 : template <typename T, typename FirstOrderT>
516 : typename std::enable_if<std::is_base_of<Edge, T>::value, bool>::type
517 6374 : TraceRay::exitsElem(const Elem * elem,
518 : const unsigned short incoming_side,
519 : Point & intersection_point,
520 : unsigned short & intersected_side,
521 : ElemExtrema & intersected_extrema,
522 : Real & intersection_distance,
523 : const Point *)
524 : {
525 6374 : ++_results[INTERSECTION_CALLS];
526 :
527 : debugRay("Called exitsElem() in 1D");
528 :
529 : // Scale the tolerance based on the element size
530 6374 : const auto tol = subdomainHmax(elem) * TRACE_TOLERANCE;
531 :
532 : // Can quickly return if we have an incoming side
533 : // There's only one other choice in 1D!
534 6374 : if (incoming_side != RayTracingCommon::invalid_side)
535 : {
536 6054 : intersected_side = (incoming_side == 1 ? 0 : 1);
537 : intersected_extrema.setVertex(intersected_side);
538 6054 : intersection_point = elem->point(intersected_side);
539 6054 : intersection_distance = (_incoming_point - intersection_point).norm();
540 :
541 : debugRay(" Incoming side is set to ", incoming_side, " so setting to other side");
542 : debugRay(" Intersected side ", intersected_side, " at ", intersection_point);
543 :
544 6054 : return true;
545 : }
546 :
547 : // End point that is for sure out of the element
548 : const Point extended_end_point =
549 320 : _incoming_point + _study.domainMaxLength() * (*_current_ray)->direction();
550 :
551 : // We're looking for the side whose point lays between the incoming point and the
552 : // extended end point
553 526 : for (MooseIndex(elem->n_sides()) side = 0; side < elem->n_sides(); ++side)
554 : {
555 514 : const Point side_point = elem->point(side);
556 : debugRay(" Checking side ", side, " at ", side_point);
557 :
558 514 : const Real incoming_to_side = (side_point - _incoming_point).norm();
559 514 : if (incoming_to_side < tol)
560 : {
561 : debugRay(" Continuing because at side");
562 134 : continue;
563 : }
564 :
565 380 : const Real incoming_to_end = (extended_end_point - _incoming_point).norm();
566 380 : const Real side_to_end = (extended_end_point - side_point).norm();
567 380 : const Real sum = incoming_to_side + side_to_end - incoming_to_end;
568 : debugRay(" Sum = ", sum);
569 :
570 380 : if (std::abs(sum) < tol)
571 : {
572 308 : intersected_side = side;
573 : intersected_extrema.setVertex(side);
574 308 : intersection_point = side_point;
575 308 : intersection_distance = incoming_to_side;
576 : debugRay(" Intersected at ", intersection_point);
577 308 : return true;
578 : }
579 : }
580 :
581 : return false;
582 : }
583 :
584 : TraceRay::ExitsElemResult
585 2586642 : TraceRay::moveThroughNeighbors(const std::vector<NeighborInfo> & neighbors,
586 : const Elem * last_elem,
587 : const Elem *& best_elem,
588 : unsigned short & best_elem_incoming_side)
589 : {
590 2586642 : ++_results[MOVED_THROUGH_NEIGHBORS];
591 :
592 : debugRay("Called moveThroughNeighbors() with ", neighbors.size(), " neighbors, and:");
593 : debugRay(" last_elem->id() = ", last_elem ? last_elem->id() : DofObject::invalid_id);
594 : debugRay(" _incoming_point = ", _incoming_point);
595 :
596 : traceAssert(!best_elem, "Best elem should be null");
597 : traceAssert(best_elem_incoming_side == RayTracingCommon::invalid_side,
598 : "Best elem side should be invalid");
599 : traceAssert(_intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
600 : traceAssert(_intersected_side == RayTracingCommon::invalid_side, "Side should be invalid");
601 : traceAssert(_intersected_extrema.isInvalid(), "Extrema should be invalid");
602 : traceAssert(_intersection_distance == RayTracingCommon::invalid_distance,
603 : "Distance should be invalid");
604 :
605 : // Marker for the longest ray segment distance we've found in a neighbor. Start with a quite
606 : // small number because at this point we're desperate and will take anything!
607 : Real longest_distance = 1.0e-12;
608 : // Temporaries for the intersection checks
609 : unsigned short current_incoming_side;
610 : Point current_intersection_point;
611 : unsigned short current_intersected_side;
612 : ElemExtrema current_intersected_extrema;
613 : Real current_intersection_distance;
614 : TraceRay::ExitsElemResult best_exit_result = NO_EXIT;
615 :
616 : // The NeighborInfo for the last_elem
617 : // (to store so we can try it later if we fail for everyone else)
618 : const NeighborInfo * last_elem_info = nullptr;
619 :
620 15453042 : for (const NeighborInfo & neighbor_info : neighbors)
621 : {
622 : // If we're at the element we want to do last, skip it and store the info so
623 : // we can get to it later in case we fail at finding anything
624 12866400 : if (neighbor_info._elem == last_elem)
625 : {
626 : debugRay("Skipping last elem ", last_elem->id());
627 : last_elem_info = &neighbor_info;
628 2588572 : continue;
629 : }
630 :
631 10277828 : const auto exit_result = moveThroughNeighbor(neighbor_info,
632 : current_incoming_side,
633 : current_intersection_point,
634 : current_intersected_side,
635 : current_intersected_extrema,
636 : current_intersection_distance);
637 :
638 : // Found one! See if it's the best way
639 10277828 : if (exit_result != NO_EXIT)
640 : {
641 : debugRay("Ray can exit through neighbor ", neighbor_info._elem->id());
642 :
643 3770989 : if (current_intersection_distance > longest_distance)
644 : {
645 2662424 : best_elem = neighbor_info._elem;
646 2662424 : best_elem_incoming_side = current_incoming_side;
647 2662424 : _intersection_point = current_intersection_point;
648 2662424 : _intersected_side = current_intersected_side;
649 : _intersected_extrema = current_intersected_extrema;
650 2662424 : _intersection_distance = current_intersection_distance;
651 : longest_distance = current_intersection_distance;
652 : best_exit_result = exit_result;
653 : }
654 : }
655 : }
656 :
657 : // Didn't find someone to exit, so try the last_elem (if any)
658 2586642 : if (!best_elem && last_elem_info)
659 : {
660 276 : const auto exit_result = moveThroughNeighbor(*last_elem_info,
661 : current_incoming_side,
662 : current_intersection_point,
663 : current_intersected_side,
664 : current_intersected_extrema,
665 : current_intersection_distance);
666 :
667 276 : if (exit_result != NO_EXIT && current_intersection_distance > longest_distance)
668 : {
669 : debugRay("Ray can exit through last_elem ", last_elem->id());
670 :
671 276 : best_elem = last_elem;
672 276 : best_elem_incoming_side = current_incoming_side;
673 276 : _intersection_point = current_intersection_point;
674 276 : _intersected_side = current_intersected_side;
675 : _intersected_extrema = current_intersected_extrema;
676 276 : _intersection_distance = current_intersection_distance;
677 : best_exit_result = exit_result;
678 : }
679 : }
680 :
681 : debugRay("moveThroughNeighbors() best result:");
682 : debugRay(" best_elem = ", best_elem ? best_elem->id() : DofObject::invalid_id);
683 : debugRay(" best_elem_incoming_side = ", best_elem_incoming_side);
684 : debugRay(" _intersection_point = ", _intersection_point);
685 : debugRay(" _intersected_side = ", _intersected_side);
686 : debugRay(" _intersected_extrema = ", _intersected_extrema);
687 : debugRay(" _intersection_distance = ", _intersection_distance);
688 : if (best_elem)
689 : {
690 : debugRay("moveThroughNeighbors() next neighbor elem info:");
691 : debugRay(best_elem->get_info());
692 : }
693 :
694 2586642 : return best_exit_result;
695 : }
696 :
697 : TraceRay::ExitsElemResult
698 10278104 : TraceRay::moveThroughNeighbor(const NeighborInfo & neighbor_info,
699 : unsigned short & incoming_side,
700 : Point & intersection_point,
701 : unsigned short & intersected_side,
702 : ElemExtrema & intersected_extrema,
703 : Real & intersection_distance)
704 : {
705 10278104 : if (!neighbor_info._valid)
706 : return NO_EXIT;
707 :
708 10277242 : const Elem * neighbor = neighbor_info._elem;
709 : debugRay("Checking neighbor ", neighbor->id(), " with centroid ", neighbor->vertex_average());
710 :
711 : // Find an entrant side (if any)
712 10277242 : incoming_side = RayTracingCommon::invalid_side;
713 13539636 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
714 13537343 : if (neighbor_info._side_normals[i] * (*_current_ray)->direction() < LOOSE_TRACE_TOLERANCE)
715 : {
716 10274949 : incoming_side = neighbor_info._sides[i];
717 10274949 : break;
718 : }
719 :
720 : // No entrant sides on this neighbor were found
721 10277242 : if (incoming_side == RayTracingCommon::invalid_side)
722 : return NO_EXIT;
723 :
724 10274949 : intersection_point = RayTracingCommon::invalid_point;
725 10274949 : intersected_side = RayTracingCommon::invalid_side;
726 : intersected_extrema.invalidate();
727 10274949 : intersection_distance = RayTracingCommon::invalid_distance;
728 :
729 : // See if there is a way through the neighbor element
730 : debugRay("Called exitsElem() from moveThroughNeighbor()");
731 : const auto exit_result =
732 10274949 : exitsElem(neighbor,
733 10274949 : neighbor->type(),
734 10274949 : incoming_side,
735 : intersection_point,
736 : intersected_side,
737 : intersected_extrema,
738 : intersection_distance,
739 10274949 : _backface_culling ? _study.getElemNormals(neighbor, _tid) : nullptr);
740 : debugRay("Done with exitsElem() from moveThroughNeighbor()");
741 :
742 10274949 : return exit_result;
743 : }
744 :
745 : void
746 972415 : TraceRay::applyOnExternalBoundary(const std::shared_ptr<Ray> & ray)
747 : {
748 : debugRay("Called applyOnExternalBoundary() with");
749 : debugRay(" _current_elem->id() = ", _current_elem->id());
750 : debugRay(" _intersection_point = ", _intersection_point);
751 : debugRay(" _intersected_side = ", _intersected_side);
752 : debugRay(" _intersected_extrema = ", _intersected_extrema);
753 :
754 : // Clear storage for the list of ConstBndElement that we need to apply RayBCs to
755 972415 : _boundary_elems.clear();
756 :
757 : // If on the elem extrema (at a vertex or within an edge), check for external sidesets on
758 : // neighbors (which will include _current_elem)
759 972415 : if (_dim != 1 && _intersected_extrema.atExtrema())
760 : {
761 587234 : const auto & neighbors = getNeighbors(_current_elem, _intersected_extrema, _intersection_point);
762 : debugRay(" Found ", neighbors.size(), " vertex/edge neighbors (including self)");
763 : traceAssert(std::count_if(neighbors.begin(),
764 : neighbors.end(),
765 : [this](const NeighborInfo & ni)
766 : { return ni._elem == _current_elem; }),
767 : "_current_elem not in neighbors");
768 :
769 2253514 : for (const auto & neighbor_info : neighbors)
770 : {
771 1666280 : if (!neighbor_info._valid)
772 0 : continue;
773 :
774 1666280 : const Elem * elem = neighbor_info._elem;
775 : const auto & sides = neighbor_info._sides;
776 : const auto & side_normals = neighbor_info._side_normals;
777 :
778 5950241 : for (MooseIndex(side_normals.size()) i = 0; i < side_normals.size(); ++i)
779 4283961 : if (!elem->neighbor_ptr(sides[i]) // is a boundary side that has our point
780 4283961 : && side_normals[i] * ray->direction() > TRACE_TOLERANCE) // and is entrant
781 : {
782 : // TODO: this could likely be optimized
783 : ElemExtrema extrema;
784 1420123 : withinExtremaOnSide(elem, _intersection_point, sides[i], _dim, extrema);
785 :
786 1420123 : _boundary_info.boundary_ids(elem, sides[i], _boundary_ids);
787 1420123 : possiblyAddToBoundaryElems(elem, sides[i], _boundary_ids, extrema);
788 : }
789 : }
790 : }
791 : // Not on the periphery (at vertex/edge), so we just need to add the external
792 : // sidesets from _current_elem on _intersected_side
793 : else
794 : {
795 385181 : _boundary_info.boundary_ids(_current_elem, _intersected_side, _boundary_ids);
796 385181 : possiblyAddToBoundaryElems(
797 385181 : _current_elem, _intersected_side, _boundary_ids, _intersected_extrema);
798 : }
799 :
800 : debugRay("Calling external onBoundary() with ", _boundary_elems.size(), " boundaries");
801 972415 : onBoundary(ray, /* external = */ true);
802 972405 : }
803 :
804 : void
805 1729 : TraceRay::applyOnInternalBoundary(const std::shared_ptr<Ray> & ray)
806 : {
807 : traceAssert(_last_elem, "Must be valid");
808 :
809 : debugRay("Called applyOnInternalBoundary() with");
810 : debugRay(" intersection_point = ", _intersection_point);
811 : debugRay(" _current_elem->id() = ", _current_elem->id());
812 : debugRay(" _incoming_side = ", _incoming_side);
813 : debugRay(" _last_elem->id() = ", _last_elem->id());
814 : debugRay(" _intersected_side = ", _intersected_side);
815 : debugRay(" _intersected_extrema = ", _intersected_extrema);
816 : traceAssert(_study.hasInternalSidesets(), "Do not have internal sidesets");
817 : traceAssert(_intersection_point.absolute_fuzzy_equals(_incoming_point, TRACE_TOLERANCE),
818 : "Intersection and incoming points should be the same");
819 :
820 : // Clear storage for the list of ConstBndElement that we need to apply RayBCs to
821 1729 : _boundary_elems.clear();
822 :
823 : ElemExtrema temp_extrema;
824 :
825 : // If on the elem extrema (at a vertex or within an edge), we need to check for internal
826 : // sidesets on neighbors (which will include _last_elem and _current_elem)
827 1729 : if (_dim != 1 && _intersected_extrema.atExtrema())
828 : {
829 : debugRay("Checking point neighbors for internal sidesets");
830 :
831 : // Get the neighbors
832 1182 : const auto & neighbors = getNeighbors(_last_elem, _intersected_extrema, _intersection_point);
833 : debugRay(" Found ", neighbors.size(), " vertex/edge neighbors");
834 : traceAssert(std::count_if(neighbors.begin(),
835 : neighbors.end(),
836 : [this](const NeighborInfo & ni)
837 : { return ni._elem == _current_elem || ni._elem == _last_elem; }) == 2,
838 : "_current_elem/_last_elem not in neighbors");
839 :
840 10362 : for (const auto & neighbor_info : neighbors)
841 : {
842 9180 : if (!neighbor_info._valid)
843 0 : continue;
844 :
845 9180 : const Elem * elem = neighbor_info._elem;
846 :
847 : // Grab the internal sidesets for this elem
848 9180 : const auto & sidesets = _study.getInternalSidesets(elem);
849 : // It has none to contribute, so we can continue
850 9180 : if (sidesets.empty())
851 : {
852 : debugRay(" Elem ", elem->id(), " has no internal sidesets");
853 7500 : continue;
854 : }
855 :
856 : // The sides on this elem that contain our point and their outward normals
857 : const auto & sides = neighbor_info._sides;
858 : const auto & side_normals = neighbor_info._side_normals;
859 :
860 : // See if any of the internal sidesets are relevant to this point
861 6036 : for (std::size_t i = 0; i < sides.size(); ++i)
862 : {
863 4356 : const auto side = sides[i];
864 : // Side has internal sidesets and is entrant
865 4356 : if (sidesets[side].size() && std::abs(side_normals[i] * ray->direction()) > TRACE_TOLERANCE)
866 : {
867 : // TODO: this could likely be optimized
868 : temp_extrema.invalidate();
869 1164 : withinExtremaOnSide(elem, _intersection_point, side, _dim, temp_extrema);
870 :
871 1164 : possiblyAddToBoundaryElems(elem, side, sidesets[side], temp_extrema);
872 : }
873 : }
874 : }
875 : }
876 : // Not on the periphery (at vertex/edge), so we just need to add the sidesets from _current_elem
877 : // on _incoming_side and _last_elem on _intersected_side
878 : else
879 : {
880 547 : const auto & current_elem_sidesets = _study.getInternalSidesets(_current_elem);
881 547 : if (current_elem_sidesets.size() && current_elem_sidesets[_incoming_side].size())
882 : {
883 : // This is possible but we need extrema checks on _incoming_elem to pass to the
884 : // boundary conditions. For now, we just pass in _intersected_extrema which could
885 : // be very wrong with adaptivity but is correct without it. I struggled with getting
886 : // the actual extrema checks and we're not using this now so I bet this error will
887 : // come back to haunt me in the future :-)
888 197 : if (!_study.hasSameLevelActiveElems())
889 0 : mooseError("Internal sidesets are not currently supported with adaptivity in tracing");
890 :
891 : // Special case for 1D
892 197 : if (_dim == 1)
893 26 : temp_extrema.setVertex(atVertex(_current_elem, _intersection_point));
894 :
895 197 : possiblyAddToBoundaryElems(_current_elem,
896 : _incoming_side,
897 197 : current_elem_sidesets[_incoming_side],
898 197 : _dim != 1 ? _intersected_extrema : temp_extrema);
899 : }
900 :
901 547 : const auto & last_elem_sidesets = _study.getInternalSidesets(_last_elem);
902 547 : if (last_elem_sidesets.size() && last_elem_sidesets[_intersected_side].size())
903 218 : possiblyAddToBoundaryElems(_last_elem,
904 : _intersected_side,
905 : last_elem_sidesets[_intersected_side],
906 218 : _intersected_extrema);
907 : }
908 :
909 1729 : if (!_boundary_elems.empty())
910 : {
911 : debugRay(" Calling internal onBoundary() with ", _boundary_elems.size(), " boundaries");
912 856 : onBoundary(ray, /* external = */ false);
913 : }
914 1725 : }
915 :
916 : void
917 1806883 : TraceRay::possiblyAddToBoundaryElems(const Elem * elem,
918 : const unsigned short side,
919 : const std::vector<BoundaryID> & bnd_ids,
920 : const ElemExtrema & extrema)
921 : {
922 : if (!_study.sideIsNonPlanar(elem, side))
923 : traceAssert(extrema.isValid(elem, _intersection_point), "Extrema not correct");
924 :
925 3613766 : for (const auto bnd_id : bnd_ids)
926 : {
927 : bool found = false;
928 1920959 : for (const auto & bnd_elem : _boundary_elems)
929 878679 : if (bnd_elem.bnd_id == bnd_id)
930 : {
931 : found = true;
932 : break;
933 : }
934 :
935 1806883 : if (!found)
936 : {
937 : debugRay(" Need to apply boundary on elem ",
938 : elem->id(),
939 : " and side ",
940 : side,
941 : " for bnd_id ",
942 : bnd_id);
943 :
944 1042280 : _boundary_elems.emplace_back(elem, side, bnd_id, extrema);
945 : }
946 : }
947 1806883 : }
948 :
949 : void
950 767413 : TraceRay::findExternalBoundarySide(unsigned short & boundary_side,
951 : ElemExtrema & boundary_extrema,
952 : const Elem *& boundary_elem)
953 : {
954 : traceAssert(_current_elem->neighbor_ptr(_intersected_side), "Already on boundary");
955 : traceAssert(boundary_side == RayTracingCommon::invalid_side, "Side should be invalid");
956 : traceAssert(boundary_extrema.isInvalid(), "Extrema should be invalid");
957 : traceAssert(!boundary_elem, "Elem should be invalid");
958 : traceAssert(_current_elem->dim() != 1, "1D traces shouldn't make it here");
959 : traceAssert(_current_elem->n_sides() == _current_elem_n_sides, "_current_elem_n_sides incorrect");
960 : traceAssert(_intersected_extrema.atExtrema(), "Should be at extrema");
961 : debugRay("Called findExternalBoundarySide() on side ", _intersected_side);
962 : debugRay(" _intersected_side = ", _intersected_side);
963 : debugRay(" _intersected_extrema", _intersected_extrema);
964 :
965 767413 : const auto & direction = (*_current_ray)->direction();
966 : const auto at_edge = _intersected_extrema.atEdge();
967 :
968 : // First, look for other sides on _current_elem that touch the intersected vertex/edge
969 : // that are on the boundary and are outgoing
970 4365852 : for (unsigned short s = 0; s < _current_elem_n_sides; ++s)
971 4664589 : if (!_current_elem->neighbor_ptr(s) && s != _intersected_side &&
972 1722049 : _current_elem->is_node_on_side(_intersected_extrema.first, s) &&
973 5045174 : (!at_edge || _current_elem->is_node_on_side(_intersected_extrema.second, s)) &&
974 732631 : !_study.sideIsIncoming(_current_elem, s, direction, _tid))
975 : {
976 : debugRay(" Side ", s, " is a boundary side and the Ray exits");
977 125034 : boundary_side = s;
978 : boundary_extrema = _intersected_extrema;
979 125034 : boundary_elem = _current_elem;
980 125034 : return;
981 : }
982 :
983 : // No luck in our element, so see if any neighbors at this vertex/edge are
984 : // on the boundary and are outgoing
985 642379 : const auto & neighbors = getNeighbors(_current_elem, _intersected_extrema, _intersection_point);
986 : debugRay("Checking current element failed, now checking neighbors");
987 :
988 : debugRay("Found ", neighbors.size(), " candidate neighbors (including self)");
989 2229898 : for (const auto & neighbor_info : neighbors)
990 : {
991 : // This will be false when we have an edge neighbor that isn't a neighbor at
992 : // _intersection_point
993 1696705 : if (!neighbor_info._valid)
994 0 : continue;
995 :
996 1696705 : const Elem * neighbor = neighbor_info._elem;
997 : // We've already checked ourself
998 1696705 : if (neighbor == _current_elem)
999 579305 : continue;
1000 :
1001 : debugRay("Checking neighbor ", neighbor->id());
1002 :
1003 : // Loop through the sides we touch and look for one that this Ray exits and is on the boundary
1004 3495090 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1005 2486876 : if (neighbor_info._side_normals[i] * direction > TRACE_TOLERANCE &&
1006 542017 : !neighbor->neighbor_ptr(neighbor_info._sides[i]))
1007 : {
1008 109186 : withinExtremaOnSide(
1009 109186 : neighbor, _intersection_point, neighbor_info._sides[i], _dim, boundary_extrema);
1010 : traceAssert(boundary_extrema.atExtrema(), "Should be at extrema");
1011 109186 : boundary_side = neighbor_info._sides[i];
1012 109186 : boundary_elem = neighbor;
1013 : return;
1014 : }
1015 : }
1016 : }
1017 :
1018 : void
1019 3267644 : TraceRay::trace(const std::shared_ptr<Ray> & ray)
1020 : {
1021 : mooseAssert(_study.currentlyPropagating(), "Should only use while propagating rays");
1022 :
1023 3267644 : _current_ray = &ray;
1024 3267644 : _current_elem = ray->currentElem();
1025 3267644 : _last_elem = nullptr;
1026 3267644 : _incoming_point = ray->currentPoint();
1027 3267644 : _incoming_side = ray->currentIncomingSide();
1028 3267644 : _should_continue = true;
1029 :
1030 3267644 : _study.preTrace(_tid, ray);
1031 :
1032 : traceAssert(_current_elem, "Current element is not set");
1033 : traceAssert(_current_elem->active(), "Current element is not active");
1034 : traceAssert(!ray->invalidCurrentPoint(), "Current point is invalid");
1035 : traceAssert(ray->shouldContinue(), "Ray should not continue");
1036 3267644 : if (_study.verifyRays() && !ray->invalidCurrentIncomingSide() && ray->maxDistance() > 0 &&
1037 5861901 : !_study.sideIsNonPlanar(_current_elem, _incoming_side) &&
1038 1286038 : !_study.sideIsIncoming(_current_elem, _incoming_side, ray->direction(), _tid))
1039 0 : failTrace("Ray incoming side is not incoming", /* warning = */ false, __LINE__);
1040 :
1041 : #ifdef DEBUG_RAY_IF
1042 : if (DEBUG_RAY_IF)
1043 : libMesh::err << "\n\n";
1044 : #endif
1045 : debugRay("At top of trace for Ray");
1046 : debugRay("Top of trace loop Ray info\n", ray->getInfo());
1047 : debugRay("Top of trace loop starting elem info\n", ray->currentElem()->get_info());
1048 :
1049 : // Invalidate this up front because it's copied immediately into _last_intersected_extrema
1050 : _intersected_extrema.invalidate();
1051 :
1052 : #ifdef DEBUG_RAY_MESH_IF
1053 : _debug_mesh = nullptr;
1054 : _debug_node_count = 0;
1055 : if (DEBUG_RAY_MESH_IF)
1056 : {
1057 : _debug_mesh = new Mesh(_debug_comm, _dim);
1058 : _debug_mesh->skip_partitioning(true);
1059 : }
1060 : #endif
1061 :
1062 : // Caching trace along the way: init for this Ray
1063 3267644 : if (_study.shouldCacheTrace(ray))
1064 : {
1065 : debugRay("Trying to init threaded cached trace");
1066 :
1067 8576 : _current_cached_trace = &_study.initThreadedCachedTrace(ray, _tid);
1068 :
1069 : // Add starting data
1070 8576 : if (_study.dataOnCacheTraces())
1071 904 : _current_cached_trace->lastPoint()._data = ray->data();
1072 8576 : if (_study.auxDataOnCacheTraces())
1073 876 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1074 : }
1075 : else
1076 3259068 : _current_cached_trace = nullptr;
1077 :
1078 : // Need to call subdomain setup
1079 3267644 : if (_current_elem->subdomain_id() != _current_subdomain_id || _study.rayDependentSubdomainSetup())
1080 3267644 : onSubdomainChanged(ray, /* same_ray = */ false);
1081 : // If we didn't change subdomains, we still need to call this on the RayKernels
1082 : else
1083 0 : for (RayKernelBase * rk : _study.currentRayKernels(_tid))
1084 0 : rk->preTrace();
1085 :
1086 : // Ray tracing loop through each segment for a single Ray on this processor
1087 : do
1088 : {
1089 : #ifdef DEBUG_RAY_IF
1090 : if (DEBUG_RAY_IF)
1091 : libMesh::err << "\n\n";
1092 : #endif
1093 : debugRay("At top of ray tracing loop");
1094 : debugRay(" ray->id() = ", ray->id());
1095 : debugRay(" _incoming_point = ", _incoming_point);
1096 : debugRay(" _incoming_side = ", _incoming_side);
1097 : debugRay(" _current_elem->id() = ", _current_elem->id());
1098 : debugRay(" _current_elem->subdomain_id() = ", _current_elem->subdomain_id());
1099 : debugRay(" _current_subdomain_id = ", _current_subdomain_id);
1100 : debugRay(" _current_elem_type = ", Utility::enum_to_string(_current_elem_type));
1101 : debugRay("Top of ray tracing loop Ray info\n", ray->getInfo());
1102 : debugRay("Top of ray tracing loop current elem info\n", _current_elem->get_info());
1103 :
1104 : traceAssert(_current_ray == &ray, "Current ray mismatch");
1105 :
1106 : // Copy over in case we need to use it
1107 : _last_intersected_extrema = _intersected_extrema;
1108 : // Invalidate all intersections as we're tracing again
1109 23718768 : _exits_elem = false;
1110 23718768 : _intersection_point = RayTracingCommon::invalid_point;
1111 23718768 : _intersected_side = RayTracingCommon::invalid_side;
1112 : _intersected_extrema.invalidate();
1113 23718768 : _intersection_distance = RayTracingCommon::invalid_distance;
1114 :
1115 : // Stationary ray
1116 23718768 : if (ray->stationary())
1117 : {
1118 : mooseAssert(ray->invalidDirection(), "Should have an invalid direction");
1119 1050 : _exits_elem = true;
1120 1050 : _intersection_point = _incoming_point;
1121 : _intersected_extrema = _last_intersected_extrema;
1122 1050 : _intersection_distance = 0;
1123 1050 : ++_results[ENDED_STATIONARY];
1124 : }
1125 : // If we haven't hit a vertex or an edge, do the normal exit algorithm first.
1126 : // In the case of a Ray that previously moved through point neighbors due to
1127 : // being at a vertex/edge, this will be true because we do not communicate the
1128 : // vertex/edge intersection. The previous processor already set us up for
1129 : // the best intersection because it already computed it and chose to
1130 : // send it our way as such
1131 23717718 : else if (!_last_intersected_extrema.atExtrema())
1132 : {
1133 : traceAssert(_current_elem->processor_id() == _pid, "Trace elem not on processor");
1134 : debugRay("Didn't hit vertex or edge: doing normal exits elem check");
1135 :
1136 21154699 : if (_backface_culling)
1137 630809 : _current_normals = _study.getElemNormals(_current_elem, _tid);
1138 :
1139 21154699 : const auto exits_elem_result = exitsElem(_current_elem,
1140 : _current_elem_type,
1141 21154699 : _incoming_side,
1142 21154699 : _intersection_point,
1143 21154699 : _intersected_side,
1144 21154699 : _intersected_extrema,
1145 21154699 : _intersection_distance,
1146 : _current_normals);
1147 :
1148 21154699 : if (exits_elem_result != NO_EXIT)
1149 : {
1150 21131076 : storeExitsElemResult(exits_elem_result);
1151 21131076 : _exits_elem = true;
1152 : ray->setCurrentPoint(_intersection_point);
1153 : }
1154 : }
1155 : else
1156 : {
1157 : debugRay("Will not do normal exits elem check because at a vertex/edge");
1158 : }
1159 :
1160 : // At this point, we either did a regular exit check and it failed, or we hit a vertex/edge
1161 : // on the previous intersection and didn't bother to do a regular exit check on
1162 : // _current_elem
1163 23718768 : if (!_exits_elem)
1164 : {
1165 : debugRay("Moving through neighbors");
1166 :
1167 : // The element we want moveThroughNeighbors() to try last. That is - if all others fail,
1168 : // this is the last resort. If we have a last intersection that tells us we're going
1169 : // through a vertex or edge, we probably won't go back to that elem. If we don't, it means
1170 : // that we failed the exitsElem() check on _current_elem and probably won't return to it
1171 : // either.
1172 : const Elem * move_through_neighbors_last = nullptr;
1173 :
1174 : // Get the neighbors
1175 : const std::vector<NeighborInfo> * neighbors = nullptr;
1176 : // If we have previous vertex/edge information, use it
1177 2586642 : if (_last_intersected_extrema.atExtrema())
1178 : {
1179 : traceAssert(_last_elem, "Should be valid");
1180 2563019 : move_through_neighbors_last = _last_elem;
1181 2563019 : neighbors = &getNeighbors(_last_elem, _last_intersected_extrema, _incoming_point);
1182 : }
1183 : // Without previous vertex/edge information, let's try to find some. We could just do
1184 : // a pure point neighbor check at the current point, but we'd rather be able to cache
1185 : // this information so that someone else can use it - this requires the more unique
1186 : // identifier of which vertex/edge we are at.
1187 : else
1188 : {
1189 : debugRay(" Searching for vertex/edge hit with incoming side ",
1190 : _incoming_side,
1191 : " on elem ",
1192 : _current_elem->id(),
1193 : " at ",
1194 : _incoming_point);
1195 :
1196 23623 : move_through_neighbors_last = _current_elem;
1197 :
1198 : // If we have side info: check for vertices on said side, otherwise, check everywhere
1199 23623 : const auto at_v = _incoming_side != RayTracingCommon::invalid_side
1200 23623 : ? atVertexOnSide(_current_elem, _incoming_point, _incoming_side)
1201 971 : : atVertex(_current_elem, _incoming_point);
1202 :
1203 23623 : if (at_v != RayTracingCommon::invalid_vertex)
1204 23601 : neighbors = &getVertexNeighbors(_current_elem, at_v);
1205 : // If still nothing and in 3D, check if we're on an edge
1206 : // TODO: Handle 2D with a similar check for if we're on a side
1207 22 : else if (_dim == 3)
1208 : {
1209 : ElemExtrema extrema;
1210 0 : if (_incoming_side != RayTracingCommon::invalid_side)
1211 0 : withinEdgeOnSide(_current_elem, _incoming_point, _incoming_side, extrema);
1212 : else
1213 0 : withinEdge(_current_elem, _incoming_point, extrema);
1214 :
1215 : if (extrema.atEdge())
1216 0 : neighbors = &getEdgeNeighbors(_current_elem, extrema, _incoming_point);
1217 : }
1218 :
1219 : // If we still haven't found anything - let's try a last-ditch effort
1220 23601 : if (!neighbors || neighbors->empty())
1221 22 : neighbors = &getPointNeighbors(_current_elem, _incoming_point);
1222 :
1223 : // Couldn't find anything
1224 23623 : if (neighbors->empty())
1225 : {
1226 0 : failTrace("Could not find neighbors to move through", _study.tolerateFailure(), __LINE__);
1227 101273 : return;
1228 : }
1229 : }
1230 :
1231 : // Move through a neighbor
1232 2586642 : const Elem * best_neighbor = nullptr;
1233 2586642 : auto best_neighbor_side = RayTracingCommon::invalid_side;
1234 2586642 : const auto exits_elem_result = moveThroughNeighbors(
1235 : *neighbors, move_through_neighbors_last, best_neighbor, best_neighbor_side);
1236 :
1237 : // If we didn't find anything... we're out of luck
1238 2586642 : if (exits_elem_result == NO_EXIT)
1239 : {
1240 0 : failTrace("Could not find intersection after trying to move through point neighbors",
1241 0 : _study.tolerateFailure(),
1242 : __LINE__);
1243 0 : return;
1244 : }
1245 :
1246 : // At this point, we've successfully made it through an element and also started
1247 : // through another element, so set that
1248 2586642 : _exits_elem = true;
1249 2586642 : _last_elem = _current_elem;
1250 2586642 : _current_elem = best_neighbor;
1251 2586642 : _incoming_side = best_neighbor_side;
1252 : ray->setCurrentElem(best_neighbor);
1253 : ray->setCurrentIncomingSide(best_neighbor_side);
1254 : ray->setCurrentPoint(_intersection_point);
1255 :
1256 : // Don't own this element - return exits trace for this Ray on this proc
1257 2586642 : if (best_neighbor->processor_id() != _pid)
1258 : {
1259 : // We've already computed the next intersection but said intersection is
1260 : // on an elem on another processor. Therefore, the next proc will re-trace
1261 : // this Ray on the perfect elem that we've picked (best_neighbor)
1262 : ray->setCurrentPoint(_incoming_point);
1263 101273 : _intersection_distance = RayTracingCommon::invalid_distance;
1264 :
1265 101273 : continueTraceOffProcessor(ray);
1266 : return;
1267 : }
1268 :
1269 : // Subdomain changed
1270 2485369 : if (_current_elem->subdomain_id() != _current_subdomain_id)
1271 146 : onSubdomainChanged(ray, /* same_ray = */ true);
1272 :
1273 : // Do own this element - tally the result as we're the ones tracing it
1274 2485369 : storeExitsElemResult(exits_elem_result);
1275 : }
1276 :
1277 : debugRay("Done with trace");
1278 : debugRay(" _exits_elem: ", _exits_elem);
1279 : debugRay(" _intersection_point: ", _intersection_point);
1280 : debugRay(" _intersected_side: ", _intersected_side);
1281 : debugRay(" _intersected_side centroid: ",
1282 : _intersected_side == RayTracingCommon::invalid_side
1283 : ? RayTracingCommon::invalid_point
1284 : : _current_elem->side_ptr(_intersected_side)->vertex_average());
1285 : debugRay(" _intersected_extrema = ", _intersected_extrema);
1286 : debugRay(" _intersection_distance = ", _intersection_distance);
1287 :
1288 : // Increment intersections
1289 23617495 : if (_intersection_distance > 0)
1290 : {
1291 : debugRay("Incrementing ray intersections by 1 to ", ray->intersections() + 1);
1292 : ray->addIntersection();
1293 23616445 : _results[INTERSECTIONS]++;
1294 : }
1295 :
1296 : // Increment distance
1297 : ray->addDistance(_intersection_distance);
1298 : debugRay("Incremented ray distance by ", _intersection_distance);
1299 :
1300 : // The effective max distance that the Ray should travel - minimum of the two
1301 23700535 : const auto max_distance = std::min(ray->maxDistance(), _study.rayMaxDistance());
1302 : debugRay("Max distance checks");
1303 : debugRay(" _study.rayMaxDistance() = ", _study.rayMaxDistance());
1304 : debugRay(" ray->maxDistance() = ", ray->maxDistance());
1305 : debugRay(" max_distance (effective) = ", max_distance);
1306 :
1307 : // At the maximum distance - distinguish between this and past the maximum
1308 : // distance because in this case we're close enough to the intersection point
1309 : // that we can keep it, its intersected side, and intersected extrema
1310 23617495 : if (MooseUtils::absoluteFuzzyEqual(ray->distance(), max_distance))
1311 : {
1312 : debugRay("At max distance");
1313 :
1314 : ray->setShouldContinue(false);
1315 15042 : _should_continue = false;
1316 : }
1317 : // Past the max distance - need to remove the additional distance we traveled,
1318 : // change the point and invalidate the intersection data (moves the Ray back)
1319 23602453 : else if (ray->distance() > max_distance)
1320 : {
1321 : debugRay("Past max distance");
1322 :
1323 : // The distance past the max distance we have traveled
1324 : const auto difference = ray->distance() - max_distance;
1325 : traceAssert(difference > 0, "Negative distance change after past_max_distance");
1326 :
1327 : debugRay("Removing distance ", difference);
1328 : ray->addDistance(-difference);
1329 : debugRay(" New ray->distance() = ", ray->distance());
1330 :
1331 1801528 : _intersection_point -= ray->direction() * difference;
1332 1801528 : _intersection_distance -= difference;
1333 1801528 : _intersected_side = RayTracingCommon::invalid_side;
1334 : _intersected_extrema.invalidate();
1335 : ray->setCurrentPoint(_intersection_point);
1336 :
1337 : ray->setShouldContinue(false);
1338 1801528 : _should_continue = false;
1339 :
1340 : traceAssert(_intersection_distance >= 0, "Negative _intersection_distance");
1341 : #ifndef NDEBUG
1342 : if (_study.verifyTraceIntersections() && !_current_elem->contains_point(_intersection_point))
1343 : failTrace("Does not contain point after past max distance",
1344 : /* warning = */ false,
1345 : __LINE__);
1346 : #endif
1347 : }
1348 :
1349 23617495 : if (!_study.currentRayKernels(_tid).empty())
1350 : {
1351 : debugRay("Calling onSegment() with");
1352 : debugRay(" current_elem->id() = ", _current_elem->id());
1353 : debugRay(" _incoming_point = ", _incoming_point);
1354 : debugRay(" _incoming_side = ", _incoming_side);
1355 : debugRay(" _intersection_point = ", _intersection_point);
1356 : debugRay(" _intersected_side = ", _intersected_side);
1357 : debugRay(" _intersected_extrema = ", _intersected_extrema);
1358 : debugRay(" _intersection_distance = ", _intersection_distance);
1359 22560690 : onSegment(ray);
1360 :
1361 22560678 : _study.postOnSegment(_tid, ray);
1362 :
1363 : // RayKernel killed a Ray or we're at the end
1364 : traceAssert(_should_continue == ray->shouldContinue(), "Should be the same");
1365 22560678 : if (!_should_continue)
1366 : {
1367 : debugRay("RayKernel killed the ray or past max distance");
1368 : traceAssert(!ray->trajectoryChanged(),
1369 : "RayKernels should not change trajectories of Rays at end");
1370 :
1371 1696000 : onCompleteTrace(ray);
1372 1696000 : return;
1373 : }
1374 :
1375 : // RayKernel moved a Ray
1376 20864678 : if (ray->trajectoryChanged())
1377 : {
1378 : debugRay("RayKernel changed the Ray's trajectory");
1379 : debugRay(" new direction = ", ray->direction());
1380 : debugRay(" old incoming point = ", _incoming_point);
1381 : debugRay(" new incoming point = ", ray->currentPoint());
1382 : possiblyAddDebugRayMeshPoint(_incoming_point, ray->currentPoint());
1383 :
1384 432 : _incoming_point = ray->currentPoint();
1385 432 : _incoming_side = RayTracingCommon::invalid_side;
1386 : _intersected_extrema.invalidate();
1387 : ray->setCurrentIncomingSide(_incoming_side);
1388 :
1389 432 : const auto new_intersection_distance = (ray->currentPoint() - _incoming_point).norm();
1390 432 : ray->addDistance(-_intersection_distance + new_intersection_distance);
1391 432 : _intersection_distance = new_intersection_distance;
1392 :
1393 432 : onTrajectoryChanged(ray);
1394 432 : onContinueTrace(ray);
1395 : continue;
1396 432 : }
1397 : else
1398 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1399 : }
1400 : else
1401 : {
1402 1056805 : _study.postOnSegment(_tid, ray);
1403 :
1404 1056805 : if (!_should_continue)
1405 : {
1406 : debugRay("Killing due to at end without RayKernels");
1407 : traceAssert(!ray->shouldContinue(), "Ray shouldn't continue");
1408 :
1409 120568 : onCompleteTrace(ray);
1410 120568 : return;
1411 : }
1412 : }
1413 :
1414 : // If at a vertex/on an edge and not on the boundary, we may actually be on the boundary,
1415 : // just not on a boundary side. Check for that here. If the domain is rectangular we will
1416 : // do a quick check against the bounding box as if we're not on the boundary of the
1417 : // bounding box for a rectangular problem, we can skip this.
1418 : // TODO: see how much the tolerance for the bounding box check can be tightened
1419 21794633 : if (_dim > 1 && _intersected_extrema.atExtrema() &&
1420 24944807 : _current_elem->neighbor_ptr(_intersected_side) &&
1421 5582620 : (!_is_rectangular_domain ||
1422 2791310 : onBoundingBoxBoundary(_study.boundingBox(),
1423 2791310 : _intersection_point,
1424 : _dim,
1425 2791310 : LOOSE_TRACE_TOLERANCE * _study.domainMaxLength())))
1426 : {
1427 767413 : auto boundary_side = RayTracingCommon::invalid_side;
1428 : ElemExtrema boundary_extrema;
1429 767413 : const Elem * boundary_elem = nullptr;
1430 :
1431 767413 : findExternalBoundarySide(boundary_side, boundary_extrema, boundary_elem);
1432 :
1433 767413 : if (boundary_elem)
1434 : {
1435 : // At this point, the new incoming side is very difficult to find
1436 : // and may require some re-tracing backwards. Let's not do that -
1437 : // we don't need it to continue the trace. This is reason why
1438 : // _incoming_side is not available in RayBCs.
1439 234220 : _last_elem = _current_elem;
1440 234220 : _current_elem = boundary_elem;
1441 234220 : _intersected_side = boundary_side;
1442 : _intersected_extrema = boundary_extrema;
1443 : ray->setCurrentElem(_current_elem);
1444 :
1445 : debugRay("Found a neighbor boundary side with:");
1446 : debugRay(" _current_elem->id() = ", _current_elem->id());
1447 : debugRay(" _intersected_side = ", _intersected_side);
1448 : debugRay(" _intersected_extrema = ", _intersected_extrema);
1449 : }
1450 : }
1451 :
1452 : // The next element
1453 : const Elem * neighbor = nullptr;
1454 :
1455 : // Set up where to go next
1456 21800483 : _incoming_point = _intersection_point;
1457 : debugRay("Set _incoming point to ", _incoming_point);
1458 :
1459 : debugRay("Looking for neighbor on side ", _intersected_side, " of elem ", _current_elem->id());
1460 :
1461 : // Get the neighbor on that side
1462 : // If the mesh has active elements of the same level, just grab the neighbor directly. In
1463 : // this case, we can avoid a call to active() on the neighbor in getActiveNeighbor(), which
1464 : // can be expensive depending on caching.
1465 21800483 : neighbor = _study.hasSameLevelActiveElems()
1466 21800483 : ? _current_elem->neighbor_ptr(_intersected_side)
1467 126255 : : getActiveNeighbor(_current_elem, _intersected_side, _intersection_point);
1468 :
1469 : // Found one - not on the boundary so set the next info
1470 21800483 : if (neighbor)
1471 : {
1472 : traceAssert(neighbor->active(), "Inactive neighbor");
1473 : traceAssert(_current_elem->subdomain_id() == _current_subdomain_id,
1474 : "_current_subdomain_id invalid");
1475 :
1476 : // If the mesh has active elements of the same level, don't use
1477 : // neighbor->which_neighbor_am_i(), which calls active() and an n_sides() virtual call
1478 20828068 : if (_study.hasSameLevelActiveElems())
1479 : {
1480 : // If the subdomain hasn't changed, we can guarantee n_sides is the same as
1481 : // _current_elem. Prefer this because neighbor->n_sides() is an expensive virtual. Even
1482 : // better, subdomain_id() isn't a virtual!
1483 20711077 : const unsigned short n_sides = neighbor->subdomain_id() == _current_subdomain_id
1484 : ? _current_elem_n_sides
1485 20711077 : : neighbor->n_sides();
1486 : traceAssert(n_sides == neighbor->n_sides(), "n_sides incorrect");
1487 :
1488 63285502 : for (_incoming_side = 0; _incoming_side < n_sides; ++_incoming_side)
1489 63285502 : if (neighbor->neighbor_ptr(_incoming_side) == _current_elem)
1490 : break;
1491 : }
1492 : else
1493 116991 : _incoming_side = neighbor->which_neighbor_am_i(_current_elem);
1494 :
1495 20828068 : _last_elem = _current_elem;
1496 20828068 : _current_elem = neighbor;
1497 : ray->setCurrentElem(neighbor);
1498 20828068 : ray->setCurrentIncomingSide(_incoming_side);
1499 :
1500 : debugRay("Next elem: ", neighbor->id(), " with centroid ", neighbor->vertex_average());
1501 : debugRay("Next _incoming_side: ",
1502 : _incoming_side,
1503 : " with centroid ",
1504 : neighbor->side_ptr(_incoming_side)->vertex_average());
1505 : traceAssert(_last_elem->subdomain_id() == _current_subdomain_id,
1506 : "_current_subdomain_id invalid");
1507 :
1508 : // Whether or not the subdomain changed
1509 20828068 : const bool subdomain_changed = neighbor->subdomain_id() != _current_subdomain_id;
1510 :
1511 : // See if we hit any internal sides leaving this element and apply
1512 : // We require that all internal RayBCs be on internal sidesets that have different
1513 : // subdomain ids on each side. Therefore, only check if at vertex/edge or if
1514 : // the subdomain changes
1515 20828068 : if (_study.hasInternalSidesets() && (subdomain_changed || _intersected_extrema.atExtrema()))
1516 : {
1517 1729 : applyOnInternalBoundary(ray);
1518 :
1519 : // Internal RayBC killed a Ray
1520 1725 : if (!_should_continue)
1521 : {
1522 : traceAssert(!ray->shouldContinue(), "Should be the same");
1523 : debugRay("Internal RayBC killed the ray");
1524 :
1525 159 : onCompleteTrace(ray);
1526 159 : return;
1527 : }
1528 :
1529 : // Internal RayBC changed the Ray
1530 1566 : if (ray->trajectoryChanged())
1531 : {
1532 : debugRay("Internal RayBC changed the trajectory:");
1533 : debugRay(" new direction = ", ray->direction());
1534 :
1535 : // Is this side still incoming?
1536 297 : const auto normal = _study.getSideNormal(_current_elem, _incoming_side, _tid);
1537 : const auto dot = normal * ray->direction();
1538 : debugRay("Dot product with new direction and side = ", dot);
1539 297 : if (dot > -TRACE_TOLERANCE)
1540 : {
1541 273 : _incoming_side = _last_elem->which_neighbor_am_i(_current_elem);
1542 273 : _current_elem = _last_elem;
1543 : neighbor = _last_elem;
1544 : ray->setCurrentElem(_current_elem);
1545 : ray->setCurrentIncomingSide(_incoming_side);
1546 : debugRay(" Dot > 0 (Ray turned around): Setting _current_elem = ",
1547 : _current_elem->id(),
1548 : " and _incoming_side = ",
1549 : _incoming_side);
1550 : }
1551 :
1552 297 : onTrajectoryChanged(ray);
1553 : }
1554 :
1555 : traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point, TRACE_TOLERANCE),
1556 : "Internal RayBC changed the Ray point");
1557 : }
1558 :
1559 : // Neighbor is off processor
1560 : // If we hit at a vertex/edge, we will continue and let the move through neighbor exit
1561 : // figure out who to send to
1562 20827905 : if (neighbor->processor_id() != _pid)
1563 : {
1564 466730 : if (_intersected_extrema.atExtrema())
1565 : {
1566 : debugRay("Neighbor is off processor but continuing to move through neighbors");
1567 : }
1568 : else
1569 : {
1570 408258 : continueTraceOffProcessor(ray);
1571 408258 : return;
1572 : }
1573 : }
1574 :
1575 : // Neighbor is on processor, call subdomain setup if needed
1576 20419647 : if (subdomain_changed)
1577 1687 : onSubdomainChanged(ray, /* same_ray = */ true);
1578 : }
1579 : // No neighbor found: on the boundary
1580 : else
1581 : {
1582 : debugRay("No neighbor found - on the boundary");
1583 :
1584 : // Apply boundary conditions
1585 972415 : applyOnExternalBoundary(ray);
1586 :
1587 972405 : if (_current_elem == ray->currentElem())
1588 : traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point, TRACE_TOLERANCE),
1589 : "RayBC changed the Ray point");
1590 :
1591 : // Quit tracing if the Ray was killed by a BC
1592 972405 : if (!_should_continue)
1593 : {
1594 : traceAssert(!ray->shouldContinue(), "Should be the same");
1595 : debugRay("Exiting due to death by BC");
1596 :
1597 938538 : onCompleteTrace(ray);
1598 938538 : return;
1599 : }
1600 : // RayBC changed the direction of the Ray
1601 33867 : if (ray->trajectoryChanged())
1602 : {
1603 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1604 :
1605 33867 : _last_elem = _current_elem;
1606 : // Direction changed
1607 33867 : if (_current_elem == ray->currentElem())
1608 : {
1609 : debugRay("RayBC reflected the ray");
1610 : debugRay(" new direction = ", ray->direction());
1611 : traceAssert(ray->direction() *
1612 : _study.getSideNormal(_current_elem, _intersected_side, _tid) <
1613 : TRACE_TOLERANCE,
1614 : "Reflected ray is not incoming");
1615 :
1616 2181 : _incoming_point = _intersection_point;
1617 2181 : _incoming_side = _intersected_side;
1618 : ray->setCurrentPoint(_incoming_point);
1619 : ray->setCurrentIncomingSide(_incoming_side);
1620 : }
1621 : // Position changed (PeriodicRayBC)
1622 : else
1623 : {
1624 : debugRay("RayBC moved the ray");
1625 : debugRay(" new point = ", ray->currentPoint());
1626 : debugRay(" new pid = ", ray->currentElem()->processor_id());
1627 : debugRay(" new elem id = ", ray->currentElem()->id());
1628 : debugRay(" new side = ", ray->currentIncomingSide());
1629 :
1630 31686 : _current_elem = ray->currentElem();
1631 31686 : _incoming_point = ray->currentPoint();
1632 31686 : _incoming_side = ray->currentIncomingSide();
1633 : _intersected_extrema.invalidate();
1634 :
1635 31686 : if (_current_elem->processor_id() != _pid)
1636 : {
1637 2822 : continueTraceOffProcessor(ray);
1638 2822 : return;
1639 : }
1640 : }
1641 31045 : onTrajectoryChanged(ray);
1642 : }
1643 : }
1644 :
1645 20450692 : onContinueTrace(ray);
1646 :
1647 : } while (true);
1648 :
1649 : // If a trace made its way down here and didn't return... it failed
1650 : failTrace("Could not find an intersection", _study.tolerateFailure(), __LINE__);
1651 : }
1652 :
1653 : void
1654 2755265 : TraceRay::onCompleteTrace(const std::shared_ptr<Ray> & ray)
1655 : {
1656 5287231 : for (RayKernelBase * rk : _study.currentRayKernels(_tid))
1657 2531966 : rk->postTrace();
1658 :
1659 : debugRay("Called onCompleteTrace()\n", (*_current_ray)->getInfo());
1660 : if (_intersection_distance > 0)
1661 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1662 : possiblySaveDebugRayMesh();
1663 :
1664 2755265 : if (_current_cached_trace)
1665 : {
1666 7571 : _current_cached_trace->_last = true;
1667 :
1668 7571 : if (_intersection_distance > 0)
1669 : {
1670 : _current_cached_trace->addPoint(ray->currentPoint());
1671 6833 : if (_study.dataOnCacheTraces())
1672 118 : _current_cached_trace->lastPoint()._data = ray->data();
1673 6833 : if (_study.auxDataOnCacheTraces())
1674 94 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1675 : }
1676 :
1677 : mooseAssert(ray->stationary() == _current_cached_trace->stationary(), "Stationary mismatch");
1678 : }
1679 2755265 : }
1680 :
1681 : void
1682 20451124 : TraceRay::onContinueTrace(const std::shared_ptr<Ray> & ray)
1683 : {
1684 : traceAssert(ray->shouldContinue(), "Ray must continue");
1685 :
1686 20451124 : if (_current_cached_trace && _study.segmentsOnCacheTraces() && _intersection_distance > 0)
1687 : {
1688 : _current_cached_trace->addPoint(ray->currentPoint());
1689 30621 : if (_study.dataOnCacheTraces())
1690 453 : _current_cached_trace->lastPoint()._data = ray->data();
1691 30621 : if (_study.auxDataOnCacheTraces())
1692 357 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1693 : }
1694 20451124 : }
1695 :
1696 : void
1697 512353 : TraceRay::continueTraceOffProcessor(const std::shared_ptr<Ray> & ray)
1698 : {
1699 : traceAssert(ray->currentElem() == _current_elem, "Ray currentElem() invalid");
1700 : traceAssert(ray->currentIncomingSide() == _incoming_side, "Ray currentIncomingSide() invalid");
1701 : traceAssert(ray->currentPoint() == _incoming_point, "Ray currentPoint() invalid");
1702 : traceAssert(_current_elem->processor_id() != _pid, "Off processor trace is not off processor");
1703 : debugRay("Ray going off processor to ", _current_elem->processor_id());
1704 :
1705 : ray->addProcessorCrossing();
1706 :
1707 512353 : if (_current_cached_trace && _intersection_distance > 0)
1708 : {
1709 894 : _current_cached_trace->addPoint(_incoming_point);
1710 894 : if (_study.dataOnCacheTraces())
1711 11 : _current_cached_trace->lastPoint()._data = ray->data();
1712 894 : if (_study.auxDataOnCacheTraces())
1713 11 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1714 : }
1715 :
1716 : if (_intersection_distance > 0)
1717 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1718 : possiblySaveDebugRayMesh();
1719 512353 : }
1720 :
1721 : void
1722 31774 : TraceRay::onTrajectoryChanged(const std::shared_ptr<Ray> & ray)
1723 : {
1724 : #ifndef NDEBUG
1725 : if (_study.verifyTraceIntersections() &&
1726 : (_intersected_extrema.atExtrema()
1727 : ? !_current_elem->close_to_point(ray->currentPoint(), LOOSE_TRACE_TOLERANCE)
1728 : : !_current_elem->contains_point(ray->currentPoint())))
1729 : failTrace("Elem does not contain point after trajectory change",
1730 : /* warning = */ false,
1731 : __LINE__);
1732 : #endif
1733 :
1734 : traceAssert(ray->shouldContinue(), "Ray should continue when trajectory is being changed");
1735 :
1736 : ray->setTrajectoryChanged(false);
1737 : ray->addTrajectoryChange();
1738 :
1739 31774 : if (_current_cached_trace && !_study.segmentsOnCacheTraces())
1740 : {
1741 40 : if (_intersection_distance > 0)
1742 : _current_cached_trace->addPoint(ray->currentPoint());
1743 40 : if (_study.dataOnCacheTraces())
1744 0 : _current_cached_trace->lastPoint()._data = ray->data();
1745 40 : if (_study.auxDataOnCacheTraces())
1746 0 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1747 : }
1748 31774 : }
1749 :
1750 : void
1751 3269477 : TraceRay::onSubdomainChanged(const std::shared_ptr<Ray> & ray, const bool same_ray)
1752 : {
1753 : debugRay("Calling onSubdomainChanged() on subdomain ", _current_elem->subdomain_id());
1754 : debugRay(" _current_subdomain_id = ", _current_subdomain_id);
1755 :
1756 3269477 : _current_subdomain_id = _current_elem->subdomain_id();
1757 3269477 : _current_subdomain_hmax = _study.subdomainHmax(_current_subdomain_id);
1758 3269477 : _current_elem_type = _current_elem->type();
1759 3269477 : _current_elem_n_sides = _current_elem->n_sides();
1760 :
1761 3269477 : if (_has_ray_kernels)
1762 : {
1763 2997087 : auto & current_ray_kernels = _study.currentRayKernels(_tid);
1764 :
1765 : // If we're still tracing the same Ray, keep track of our old RayKernels
1766 : // so that we don't call preTrace() on them again
1767 2997087 : if (same_ray)
1768 1130 : _old_ray_kernels.insert(current_ray_kernels.begin(), current_ray_kernels.end());
1769 : // If we're not tracing the same Ray, we need to call preTrace() on everything
1770 : else
1771 : _old_ray_kernels.clear();
1772 :
1773 : // Call segmentSubdomainSetup to get new kernels etc
1774 2997087 : _study.segmentSubdomainSetup(_current_subdomain_id, _tid, ray->id());
1775 :
1776 : // Call preTrace() on all of the RayKernels that need it
1777 5998762 : for (RayKernelBase * rk : current_ray_kernels)
1778 : // Haven't called preTrace() for this Ray on this RayKernel yet
1779 : if (!_old_ray_kernels.count(rk))
1780 3001629 : rk->preTrace();
1781 : }
1782 3269477 : }
1783 :
1784 : std::string
1785 8 : TraceRay::failTraceMessage(const std::string & reason, const int line)
1786 : {
1787 8 : std::stringstream oss;
1788 8 : oss << "Ray on processor " << _pid << " and thread " << _tid << " failed to trace";
1789 8 : if (line != -1)
1790 8 : oss << " at line " << line;
1791 16 : oss << "\n\n" << reason << "\n\n";
1792 16 : oss << ((*_current_ray))->getInfo() << "\n";
1793 8 : oss << "Current trace information\n";
1794 8 : oss << " _current_subdomain_id = ";
1795 8 : if (_current_subdomain_id == Elem::invalid_subdomain_id)
1796 0 : oss << "invalid subdomain id\n";
1797 : else
1798 8 : oss << _current_subdomain_id << "\n";
1799 16 : oss << " _current_elem_type = " << Utility::enum_to_string(_current_elem_type) << "\n";
1800 8 : oss << " _current_elem_n_sides = " << _current_elem_n_sides << "\n";
1801 8 : oss << " _incoming_point = ";
1802 : if (_incoming_point == RayTracingCommon::invalid_point)
1803 0 : oss << "invalid point\n";
1804 : else
1805 8 : oss << _incoming_point << "\n";
1806 8 : oss << " _incoming_side = ";
1807 8 : if (_incoming_side == RayTracingCommon::invalid_side)
1808 0 : oss << "invalid side\n";
1809 : else
1810 8 : oss << _incoming_side << "\n";
1811 8 : oss << " _intersection_point = ";
1812 : if (_intersection_point == RayTracingCommon::invalid_point)
1813 0 : oss << "invalid point\n";
1814 : else
1815 8 : oss << _intersection_point << "\n";
1816 8 : oss << " _intersected_side = ";
1817 8 : if (_intersected_side == RayTracingCommon::invalid_side)
1818 0 : oss << "invalid side\n";
1819 : else
1820 8 : oss << _intersected_side << "\n";
1821 8 : oss << " _intersected_extrema = " << _intersected_extrema << "\n";
1822 8 : oss << " _exits_elem = " << _exits_elem << "\n";
1823 8 : if (_current_elem)
1824 16 : oss << _current_elem->get_info();
1825 : else
1826 0 : oss << "_current_elem = invalid\n";
1827 :
1828 : possiblySaveDebugRayMesh();
1829 8 : return oss.str();
1830 8 : }
1831 :
1832 : void
1833 8 : TraceRay::failTrace(const std::string & reason, const bool warning, const int line)
1834 : {
1835 8 : const auto message = failTraceMessage(reason, line);
1836 :
1837 8 : if (warning)
1838 : {
1839 4 : ++_results[FAILED_TRACES];
1840 : mooseWarning(message);
1841 4 : (*_current_ray)->setShouldContinue(false);
1842 4 : _should_continue = false;
1843 : }
1844 : else
1845 4 : mooseError(message);
1846 4 : }
1847 :
1848 : const std::vector<NeighborInfo> &
1849 933020 : TraceRay::getVertexNeighbors(const Elem * elem, const Node * vertex)
1850 : {
1851 : traceAssert(elem, "Elem must be valid");
1852 : traceAssert(vertex, "Vertex must be valid");
1853 :
1854 : debugRay("Called getVertexNeighbors() with:");
1855 : debugRay(" elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
1856 : debugRay(" vertex->id() = ", vertex->id(), ", at ", (Point)*vertex);
1857 :
1858 : traceAssert(elem->get_node_index(vertex) != libMesh::invalid_uint, "Doesn't contain node");
1859 : traceAssert(elem->is_vertex(elem->get_node_index(vertex)), "Node is not a vertex");
1860 :
1861 933020 : ++_results[VERTEX_NEIGHBOR_LOOKUPS];
1862 :
1863 : // Return the entry if we have it cached
1864 : auto search = _vertex_neighbors.find(vertex);
1865 933020 : if (search != _vertex_neighbors.end())
1866 866504 : return search->second;
1867 :
1868 66516 : ++_results[VERTEX_NEIGHBOR_BUILDS];
1869 :
1870 : // Make a new entry
1871 : debugRay("Building vertex neighbors");
1872 : std::vector<NeighborInfo> & entry =
1873 66516 : _vertex_neighbors.emplace(vertex, std::vector<NeighborInfo>()).first->second;
1874 :
1875 66516 : findNodeNeighbors(elem,
1876 : vertex,
1877 66516 : _neighbor_set,
1878 66516 : _neighbor_untested_set,
1879 66516 : _neighbor_next_untested_set,
1880 66516 : _neighbor_active_neighbor_children,
1881 : entry);
1882 :
1883 : // Fill the side normals
1884 505463 : for (auto & neighbor_info : entry)
1885 1759784 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1886 1320837 : neighbor_info._side_normals[i] =
1887 1320837 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
1888 :
1889 : return entry;
1890 : }
1891 :
1892 : const std::vector<NeighborInfo> &
1893 933020 : TraceRay::getVertexNeighbors(const Elem * elem, const unsigned short vertex)
1894 : {
1895 : traceAssert(vertex < elem->n_vertices(), "Invalid vertex");
1896 :
1897 933020 : return getVertexNeighbors(elem, elem->node_ptr(vertex));
1898 : }
1899 :
1900 : const std::vector<NeighborInfo> &
1901 2884395 : TraceRay::getEdgeNeighbors(const Elem * elem,
1902 : const std::pair<const Node *, const Node *> & vertices,
1903 : const Point & point)
1904 : {
1905 : traceAssert(elem, "Invalid elem");
1906 : traceAssert(vertices.first, "Must be valid");
1907 : traceAssert(vertices.second, "Must be valid");
1908 :
1909 : debugRay("Called getEdgeNeighbors() with:");
1910 : debugRay(" elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
1911 : debugRay(" vertices.first = ", vertices.first->id(), " at ", (Point)*vertices.first);
1912 : debugRay(" vertices.second = ", vertices.second->id(), " at ", (Point)*vertices.second);
1913 : debugRay(" point = ", point);
1914 :
1915 : traceAssert(elem->get_node_index(vertices.first) != libMesh::invalid_uint,
1916 : "Doesn't contain vertex");
1917 : traceAssert(elem->get_node_index(vertices.second) != libMesh::invalid_uint,
1918 : "Doesn't contain vertex");
1919 : traceAssert(isWithinSegment(
1920 : (Point)*vertices.first, (Point)*vertices.second, point, LOOSE_TRACE_TOLERANCE),
1921 : "Point not within edge");
1922 :
1923 2884395 : ++_results[EDGE_NEIGHBOR_LOOKUPS];
1924 :
1925 2884395 : const auto ordered_vertices = vertices.first->id() < vertices.second->id()
1926 2884395 : ? vertices
1927 2884395 : : std::make_pair(vertices.second, vertices.first);
1928 :
1929 : // Look for the entry and build if necessary
1930 : std::pair<bool, std::vector<NeighborInfo>> * entry;
1931 : auto search = _edge_neighbors.find(ordered_vertices);
1932 2884395 : if (search != _edge_neighbors.end())
1933 2752570 : entry = &search->second;
1934 : else
1935 : {
1936 : debugRay("Building edge neighbors");
1937 131825 : ++_results[EDGE_NEIGHBOR_BUILDS];
1938 131825 : entry = &_edge_neighbors
1939 131825 : .emplace(ordered_vertices, std::make_pair(true, std::vector<NeighborInfo>()))
1940 : .first->second;
1941 131825 : findEdgeNeighbors(elem,
1942 131825 : ordered_vertices.first,
1943 131825 : ordered_vertices.second,
1944 131825 : _neighbor_set,
1945 131825 : _neighbor_untested_set,
1946 131825 : _neighbor_next_untested_set,
1947 131825 : _neighbor_active_neighbor_children,
1948 131825 : entry->second);
1949 :
1950 : bool all_same_edge = true;
1951 646538 : for (auto & neighbor_info : entry->second)
1952 : {
1953 : traceAssert(neighbor_info._lower_bound <= neighbor_info._upper_bound,
1954 : "Bound order incorrect");
1955 :
1956 : // Fill the side normals
1957 1544095 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1958 1029382 : neighbor_info._side_normals[i] =
1959 1029382 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
1960 :
1961 : // See if the bounds are the same as the target edge
1962 514713 : if (neighbor_info._lower_bound != 0 || neighbor_info._upper_bound != 1)
1963 : all_same_edge = false;
1964 : }
1965 131825 : entry->first = all_same_edge;
1966 : }
1967 :
1968 : // Means that all neighbors are not on the exact same edge, so we must
1969 : // validate/invalidate based on if the neighbor's edge contains our point
1970 2884395 : if (!entry->first)
1971 : {
1972 : const auto edge_length =
1973 596 : ((Point)*ordered_vertices.first - (Point)*ordered_vertices.second).norm();
1974 596 : const auto point_location = ((Point)*ordered_vertices.first - point).norm() / edge_length;
1975 4282 : for (auto & info : entry->second)
1976 7372 : info._valid = (info._lower_bound - TRACE_TOLERANCE) < point_location &&
1977 3074 : point_location < (info._upper_bound + TRACE_TOLERANCE);
1978 : }
1979 :
1980 2884395 : return entry->second;
1981 : }
1982 :
1983 : const std::vector<NeighborInfo> &
1984 2884395 : TraceRay::getEdgeNeighbors(const Elem * elem,
1985 : const std::pair<unsigned short, unsigned short> & vertices,
1986 : const Point & point)
1987 : {
1988 : debugRay("Called getEdgeNeighbors(), local index version with:");
1989 : debugRay(" vertices.first = ", vertices.first);
1990 : debugRay(" vertices.second = ", vertices.second);
1991 : traceAssert(vertices.first < elem->n_vertices(),
1992 : "Invalid vertex with ray " + std::to_string((*_current_ray)->id()));
1993 : traceAssert(vertices.second < elem->n_vertices(), "Invalid vertex");
1994 :
1995 2884395 : return getEdgeNeighbors(
1996 2884395 : elem, std::make_pair(elem->node_ptr(vertices.first), elem->node_ptr(vertices.second)), point);
1997 : }
1998 :
1999 : const std::vector<NeighborInfo> &
2000 3793814 : TraceRay::getNeighbors(const Elem * elem, const ElemExtrema & extrema, const Point & point)
2001 : {
2002 3793814 : if (!extrema.atExtrema())
2003 0 : return getPointNeighbors(elem, point);
2004 : if (extrema.atVertex())
2005 909419 : return getVertexNeighbors(elem, extrema.vertex());
2006 2884395 : return getEdgeNeighbors(elem, extrema.edgeVertices(), point);
2007 : }
2008 :
2009 : const std::vector<NeighborInfo> &
2010 22 : TraceRay::getPointNeighbors(const Elem * elem, const Point & point)
2011 : {
2012 : traceAssert(elem, "Invalid elem");
2013 :
2014 : debugRay("Called getPointNeighbors()");
2015 : debugRay(" elem = ", elem->id());
2016 : debugRay(" point = ", point);
2017 :
2018 22 : ++_results[POINT_NEIGHBOR_BUILDS];
2019 22 : _point_neighbor_helper.clear();
2020 :
2021 22 : findPointNeighbors(elem,
2022 : point,
2023 22 : _neighbor_set,
2024 22 : _neighbor_untested_set,
2025 22 : _neighbor_next_untested_set,
2026 22 : _neighbor_active_neighbor_children,
2027 : _point_neighbor_helper);
2028 :
2029 : // Fill the side normals
2030 44 : for (auto & neighbor_info : _point_neighbor_helper)
2031 44 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
2032 22 : neighbor_info._side_normals[i] =
2033 22 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
2034 :
2035 22 : return _point_neighbor_helper;
2036 : }
2037 :
2038 : void
2039 23616445 : TraceRay::storeExitsElemResult(const TraceRay::ExitsElemResult result)
2040 : {
2041 23616445 : if (result == HIT_FACE)
2042 20350947 : ++_results[FACE_HITS];
2043 3265498 : else if (result == HIT_VERTEX)
2044 812794 : ++_results[VERTEX_HITS];
2045 2452704 : else if (result == HIT_EDGE)
2046 2452704 : ++_results[EDGE_HITS];
2047 : else
2048 0 : mooseError("Should not call storeExitsElemResult() with result ", result);
2049 23616445 : }
2050 :
2051 : void
2052 22560690 : TraceRay::onSegment(const std::shared_ptr<Ray> & ray)
2053 : {
2054 : traceAssert((*_current_ray)->currentElem() == _current_elem, "Ray currentElem() incorrect");
2055 : traceAssert((*_current_ray)->currentPoint() == _intersection_point,
2056 : "Ray currentPoint() incorrect");
2057 : traceAssert((*_current_ray)->currentIncomingSide() == _incoming_side,
2058 : "Ray currentIncomingSide() incorrect");
2059 : #ifndef NDEBUG
2060 : if (_study.verifyTraceIntersections())
2061 : {
2062 : if (_current_elem->has_affine_map())
2063 : traceAssert(_current_elem->contains_point(_incoming_point),
2064 : "_current_elem does not contain incoming point");
2065 :
2066 : if (_intersected_side != RayTracingCommon::invalid_side &&
2067 : !_study.sideIsNonPlanar(_current_elem, _intersected_side))
2068 : {
2069 : traceAssert(_elem_side_builder(*_current_elem, _intersected_side)
2070 : .close_to_point(_intersection_point, LOOSE_TRACE_TOLERANCE),
2071 : "Intersected point is not on intersected side");
2072 : traceAssert(!_study.sideIsIncoming(
2073 : _current_elem, _intersected_side, (*_current_ray)->direction(), _tid),
2074 : "Intersected side is not outgoing");
2075 : }
2076 : if (_incoming_side != RayTracingCommon::invalid_side &&
2077 : !_study.sideIsNonPlanar(_current_elem, _incoming_side))
2078 : {
2079 : traceAssert(_elem_side_builder(*_current_elem, _incoming_side)
2080 : .close_to_point(_incoming_point, LOOSE_TRACE_TOLERANCE),
2081 : "Incoming point is not on incoming side");
2082 : if (ray->intersections() != 0 && ray->maxDistance() != 0)
2083 : traceAssert(_study.sideIsIncoming(
2084 : _current_elem, _incoming_side, (*_current_ray)->direction(), _tid),
2085 : "Incoming side is not incoming");
2086 : }
2087 : }
2088 : #endif
2089 : traceAssert(MooseUtils::absoluteFuzzyEqual(_intersection_distance,
2090 : (_intersection_point - _incoming_point).norm()),
2091 : "_intersection_distance is incorrect");
2092 : traceAssert(_current_subdomain_id == _current_elem->subdomain_id(), "Subdomain incorrect");
2093 : traceAssert(MooseUtils::absoluteFuzzyEqual((_incoming_point - _intersection_point).norm(),
2094 : _intersection_distance),
2095 : "Invalid intersection distance");
2096 :
2097 22560690 : _study.reinitSegment(
2098 22560690 : _current_elem, _incoming_point, _intersection_point, _intersection_distance, _tid);
2099 :
2100 22560690 : const auto & rks = _study.currentRayKernels(_tid);
2101 45130632 : for (auto & rk : rks)
2102 : {
2103 22569954 : rk->onSegment();
2104 22569942 : postRayTracingObject(ray, rk);
2105 : }
2106 22560678 : }
2107 :
2108 : void
2109 973271 : TraceRay::onBoundary(const std::shared_ptr<Ray> & ray, const bool external)
2110 : {
2111 : traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point),
2112 : "Ray currentPoint() not set before onBoundary()");
2113 :
2114 : // Get the RayBCs on bnd_elems
2115 973271 : _study.getRayBCs(_on_boundary_ray_bcs, _boundary_elems, _tid, ray->id());
2116 :
2117 : // Store this information temprorarily because we are going to change it as we
2118 : // apply each boundary condition
2119 973271 : const auto old_current_elem = _current_elem;
2120 973271 : const auto old_intersected_side = _intersected_side;
2121 973271 : const auto old_intersected_extrema = _intersected_extrema;
2122 973271 : const auto old_subdomain_id = _current_subdomain_id;
2123 :
2124 : // For each RayBC we found, apply it on the boundaries that we need
2125 1962950 : for (RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2126 : {
2127 : // First, find the boundaries this RayBC is valid on that are also in _boundary_elems.
2128 : // We do this ahead of time so that we can pass in to RayBC::apply if the same
2129 : // boundary condition is being appled multiple times on different boundaries.
2130 : // This is useful in situations like reflection where multiple reflections are
2131 : // necessary at a corner to perfectly reflect.
2132 989689 : _on_boundary_apply_index.clear();
2133 2053953 : for (MooseIndex(_boundary_elems.size()) bnd_elems_i = 0; bnd_elems_i < _boundary_elems.size();
2134 : ++bnd_elems_i)
2135 1064264 : if (rbc->hasBoundary(_boundary_elems[bnd_elems_i].bnd_id))
2136 1064214 : _on_boundary_apply_index.push_back(bnd_elems_i);
2137 :
2138 : traceAssert(!_on_boundary_apply_index.empty(), "Must not be empty");
2139 :
2140 : // Apply the RayBC on each of the relevant boundary elements
2141 2053893 : for (const auto bnd_elems_i : _on_boundary_apply_index)
2142 : {
2143 : auto & bnd_elem = _boundary_elems[bnd_elems_i];
2144 :
2145 : debugRay("Calling ",
2146 : rbc->type(),
2147 : "::onBoundary for \"",
2148 : rbc->name(),
2149 : "\" on elem ",
2150 : bnd_elem.elem->id(),
2151 : " and side ",
2152 : bnd_elem.side,
2153 : " for bnd_id ",
2154 : bnd_elem.bnd_id);
2155 :
2156 1064214 : _current_elem = bnd_elem.elem;
2157 1064214 : _current_bnd_id = bnd_elem.bnd_id;
2158 1064214 : _intersected_side = bnd_elem.side;
2159 : _intersected_extrema = bnd_elem.extrema;
2160 1064214 : _current_subdomain_id = _current_elem->subdomain_id();
2161 :
2162 1064214 : rbc->onBoundary(_on_boundary_apply_index.size());
2163 1064204 : postRayTracingObject(ray, rbc);
2164 : }
2165 : }
2166 :
2167 : // Set this info back now that we're done applying BCs
2168 973261 : _current_elem = old_current_elem;
2169 973261 : _intersected_side = old_intersected_side;
2170 : _intersected_extrema = old_intersected_extrema;
2171 973261 : _current_bnd_id = BoundaryInfo::invalid_id;
2172 973261 : _current_subdomain_id = old_subdomain_id;
2173 :
2174 : // When on an external boundary, the Ray must have been changed or killed.
2175 : // Otherwise, we don't know what to do with it now! If this didn't happen,
2176 : // output a detailed error message.
2177 973261 : if (external && !ray->trajectoryChanged() && ray->shouldContinue())
2178 : {
2179 8 : std::stringstream oss;
2180 8 : oss << "Don't know what to do with a Ray after it hit an external\n";
2181 8 : oss << "boundary at point " << _intersection_point << "!\n\n";
2182 8 : oss << "When hitting an external RayBC, a Ray must either:\n";
2183 8 : oss << " Be killed by a RayBC\n";
2184 8 : oss << " Have its trajectory changed by the RayBC\n";
2185 8 : oss << "by at least one of the executed RayBCs.\n\n";
2186 8 : oss << "You need to either:\n";
2187 8 : oss << " Kill/change the Ray sooner with RayKernels, internal RayBCs, or a max distance\n";
2188 8 : oss << " Kill/change the Ray on the boundary with a RayBC\n\n";
2189 8 : if (!_on_boundary_ray_bcs.empty())
2190 : {
2191 6 : oss << "RayBCs executed that did not kill or change the Ray:\n";
2192 12 : for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2193 14 : for (const auto & bnd_elem : _boundary_elems)
2194 8 : if (rbc->hasBoundary(bnd_elem.bnd_id))
2195 18 : oss << " " << rbc->typeAndName() << " on boundary " << bnd_elem.bnd_id << " ("
2196 12 : << _mesh.getBoundaryName(bnd_elem.bnd_id) << ")\n";
2197 6 : oss << "\n";
2198 : }
2199 : bool output_header = false;
2200 20 : for (std::size_t i = 0; i < _boundary_elems.size(); ++i)
2201 : {
2202 12 : const auto bnd_id = _boundary_elems[i].bnd_id;
2203 : bool found = false;
2204 14 : for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2205 8 : if (rbc->hasBoundary(bnd_id))
2206 : {
2207 : found = true;
2208 : break;
2209 : }
2210 :
2211 12 : if (!found)
2212 : {
2213 6 : if (!output_header)
2214 : {
2215 4 : oss << "Boundaries that did not have any RayBCs:\n";
2216 : output_header = true;
2217 : }
2218 12 : oss << " " << bnd_id << " (" << _mesh.getBoundaryName(bnd_id) << ")\n";
2219 : }
2220 : }
2221 :
2222 8 : failTrace(oss.str(), _study.tolerateFailure(), __LINE__);
2223 4 : }
2224 973257 : }
2225 :
2226 : Real
2227 31429648 : TraceRay::subdomainHmax(const Elem * elem) const
2228 : {
2229 : const auto subdomain_id = elem->subdomain_id();
2230 31429648 : return subdomain_id == _current_subdomain_id ? _current_subdomain_hmax
2231 1449 : : _study.subdomainHmax(subdomain_id);
2232 : }
2233 :
2234 : void
2235 23634146 : TraceRay::postRayTracingObject(const std::shared_ptr<Ray> & ray, const RayTracingObject * rto)
2236 : {
2237 23634146 : if (!ray->shouldContinue())
2238 : {
2239 2719014 : if (_should_continue)
2240 938697 : _should_continue = false;
2241 : }
2242 20915132 : else if (!_should_continue)
2243 0 : failTrace(rto->typeAndName() +
2244 0 : " set a Ray to continue that was previously set to not continue.\n\n" +
2245 : "Once a Ray has been set to not continue, its continue status cannot change.",
2246 : /* warning = */ false,
2247 : __LINE__);
2248 :
2249 23634146 : if (!_should_continue && ray->trajectoryChanged())
2250 0 : failTrace(rto->typeAndName() +
2251 0 : " changed the trajectory of a Ray that was set to not continue,\n" +
2252 : "or set a Ray whose trajectory was changed to not continue.",
2253 : /* warning = */ false,
2254 : __LINE__);
2255 23634146 : }
|
