Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : // Local includes
13 : #include "DebugRay.h"
14 : #include "ElemExtrema.h"
15 : #include "NeighborInfo.h"
16 : #include "RayTracingCommon.h"
17 : #include "TraceRayBndElement.h"
18 :
19 : // MOOSE Includes
20 : #include "MooseMesh.h"
21 : #include "MooseHashing.h"
22 : #include "MooseTypes.h"
23 : #include "StaticallyAllocatedSet.h"
24 :
25 : // libMesh includes
26 : #include "libmesh/enum_elem_type.h"
27 : #include "libmesh/elem_side_builder.h"
28 :
29 : // Forward declarations
30 : class Ray;
31 : class RayBoundaryConditionBase;
32 : class RayKernelBase;
33 : class RayTracingStudy;
34 : class RayTracingObject;
35 : struct TraceData;
36 : namespace libMesh
37 : {
38 : class BoundaryInfo;
39 : class Edge;
40 : class Mesh;
41 : }
42 :
43 : /**
44 : * Traces Rays through the mesh on a single processor
45 : */
46 : class TraceRay
47 : {
48 : public:
49 : TraceRay(RayTracingStudy & study, const THREAD_ID tid);
50 :
51 12483 : virtual ~TraceRay() {}
52 :
53 : /**
54 : * Enum for the various results reported by this object
55 : */
56 : enum TraceRayResult
57 : {
58 : INTERSECTIONS = 0, // intersections
59 : FACE_HITS = 1, // intersections at an element face
60 : VERTEX_HITS = 2, // intersections at an element vertex
61 : EDGE_HITS = 3, // intersections at an element edge
62 : MOVED_THROUGH_NEIGHBORS = 4, // intersections through point/vertex/edge neighbors
63 : BACKFACE_CULLING_SUCCESSES = 5, // face intersection calls that backface culling worked on
64 : BACKFACE_CULLING_FAILURES = 6, // face interesction calls that backface culling failed on
65 : INTERSECTION_CALLS = 7, // face intersection calls made
66 : VERTEX_NEIGHBOR_BUILDS = 8, // builds for neighbors at a vertex
67 : VERTEX_NEIGHBOR_LOOKUPS = 9, // lookups for neighbors at a vertex
68 : EDGE_NEIGHBOR_BUILDS = 10, // builds for neighbors on an edge
69 : EDGE_NEIGHBOR_LOOKUPS = 11, // lookups for neighborson an edge
70 : POINT_NEIGHBOR_BUILDS = 12, // builds for point neighbors (non edge/vertex)
71 : FAILED_TRACES = 13, // rays that fail (allowed when _tolerate_failure == true),
72 : ENDED_STATIONARY = 14 // rays the end because they are stationary
73 : };
74 :
75 : /**
76 : * Should be called immediately before calling any traces.
77 : */
78 : void preExecute();
79 :
80 : /**
81 : * Called on mesh change
82 : */
83 : void meshChanged();
84 :
85 : /**
86 : * Traces a ray.
87 : */
88 : void trace(const std::shared_ptr<Ray> & ray);
89 :
90 : /**
91 : * Get the various results reported by this object, indexed by TraceRayResult
92 : */
93 143500 : std::vector<unsigned long long int> results() const { return _results; }
94 :
95 : /**
96 : * Enable or disable the use of element normals for backface culling through the getElemNormals()
97 : * method
98 : *
99 : * By default, the _study.getElemNormals() method is not implemented. This is because in order to
100 : * be able to obtain the element normals often enough during a trace, one needs to implement a
101 : * method to cache the normals. Otherwise, it is not worthwhile to use the normals during the
102 : * trace because the time spent computing the normals on the fly will be greater than the time
103 : * gain by using the normals in the tracing process.
104 : */
105 44 : void setBackfaceCulling(const bool backface_culling) { _backface_culling = backface_culling; }
106 :
107 : /**
108 : * Gets the current Ray that is being traced
109 : */
110 8443 : const std::shared_ptr<Ray> * const & currentRay() const { return _current_ray; }
111 : /**
112 : * Gets the element that the current Ray is being traced in
113 : */
114 8443 : const Elem * const & currentElem() const { return _current_elem; }
115 : /**
116 : * Gets the current intersection point for the Ray that is being traced.
117 : *
118 : * If the Ray ends within the current element, this will be the point at which
119 : * the Ray ends within the current element.
120 : */
121 8443 : const Point & currentIntersectionPoint() const { return _intersection_point; }
122 : /**
123 : * Gets the current incoming point for the Ray that is being traced.
124 : *
125 : * This is ONLY guaranteed to be valid during execution of RayKernels!
126 : */
127 3825 : const Point & currentIncomingPoint() const { return _incoming_point; }
128 : /**
129 : * Gets the current incoming side for the Ray that is being traced.
130 : *
131 : * This is only valid during onSegment() (when RayKernels are being called)
132 : */
133 3825 : const unsigned short & currentIncomingSide() const { return _incoming_side; }
134 : /**
135 : * Gets the side that the current Ray intersected.
136 : *
137 : * If the ray ends within the current element, this will be
138 : * RayTracingCommon::invalid_side.
139 : */
140 8443 : const unsigned short & currentIntersectedSide() const { return _intersected_side; }
141 : /**
142 : * Gets the element extrema (vertex/edge) that the current Ray intersected
143 : */
144 8443 : const ElemExtrema & currentIntersectedExtrema() const { return _intersected_extrema; }
145 : /**
146 : * Gets the current intersection distance for the Ray that is being traced.
147 : *
148 : * When within an element (executing with a RayKernel), this is the segment distance.
149 : */
150 3825 : const Real & currentIntersectionDistance() const { return _intersection_distance; }
151 : /**
152 : * Gets the BoundaryID of the boundary that the Ray intersected and is being applied
153 : * a boundary condition.
154 : */
155 4618 : const BoundaryID & currentBoundaryID() const { return _current_bnd_id; }
156 : /**
157 : * Gets the subdomain of the current element that the Ray is being traced in.
158 : */
159 8443 : const SubdomainID & currentSubdomainID() const { return _current_subdomain_id; }
160 :
161 : /**
162 : * Gets the neighbors at a vertex.
163 : * @param elem An elem that contains the vertex
164 : * @param vertex The Node that is the vertex
165 : */
166 : const std::vector<NeighborInfo> & getVertexNeighbors(const Elem * elem, const Node * vertex);
167 : /**
168 : * Gets the neighbors at a vertex.
169 : * @param elem An elem that contains the vertex
170 : * @param vertex The local ID of the vertex on elem
171 : */
172 : const std::vector<NeighborInfo> & getVertexNeighbors(const Elem * elem,
173 : const unsigned short vertex);
174 : /**
175 : * Get the neighbors at an edge.
176 : * @param elem An elem that contains the edge
177 : * @param vertices A pair of Nodes that are the vertices that contain the edge fully
178 : * @param point The point on the edge where neighbors are desired
179 : */
180 : const std::vector<NeighborInfo> &
181 : getEdgeNeighbors(const Elem * elem,
182 : const std::pair<const Node *, const Node *> & vertices,
183 : const Point & point);
184 : /**
185 : * Get the neighbors at an edge.
186 : * @param elem An elem that contains the edge
187 : * @param vertices A pair of local vertex IDs of the vertices that contain the edge fully
188 : * @param point The point on the edge where neighbors are desired
189 : */
190 : const std::vector<NeighborInfo> &
191 : getEdgeNeighbors(const Elem * elem,
192 : const std::pair<unsigned short, unsigned short> & vertices,
193 : const Point & point);
194 : /**
195 : * Get the point neighbors
196 : * @param elem The elem
197 : * @param point The point
198 : */
199 : const std::vector<NeighborInfo> & getPointNeighbors(const Elem * elem, const Point & point);
200 : /**
201 : * Get the point/vertex/edge neighbors depending on extrema
202 : */
203 : const std::vector<NeighborInfo> &
204 : getNeighbors(const Elem * elem, const ElemExtrema & extrema, const Point & point);
205 :
206 : private:
207 : /**
208 : * Called on a single segment traced by a Ray.
209 : */
210 : void onSegment(const std::shared_ptr<Ray> & ray);
211 :
212 : /**
213 : * Called when a Ray hits a boundary.
214 : */
215 : void onBoundary(const std::shared_ptr<Ray> & ray, const bool external);
216 :
217 : /**
218 : * Called when a Ray is finished tracing (whenever !ray->shouldContinue())
219 : * @param ray The ray that is finished tracing
220 : */
221 : void onCompleteTrace(const std::shared_ptr<Ray> & ray);
222 : /**
223 : * Called when a Ray is continuing to trace after segment
224 : * @param ray The ray
225 : */
226 : void onContinueTrace(const std::shared_ptr<Ray> & /* ray */);
227 : /**
228 : * Called when a Ray's trajectory changes
229 : * @param ray The ray
230 : */
231 : void onTrajectoryChanged(const std::shared_ptr<Ray> & ray);
232 : /**
233 : * Called when the subdomain changes.
234 : * @param ray The current Ray
235 : * @param same_ray Whether or not this is being called on the same Ray as previously
236 : */
237 : void onSubdomainChanged(const std::shared_ptr<Ray> & ray, const bool same_ray);
238 :
239 : /**
240 : * Creates a useful error string with current tracing information.
241 : *
242 : * If line is provided, will include the line number.
243 : */
244 : std::string failTraceMessage(const std::string & reason, const int line = -1);
245 :
246 : /**
247 : * Specialized mooseError for a failed Ray trace with detailed information regarding the trace.
248 : *
249 : * If warning = true, use mooseWarning instead.
250 : * If line is provided, output will include the line number.
251 : */
252 : void failTrace(const std::string & reason, const bool warning, const int line = -1);
253 :
254 : /**
255 : * Enum for the different exit results for exitElem()
256 : */
257 : enum ExitsElemResult
258 : {
259 : NO_EXIT = 0, // doesn't exit
260 : HIT_FACE = 1, // exits through a face
261 : HIT_VERTEX = 2, // exits through a vertex
262 : HIT_EDGE = 3 // exits through an edge
263 : };
264 :
265 : /**
266 : * Determines if _current_ray moving in direction _direction exits elem.
267 : * @param elem The element to check
268 : * @param elem_type The type of elem (used to avoid virtual function calls on elem)
269 : * @param incoming_side The incoming side (if any)
270 : * @param intersection_point To be modified with the intersection point, if any
271 : * @param intersected_side To be modified with the intersected side, if any
272 : * @param intersected_extrema To be modified with the intersected extremum
273 : * (vertex/edge), if any
274 : * @param intersection_distance To be modified with the intersection distance, if any
275 : * @param normals Pointer to the normals for the current element, if any
276 : * @return Whether or not _current_ray exits elem
277 : */
278 : ExitsElemResult exitsElem(const Elem * elem,
279 : const ElemType elem_type,
280 : const unsigned short incoming_side,
281 : Point & intersection_point,
282 : unsigned short & intersected_side,
283 : ElemExtrema & intersected_extrema,
284 : Real & intersection_distance,
285 : const Point * normals);
286 :
287 : /**
288 : * Determines if _current_ray moving in direction _direction exits \p elem.
289 : *
290 : * This function is templated to work only with elements that are faces or cells (2D or 3D)
291 : *
292 : * @param elem The element to check
293 : * @param incoming_side The incoming side (if any)
294 : * @param intersection_point To be modified with the intersection point, if any
295 : * @param intersected_side To be modified with the intersected side, if any
296 : * @param intersected_extrema To be modified with the intersected extremum
297 : * (vertex/edge), if any
298 : * @param intersection_distance To be modified with the intersection distance, if any
299 : * @param normals Pointer to the normals for the current element, if any
300 : * @return Whether or not _current_ray exits elem through a face
301 : */
302 : template <typename T, typename FirstOrderT>
303 : typename std::enable_if<!std::is_base_of<Edge, T>::value, bool>::type
304 : exitsElem(const Elem * elem,
305 : const unsigned short incoming_side,
306 : Point & intersection_point,
307 : unsigned short & intersected_side,
308 : ElemExtrema & intersected_extrema,
309 : Real & intersection_distance,
310 : const Point * normals);
311 :
312 : /**
313 : * Determines if _current_ray moving in direction _direction exits \p elem.
314 : *
315 : * This function is templated to work only with elements that are edges (1D).
316 : *
317 : * @param elem The element to check
318 : * @param incoming_side The incoming side (if any)
319 : * @param intersection_point To be modified with the intersection point, if any
320 : * @param intersected_side To be modified with the intersected side, if any
321 : * @param intersected_extrema To be modified with the intersected extremum
322 : * (vertex/edge), if any
323 : * @param intersection_distance To be modified with the intersection distance, if any
324 : * @param normals Pointer to the normals for the current element, if any
325 : * @return Whether or not _current_ray exits elem through a face
326 : */
327 : template <typename T, typename FirstOrderT>
328 : typename std::enable_if<std::is_base_of<Edge, T>::value, bool>::type
329 : exitsElem(const Elem * elem,
330 : const unsigned short incoming_side,
331 : Point & intersection_point,
332 : unsigned short & intersected_side,
333 : ElemExtrema & intersected_extrema,
334 : Real & intersection_distance,
335 : const Point * normals);
336 :
337 : /**
338 : * Moves the Ray though neighbors (vertex/edge/point)
339 : * @param neighbors The neighbors to try to move through
340 : * @param last_elem The last element that was moved through - this one will be tried last
341 : * @param best_elem The resulting element for which an intersection was found
342 : * @param best_elem_incoming_side The incoming side on the resulting element for which
343 : * @return The trace result
344 : */
345 : ExitsElemResult moveThroughNeighbors(const std::vector<NeighborInfo> & neighbors,
346 : const Elem * last_elem,
347 : const Elem *& best_elem,
348 : unsigned short & best_elem_incoming_side);
349 :
350 : /**
351 : * Sees if a Ray can move through a neighbor (vertex/edge/point)
352 : * @param neighbor_info The NeighborInfo for the neighbor
353 : * @param incoming_side To be filled with the incoming side on the neighbor, if any
354 : * @param intersection_point To be filled with the intersection point on the neighbor, if any
355 : * @param intersected_side To be filled with the intersected side on the neighbor, if any
356 : * @param intersected_extrema To be filled with the intersected extrema on the neighbor, if any
357 : * @param intersection_distance To be filled with the intersection distance, if any
358 : * @return The trace result
359 : */
360 : ExitsElemResult moveThroughNeighbor(const NeighborInfo & neighbor_info,
361 : unsigned short & incoming_side,
362 : Point & intersection_point,
363 : unsigned short & intersected_side,
364 : ElemExtrema & intersected_extrema,
365 : Real & intersection_distance);
366 :
367 : /**
368 : * Gets and applies external boundary conditions in _current_elem on side _intersected_side at
369 : * _intersection_point.
370 : */
371 : void applyOnExternalBoundary(const std::shared_ptr<Ray> & ray);
372 : /**
373 : * Gets and applies internal boundary conditions (if any) from _current_elem, _last_elem, and
374 : * any other point neighbors that have internal sidesets at _intersection_point.
375 : */
376 : void applyOnInternalBoundary(const std::shared_ptr<Ray> & ray);
377 : /**
378 : * Helper for possibly storing boundary information in _boundary_elems, which is storage for
379 : * boundary elements (elem, side, bnd_id) that need to have RayBCs applied to them.
380 : *
381 : * For each BoundaryID in bnd_ids, a ConstBndElement will be added if said BoundaryID
382 : * does not already exist in _boundary_elems.
383 : *
384 : * @param elem The element to possibly add
385 : * @param side The side to possibly add
386 : * @param bnd_ids The BoundaryIDs to possibly add
387 : * @param extrema The intersected elem extrema (vertex/edge), if any
388 : */
389 : void possiblyAddToBoundaryElems(const Elem * elem,
390 : const unsigned short side,
391 : const std::vector<BoundaryID> & bnd_ids,
392 : const ElemExtrema & extrema);
393 : /**
394 : * Sets up a ray to continue tracing off processor
395 : * @param ray The ray
396 : * @param elem The element that is owned by another processor
397 : * @param incoming_side The incoming side elem
398 : * @param point The point on elem at which ray continues
399 : */
400 : void continueTraceOffProcessor(const std::shared_ptr<Ray> & ray);
401 :
402 : /**
403 : * Finds (if any) an element side that is on the boundary and is outgoing
404 : * at _intersection_point that is on the extrema _intersected_extrema for _current_elem.
405 : *
406 : * This is necessary when a Ray hits a point that is on the boundary but
407 : * is not actually on a boundary side
408 : *
409 : * Fills the found boundary side into \p boundary_side, the element said side is on into
410 : * \p boundary_elem and the extrema on said elem into \p boundary_extrema
411 : */
412 : void findExternalBoundarySide(unsigned short & boundary_side,
413 : ElemExtrema & boundary_extrema,
414 : const Elem *& boundary_elem);
415 :
416 : /**
417 : * Stores the result given by an intersection in _results as necessary
418 : */
419 : void storeExitsElemResult(const ExitsElemResult result);
420 :
421 : /**
422 : * Get the approximate subdomain hmax for an element.
423 : *
424 : * Uses the cached value of _current_subdomain_hmax if the element is in the cached subdomain.
425 : */
426 : Real subdomainHmax(const Elem * elem) const;
427 :
428 : /**
429 : * Called after executing a RayTracingObject (RayBCs and RayKernels)
430 : *
431 : * Verifies the configuration of ray->shouldContinue() and ray->trajectoryChanged()
432 : */
433 : void postRayTracingObject(const std::shared_ptr<Ray> & ray, const RayTracingObject * rto);
434 :
435 : /// The RayTracingStudy
436 : RayTracingStudy & _study;
437 : /// The MooseMesh
438 : MooseMesh & _mesh;
439 : /// The mesh dimension
440 : const unsigned int _dim;
441 : /// The BoundaryInfo for the mesh
442 : const BoundaryInfo & _boundary_info;
443 : /// The processor id
444 : const processor_id_type _pid;
445 : /// The thread id
446 : const THREAD_ID _tid;
447 :
448 : /// Whether or not to use element normals for backface culling
449 : bool _backface_culling;
450 :
451 : /// The TraceData for the current cached trace (if any)
452 : TraceData * _current_cached_trace;
453 :
454 : /// The current ray being traced
455 : const std::shared_ptr<Ray> * _current_ray;
456 : /// The element the current Ray is being traced in
457 : const Elem * _current_elem;
458 : /// The last element the current Ray was traced in
459 : const Elem * _last_elem;
460 : /// The current SubdomainID
461 : SubdomainID _current_subdomain_id;
462 : /// The current subdomain hmax
463 : Real _current_subdomain_hmax;
464 : /// The current elem type (constant on subdomain), used to avoid elem->type() calls
465 : libMesh::ElemType _current_elem_type;
466 : /// The number of sides on the current elem, used to avoid elem->n_sides() virtual calls
467 : unsigned short _current_elem_n_sides;
468 : /// The incoming point of the current Ray
469 : Point _incoming_point;
470 : /// The incoming side of the current Ray
471 : unsigned short _incoming_side;
472 : /// Whether or not the current Ray should continue
473 : bool _should_continue;
474 :
475 : /// The work point for the intersection of the current Ray
476 : Point _intersection_point;
477 : /// The work point for the intersected side of the current Ray
478 : unsigned short _intersected_side;
479 : /// The work point for the intersected vertex/edge vertices of the current Ray, if any
480 : ElemExtrema _intersected_extrema;
481 : /// The intersected vertex/edge vertices for the previous intersection, if any
482 : ElemExtrema _last_intersected_extrema;
483 : /// The work point for the intersection distance of the current Ray
484 : Real _intersection_distance;
485 : /// The current BoundaryID (used when calling RayBoundaryConditionBase::onBoundary())
486 : BoundaryID _current_bnd_id;
487 :
488 : /// Whether or not the current trace exits an element
489 : bool _exits_elem;
490 :
491 : /// The normals for the current element for backface culling (pointer to the first normal - optional)
492 : const Point * _current_normals;
493 :
494 : /// Reusable vector for calling _boundary_info.boundary_ids()
495 : std::vector<BoundaryID> _boundary_ids;
496 : /// Boundary elements that need RayBCs to be applied
497 : std::vector<TraceRayBndElement> _boundary_elems;
498 :
499 : /// The cached vertex neighbors
500 : std::unordered_map<const Node *, std::vector<NeighborInfo>> _vertex_neighbors;
501 : /// The cached edge neighbors
502 : std::unordered_map<std::pair<const Node *, const Node *>,
503 : std::pair<bool, std::vector<NeighborInfo>>>
504 : _edge_neighbors;
505 :
506 : /// Reusable for building neighbors
507 : std::vector<NeighborInfo> _point_neighbor_helper;
508 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> _neighbor_set;
509 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> _neighbor_untested_set;
510 : MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> _neighbor_next_untested_set;
511 : std::vector<const Elem *> _neighbor_active_neighbor_children;
512 :
513 : /// Results over all of the local traces, indexed by TraceRayResult
514 : std::vector<unsigned long long int> _results;
515 :
516 : /// Helper for building element sides without excessive allocation
517 : libMesh::ElemSideBuilder _elem_side_builder;
518 :
519 : /// Reusable for getting the RayBCs in onBoundary()
520 : std::vector<RayBoundaryConditionBase *> _on_boundary_ray_bcs;
521 : /// Reusable for which boundary elements to apply for a specific RayBC in onBoundary()
522 : std::vector<std::size_t> _on_boundary_apply_index;
523 :
524 : /// Whether or not the RayTracingStudy has any RayKernels
525 : bool _has_ray_kernels;
526 : /// Whether or not the domain is rectangular (defined perfectly by its bounding box)
527 : bool _is_rectangular_domain;
528 :
529 : /// Helper for avoiding calling preTrace() on the same RayKernel multiple times
530 : std::set<RayKernelBase *> _old_ray_kernels;
531 :
532 : #ifdef DEBUG_RAY_MESH_IF
533 : Mesh * _debug_mesh;
534 : Parallel::Communicator _debug_comm;
535 : unsigned int _debug_node_count;
536 : #endif
537 : };
538 :
539 : // Helper for a mooseAssert() with useful trace information
540 : #ifdef NDEBUG
541 : #define traceAssert(asserted, msg) ((void)0)
542 : #else
543 : #define traceAssert(asserted, msg) \
544 : do \
545 : { \
546 : if (!(asserted)) \
547 : mooseAssert(asserted, failTraceMessage(msg, __LINE__)); \
548 : } while (0)
549 : #endif
|
