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 4551 : TraceRay::TraceRay(RayTracingStudy & study, const THREAD_ID tid)
48 4551 : : _study(study),
49 4551 : _mesh(study.getSubProblem().mesh()),
50 4551 : _dim(_mesh.dimension()),
51 4551 : _boundary_info(_mesh.getMesh().get_boundary_info()),
52 4551 : _pid(_study.comm().rank()),
53 4551 : _tid(tid),
54 4551 : _backface_culling(false),
55 4551 : _current_normals(nullptr),
56 9102 : _results(ENDED_STATIONARY + 1)
57 : {
58 4551 : }
59 :
60 : void
61 10394 : TraceRay::preExecute()
62 : {
63 10394 : _current_subdomain_id = Elem::invalid_subdomain_id;
64 10394 : _current_elem_type = INVALID_ELEM;
65 :
66 : // Zero out all results
67 166304 : for (auto & val : _results)
68 155910 : val = 0;
69 :
70 10394 : _has_ray_kernels = _study.hasRayKernels(_tid);
71 10394 : _is_rectangular_domain = _study.isRectangularDomain();
72 10394 : }
73 :
74 : void
75 278 : TraceRay::meshChanged()
76 : {
77 : // Invalidate the vertex and edge neighbor caches
78 : _vertex_neighbors.clear();
79 : _edge_neighbors.clear();
80 278 : }
81 :
82 : TraceRay::ExitsElemResult
83 46867642 : 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 46867642 : switch (elem_type)
107 : {
108 5551245 : case HEX8:
109 5551245 : intersected = exitsElem<Hex8, Hex8>(elem,
110 : incoming_side,
111 : intersection_point,
112 : intersected_side,
113 : intersected_extrema,
114 : intersection_distance,
115 : normals);
116 5551245 : break;
117 7781387 : case TET4:
118 7781387 : intersected = exitsElem<Tet4, Tet4>(elem,
119 : incoming_side,
120 : intersection_point,
121 : intersected_side,
122 : intersected_extrema,
123 : intersection_distance,
124 : normals);
125 7781387 : break;
126 2531714 : case PYRAMID5:
127 2531714 : intersected = exitsElem<Pyramid5, Pyramid5>(elem,
128 : incoming_side,
129 : intersection_point,
130 : intersected_side,
131 : intersected_extrema,
132 : intersection_distance,
133 : normals);
134 2531714 : break;
135 1921082 : case PRISM6:
136 1921082 : intersected = exitsElem<Prism6, Prism6>(elem,
137 : incoming_side,
138 : intersection_point,
139 : intersected_side,
140 : intersected_extrema,
141 : intersection_distance,
142 : normals);
143 1921082 : break;
144 830453 : case QUAD4:
145 830453 : intersected = exitsElem<Quad4, Quad4>(elem,
146 : incoming_side,
147 : intersection_point,
148 : intersected_side,
149 : intersected_extrema,
150 : intersection_distance,
151 : normals);
152 830453 : break;
153 54055 : case TRI3:
154 54055 : intersected = exitsElem<Tri3, Tri3>(elem,
155 : incoming_side,
156 : intersection_point,
157 : intersected_side,
158 : intersected_extrema,
159 : intersection_distance,
160 : normals);
161 54055 : break;
162 1796521 : case HEX20:
163 1796521 : intersected = exitsElem<Hex20, Hex8>(elem,
164 : incoming_side,
165 : intersection_point,
166 : intersected_side,
167 : intersected_extrema,
168 : intersection_distance,
169 : normals);
170 1796521 : break;
171 1796510 : case HEX27:
172 1796510 : intersected = exitsElem<Hex27, Hex8>(elem,
173 : incoming_side,
174 : intersection_point,
175 : intersected_side,
176 : intersected_extrema,
177 : intersection_distance,
178 : normals);
179 1796510 : break;
180 29306 : case QUAD8:
181 29306 : intersected = exitsElem<Quad8, Quad4>(elem,
182 : incoming_side,
183 : intersection_point,
184 : intersected_side,
185 : intersected_extrema,
186 : intersection_distance,
187 : normals);
188 29306 : break;
189 29306 : case QUAD9:
190 29306 : intersected = exitsElem<Quad9, Quad4>(elem,
191 : incoming_side,
192 : intersection_point,
193 : intersected_side,
194 : intersected_extrema,
195 : intersection_distance,
196 : normals);
197 29306 : break;
198 44329 : case TRI6:
199 44329 : intersected = exitsElem<Tri6, Tri3>(elem,
200 : incoming_side,
201 : intersection_point,
202 : intersected_side,
203 : intersected_extrema,
204 : intersection_distance,
205 : normals);
206 44329 : break;
207 44329 : case TRI7:
208 44329 : intersected = exitsElem<Tri7, Tri3>(elem,
209 : incoming_side,
210 : intersection_point,
211 : intersected_side,
212 : intersected_extrema,
213 : intersection_distance,
214 : normals);
215 44329 : break;
216 7773339 : case TET10:
217 7773339 : intersected = exitsElem<Tet10, Tet4>(elem,
218 : incoming_side,
219 : intersection_point,
220 : intersected_side,
221 : intersected_extrema,
222 : intersection_distance,
223 : normals);
224 7773339 : break;
225 7773342 : case TET14:
226 7773342 : intersected = exitsElem<Tet14, Tet4>(elem,
227 : incoming_side,
228 : intersection_point,
229 : intersected_side,
230 : intersected_extrema,
231 : intersection_distance,
232 : normals);
233 7773342 : break;
234 2531456 : case PYRAMID13:
235 2531456 : intersected = exitsElem<Pyramid13, Pyramid5>(elem,
236 : incoming_side,
237 : intersection_point,
238 : intersected_side,
239 : intersected_extrema,
240 : intersection_distance,
241 : normals);
242 2531456 : break;
243 2531456 : case PYRAMID14:
244 2531456 : intersected = exitsElem<Pyramid14, Pyramid5>(elem,
245 : incoming_side,
246 : intersection_point,
247 : intersected_side,
248 : intersected_extrema,
249 : intersection_distance,
250 : normals);
251 2531456 : break;
252 1921074 : case PRISM15:
253 1921074 : intersected = exitsElem<Prism15, Prism6>(elem,
254 : incoming_side,
255 : intersection_point,
256 : intersected_side,
257 : intersected_extrema,
258 : intersection_distance,
259 : normals);
260 1921074 : break;
261 1921105 : case PRISM18:
262 1921105 : intersected = exitsElem<Prism18, Prism6>(elem,
263 : incoming_side,
264 : intersection_point,
265 : intersected_side,
266 : intersected_extrema,
267 : intersection_distance,
268 : normals);
269 1921105 : break;
270 4997 : case EDGE2:
271 4997 : intersected = exitsElem<Edge2, Edge2>(elem,
272 : incoming_side,
273 : intersection_point,
274 : intersected_side,
275 : intersected_extrema,
276 : intersection_distance,
277 : normals);
278 4997 : break;
279 318 : case EDGE3:
280 318 : intersected = exitsElem<Edge3, Edge2>(elem,
281 : incoming_side,
282 : intersection_point,
283 : intersected_side,
284 : intersected_extrema,
285 : intersection_distance,
286 : normals);
287 318 : break;
288 318 : case EDGE4:
289 318 : intersected = exitsElem<Edge4, Edge2>(elem,
290 : incoming_side,
291 : intersection_point,
292 : intersected_side,
293 : intersected_extrema,
294 : intersection_distance,
295 : normals);
296 318 : break;
297 0 : default:
298 0 : mooseError(
299 0 : "Element type ", Utility::enum_to_string(elem->type()), " not supported by TraceRay");
300 : }
301 :
302 46867642 : if (intersected)
303 : {
304 37267646 : 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 46862009 : 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 46862009 : ++_results[INTERSECTION_CALLS];
323 :
324 : debugRay("Called exitsElem() in 2D or 3D");
325 :
326 46862009 : 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 46862009 : 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 46862009 : 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 46862009 : 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 221942212 : 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 221835506 : if (s == incoming_side)
369 : {
370 : debugRay(" Skipping due to incoming side");
371 44407307 : if (++s == T::num_sides)
372 : break;
373 : else
374 32743811 : continue;
375 : }
376 :
377 : // Using backface culling on this run through the sides
378 : // If the direction is non-entrant, skip this side
379 177428199 : if (use_backface_culling)
380 : {
381 : // Side is non-entrant per the culling, so skip
382 9214908 : if (normals[s] * direction < -LOOSE_TRACE_TOLERANCE)
383 : {
384 : debugRay(" Skipping due to backface culling dot = ", normals[s] * direction);
385 :
386 2973349 : if (++s == T::num_sides)
387 : break;
388 : else
389 2422158 : 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 168213291 : else if (_backface_culling)
399 : {
400 : // Side was non-entrant per the culling, so try again
401 2401920 : if (normals[s] * direction >= -LOOSE_TRACE_TOLERANCE)
402 : {
403 : debugRay(" Skipping because we already checked this side with culling enabled");
404 1568016 : if (++s == T::num_sides)
405 : break;
406 : else
407 1283281 : continue;
408 : }
409 : else
410 : debugRay(" Side that was skipped due to culling");
411 :
412 833904 : ++_results[BACKFACE_CULLING_FAILURES];
413 : }
414 : }
415 :
416 : // Look for an intersection!
417 172993540 : current_intersection_point = RayTracingCommon::invalid_point;
418 : current_intersected_extrema.invalidate();
419 172993540 : const bool intersected = sideIntersectedByLine<FirstOrderT>(elem,
420 172993540 : _incoming_point,
421 : direction,
422 : s,
423 172993540 : _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 172993540 : 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 40772661 : if (current_intersection_distance > best_intersection_distance)
456 : {
457 : debugRay(" Best intersection so far");
458 :
459 37909849 : intersected_side = s;
460 37909849 : intersection_distance = current_intersection_distance;
461 37909849 : intersection_point = current_intersection_point;
462 : intersected_extrema = current_intersected_extrema;
463 : best_intersection_distance = current_intersection_distance;
464 : }
465 : }
466 :
467 172993540 : if (++s == T::num_sides || try_nonplanar_incoming_side)
468 : break;
469 : else
470 138042423 : continue;
471 : } // while(true)
472 :
473 : // Found an intersection
474 47450539 : 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 10188508 : 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 10081802 : 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 481824 : 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 9599978 : if (_incoming_side != RayTracingCommon::invalid_side &&
502 9595583 : _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 106706 : 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 5633 : 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 5633 : ++_results[INTERSECTION_CALLS];
526 :
527 : debugRay("Called exitsElem() in 1D");
528 :
529 : // Scale the tolerance based on the element size
530 5633 : 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 5633 : if (incoming_side != RayTracingCommon::invalid_side)
535 : {
536 5157 : intersected_side = (incoming_side == 1 ? 0 : 1);
537 : intersected_extrema.setVertex(intersected_side);
538 5157 : intersection_point = elem->point(intersected_side);
539 5157 : 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 5157 : return true;
545 : }
546 :
547 : // End point that is for sure out of the element
548 : const Point extended_end_point =
549 476 : _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 784 : for (MooseIndex(elem->n_sides()) side = 0; side < elem->n_sides(); ++side)
554 : {
555 766 : const Point side_point = elem->point(side);
556 : debugRay(" Checking side ", side, " at ", side_point);
557 :
558 766 : const Real incoming_to_side = (side_point - _incoming_point).norm();
559 766 : if (incoming_to_side < tol)
560 : {
561 : debugRay(" Continuing because at side");
562 200 : continue;
563 : }
564 :
565 566 : const Real incoming_to_end = (extended_end_point - _incoming_point).norm();
566 566 : const Real side_to_end = (extended_end_point - side_point).norm();
567 566 : const Real sum = incoming_to_side + side_to_end - incoming_to_end;
568 : debugRay(" Sum = ", sum);
569 :
570 566 : if (std::abs(sum) < tol)
571 : {
572 458 : intersected_side = side;
573 : intersected_extrema.setVertex(side);
574 458 : intersection_point = side_point;
575 458 : intersection_distance = incoming_to_side;
576 : debugRay(" Intersected at ", intersection_point);
577 458 : return true;
578 : }
579 : }
580 :
581 : return false;
582 : }
583 :
584 : TraceRay::ExitsElemResult
585 3803893 : 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 3803893 : ++_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 22729628 : 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 18925735 : if (neighbor_info._elem == last_elem)
625 : {
626 : debugRay("Skipping last elem ", last_elem->id());
627 : last_elem_info = &neighbor_info;
628 3806736 : continue;
629 : }
630 :
631 15118999 : 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 15118999 : if (exit_result != NO_EXIT)
640 : {
641 : debugRay("Ray can exit through neighbor ", neighbor_info._elem->id());
642 :
643 5530808 : if (current_intersection_distance > longest_distance)
644 : {
645 3914385 : best_elem = neighbor_info._elem;
646 3914385 : best_elem_incoming_side = current_incoming_side;
647 3914385 : _intersection_point = current_intersection_point;
648 3914385 : _intersected_side = current_intersected_side;
649 : _intersected_extrema = current_intersected_extrema;
650 3914385 : _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 3803893 : if (!best_elem && last_elem_info)
659 : {
660 414 : 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 414 : if (exit_result != NO_EXIT && current_intersection_distance > longest_distance)
668 : {
669 : debugRay("Ray can exit through last_elem ", last_elem->id());
670 :
671 414 : best_elem = last_elem;
672 414 : best_elem_incoming_side = current_incoming_side;
673 414 : _intersection_point = current_intersection_point;
674 414 : _intersected_side = current_intersected_side;
675 : _intersected_extrema = current_intersected_extrema;
676 414 : _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 3803893 : return best_exit_result;
695 : }
696 :
697 : TraceRay::ExitsElemResult
698 15119413 : 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 15119413 : if (!neighbor_info._valid)
706 : return NO_EXIT;
707 :
708 15118121 : 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 15118121 : incoming_side = RayTracingCommon::invalid_side;
713 19929218 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
714 19925777 : if (neighbor_info._side_normals[i] * (*_current_ray)->direction() < LOOSE_TRACE_TOLERANCE)
715 : {
716 15114680 : incoming_side = neighbor_info._sides[i];
717 15114680 : break;
718 : }
719 :
720 : // No entrant sides on this neighbor were found
721 15118121 : if (incoming_side == RayTracingCommon::invalid_side)
722 : return NO_EXIT;
723 :
724 15114680 : intersection_point = RayTracingCommon::invalid_point;
725 15114680 : intersected_side = RayTracingCommon::invalid_side;
726 : intersected_extrema.invalidate();
727 15114680 : 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 15114680 : exitsElem(neighbor,
733 15114680 : neighbor->type(),
734 15114680 : incoming_side,
735 : intersection_point,
736 : intersected_side,
737 : intersected_extrema,
738 : intersection_distance,
739 15114680 : _backface_culling ? _study.getElemNormals(neighbor, _tid) : nullptr);
740 : debugRay("Done with exitsElem() from moveThroughNeighbor()");
741 :
742 15114680 : return exit_result;
743 : }
744 :
745 : void
746 1410360 : 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 1410360 : _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 1410360 : if (_dim != 1 && _intersected_extrema.atExtrema())
760 : {
761 836807 : 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 3206576 : for (const auto & neighbor_info : neighbors)
770 : {
771 2369769 : if (!neighbor_info._valid)
772 0 : continue;
773 :
774 2369769 : const Elem * elem = neighbor_info._elem;
775 : const auto & sides = neighbor_info._sides;
776 : const auto & side_normals = neighbor_info._side_normals;
777 :
778 8415173 : for (MooseIndex(side_normals.size()) i = 0; i < side_normals.size(); ++i)
779 6045404 : if (!elem->neighbor_ptr(sides[i]) // is a boundary side that has our point
780 6045404 : && side_normals[i] * ray->direction() > TRACE_TOLERANCE) // and is entrant
781 : {
782 : // TODO: this could likely be optimized
783 : ElemExtrema extrema;
784 1983839 : withinExtremaOnSide(elem, _intersection_point, sides[i], _dim, extrema);
785 :
786 1983839 : _boundary_info.boundary_ids(elem, sides[i], _boundary_ids);
787 1983839 : 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 573553 : _boundary_info.boundary_ids(_current_elem, _intersected_side, _boundary_ids);
796 573553 : possiblyAddToBoundaryElems(
797 573553 : _current_elem, _intersected_side, _boundary_ids, _intersected_extrema);
798 : }
799 :
800 : debugRay("Calling external onBoundary() with ", _boundary_elems.size(), " boundaries");
801 1410360 : onBoundary(ray, /* external = */ true);
802 1410350 : }
803 :
804 : void
805 2506 : 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 2506 : _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 2506 : if (_dim != 1 && _intersected_extrema.atExtrema())
828 : {
829 : debugRay("Checking point neighbors for internal sidesets");
830 :
831 : // Get the neighbors
832 1728 : 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 15354 : for (const auto & neighbor_info : neighbors)
841 : {
842 13626 : if (!neighbor_info._valid)
843 0 : continue;
844 :
845 13626 : const Elem * elem = neighbor_info._elem;
846 :
847 : // Grab the internal sidesets for this elem
848 13626 : const auto & sidesets = _study.getInternalSidesets(elem);
849 : // It has none to contribute, so we can continue
850 13626 : if (sidesets.empty())
851 : {
852 : debugRay(" Elem ", elem->id(), " has no internal sidesets");
853 11142 : 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 8946 : for (std::size_t i = 0; i < sides.size(); ++i)
862 : {
863 6462 : const auto side = sides[i];
864 : // Side has internal sidesets and is entrant
865 6462 : 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 1710 : withinExtremaOnSide(elem, _intersection_point, side, _dim, temp_extrema);
870 :
871 1710 : 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 778 : const auto & current_elem_sidesets = _study.getInternalSidesets(_current_elem);
881 778 : 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 272 : if (!_study.hasSameLevelActiveElems())
889 0 : mooseError("Internal sidesets are not currently supported with adaptivity in tracing");
890 :
891 : // Special case for 1D
892 272 : if (_dim == 1)
893 38 : temp_extrema.setVertex(atVertex(_current_elem, _intersection_point));
894 :
895 272 : possiblyAddToBoundaryElems(_current_elem,
896 : _incoming_side,
897 272 : current_elem_sidesets[_incoming_side],
898 272 : _dim != 1 ? _intersected_extrema : temp_extrema);
899 : }
900 :
901 778 : const auto & last_elem_sidesets = _study.getInternalSidesets(_last_elem);
902 778 : if (last_elem_sidesets.size() && last_elem_sidesets[_intersected_side].size())
903 308 : possiblyAddToBoundaryElems(_last_elem,
904 : _intersected_side,
905 : last_elem_sidesets[_intersected_side],
906 308 : _intersected_extrema);
907 : }
908 :
909 2506 : if (!_boundary_elems.empty())
910 : {
911 : debugRay(" Calling internal onBoundary() with ", _boundary_elems.size(), " boundaries");
912 1219 : onBoundary(ray, /* external = */ false);
913 : }
914 2502 : }
915 :
916 : void
917 2559682 : 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 5119364 : for (const auto bnd_id : bnd_ids)
926 : {
927 : bool found = false;
928 2713309 : for (const auto & bnd_elem : _boundary_elems)
929 1208041 : if (bnd_elem.bnd_id == bnd_id)
930 : {
931 : found = true;
932 : break;
933 : }
934 :
935 2559682 : 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 1505268 : _boundary_elems.emplace_back(elem, side, bnd_id, extrema);
945 : }
946 : }
947 2559682 : }
948 :
949 : void
950 1117190 : 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 1117190 : 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 6376084 : for (unsigned short s = 0; s < _current_elem_n_sides; ++s)
971 6784724 : if (!_current_elem->neighbor_ptr(s) && s != _intersected_side &&
972 2489073 : _current_elem->is_node_on_side(_intersected_extrema.first, s) &&
973 7377325 : (!at_edge || _current_elem->is_node_on_side(_intersected_extrema.second, s)) &&
974 1062808 : !_study.sideIsIncoming(_current_elem, s, direction, _tid))
975 : {
976 : debugRay(" Side ", s, " is a boundary side and the Ray exits");
977 172018 : boundary_side = s;
978 : boundary_extrema = _intersected_extrema;
979 172018 : boundary_elem = _current_elem;
980 172018 : 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 945172 : 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 3267297 : 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 2484948 : if (!neighbor_info._valid)
994 0 : continue;
995 :
996 2484948 : const Elem * neighbor = neighbor_info._elem;
997 : // We've already checked ourself
998 2484948 : if (neighbor == _current_elem)
999 851710 : 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 5076268 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1005 3605853 : if (neighbor_info._side_normals[i] * direction > TRACE_TOLERANCE &&
1006 795137 : !neighbor->neighbor_ptr(neighbor_info._sides[i]))
1007 : {
1008 162823 : withinExtremaOnSide(
1009 162823 : neighbor, _intersection_point, neighbor_info._sides[i], _dim, boundary_extrema);
1010 : traceAssert(boundary_extrema.atExtrema(), "Should be at extrema");
1011 162823 : boundary_side = neighbor_info._sides[i];
1012 162823 : boundary_elem = neighbor;
1013 : return;
1014 : }
1015 : }
1016 : }
1017 :
1018 : void
1019 5380196 : TraceRay::trace(const std::shared_ptr<Ray> & ray)
1020 : {
1021 : mooseAssert(_study.currentlyPropagating(), "Should only use while propagating rays");
1022 :
1023 5380196 : _current_ray = &ray;
1024 5380196 : _current_elem = ray->currentElem();
1025 5380196 : _last_elem = nullptr;
1026 5380196 : _incoming_point = ray->currentPoint();
1027 5380196 : _incoming_side = ray->currentIncomingSide();
1028 5380196 : _should_continue = true;
1029 :
1030 5380196 : _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 5380196 : if (_study.verifyRays() && !ray->invalidCurrentIncomingSide() && ray->maxDistance() > 0 &&
1037 10209177 : !_study.sideIsNonPlanar(_current_elem, _incoming_side) &&
1038 2386714 : !_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 immedtiately 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 5380196 : if (_study.shouldCacheTrace(ray))
1064 : {
1065 : debugRay("Trying to init threaded cached trace");
1066 :
1067 13787 : _current_cached_trace = &_study.initThreadedCachedTrace(ray, _tid);
1068 :
1069 : // Add starting data
1070 13787 : if (_study.dataOnCacheTraces())
1071 1371 : _current_cached_trace->lastPoint()._data = ray->data();
1072 13787 : if (_study.auxDataOnCacheTraces())
1073 1325 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1074 : }
1075 : else
1076 5366409 : _current_cached_trace = nullptr;
1077 :
1078 : // Need to call subdomain setup
1079 5380196 : if (_current_elem->subdomain_id() != _current_subdomain_id || _study.rayDependentSubdomainSetup())
1080 5380196 : 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 35541892 : _exits_elem = false;
1110 35541892 : _intersection_point = RayTracingCommon::invalid_point;
1111 35541892 : _intersected_side = RayTracingCommon::invalid_side;
1112 : _intersected_extrema.invalidate();
1113 35541892 : _intersection_distance = RayTracingCommon::invalid_distance;
1114 :
1115 : // Stationary ray
1116 35541892 : if (ray->stationary())
1117 : {
1118 : mooseAssert(ray->invalidDirection(), "Should have an invalid direction");
1119 1575 : _exits_elem = true;
1120 1575 : _intersection_point = _incoming_point;
1121 : _intersected_extrema = _last_intersected_extrema;
1122 1575 : _intersection_distance = 0;
1123 1575 : ++_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 35540317 : 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 31752962 : if (_backface_culling)
1137 955942 : _current_normals = _study.getElemNormals(_current_elem, _tid);
1138 :
1139 31752962 : const auto exits_elem_result = exitsElem(_current_elem,
1140 : _current_elem_type,
1141 31752962 : _incoming_side,
1142 31752962 : _intersection_point,
1143 31752962 : _intersected_side,
1144 31752962 : _intersected_extrema,
1145 31752962 : _intersection_distance,
1146 : _current_normals);
1147 :
1148 31752962 : if (exits_elem_result != NO_EXIT)
1149 : {
1150 31736424 : storeExitsElemResult(exits_elem_result);
1151 31736424 : _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 35541892 : 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 3803893 : if (_last_intersected_extrema.atExtrema())
1178 : {
1179 : traceAssert(_last_elem, "Should be valid");
1180 3787355 : move_through_neighbors_last = _last_elem;
1181 3787355 : 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 16538 : move_through_neighbors_last = _current_elem;
1197 :
1198 : // If we have side info: check for vertices on said side, otherwise, check everywhere
1199 16538 : const auto at_v = _incoming_side != RayTracingCommon::invalid_side
1200 16538 : ? atVertexOnSide(_current_elem, _incoming_point, _incoming_side)
1201 1479 : : atVertex(_current_elem, _incoming_point);
1202 :
1203 16538 : if (at_v != RayTracingCommon::invalid_vertex)
1204 16508 : 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 30 : 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 16508 : if (!neighbors || neighbors->empty())
1221 30 : neighbors = &getPointNeighbors(_current_elem, _incoming_point);
1222 :
1223 : // Couldn't find anything
1224 16538 : if (neighbors->empty())
1225 : {
1226 0 : failTrace("Could not find neighbors to move through", _study.tolerateFailure(), __LINE__);
1227 247985 : return;
1228 : }
1229 : }
1230 :
1231 : // Move through a neighbor
1232 3803893 : const Elem * best_neighbor = nullptr;
1233 3803893 : auto best_neighbor_side = RayTracingCommon::invalid_side;
1234 3803893 : 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 3803893 : 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 3803893 : _exits_elem = true;
1249 3803893 : _last_elem = _current_elem;
1250 3803893 : _current_elem = best_neighbor;
1251 3803893 : _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 3803893 : 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 247985 : _intersection_distance = RayTracingCommon::invalid_distance;
1264 :
1265 247985 : continueTraceOffProcessor(ray);
1266 : return;
1267 : }
1268 :
1269 : // Subdomain changed
1270 3555908 : if (_current_elem->subdomain_id() != _current_subdomain_id)
1271 202 : onSubdomainChanged(ray, /* same_ray = */ true);
1272 :
1273 : // Do own this element - tally the result as we're the ones tracing it
1274 3555908 : 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 35293907 : if (_intersection_distance > 0)
1290 : {
1291 : debugRay("Incrementing ray intersections by 1 to ", ray->intersections() + 1);
1292 : ray->addIntersection();
1293 35292332 : _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 35303087 : 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 35293907 : if (MooseUtils::absoluteFuzzyEqual(ray->distance(), max_distance))
1311 : {
1312 : debugRay("At max distance");
1313 :
1314 : ray->setShouldContinue(false);
1315 20048 : _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 35273859 : 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 2693517 : _intersection_point -= ray->direction() * difference;
1332 2693517 : _intersection_distance -= difference;
1333 2693517 : _intersected_side = RayTracingCommon::invalid_side;
1334 : _intersected_extrema.invalidate();
1335 : ray->setCurrentPoint(_intersection_point);
1336 :
1337 : ray->setShouldContinue(false);
1338 2693517 : _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 35293907 : 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 33831036 : onSegment(ray);
1360 :
1361 33831024 : _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 33831024 : 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 2542122 : onCompleteTrace(ray);
1372 2542122 : return;
1373 : }
1374 :
1375 : // RayKernel moved a Ray
1376 31288902 : 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 648 : _incoming_point = ray->currentPoint();
1385 648 : _incoming_side = RayTracingCommon::invalid_side;
1386 : _intersected_extrema.invalidate();
1387 : ray->setCurrentIncomingSide(_incoming_side);
1388 :
1389 648 : const auto new_intersection_distance = (ray->currentPoint() - _incoming_point).norm();
1390 648 : ray->addDistance(-_intersection_distance + new_intersection_distance);
1391 648 : _intersection_distance = new_intersection_distance;
1392 :
1393 648 : onTrajectoryChanged(ray);
1394 648 : onContinueTrace(ray);
1395 : continue;
1396 648 : }
1397 : else
1398 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1399 : }
1400 : else
1401 : {
1402 1462871 : _study.postOnSegment(_tid, ray);
1403 :
1404 1462871 : if (!_should_continue)
1405 : {
1406 : debugRay("Killing due to at end without RayKernels");
1407 : traceAssert(!ray->shouldContinue(), "Ray shouldn't continue");
1408 :
1409 171441 : onCompleteTrace(ray);
1410 171441 : 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 32574837 : if (_dim > 1 && _intersected_extrema.atExtrema() &&
1420 37198493 : _current_elem->neighbor_ptr(_intersected_side) &&
1421 8233686 : (!_is_rectangular_domain ||
1422 4116843 : onBoundingBoxBoundary(_study.boundingBox(),
1423 4116843 : _intersection_point,
1424 : _dim,
1425 4116843 : LOOSE_TRACE_TOLERANCE * _study.domainMaxLength())))
1426 : {
1427 1117190 : auto boundary_side = RayTracingCommon::invalid_side;
1428 : ElemExtrema boundary_extrema;
1429 1117190 : const Elem * boundary_elem = nullptr;
1430 :
1431 1117190 : findExternalBoundarySide(boundary_side, boundary_extrema, boundary_elem);
1432 :
1433 1117190 : 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 334841 : _last_elem = _current_elem;
1440 334841 : _current_elem = boundary_elem;
1441 334841 : _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 32579684 : _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 32579684 : neighbor = _study.hasSameLevelActiveElems()
1466 32579684 : ? _current_elem->neighbor_ptr(_intersected_side)
1467 189283 : : getActiveNeighbor(_current_elem, _intersected_side, _intersection_point);
1468 :
1469 : // Found one - not on the boundary so set the next info
1470 32579684 : 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 31169324 : 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 30993925 : const unsigned short n_sides = neighbor->subdomain_id() == _current_subdomain_id
1484 : ? _current_elem_n_sides
1485 30993925 : : neighbor->n_sides();
1486 : traceAssert(n_sides == neighbor->n_sides(), "n_sides incorrect");
1487 :
1488 94688113 : for (_incoming_side = 0; _incoming_side < n_sides; ++_incoming_side)
1489 94688113 : if (neighbor->neighbor_ptr(_incoming_side) == _current_elem)
1490 : break;
1491 : }
1492 : else
1493 175399 : _incoming_side = neighbor->which_neighbor_am_i(_current_elem);
1494 :
1495 31169324 : _last_elem = _current_elem;
1496 31169324 : _current_elem = neighbor;
1497 : ray->setCurrentElem(neighbor);
1498 31169324 : 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 31169324 : 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 31169324 : if (_study.hasInternalSidesets() && (subdomain_changed || _intersected_extrema.atExtrema()))
1516 : {
1517 2506 : applyOnInternalBoundary(ray);
1518 :
1519 : // Internal RayBC killed a Ray
1520 2502 : if (!_should_continue)
1521 : {
1522 : traceAssert(!ray->shouldContinue(), "Should be the same");
1523 : debugRay("Internal RayBC killed the ray");
1524 :
1525 234 : onCompleteTrace(ray);
1526 234 : return;
1527 : }
1528 :
1529 : // Internal RayBC changed the Ray
1530 2268 : 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 441 : 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 441 : if (dot > -TRACE_TOLERANCE)
1540 : {
1541 405 : _incoming_side = _last_elem->which_neighbor_am_i(_current_elem);
1542 405 : _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 441 : 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 31169086 : if (neighbor->processor_id() != _pid)
1563 : {
1564 1156350 : if (_intersected_extrema.atExtrema())
1565 : {
1566 : debugRay("Neighbor is off processor but continuing to move through neighbors");
1567 : }
1568 : else
1569 : {
1570 1011341 : continueTraceOffProcessor(ray);
1571 1011341 : return;
1572 : }
1573 : }
1574 :
1575 : // Neighbor is on processor, call subdomain setup if needed
1576 30157745 : if (subdomain_changed)
1577 2422 : 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 1410360 : applyOnExternalBoundary(ray);
1586 :
1587 : traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point, TRACE_TOLERANCE),
1588 : "RayBC changed the Ray point");
1589 :
1590 : // Quit tracing if the Ray was killed by a BC
1591 1410350 : if (!_should_continue)
1592 : {
1593 : traceAssert(!ray->shouldContinue(), "Should be the same");
1594 : debugRay("Exiting due to death by BC");
1595 :
1596 1407047 : onCompleteTrace(ray);
1597 1407047 : return;
1598 : }
1599 : // RayBC changed the direction of the Ray
1600 3303 : if (ray->trajectoryChanged())
1601 : {
1602 : debugRay("RayBC changed the trajectory");
1603 : debugRay(" new direction = ", ray->direction());
1604 : traceAssert(ray->direction() *
1605 : _study.getSideNormal(_current_elem, _intersected_side, _tid) <
1606 : TRACE_TOLERANCE,
1607 : "Reflected ray is not incoming");
1608 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1609 :
1610 3303 : _last_elem = _current_elem;
1611 3303 : _incoming_point = _intersection_point;
1612 3303 : _incoming_side = _intersected_side;
1613 : ray->setCurrentPoint(_incoming_point);
1614 : ray->setCurrentIncomingSide(_incoming_side);
1615 3303 : onTrajectoryChanged(ray);
1616 : }
1617 : }
1618 :
1619 30161048 : onContinueTrace(ray);
1620 : } while (true);
1621 :
1622 : // If a trace made its way down here and didn't return... it failed
1623 : failTrace("Could not find an intersection", _study.tolerateFailure(), __LINE__);
1624 : }
1625 :
1626 : void
1627 4120844 : TraceRay::onCompleteTrace(const std::shared_ptr<Ray> & ray)
1628 : {
1629 7916772 : for (RayKernelBase * rk : _study.currentRayKernels(_tid))
1630 3795928 : rk->postTrace();
1631 :
1632 : debugRay("Called onCompleteTrace()\n", (*_current_ray)->getInfo());
1633 : if (_intersection_distance > 0)
1634 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1635 : possiblySaveDebugRayMesh();
1636 :
1637 4120844 : if (_current_cached_trace)
1638 : {
1639 11392 : _current_cached_trace->_last = true;
1640 :
1641 11392 : if (_intersection_distance > 0)
1642 : {
1643 : _current_cached_trace->addPoint(ray->currentPoint());
1644 10285 : if (_study.dataOnCacheTraces())
1645 170 : _current_cached_trace->lastPoint()._data = ray->data();
1646 10285 : if (_study.auxDataOnCacheTraces())
1647 134 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1648 : }
1649 :
1650 : mooseAssert(ray->stationary() == _current_cached_trace->stationary(), "Stationary mismatch");
1651 : }
1652 4120844 : }
1653 :
1654 : void
1655 30161696 : TraceRay::onContinueTrace(const std::shared_ptr<Ray> & ray)
1656 : {
1657 : traceAssert(ray->shouldContinue(), "Ray must continue");
1658 :
1659 30161696 : if (_current_cached_trace && _study.segmentsOnCacheTraces() && _intersection_distance > 0)
1660 : {
1661 : _current_cached_trace->addPoint(ray->currentPoint());
1662 45860 : if (_study.dataOnCacheTraces())
1663 651 : _current_cached_trace->lastPoint()._data = ray->data();
1664 45860 : if (_study.auxDataOnCacheTraces())
1665 507 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1666 : }
1667 30161696 : }
1668 :
1669 : void
1670 1259326 : TraceRay::continueTraceOffProcessor(const std::shared_ptr<Ray> & ray)
1671 : {
1672 : traceAssert(ray->currentElem() == _current_elem, "Ray currentElem() invalid");
1673 : traceAssert(ray->currentIncomingSide() == _incoming_side, "Ray currentIncomingSide() invalid");
1674 : traceAssert(ray->currentPoint() == _incoming_point, "Ray currentPoint() invalid");
1675 : traceAssert(_current_elem->processor_id() != _pid, "Off processor trace is not off processor");
1676 : debugRay("Ray going off processor to ", _current_elem->processor_id());
1677 :
1678 : ray->addProcessorCrossing();
1679 :
1680 1259326 : if (_current_cached_trace && _intersection_distance > 0)
1681 : {
1682 2194 : _current_cached_trace->addPoint(_incoming_point);
1683 2194 : if (_study.dataOnCacheTraces())
1684 21 : _current_cached_trace->lastPoint()._data = ray->data();
1685 2194 : if (_study.auxDataOnCacheTraces())
1686 21 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1687 : }
1688 :
1689 : if (_intersection_distance > 0)
1690 : possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
1691 : possiblySaveDebugRayMesh();
1692 1259326 : }
1693 :
1694 : void
1695 4392 : TraceRay::onTrajectoryChanged(const std::shared_ptr<Ray> & ray)
1696 : {
1697 : #ifndef NDEBUG
1698 : if (_study.verifyTraceIntersections() &&
1699 : (_intersected_extrema.atExtrema()
1700 : ? !_current_elem->close_to_point(ray->currentPoint(), LOOSE_TRACE_TOLERANCE)
1701 : : !_current_elem->contains_point(ray->currentPoint())))
1702 : failTrace("Elem does not contain point after trajectory change",
1703 : /* warning = */ false,
1704 : __LINE__);
1705 : #endif
1706 :
1707 : traceAssert(ray->shouldContinue(), "Ray should continue when trajectory is being changed");
1708 :
1709 : ray->setTrajectoryChanged(false);
1710 : ray->addTrajectoryChange();
1711 :
1712 4392 : if (_current_cached_trace && !_study.segmentsOnCacheTraces())
1713 : {
1714 48 : if (_intersection_distance > 0)
1715 : _current_cached_trace->addPoint(ray->currentPoint());
1716 48 : if (_study.dataOnCacheTraces())
1717 0 : _current_cached_trace->lastPoint()._data = ray->data();
1718 48 : if (_study.auxDataOnCacheTraces())
1719 0 : _current_cached_trace->lastPoint()._aux_data = ray->auxData();
1720 : }
1721 4392 : }
1722 :
1723 : void
1724 5382820 : TraceRay::onSubdomainChanged(const std::shared_ptr<Ray> & ray, const bool same_ray)
1725 : {
1726 : debugRay("Calling onSubdomainChanged() on subdomain ", _current_elem->subdomain_id());
1727 : debugRay(" _current_subdomain_id = ", _current_subdomain_id);
1728 :
1729 5382820 : _current_subdomain_id = _current_elem->subdomain_id();
1730 5382820 : _current_subdomain_hmax = _study.subdomainHmax(_current_subdomain_id);
1731 5382820 : _current_elem_type = _current_elem->type();
1732 5382820 : _current_elem_n_sides = _current_elem->n_sides();
1733 :
1734 5382820 : if (_has_ray_kernels)
1735 : {
1736 4955206 : auto & current_ray_kernels = _study.currentRayKernels(_tid);
1737 :
1738 : // If we're still tracing the same Ray, keep track of our old RayKernels
1739 : // so that we don't call preTrace() on them again
1740 4955206 : if (same_ray)
1741 1650 : _old_ray_kernels.insert(current_ray_kernels.begin(), current_ray_kernels.end());
1742 : // If we're not tracing the same Ray, we need to call preTrace() on everything
1743 : else
1744 : _old_ray_kernels.clear();
1745 :
1746 : // Call segmentSubdomainSetup to get new kernels etc
1747 4955206 : _study.segmentSubdomainSetup(_current_subdomain_id, _tid, ray->id());
1748 :
1749 : // Call preTrace() on all of the RayKernels that need it
1750 9917397 : for (RayKernelBase * rk : current_ray_kernels)
1751 : // Haven't called preTrace() for this Ray on this RayKernel yet
1752 : if (!_old_ray_kernels.count(rk))
1753 4962123 : rk->preTrace();
1754 : }
1755 5382820 : }
1756 :
1757 : std::string
1758 9 : TraceRay::failTraceMessage(const std::string & reason, const int line)
1759 : {
1760 9 : std::stringstream oss;
1761 9 : oss << "Ray on processor " << _pid << " and thread " << _tid << " failed to trace";
1762 9 : if (line != -1)
1763 9 : oss << " at line " << line;
1764 18 : oss << "\n\n" << reason << "\n\n";
1765 18 : oss << ((*_current_ray))->getInfo() << "\n";
1766 9 : oss << "Current trace information\n";
1767 9 : oss << " _current_subdomain_id = ";
1768 9 : if (_current_subdomain_id == Elem::invalid_subdomain_id)
1769 0 : oss << "invalid subdomain id\n";
1770 : else
1771 9 : oss << _current_subdomain_id << "\n";
1772 18 : oss << " _current_elem_type = " << Utility::enum_to_string(_current_elem_type) << "\n";
1773 9 : oss << " _current_elem_n_sides = " << _current_elem_n_sides << "\n";
1774 9 : oss << " _incoming_point = ";
1775 : if (_incoming_point == RayTracingCommon::invalid_point)
1776 0 : oss << "invalid point\n";
1777 : else
1778 9 : oss << _incoming_point << "\n";
1779 9 : oss << " _incoming_side = ";
1780 9 : if (_incoming_side == RayTracingCommon::invalid_side)
1781 0 : oss << "invalid side\n";
1782 : else
1783 9 : oss << _incoming_side << "\n";
1784 9 : oss << " _intersection_point = ";
1785 : if (_intersection_point == RayTracingCommon::invalid_point)
1786 0 : oss << "invalid point\n";
1787 : else
1788 9 : oss << _intersection_point << "\n";
1789 9 : oss << " _intersected_side = ";
1790 9 : if (_intersected_side == RayTracingCommon::invalid_side)
1791 0 : oss << "invalid side\n";
1792 : else
1793 9 : oss << _intersected_side << "\n";
1794 9 : oss << " _intersected_extrema = " << _intersected_extrema << "\n";
1795 9 : oss << " _exits_elem = " << _exits_elem << "\n";
1796 9 : if (_current_elem)
1797 18 : oss << _current_elem->get_info();
1798 : else
1799 0 : oss << "_current_elem = invalid\n";
1800 :
1801 : possiblySaveDebugRayMesh();
1802 9 : return oss.str();
1803 9 : }
1804 :
1805 : void
1806 9 : TraceRay::failTrace(const std::string & reason, const bool warning, const int line)
1807 : {
1808 9 : const auto message = failTraceMessage(reason, line);
1809 :
1810 9 : if (warning)
1811 : {
1812 5 : ++_results[FAILED_TRACES];
1813 : mooseWarning(message);
1814 5 : (*_current_ray)->setShouldContinue(false);
1815 5 : _should_continue = false;
1816 : }
1817 : else
1818 4 : mooseError(message);
1819 5 : }
1820 :
1821 : const std::vector<NeighborInfo> &
1822 1261716 : TraceRay::getVertexNeighbors(const Elem * elem, const Node * vertex)
1823 : {
1824 : traceAssert(elem, "Elem must be valid");
1825 : traceAssert(vertex, "Vertex must be valid");
1826 :
1827 : debugRay("Called getVertexNeighbors() with:");
1828 : debugRay(" elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
1829 : debugRay(" vertex->id() = ", vertex->id(), ", at ", (Point)*vertex);
1830 :
1831 : traceAssert(elem->get_node_index(vertex) != libMesh::invalid_uint, "Doesn't contain node");
1832 : traceAssert(elem->is_vertex(elem->get_node_index(vertex)), "Node is not a vertex");
1833 :
1834 1261716 : ++_results[VERTEX_NEIGHBOR_LOOKUPS];
1835 :
1836 : // Return the entry if we have it cached
1837 : auto search = _vertex_neighbors.find(vertex);
1838 1261716 : if (search != _vertex_neighbors.end())
1839 1155977 : return search->second;
1840 :
1841 105739 : ++_results[VERTEX_NEIGHBOR_BUILDS];
1842 :
1843 : // Make a new entry
1844 : debugRay("Building vertex neighbors");
1845 : std::vector<NeighborInfo> & entry =
1846 105739 : _vertex_neighbors.emplace(vertex, std::vector<NeighborInfo>()).first->second;
1847 :
1848 105739 : findNodeNeighbors(elem,
1849 : vertex,
1850 105739 : _neighbor_set,
1851 105739 : _neighbor_untested_set,
1852 105739 : _neighbor_next_untested_set,
1853 105739 : _neighbor_active_neighbor_children,
1854 : entry);
1855 :
1856 : // Fill the side normals
1857 815937 : for (auto & neighbor_info : entry)
1858 2847464 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1859 2137266 : neighbor_info._side_normals[i] =
1860 2137266 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
1861 :
1862 : return entry;
1863 : }
1864 :
1865 : const std::vector<NeighborInfo> &
1866 1261716 : TraceRay::getVertexNeighbors(const Elem * elem, const unsigned short vertex)
1867 : {
1868 : traceAssert(vertex < elem->n_vertices(), "Invalid vertex");
1869 :
1870 1261716 : return getVertexNeighbors(elem, elem->node_ptr(vertex));
1871 : }
1872 :
1873 : const std::vector<NeighborInfo> &
1874 4325854 : TraceRay::getEdgeNeighbors(const Elem * elem,
1875 : const std::pair<const Node *, const Node *> & vertices,
1876 : const Point & point)
1877 : {
1878 : traceAssert(elem, "Invalid elem");
1879 : traceAssert(vertices.first, "Must be valid");
1880 : traceAssert(vertices.second, "Must be valid");
1881 :
1882 : debugRay("Called getEdgeNeighbors() with:");
1883 : debugRay(" elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
1884 : debugRay(" vertices.first = ", vertices.first->id(), " at ", (Point)*vertices.first);
1885 : debugRay(" vertices.second = ", vertices.second->id(), " at ", (Point)*vertices.second);
1886 : debugRay(" point = ", point);
1887 :
1888 : traceAssert(elem->get_node_index(vertices.first) != libMesh::invalid_uint,
1889 : "Doesn't contain vertex");
1890 : traceAssert(elem->get_node_index(vertices.second) != libMesh::invalid_uint,
1891 : "Doesn't contain vertex");
1892 : traceAssert(isWithinSegment(
1893 : (Point)*vertices.first, (Point)*vertices.second, point, LOOSE_TRACE_TOLERANCE),
1894 : "Point not within edge");
1895 :
1896 4325854 : ++_results[EDGE_NEIGHBOR_LOOKUPS];
1897 :
1898 4325854 : const auto ordered_vertices = vertices.first->id() < vertices.second->id()
1899 4325854 : ? vertices
1900 4325854 : : std::make_pair(vertices.second, vertices.first);
1901 :
1902 : // Look for the entry and build if necessary
1903 : std::pair<bool, std::vector<NeighborInfo>> * entry;
1904 : auto search = _edge_neighbors.find(ordered_vertices);
1905 4325854 : if (search != _edge_neighbors.end())
1906 4114249 : entry = &search->second;
1907 : else
1908 : {
1909 : debugRay("Building edge neighbors");
1910 211605 : ++_results[EDGE_NEIGHBOR_BUILDS];
1911 211605 : entry = &_edge_neighbors
1912 211605 : .emplace(ordered_vertices, std::make_pair(true, std::vector<NeighborInfo>()))
1913 : .first->second;
1914 211605 : findEdgeNeighbors(elem,
1915 211605 : ordered_vertices.first,
1916 211605 : ordered_vertices.second,
1917 211605 : _neighbor_set,
1918 211605 : _neighbor_untested_set,
1919 211605 : _neighbor_next_untested_set,
1920 211605 : _neighbor_active_neighbor_children,
1921 211605 : entry->second);
1922 :
1923 : bool all_same_edge = true;
1924 1042460 : for (auto & neighbor_info : entry->second)
1925 : {
1926 : traceAssert(neighbor_info._lower_bound <= neighbor_info._upper_bound,
1927 : "Bound order incorrect");
1928 :
1929 : // Fill the side normals
1930 2492501 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
1931 1661646 : neighbor_info._side_normals[i] =
1932 1661646 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
1933 :
1934 : // See if the bounds are the same as the target edge
1935 830855 : if (neighbor_info._lower_bound != 0 || neighbor_info._upper_bound != 1)
1936 : all_same_edge = false;
1937 : }
1938 211605 : entry->first = all_same_edge;
1939 : }
1940 :
1941 : // Means that all neighbors are not on the exact same edge, so we must
1942 : // validate/invalidate based on if the neighbor's edge contains our point
1943 4325854 : if (!entry->first)
1944 : {
1945 : const auto edge_length =
1946 892 : ((Point)*ordered_vertices.first - (Point)*ordered_vertices.second).norm();
1947 892 : const auto point_location = ((Point)*ordered_vertices.first - point).norm() / edge_length;
1948 6410 : for (auto & info : entry->second)
1949 11036 : info._valid = (info._lower_bound - TRACE_TOLERANCE) < point_location &&
1950 4600 : point_location < (info._upper_bound + TRACE_TOLERANCE);
1951 : }
1952 :
1953 4325854 : return entry->second;
1954 : }
1955 :
1956 : const std::vector<NeighborInfo> &
1957 4325854 : TraceRay::getEdgeNeighbors(const Elem * elem,
1958 : const std::pair<unsigned short, unsigned short> & vertices,
1959 : const Point & point)
1960 : {
1961 : debugRay("Called getEdgeNeighbors(), local index version with:");
1962 : debugRay(" vertices.first = ", vertices.first);
1963 : debugRay(" vertices.second = ", vertices.second);
1964 : traceAssert(vertices.first < elem->n_vertices(),
1965 : "Invalid vertex with ray " + std::to_string((*_current_ray)->id()));
1966 : traceAssert(vertices.second < elem->n_vertices(), "Invalid vertex");
1967 :
1968 4325854 : return getEdgeNeighbors(
1969 4325854 : elem, std::make_pair(elem->node_ptr(vertices.first), elem->node_ptr(vertices.second)), point);
1970 : }
1971 :
1972 : const std::vector<NeighborInfo> &
1973 5571062 : TraceRay::getNeighbors(const Elem * elem, const ElemExtrema & extrema, const Point & point)
1974 : {
1975 5571062 : if (!extrema.atExtrema())
1976 0 : return getPointNeighbors(elem, point);
1977 : if (extrema.atVertex())
1978 1245208 : return getVertexNeighbors(elem, extrema.vertex());
1979 4325854 : return getEdgeNeighbors(elem, extrema.edgeVertices(), point);
1980 : }
1981 :
1982 : const std::vector<NeighborInfo> &
1983 30 : TraceRay::getPointNeighbors(const Elem * elem, const Point & point)
1984 : {
1985 : traceAssert(elem, "Invalid elem");
1986 :
1987 : debugRay("Called getPointNeighbors()");
1988 : debugRay(" elem = ", elem->id());
1989 : debugRay(" point = ", point);
1990 :
1991 30 : ++_results[POINT_NEIGHBOR_BUILDS];
1992 30 : _point_neighbor_helper.clear();
1993 :
1994 30 : findPointNeighbors(elem,
1995 : point,
1996 30 : _neighbor_set,
1997 30 : _neighbor_untested_set,
1998 30 : _neighbor_next_untested_set,
1999 30 : _neighbor_active_neighbor_children,
2000 : _point_neighbor_helper);
2001 :
2002 : // Fill the side normals
2003 60 : for (auto & neighbor_info : _point_neighbor_helper)
2004 60 : for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
2005 30 : neighbor_info._side_normals[i] =
2006 30 : _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
2007 :
2008 30 : return _point_neighbor_helper;
2009 : }
2010 :
2011 : void
2012 35292332 : TraceRay::storeExitsElemResult(const TraceRay::ExitsElemResult result)
2013 : {
2014 35292332 : if (result == HIT_FACE)
2015 30505295 : ++_results[FACE_HITS];
2016 4787037 : else if (result == HIT_VERTEX)
2017 1108104 : ++_results[VERTEX_HITS];
2018 3678933 : else if (result == HIT_EDGE)
2019 3678933 : ++_results[EDGE_HITS];
2020 : else
2021 0 : mooseError("Should not call storeExitsElemResult() with result ", result);
2022 35292332 : }
2023 :
2024 : void
2025 33831036 : TraceRay::onSegment(const std::shared_ptr<Ray> & ray)
2026 : {
2027 : traceAssert((*_current_ray)->currentElem() == _current_elem, "Ray currentElem() incorrect");
2028 : traceAssert((*_current_ray)->currentPoint() == _intersection_point,
2029 : "Ray currentPoint() incorrect");
2030 : traceAssert((*_current_ray)->currentIncomingSide() == _incoming_side,
2031 : "Ray currentIncomingSide() incorrect");
2032 : #ifndef NDEBUG
2033 : if (_study.verifyTraceIntersections())
2034 : {
2035 : if (_current_elem->has_affine_map())
2036 : traceAssert(_current_elem->contains_point(_incoming_point),
2037 : "_current_elem does not contain incoming point");
2038 :
2039 : if (_intersected_side != RayTracingCommon::invalid_side &&
2040 : !_study.sideIsNonPlanar(_current_elem, _intersected_side))
2041 : {
2042 : traceAssert(_elem_side_builder(*_current_elem, _intersected_side)
2043 : .close_to_point(_intersection_point, LOOSE_TRACE_TOLERANCE),
2044 : "Intersected point is not on intersected side");
2045 : traceAssert(!_study.sideIsIncoming(
2046 : _current_elem, _intersected_side, (*_current_ray)->direction(), _tid),
2047 : "Intersected side is not outgoing");
2048 : }
2049 : if (_incoming_side != RayTracingCommon::invalid_side &&
2050 : !_study.sideIsNonPlanar(_current_elem, _incoming_side))
2051 : {
2052 : traceAssert(_elem_side_builder(*_current_elem, _incoming_side)
2053 : .close_to_point(_incoming_point, LOOSE_TRACE_TOLERANCE),
2054 : "Incoming point is not on incoming side");
2055 : if (ray->intersections() != 0 && ray->maxDistance() != 0)
2056 : traceAssert(_study.sideIsIncoming(
2057 : _current_elem, _incoming_side, (*_current_ray)->direction(), _tid),
2058 : "Incoming side is not incoming");
2059 : }
2060 : }
2061 : #endif
2062 : traceAssert(MooseUtils::absoluteFuzzyEqual(_intersection_distance,
2063 : (_intersection_point - _incoming_point).norm()),
2064 : "_intersection_distance is incorrect");
2065 : traceAssert(_current_subdomain_id == _current_elem->subdomain_id(), "Subdomain incorrect");
2066 : traceAssert(MooseUtils::absoluteFuzzyEqual((_incoming_point - _intersection_point).norm(),
2067 : _intersection_distance),
2068 : "Invalid intersection distance");
2069 :
2070 33831036 : _study.reinitSegment(
2071 33831036 : _current_elem, _incoming_point, _intersection_point, _intersection_distance, _tid);
2072 :
2073 33831036 : const auto & rks = _study.currentRayKernels(_tid);
2074 67676133 : for (auto & rk : rks)
2075 : {
2076 33845109 : rk->onSegment();
2077 33845097 : postRayTracingObject(ray, rk);
2078 : }
2079 33831024 : }
2080 :
2081 : void
2082 1411579 : TraceRay::onBoundary(const std::shared_ptr<Ray> & ray, const bool external)
2083 : {
2084 : traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point),
2085 : "Ray currentPoint() not set before onBoundary()");
2086 :
2087 : // Get the RayBCs on bnd_elems
2088 1411579 : _study.getRayBCs(_on_boundary_ray_bcs, _boundary_elems, _tid, ray->id());
2089 :
2090 : // Store this information temprorarily because we are going to change it as we
2091 : // apply each boundary condition
2092 1411579 : const auto old_current_elem = _current_elem;
2093 1411579 : const auto old_intersected_side = _intersected_side;
2094 1411579 : const auto old_intersected_extrema = _intersected_extrema;
2095 1411579 : const auto old_subdomain_id = _current_subdomain_id;
2096 :
2097 : // For each RayBC we found, apply it on the boundaries that we need
2098 2847774 : for (RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2099 : {
2100 : // First, find the boundaries this RayBC is valid on that are also in _boundary_elems.
2101 : // We do this ahead of time so that we can pass in to RayBC::apply if the same
2102 : // boundary condition is being appled multiple times on different boundaries.
2103 : // This is useful in situations like reflection where multiple reflections are
2104 : // necessary at a corner to perfectly reflect.
2105 1436205 : _on_boundary_apply_index.clear();
2106 2974449 : for (MooseIndex(_boundary_elems.size()) bnd_elems_i = 0; bnd_elems_i < _boundary_elems.size();
2107 : ++bnd_elems_i)
2108 1538244 : if (rbc->hasBoundary(_boundary_elems[bnd_elems_i].bnd_id))
2109 1538170 : _on_boundary_apply_index.push_back(bnd_elems_i);
2110 :
2111 : traceAssert(!_on_boundary_apply_index.empty(), "Must not be empty");
2112 :
2113 : // Apply the RayBC on each of the relevant boundary elements
2114 2974365 : for (const auto bnd_elems_i : _on_boundary_apply_index)
2115 : {
2116 : auto & bnd_elem = _boundary_elems[bnd_elems_i];
2117 1538170 : _current_elem = bnd_elem.elem;
2118 1538170 : _current_bnd_id = bnd_elem.bnd_id;
2119 1538170 : _intersected_side = bnd_elem.side;
2120 : _intersected_extrema = bnd_elem.extrema;
2121 1538170 : _current_subdomain_id = _current_elem->subdomain_id();
2122 :
2123 1538170 : rbc->onBoundary(_on_boundary_apply_index.size());
2124 1538160 : postRayTracingObject(ray, rbc);
2125 : }
2126 : }
2127 :
2128 : // Set this info back now that we're done applying BCs
2129 1411569 : _current_elem = old_current_elem;
2130 1411569 : _intersected_side = old_intersected_side;
2131 : _intersected_extrema = old_intersected_extrema;
2132 1411569 : _current_bnd_id = BoundaryInfo::invalid_id;
2133 1411569 : _current_subdomain_id = old_subdomain_id;
2134 :
2135 : // When on an external boundary, the Ray must have been changed or killed.
2136 : // Otherwise, we don't know what to do with it now! If this didn't happen,
2137 : // output a detailed error message.
2138 1411569 : if (external && !ray->trajectoryChanged() && ray->shouldContinue())
2139 : {
2140 9 : std::stringstream oss;
2141 9 : oss << "Don't know what to do with a Ray after it hit an external\n";
2142 9 : oss << "boundary at point " << _intersection_point << "!\n\n";
2143 9 : oss << "When hitting an external RayBC, a Ray must either:\n";
2144 9 : oss << " Be killed by a RayBC\n";
2145 9 : oss << " Have its trajectory changed by the RayBC\n";
2146 9 : oss << "by at least one of the executed RayBCs.\n\n";
2147 9 : oss << "You need to either:\n";
2148 9 : oss << " Kill/change the Ray sooner with RayKernels, internal RayBCs, or a max distance\n";
2149 9 : oss << " Kill/change the Ray on the boundary with a RayBC\n\n";
2150 9 : if (!_on_boundary_ray_bcs.empty())
2151 : {
2152 7 : oss << "RayBCs executed that did not kill or change the Ray:\n";
2153 14 : for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2154 16 : for (const auto & bnd_elem : _boundary_elems)
2155 9 : if (rbc->hasBoundary(bnd_elem.bnd_id))
2156 21 : oss << " " << rbc->typeAndName() << " on boundary " << bnd_elem.bnd_id << " ("
2157 14 : << _mesh.getBoundaryName(bnd_elem.bnd_id) << ")\n";
2158 7 : oss << "\n";
2159 : }
2160 : bool output_header = false;
2161 22 : for (std::size_t i = 0; i < _boundary_elems.size(); ++i)
2162 : {
2163 13 : const auto bnd_id = _boundary_elems[i].bnd_id;
2164 : bool found = false;
2165 15 : for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
2166 9 : if (rbc->hasBoundary(bnd_id))
2167 : {
2168 : found = true;
2169 : break;
2170 : }
2171 :
2172 13 : if (!found)
2173 : {
2174 6 : if (!output_header)
2175 : {
2176 4 : oss << "Boundaries that did not have any RayBCs:\n";
2177 : output_header = true;
2178 : }
2179 12 : oss << " " << bnd_id << " (" << _mesh.getBoundaryName(bnd_id) << ")\n";
2180 : }
2181 : }
2182 :
2183 9 : failTrace(oss.str(), _study.tolerateFailure(), __LINE__);
2184 5 : }
2185 1411565 : }
2186 :
2187 : Real
2188 46867642 : TraceRay::subdomainHmax(const Elem * elem) const
2189 : {
2190 : const auto subdomain_id = elem->subdomain_id();
2191 46867642 : return subdomain_id == _current_subdomain_id ? _current_subdomain_hmax
2192 2151 : : _study.subdomainHmax(subdomain_id);
2193 : }
2194 :
2195 : void
2196 35383257 : TraceRay::postRayTracingObject(const std::shared_ptr<Ray> & ray, const RayTracingObject * rto)
2197 : {
2198 35383257 : if (!ray->shouldContinue())
2199 : {
2200 4075714 : if (_should_continue)
2201 1407280 : _should_continue = false;
2202 : }
2203 31307543 : else if (!_should_continue)
2204 0 : failTrace(rto->typeAndName() +
2205 0 : " set a Ray to continue that was previously set to not continue.\n\n" +
2206 : "Once a Ray has been set to not continue, its continue status cannot change.",
2207 : /* warning = */ false,
2208 : __LINE__);
2209 :
2210 35383257 : if (!_should_continue && ray->trajectoryChanged())
2211 0 : failTrace(rto->typeAndName() +
2212 0 : " changed the trajectory of a Ray that was set to not continue,\n" +
2213 : "or set a Ray whose trajectory was changed to not continue.",
2214 : /* warning = */ false,
2215 : __LINE__);
2216 35383257 : }
|
