TraceRay Class Reference

Traces Rays through the mesh on a single processor. More...

#include <TraceRay.h>

Public Types
enum	TraceRayResult {   INTERSECTIONS = 0 , FACE_HITS = 1 , VERTEX_HITS = 2 , EDGE_HITS = 3 ,   MOVED_THROUGH_NEIGHBORS = 4 , BACKFACE_CULLING_SUCCESSES = 5 , BACKFACE_CULLING_FAILURES = 6 , INTERSECTION_CALLS = 7 ,   VERTEX_NEIGHBOR_BUILDS = 8 , VERTEX_NEIGHBOR_LOOKUPS = 9 , EDGE_NEIGHBOR_BUILDS = 10 , EDGE_NEIGHBOR_LOOKUPS = 11 ,   POINT_NEIGHBOR_BUILDS = 12 , FAILED_TRACES = 13 , ENDED_STATIONARY = 14 }
	Enum for the various results reported by this object. More...

Public Member Functions
	TraceRay (RayTracingStudy &study, const THREAD_ID tid)
virtual	~TraceRay ()
void	preExecute ()
	Should be called immediately before calling any traces.
void	meshChanged ()
	Called on mesh change.
void	trace (const std::shared_ptr< Ray > &ray)
	Traces a ray.
std::vector< unsigned long long int >	results () const
	Get the various results reported by this object, indexed by TraceRayResult.
void	setBackfaceCulling (const bool backface_culling)
	Enable or disable the use of element normals for backface culling through the getElemNormals() method.
const std::shared_ptr< Ray > *const &	currentRay () const
	Gets the current Ray that is being traced.
const Elem *const &	currentElem () const
	Gets the element that the current Ray is being traced in.
const Point &	currentIntersectionPoint () const
	Gets the current intersection point for the Ray that is being traced.
const Point &	currentIncomingPoint () const
	Gets the current incoming point for the Ray that is being traced.
const unsigned short &	currentIncomingSide () const
	Gets the current incoming side for the Ray that is being traced.
const unsigned short &	currentIntersectedSide () const
	Gets the side that the current Ray intersected.
const ElemExtrema &	currentIntersectedExtrema () const
	Gets the element extrema (vertex/edge) that the current Ray intersected.
const Real &	currentIntersectionDistance () const
	Gets the current intersection distance for the Ray that is being traced.
const BoundaryID &	currentBoundaryID () const
	Gets the BoundaryID of the boundary that the Ray intersected and is being applied a boundary condition.
const SubdomainID &	currentSubdomainID () const
	Gets the subdomain of the current element that the Ray is being traced in.
const std::vector< NeighborInfo > &	getVertexNeighbors (const Elem *elem, const Node *vertex)
	Gets the neighbors at a vertex.
const std::vector< NeighborInfo > &	getVertexNeighbors (const Elem *elem, const unsigned short vertex)
	Gets the neighbors at a vertex.
const std::vector< NeighborInfo > &	getEdgeNeighbors (const Elem *elem, const std::pair< const Node *, const Node * > &vertices, const Point &point)
	Get the neighbors at an edge.
const std::vector< NeighborInfo > &	getEdgeNeighbors (const Elem *elem, const std::pair< unsigned short, unsigned short > &vertices, const Point &point)
	Get the neighbors at an edge.
const std::vector< NeighborInfo > &	getPointNeighbors (const Elem *elem, const Point &point)
	Get the point neighbors.
const std::vector< NeighborInfo > &	getNeighbors (const Elem *elem, const ElemExtrema &extrema, const Point &point)
	Get the point/vertex/edge neighbors depending on extrema.

Private Types
enum	ExitsElemResult { NO_EXIT = 0 , HIT_FACE = 1 , HIT_VERTEX = 2 , HIT_EDGE = 3 }
	Enum for the different exit results for exitElem() More...

Private Member Functions
void	onSegment (const std::shared_ptr< Ray > &ray)
	Called on a single segment traced by a Ray.
void	onBoundary (const std::shared_ptr< Ray > &ray, const bool external)
	Called when a Ray hits a boundary.
void	onCompleteTrace (const std::shared_ptr< Ray > &ray)
	Called when a Ray is finished tracing (whenever !ray->shouldContinue())
void	onContinueTrace (const std::shared_ptr< Ray > &)
	Called when a Ray is continuing to trace after segment.
void	onTrajectoryChanged (const std::shared_ptr< Ray > &ray)
	Called when a Ray's trajectory changes.
void	onSubdomainChanged (const std::shared_ptr< Ray > &ray, const bool same_ray)
	Called when the subdomain changes.
std::string	failTraceMessage (const std::string &reason, const int line=-1)
	Creates a useful error string with current tracing information.
void	failTrace (const std::string &reason, const bool warning, const int line=-1)
	Specialized mooseError for a failed Ray trace with detailed information regarding the trace.
ExitsElemResult	exitsElem (const Elem *elem, const ElemType elem_type, const unsigned short incoming_side, Point &intersection_point, unsigned short &intersected_side, ElemExtrema &intersected_extrema, Real &intersection_distance, const Point *normals)
	Determines if _current_ray moving in direction _direction exits elem.
template<typename T , typename FirstOrderT >
std::enable_if<!std::is_base_of< Edge, T >::value, bool >::type	exitsElem (const Elem *elem, const unsigned short incoming_side, Point &intersection_point, unsigned short &intersected_side, ElemExtrema &intersected_extrema, Real &intersection_distance, const Point *normals)
	Determines if _current_ray moving in direction _direction exits elem.
template<typename T , typename FirstOrderT >
std::enable_if< std::is_base_of< Edge, T >::value, bool >::type	exitsElem (const Elem *elem, const unsigned short incoming_side, Point &intersection_point, unsigned short &intersected_side, ElemExtrema &intersected_extrema, Real &intersection_distance, const Point *normals)
	Determines if _current_ray moving in direction _direction exits elem.
ExitsElemResult	moveThroughNeighbors (const std::vector< NeighborInfo > &neighbors, const Elem *last_elem, const Elem *&best_elem, unsigned short &best_elem_incoming_side)
	Moves the Ray though neighbors (vertex/edge/point)
ExitsElemResult	moveThroughNeighbor (const NeighborInfo &neighbor_info, unsigned short &incoming_side, Point &intersection_point, unsigned short &intersected_side, ElemExtrema &intersected_extrema, Real &intersection_distance)
	Sees if a Ray can move through a neighbor (vertex/edge/point)
void	applyOnExternalBoundary (const std::shared_ptr< Ray > &ray)
	Gets and applies external boundary conditions in _current_elem on side _intersected_side at _intersection_point.
void	applyOnInternalBoundary (const std::shared_ptr< Ray > &ray)
	Gets and applies internal boundary conditions (if any) from _current_elem, _last_elem, and any other point neighbors that have internal sidesets at _intersection_point.
void	possiblyAddToBoundaryElems (const Elem *elem, const unsigned short side, const std::vector< BoundaryID > &bnd_ids, const ElemExtrema &extrema)
	Helper for possibly storing boundary information in _boundary_elems, which is storage for boundary elements (elem, side, bnd_id) that need to have RayBCs applied to them.
void	continueTraceOffProcessor (const std::shared_ptr< Ray > &ray)
	Sets up a ray to continue tracing off processor.
void	findExternalBoundarySide (unsigned short &boundary_side, ElemExtrema &boundary_extrema, const Elem *&boundary_elem)
	Finds (if any) an element side that is on the boundary and is outgoing at _intersection_point that is on the extrema _intersected_extrema for _current_elem.
void	storeExitsElemResult (const ExitsElemResult result)
	Stores the result given by an intersection in _results as necessary.
Real	subdomainHmax (const Elem *elem) const
	Get the approximate subdomain hmax for an element.
void	postRayTracingObject (const std::shared_ptr< Ray > &ray, const RayTracingObject *rto)
	Called after executing a RayTracingObject (RayBCs and RayKernels)

Private Attributes
RayTracingStudy &	_study
	The RayTracingStudy.
MooseMesh &	_mesh
	The MooseMesh.
const unsigned int	_dim
	The mesh dimension.
const BoundaryInfo &	_boundary_info
	The BoundaryInfo for the mesh.
const processor_id_type	_pid
	The processor id.
const THREAD_ID	_tid
	The thread id.
bool	_backface_culling
	Whether or not to use element normals for backface culling.
TraceData *	_current_cached_trace
	The TraceData for the current cached trace (if any)
const std::shared_ptr< Ray > *	_current_ray
	The current ray being traced.
const Elem *	_current_elem
	The element the current Ray is being traced in.
const Elem *	_last_elem
	The last element the current Ray was traced in.
SubdomainID	_current_subdomain_id
	The current SubdomainID.
Real	_current_subdomain_hmax
	The current subdomain hmax.
libMesh::ElemType	_current_elem_type
	The current elem type (constant on subdomain), used to avoid elem->type() calls.
unsigned short	_current_elem_n_sides
	The number of sides on the current elem, used to avoid elem->n_sides() virtual calls.
Point	_incoming_point
	The incoming point of the current Ray.
unsigned short	_incoming_side
	The incoming side of the current Ray.
bool	_should_continue
	Whether or not the current Ray should continue.
Point	_intersection_point
	The work point for the intersection of the current Ray.
unsigned short	_intersected_side
	The work point for the intersected side of the current Ray.
ElemExtrema	_intersected_extrema
	The work point for the intersected vertex/edge vertices of the current Ray, if any.
ElemExtrema	_last_intersected_extrema
	The intersected vertex/edge vertices for the previous intersection, if any.
Real	_intersection_distance
	The work point for the intersection distance of the current Ray.
BoundaryID	_current_bnd_id
	The current BoundaryID (used when calling RayBoundaryConditionBase::onBoundary())
bool	_exits_elem
	Whether or not the current trace exits an element.
const Point *	_current_normals
	The normals for the current element for backface culling (pointer to the first normal - optional)
std::vector< BoundaryID >	_boundary_ids
	Reusable vector for calling _boundary_info.boundary_ids()
std::vector< TraceRayBndElement >	_boundary_elems
	Boundary elements that need RayBCs to be applied.
std::unordered_map< const Node *, std::vector< NeighborInfo > >	_vertex_neighbors
	The cached vertex neighbors.
std::unordered_map< std::pair< const Node *, const Node * >, std::pair< bool, std::vector< NeighborInfo > > >	_edge_neighbors
	The cached edge neighbors.
std::vector< NeighborInfo >	_point_neighbor_helper
	Reusable for building neighbors.
MooseUtils::StaticallyAllocatedSet< const Elem *, MAX_POINT_NEIGHBORS >	_neighbor_set
MooseUtils::StaticallyAllocatedSet< const Elem *, MAX_POINT_NEIGHBORS >	_neighbor_untested_set
MooseUtils::StaticallyAllocatedSet< const Elem *, MAX_POINT_NEIGHBORS >	_neighbor_next_untested_set
std::vector< const Elem * >	_neighbor_active_neighbor_children
std::vector< unsigned long long int >	_results
	Results over all of the local traces, indexed by TraceRayResult.
libMesh::ElemSideBuilder	_elem_side_builder
	Helper for building element sides without excessive allocation.
std::vector< RayBoundaryConditionBase * >	_on_boundary_ray_bcs
	Reusable for getting the RayBCs in onBoundary()
std::vector< std::size_t >	_on_boundary_apply_index
	Reusable for which boundary elements to apply for a specific RayBC in onBoundary()
bool	_has_ray_kernels
	Whether or not the RayTracingStudy has any RayKernels.
bool	_is_rectangular_domain
	Whether or not the domain is rectangular (defined perfectly by its bounding box)
std::set< RayKernelBase * >	_old_ray_kernels
	Helper for avoiding calling preTrace() on the same RayKernel multiple times.
Mesh *	_debug_mesh
Parallel::Communicator	_debug_comm
unsigned int	_debug_node_count

Detailed Description

Traces Rays through the mesh on a single processor.

Definition at line 46 of file TraceRay.h.

Member Enumeration Documentation

ExitsElemResult

enum TraceRay::ExitsElemResult

private

Enum for the different exit results for exitElem()

Enumerator
NO_EXIT
HIT_FACE
HIT_VERTEX
HIT_EDGE

Definition at line 257 of file TraceRay.h.

  {
    NO_EXIT = 0,    // doesn't exit
    HIT_FACE = 1,   // exits through a face
    HIT_VERTEX = 2, // exits through a vertex
    HIT_EDGE = 3    // exits through an edge
  };

TraceRayResult

enum TraceRay::TraceRayResult

Enum for the various results reported by this object.

Enumerator
INTERSECTIONS
FACE_HITS
VERTEX_HITS
EDGE_HITS
MOVED_THROUGH_NEIGHBORS
BACKFACE_CULLING_SUCCESSES
BACKFACE_CULLING_FAILURES
INTERSECTION_CALLS
VERTEX_NEIGHBOR_BUILDS
VERTEX_NEIGHBOR_LOOKUPS
EDGE_NEIGHBOR_BUILDS
EDGE_NEIGHBOR_LOOKUPS
POINT_NEIGHBOR_BUILDS
FAILED_TRACES
ENDED_STATIONARY

Definition at line 56 of file TraceRay.h.

  {
    INTERSECTIONS = 0,              // intersections
    FACE_HITS = 1,                  // intersections at an element face
    VERTEX_HITS = 2,                // intersections at an element vertex
    EDGE_HITS = 3,                  // intersections at an element edge
    MOVED_THROUGH_NEIGHBORS = 4,    // intersections through point/vertex/edge neighbors
    BACKFACE_CULLING_SUCCESSES = 5, // face intersection calls that backface culling worked on
    BACKFACE_CULLING_FAILURES = 6,  // face interesction calls that backface culling failed on
    INTERSECTION_CALLS = 7,         // face intersection calls made
    VERTEX_NEIGHBOR_BUILDS = 8,     // builds for neighbors at a vertex
    VERTEX_NEIGHBOR_LOOKUPS = 9,    // lookups for neighbors at a vertex
    EDGE_NEIGHBOR_BUILDS = 10,      // builds for neighbors on an edge
    EDGE_NEIGHBOR_LOOKUPS = 11,     // lookups for neighborson an edge
    POINT_NEIGHBOR_BUILDS = 12,     // builds for point neighbors (non edge/vertex)
    FAILED_TRACES = 13,             // rays that fail (allowed when _tolerate_failure == true),
    ENDED_STATIONARY = 14           // rays the end because they are stationary
  };

Constructor & Destructor Documentation

TraceRay()

TraceRay::TraceRay	(	RayTracingStudy &	study,
		const THREAD_ID	tid
	)

Definition at line 46 of file TraceRay.C.

  : _study(study),
    _mesh(study.getSubProblem().mesh()),
    _dim(_mesh.dimension()),
    _boundary_info(_mesh.getMesh().get_boundary_info()),
    _pid(_study.comm().rank()),
    _tid(tid),
    _backface_culling(false),
    _current_normals(nullptr),
    _results(ENDED_STATIONARY + 1)
{
}

~TraceRay()

virtual TraceRay::~TraceRay ( )

inlinevirtual

Definition at line 51 of file TraceRay.h.

{}

Member Function Documentation

applyOnExternalBoundary()

void TraceRay::applyOnExternalBoundary ( const std::shared_ptr< Ray > &  ray )

private

Gets and applies external boundary conditions in _current_elem on side _intersected_side at _intersection_point.

Definition at line 745 of file TraceRay.C.

{
  debugRay("Called applyOnExternalBoundary() with");
  debugRay("  _current_elem->id() = ", _current_elem->id());
  debugRay("  _intersection_point = ", _intersection_point);
  debugRay("  _intersected_side = ", _intersected_side);
  debugRay("  _intersected_extrema = ", _intersected_extrema);
 
  // Clear storage for the list of ConstBndElement that we need to apply RayBCs to
  _boundary_elems.clear();
 
  // If on the elem extrema (at a vertex or within an edge), check for external sidesets on
  // neighbors (which will include _current_elem)
  if (_dim != 1 && _intersected_extrema.atExtrema())
  {
    const auto & neighbors = getNeighbors(_current_elem, _intersected_extrema, _intersection_point);
    debugRay("  Found ", neighbors.size(), " vertex/edge neighbors (including self)");
    traceAssert(std::count_if(neighbors.begin(),
                              neighbors.end(),
                              [this](const NeighborInfo & ni)
                              { return ni._elem == _current_elem; }),
                "_current_elem not in neighbors");
 
    for (const auto & neighbor_info : neighbors)
    {
      if (!neighbor_info._valid)
        continue;
 
      const Elem * elem = neighbor_info._elem;
      const auto & sides = neighbor_info._sides;
      const auto & side_normals = neighbor_info._side_normals;
 
      for (MooseIndex(side_normals.size()) i = 0; i < side_normals.size(); ++i)
        if (!elem->neighbor_ptr(sides[i]) // is a boundary side that has our point
            && side_normals[i] * ray->direction() > TRACE_TOLERANCE) // and is entrant
        {
          // TODO: this could likely be optimized
          ElemExtrema extrema;
          withinExtremaOnSide(elem, _intersection_point, sides[i], _dim, extrema);
 
          _boundary_info.boundary_ids(elem, sides[i], _boundary_ids);
          possiblyAddToBoundaryElems(elem, sides[i], _boundary_ids, extrema);
        }
    }
  }
  // Not on the periphery (at vertex/edge), so we just need to add the external
  // sidesets from _current_elem on _intersected_side
  else
  {
    _boundary_info.boundary_ids(_current_elem, _intersected_side, _boundary_ids);
    possiblyAddToBoundaryElems(
        _current_elem, _intersected_side, _boundary_ids, _intersected_extrema);
  }
 
  debugRay("Calling external onBoundary() with ", _boundary_elems.size(), " boundaries");
  onBoundary(ray, /* external = */ true);
}

Referenced by trace().

applyOnInternalBoundary()

void TraceRay::applyOnInternalBoundary ( const std::shared_ptr< Ray > &  ray )

private

Gets and applies internal boundary conditions (if any) from _current_elem, _last_elem, and any other point neighbors that have internal sidesets at _intersection_point.

Definition at line 804 of file TraceRay.C.

{
  traceAssert(_last_elem, "Must be valid");
 
  debugRay("Called applyOnInternalBoundary() with");
  debugRay("  intersection_point = ", _intersection_point);
  debugRay("  _current_elem->id() = ", _current_elem->id());
  debugRay("  _incoming_side = ", _incoming_side);
  debugRay("  _last_elem->id() = ", _last_elem->id());
  debugRay("  _intersected_side = ", _intersected_side);
  debugRay("  _intersected_extrema = ", _intersected_extrema);
  traceAssert(_study.hasInternalSidesets(), "Do not have internal sidesets");
  traceAssert(_intersection_point.absolute_fuzzy_equals(_incoming_point, TRACE_TOLERANCE),
              "Intersection and incoming points should be the same");
 
  // Clear storage for the list of ConstBndElement that we need to apply RayBCs to
  _boundary_elems.clear();
 
  ElemExtrema temp_extrema;
 
  // If on the elem extrema (at a vertex or within an edge), we need to check for internal
  // sidesets on neighbors (which will include _last_elem and _current_elem)
  if (_dim != 1 && _intersected_extrema.atExtrema())
  {
    debugRay("Checking point neighbors for internal sidesets");
 
    // Get the neighbors
    const auto & neighbors = getNeighbors(_last_elem, _intersected_extrema, _intersection_point);
    debugRay("  Found ", neighbors.size(), " vertex/edge neighbors");
    traceAssert(std::count_if(neighbors.begin(),
                              neighbors.end(),
                              [this](const NeighborInfo & ni)
                              { return ni._elem == _current_elem || ni._elem == _last_elem; }) == 2,
                "_current_elem/_last_elem not in neighbors");
 
    for (const auto & neighbor_info : neighbors)
    {
      if (!neighbor_info._valid)
        continue;
 
      const Elem * elem = neighbor_info._elem;
 
      // Grab the internal sidesets for this elem
      const auto & sidesets = _study.getInternalSidesets(elem);
      // It has none to contribute, so we can continue
      if (sidesets.empty())
      {
        debugRay("    Elem ", elem->id(), " has no internal sidesets");
        continue;
      }
 
      // The sides on this elem that contain our point and their outward normals
      const auto & sides = neighbor_info._sides;
      const auto & side_normals = neighbor_info._side_normals;
 
      // See if any of the internal sidesets are relevant to this point
      for (std::size_t i = 0; i < sides.size(); ++i)
      {
        const auto side = sides[i];
        // Side has internal sidesets and is entrant
        if (sidesets[side].size() && std::abs(side_normals[i] * ray->direction()) > TRACE_TOLERANCE)
        {
          // TODO: this could likely be optimized
          temp_extrema.invalidate();
          withinExtremaOnSide(elem, _intersection_point, side, _dim, temp_extrema);
 
          possiblyAddToBoundaryElems(elem, side, sidesets[side], temp_extrema);
        }
      }
    }
  }
  // Not on the periphery (at vertex/edge), so we just need to add the sidesets from _current_elem
  // on _incoming_side and _last_elem on _intersected_side
  else
  {
    const auto & current_elem_sidesets = _study.getInternalSidesets(_current_elem);
    if (current_elem_sidesets.size() && current_elem_sidesets[_incoming_side].size())
    {
      // This is possible but we need extrema checks on _incoming_elem to pass to the
      // boundary conditions. For now, we just pass in _intersected_extrema which could
      // be very wrong with adaptivity but is correct without it. I struggled with getting
      // the actual extrema checks and we're not using this now so I bet this error will
      // come back to haunt me in the future :-)
      if (!_study.hasSameLevelActiveElems())
        mooseError("Internal sidesets are not currently supported with adaptivity in tracing");
 
      // Special case for 1D
      if (_dim == 1)
        temp_extrema.setVertex(atVertex(_current_elem, _intersection_point));
 
      possiblyAddToBoundaryElems(_current_elem,
                                 _incoming_side,
                                 current_elem_sidesets[_incoming_side],
                                 _dim != 1 ? _intersected_extrema : temp_extrema);
    }
 
    const auto & last_elem_sidesets = _study.getInternalSidesets(_last_elem);
    if (last_elem_sidesets.size() && last_elem_sidesets[_intersected_side].size())
      possiblyAddToBoundaryElems(_last_elem,
                                 _intersected_side,
                                 last_elem_sidesets[_intersected_side],
                                 _intersected_extrema);
  }
 
  if (!_boundary_elems.empty())
  {
    debugRay("  Calling internal onBoundary() with ", _boundary_elems.size(), " boundaries");
    onBoundary(ray, /* external = */ false);
  }
}

Referenced by trace().

continueTraceOffProcessor()

void TraceRay::continueTraceOffProcessor ( const std::shared_ptr< Ray > &  ray )

private

Sets up a ray to continue tracing off processor.

Parameters

ray	The ray
elem	The element that is owned by another processor
incoming_side	The incoming side elem
point	The point on elem at which ray continues

Definition at line 1696 of file TraceRay.C.

{
  traceAssert(ray->currentElem() == _current_elem, "Ray currentElem() invalid");
  traceAssert(ray->currentIncomingSide() == _incoming_side, "Ray currentIncomingSide() invalid");
  traceAssert(ray->currentPoint() == _incoming_point, "Ray currentPoint() invalid");
  traceAssert(_current_elem->processor_id() != _pid, "Off processor trace is not off processor");
  debugRay("Ray going off processor to ", _current_elem->processor_id());
 
  ray->addProcessorCrossing();
 
  if (_current_cached_trace && _intersection_distance > 0)
  {
    _current_cached_trace->addPoint(_incoming_point);
    if (_study.dataOnCacheTraces())
      _current_cached_trace->lastPoint()._data = ray->data();
    if (_study.auxDataOnCacheTraces())
      _current_cached_trace->lastPoint()._aux_data = ray->auxData();
  }
 
  if (_intersection_distance > 0)
    possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
  possiblySaveDebugRayMesh();
}

Referenced by trace().

currentBoundaryID()

const BoundaryID & TraceRay::currentBoundaryID ( ) const

inline

Gets the BoundaryID of the boundary that the Ray intersected and is being applied a boundary condition.

Definition at line 155 of file TraceRay.h.

{ return _current_bnd_id; }

currentElem()

const Elem *const & TraceRay::currentElem ( ) const

inline

Gets the element that the current Ray is being traced in.

Definition at line 114 of file TraceRay.h.

{ return _current_elem; }

currentIncomingPoint()

const Point & TraceRay::currentIncomingPoint ( ) const

inline

Gets the current incoming point for the Ray that is being traced.

This is ONLY guaranteed to be valid during execution of RayKernels!

Definition at line 127 of file TraceRay.h.

{ return _incoming_point; }

currentIncomingSide()

const unsigned short & TraceRay::currentIncomingSide ( ) const

inline

Gets the current incoming side for the Ray that is being traced.

This is only valid during onSegment() (when RayKernels are being called)

Definition at line 133 of file TraceRay.h.

{ return _incoming_side; }

currentIntersectedExtrema()

const ElemExtrema & TraceRay::currentIntersectedExtrema ( ) const

inline

Gets the element extrema (vertex/edge) that the current Ray intersected.

Definition at line 144 of file TraceRay.h.

{ return _intersected_extrema; }

currentIntersectedSide()

const unsigned short & TraceRay::currentIntersectedSide ( ) const

inline

Gets the side that the current Ray intersected.

If the ray ends within the current element, this will be RayTracingCommon::invalid_side.

Definition at line 140 of file TraceRay.h.

{ return _intersected_side; }

currentIntersectionDistance()

const Real & TraceRay::currentIntersectionDistance ( ) const

inline

Gets the current intersection distance for the Ray that is being traced.

When within an element (executing with a RayKernel), this is the segment distance.

Definition at line 150 of file TraceRay.h.

{ return _intersection_distance; }

currentIntersectionPoint()

const Point & TraceRay::currentIntersectionPoint ( ) const

inline

Gets the current intersection point for the Ray that is being traced.

If the Ray ends within the current element, this will be the point at which the Ray ends within the current element.

Definition at line 121 of file TraceRay.h.

{ return _intersection_point; }

currentRay()

const std::shared_ptr< Ray > *const & TraceRay::currentRay ( ) const

inline

Gets the current Ray that is being traced.

Definition at line 110 of file TraceRay.h.

{ return _current_ray; }

currentSubdomainID()

const SubdomainID & TraceRay::currentSubdomainID ( ) const

inline

Gets the subdomain of the current element that the Ray is being traced in.

Definition at line 159 of file TraceRay.h.

{ return _current_subdomain_id; }

exitsElem() [1/3]

TraceRay::ExitsElemResult TraceRay::exitsElem ( const Elem *  elem, const ElemType  elem_type, const unsigned short  incoming_side, Point &  intersection_point, unsigned short &  intersected_side, ElemExtrema &  intersected_extrema, Real &  intersection_distance, const Point *  normals  )

private

Determines if _current_ray moving in direction _direction exits elem.

Parameters

elem	The element to check
elem_type	The type of elem (used to avoid virtual function calls on elem)
incoming_side	The incoming side (if any)
intersection_point	To be modified with the intersection point, if any
intersected_side	To be modified with the intersected side, if any
intersected_extrema	To be modified with the intersected extremum (vertex/edge), if any
intersection_distance	To be modified with the intersection distance, if any
normals	Pointer to the normals for the current element, if any

Returns

Whether or not _current_ray exits elem

Definition at line 82 of file TraceRay.C.

{
  debugRay("Called exitsElem()");
 
  traceAssert(elem_type == elem->type(), "elem_type incorrect");
  traceAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  traceAssert(intersected_side == RayTracingCommon::invalid_side, "Side should be invalid");
  traceAssert(intersected_extrema.isInvalid(), "Extrema should be invalid");
  traceAssert(intersection_distance == RayTracingCommon::invalid_distance,
              "Distance should be invalid");
  if (_study.verifyRays() && incoming_side != RayTracingCommon::invalid_side &&
      !_study.sideIsNonPlanar(elem, incoming_side))
    traceAssert(_study.sideIsIncoming(elem, incoming_side, (*_current_ray)->direction(), _tid),
                "Incoming side is non-entrant");
 
  bool intersected;
  switch (elem_type)
  {
    case HEX8:
      intersected = exitsElem<Hex8, Hex8>(elem,
                                          incoming_side,
                                          intersection_point,
                                          intersected_side,
                                          intersected_extrema,
                                          intersection_distance,
                                          normals);
      break;
    case TET4:
      intersected = exitsElem<Tet4, Tet4>(elem,
                                          incoming_side,
                                          intersection_point,
                                          intersected_side,
                                          intersected_extrema,
                                          intersection_distance,
                                          normals);
      break;
    case PYRAMID5:
      intersected = exitsElem<Pyramid5, Pyramid5>(elem,
                                                  incoming_side,
                                                  intersection_point,
                                                  intersected_side,
                                                  intersected_extrema,
                                                  intersection_distance,
                                                  normals);
      break;
    case PRISM6:
      intersected = exitsElem<Prism6, Prism6>(elem,
                                              incoming_side,
                                              intersection_point,
                                              intersected_side,
                                              intersected_extrema,
                                              intersection_distance,
                                              normals);
      break;
    case QUAD4:
      intersected = exitsElem<Quad4, Quad4>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    case TRI3:
      intersected = exitsElem<Tri3, Tri3>(elem,
                                          incoming_side,
                                          intersection_point,
                                          intersected_side,
                                          intersected_extrema,
                                          intersection_distance,
                                          normals);
      break;
    case HEX20:
      intersected = exitsElem<Hex20, Hex8>(elem,
                                           incoming_side,
                                           intersection_point,
                                           intersected_side,
                                           intersected_extrema,
                                           intersection_distance,
                                           normals);
      break;
    case HEX27:
      intersected = exitsElem<Hex27, Hex8>(elem,
                                           incoming_side,
                                           intersection_point,
                                           intersected_side,
                                           intersected_extrema,
                                           intersection_distance,
                                           normals);
      break;
    case QUAD8:
      intersected = exitsElem<Quad8, Quad4>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    case QUAD9:
      intersected = exitsElem<Quad9, Quad4>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    case TRI6:
      intersected = exitsElem<Tri6, Tri3>(elem,
                                          incoming_side,
                                          intersection_point,
                                          intersected_side,
                                          intersected_extrema,
                                          intersection_distance,
                                          normals);
      break;
    case TRI7:
      intersected = exitsElem<Tri7, Tri3>(elem,
                                          incoming_side,
                                          intersection_point,
                                          intersected_side,
                                          intersected_extrema,
                                          intersection_distance,
                                          normals);
      break;
    case TET10:
      intersected = exitsElem<Tet10, Tet4>(elem,
                                           incoming_side,
                                           intersection_point,
                                           intersected_side,
                                           intersected_extrema,
                                           intersection_distance,
                                           normals);
      break;
    case TET14:
      intersected = exitsElem<Tet14, Tet4>(elem,
                                           incoming_side,
                                           intersection_point,
                                           intersected_side,
                                           intersected_extrema,
                                           intersection_distance,
                                           normals);
      break;
    case PYRAMID13:
      intersected = exitsElem<Pyramid13, Pyramid5>(elem,
                                                   incoming_side,
                                                   intersection_point,
                                                   intersected_side,
                                                   intersected_extrema,
                                                   intersection_distance,
                                                   normals);
      break;
    case PYRAMID14:
      intersected = exitsElem<Pyramid14, Pyramid5>(elem,
                                                   incoming_side,
                                                   intersection_point,
                                                   intersected_side,
                                                   intersected_extrema,
                                                   intersection_distance,
                                                   normals);
      break;
    case PRISM15:
      intersected = exitsElem<Prism15, Prism6>(elem,
                                               incoming_side,
                                               intersection_point,
                                               intersected_side,
                                               intersected_extrema,
                                               intersection_distance,
                                               normals);
      break;
    case PRISM18:
      intersected = exitsElem<Prism18, Prism6>(elem,
                                               incoming_side,
                                               intersection_point,
                                               intersected_side,
                                               intersected_extrema,
                                               intersection_distance,
                                               normals);
      break;
    case EDGE2:
      intersected = exitsElem<Edge2, Edge2>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    case EDGE3:
      intersected = exitsElem<Edge3, Edge2>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    case EDGE4:
      intersected = exitsElem<Edge4, Edge2>(elem,
                                            incoming_side,
                                            intersection_point,
                                            intersected_side,
                                            intersected_extrema,
                                            intersection_distance,
                                            normals);
      break;
    default:
      mooseError(
          "Element type ", Utility::enum_to_string(elem->type()), " not supported by TraceRay");
  }
 
  if (intersected)
  {
    if (intersected_extrema.atExtrema())
      return intersected_extrema.atVertex() ? HIT_VERTEX : HIT_EDGE;
    return HIT_FACE;
  }
 
  return NO_EXIT;
}

Referenced by moveThroughNeighbor(), and trace().

exitsElem() [2/3]

template<typename T , typename FirstOrderT >

std::enable_if< std::is_base_of< Edge, T >::value, bool >::type TraceRay::exitsElem ( const Elem *  elem, const unsigned short  incoming_side, Point &  intersection_point, unsigned short &  intersected_side, ElemExtrema &  intersected_extrema, Real &  intersection_distance, const Point *  normals  )

private

Determines if _current_ray moving in direction _direction exits elem.

This function is templated to work only with elements that are faces or cells (2D or 3D)

Parameters

elem	The element to check
incoming_side	The incoming side (if any)
intersection_point	To be modified with the intersection point, if any
intersected_side	To be modified with the intersected side, if any
intersected_extrema	To be modified with the intersected extremum (vertex/edge), if any
intersection_distance	To be modified with the intersection distance, if any
normals	Pointer to the normals for the current element, if any

Returns

Whether or not _current_ray exits elem through a face

Definition at line 313 of file TraceRay.C.

{
  ++_results[INTERSECTION_CALLS];
 
  debugRay("Called exitsElem() in 2D or 3D");
 
  const auto hmax = subdomainHmax(elem);
 
  // Current point and distance (maybe not the best!)
  Point current_intersection_point;
  Real current_intersection_distance;
  ElemExtrema current_intersected_extrema;
  // Scale the minimum intersection distance required by hmax
  Real best_intersection_distance = TRACE_TOLERANCE * hmax;
  // Whether or not we are going to try backface culling the first time around
  // This depends on if the user set it to use culling
  bool use_backface_culling = _backface_culling;
  // If all sides have failed to find an intersection, whether or not we
  // are going to see if the nonplanar skip side (only if it IS nonplanar)
  // is also an exiting point
  bool try_nonplanar_incoming_side = false;
  // Easy access into the current direction
  const auto & direction = (*_current_ray)->direction();
  // The current side that we're checking
  unsigned short s = 0;
 
  // Loop over the side loop. We need this in the cases that:
  // - Backface culling is enabled and we were not able to find an intersection
  //   so we will go through all culled sides and see if we can find one
  // - The incoming_side is non-planar and the ray also exits said side,
  //   which will be checked after all other sides fail
  while (true)
  {
    debugRay("  use_backface_culling = ", use_backface_culling);
    debugRay("  try_nonplanar_incoming_side = ", try_nonplanar_incoming_side);
    // Loop over all of the sides
    while (true)
    {
      debugRay("  Side ", s, " with centroid ", _elem_side_builder(*elem, s).vertex_average());
 
      // All of the checks that follow are done while we're still searching through
      // all of the sides. try_nonplanar_incoming_side is our last possible check
      // and does not involve looping through all of the sides, so skip these
      // checks if we're at that point.
      if (!try_nonplanar_incoming_side)
      {
        // Don't search backwards. If we have a non-planar incoming side, we will
        // check it if all other sides have failed
        if (s == incoming_side)
        {
          debugRay("    Skipping due to incoming side");
          if (++s == T::num_sides)
            break;
          else
            continue;
        }
 
        // Using backface culling on this run through the sides
        // If the direction is non-entrant, skip this side
        if (use_backface_culling)
        {
          // Side is non-entrant per the culling, so skip
          if (normals[s] * direction < -LOOSE_TRACE_TOLERANCE)
          {
            debugRay("    Skipping due to backface culling dot = ", normals[s] * direction);
 
            if (++s == T::num_sides)
              break;
            else
              continue;
 
            ++_results[BACKFACE_CULLING_SUCCESSES];
          }
        }
        // We're not using backface culling but it is enabled, which means
        // we're on our second run through the sides because we could not
        // find any intersections while culling the sides. Try again and
        // check the sides that we previously skipped for intersections
        else if (_backface_culling)
        {
          // Side was non-entrant per the culling, so try again
          if (normals[s] * direction >= -LOOSE_TRACE_TOLERANCE)
          {
            debugRay("    Skipping because we already checked this side with culling enabled");
            if (++s == T::num_sides)
              break;
            else
              continue;
          }
          else
            debugRay("    Side that was skipped due to culling");
 
          ++_results[BACKFACE_CULLING_FAILURES];
        }
      }
 
      // Look for an intersection!
      current_intersection_point = RayTracingCommon::invalid_point;
      current_intersected_extrema.invalidate();
      const bool intersected = sideIntersectedByLine<FirstOrderT>(elem,
                                                                  _incoming_point,
                                                                  direction,
                                                                  s,
                                                                  _study.domainMaxLength(),
                                                                  current_intersection_point,
                                                                  current_intersection_distance,
                                                                  current_intersected_extrema,
                                                                  hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                                                  ,
                                                                  DEBUG_RAY_IF
#endif
      );
 
      // Do they intersect and is it further down the path than any other intersection?
      // If so, keep track of the intersection with the furthest distance
      if (intersected)
      {
        debugRay("    Intersected at point ",
                 current_intersection_point,
                 " with distance ",
                 current_intersection_distance);
        debugRay("    Best intersection distance = ", best_intersection_distance);
 
#ifndef NDEBUG
        // Only validate intersections if the side is planar
        if (_study.verifyTraceIntersections() && !_study.sideIsNonPlanar(elem, s) &&
            !_elem_side_builder(*elem, s).contains_point(current_intersection_point))
          failTrace("Intersected side does not contain intersection point",
                    _study.tolerateFailure(),
                    __LINE__);
#endif
 
        // The intersection we just computed is further than any other intersection
        // that was found so far - mark it as the best
        if (current_intersection_distance > best_intersection_distance)
        {
          debugRay("    Best intersection so far");
 
          intersected_side = s;
          intersection_distance = current_intersection_distance;
          intersection_point = current_intersection_point;
          intersected_extrema = current_intersected_extrema;
          best_intersection_distance = current_intersection_distance;
        }
      }
 
      if (++s == T::num_sides || try_nonplanar_incoming_side)
        break;
      else
        continue;
    } // while(true)
 
    // Found an intersection
    if (intersected_side != RayTracingCommon::invalid_side)
    {
      debugRay("  Exiting with intersection");
      debugRay("    intersected_side = ", intersected_side);
      debugRay("    intersection_distance = ", intersection_distance);
      debugRay("    intersection_point = ", intersection_point);
      debugRay("    intersected_extrema = ", intersected_extrema);
 
      return true;
    }
 
    // This was our last possible check - no dice!
    if (try_nonplanar_incoming_side)
      return false;
 
    // We didn't find an intersection but we used backface culling. Try again without
    // it, checking only the sides that we skipped due to culling
    if (use_backface_culling)
    {
      debugRay("  Didn't find an intersection, retrying without backface culling");
      use_backface_culling = false;
      s = 0;
      continue;
    }
 
    // Have tried all sides (potentially with and without culling). If the incoming
    // side is valid is non planar, see if we also exit out of it
    if (_incoming_side != RayTracingCommon::invalid_side &&
        _study.sideIsNonPlanar(elem, incoming_side))
    {
      debugRay("  Didn't find an intersection, trying non-planar incoming_side");
      try_nonplanar_incoming_side = true;
      s = incoming_side;
      continue;
    }
 
    // No intersection found!
    return false;
  } // while(true)
}

exitsElem() [3/3]

template<typename T , typename FirstOrderT >

std::enable_if< std::is_base_of< Edge, T >::value, bool >::type TraceRay::exitsElem ( const Elem *  elem, const unsigned short  incoming_side, Point &  intersection_point, unsigned short &  intersected_side, ElemExtrema &  intersected_extrema, Real &  intersection_distance, const Point *  normals  )

private

Determines if _current_ray moving in direction _direction exits elem.

This function is templated to work only with elements that are edges (1D).

Parameters

elem	The element to check
incoming_side	The incoming side (if any)
intersection_point	To be modified with the intersection point, if any
intersected_side	To be modified with the intersected side, if any
intersected_extrema	To be modified with the intersected extremum (vertex/edge), if any
intersection_distance	To be modified with the intersection distance, if any
normals	Pointer to the normals for the current element, if any

Returns

Whether or not _current_ray exits elem through a face

failTrace()

void TraceRay::failTrace ( const std::string &  reason, const bool  warning, const int  line = -1  )

private

Specialized mooseError for a failed Ray trace with detailed information regarding the trace.

If warning = true, use mooseWarning instead. If line is provided, output will include the line number.

Definition at line 1832 of file TraceRay.C.

{
  const auto message = failTraceMessage(reason, line);
 
  if (warning)
  {
    ++_results[FAILED_TRACES];
    mooseWarning(message);
    (*_current_ray)->setShouldContinue(false);
    _should_continue = false;
  }
  else
    mooseError(message);
}

Referenced by exitsElem(), onBoundary(), onTrajectoryChanged(), postRayTracingObject(), and trace().

failTraceMessage()

std::string TraceRay::failTraceMessage ( const std::string &  reason, const int  line = -1  )

private

Creates a useful error string with current tracing information.

If line is provided, will include the line number.

Definition at line 1784 of file TraceRay.C.

{
  std::stringstream oss;
  oss << "Ray on processor " << _pid << " and thread " << _tid << " failed to trace";
  if (line != -1)
    oss << " at line " << line;
  oss << "\n\n" << reason << "\n\n";
  oss << ((*_current_ray))->getInfo() << "\n";
  oss << "Current trace information\n";
  oss << "  _current_subdomain_id = ";
  if (_current_subdomain_id == Elem::invalid_subdomain_id)
    oss << "invalid subdomain id\n";
  else
    oss << _current_subdomain_id << "\n";
  oss << "  _current_elem_type = " << Utility::enum_to_string(_current_elem_type) << "\n";
  oss << "  _current_elem_n_sides = " << _current_elem_n_sides << "\n";
  oss << "  _incoming_point = ";
  if (_incoming_point == RayTracingCommon::invalid_point)
    oss << "invalid point\n";
  else
    oss << _incoming_point << "\n";
  oss << "  _incoming_side = ";
  if (_incoming_side == RayTracingCommon::invalid_side)
    oss << "invalid side\n";
  else
    oss << _incoming_side << "\n";
  oss << "  _intersection_point = ";
  if (_intersection_point == RayTracingCommon::invalid_point)
    oss << "invalid point\n";
  else
    oss << _intersection_point << "\n";
  oss << "  _intersected_side = ";
  if (_intersected_side == RayTracingCommon::invalid_side)
    oss << "invalid side\n";
  else
    oss << _intersected_side << "\n";
  oss << "  _intersected_extrema = " << _intersected_extrema << "\n";
  oss << "  _exits_elem = " << _exits_elem << "\n";
  if (_current_elem)
    oss << _current_elem->get_info();
  else
    oss << "_current_elem = invalid\n";
 
  possiblySaveDebugRayMesh();
  return oss.str();
}

Referenced by failTrace().

findExternalBoundarySide()

void TraceRay::findExternalBoundarySide ( unsigned short &  boundary_side, ElemExtrema &  boundary_extrema, const Elem *&  boundary_elem  )

private

Finds (if any) an element side that is on the boundary and is outgoing at _intersection_point that is on the extrema _intersected_extrema for _current_elem.

This is necessary when a Ray hits a point that is on the boundary but is not actually on a boundary side

Fills the found boundary side into boundary_side, the element said side is on into boundary_elem and the extrema on said elem into boundary_extrema

Definition at line 949 of file TraceRay.C.

{
  traceAssert(_current_elem->neighbor_ptr(_intersected_side), "Already on boundary");
  traceAssert(boundary_side == RayTracingCommon::invalid_side, "Side should be invalid");
  traceAssert(boundary_extrema.isInvalid(), "Extrema should be invalid");
  traceAssert(!boundary_elem, "Elem should be invalid");
  traceAssert(_current_elem->dim() != 1, "1D traces shouldn't make it here");
  traceAssert(_current_elem->n_sides() == _current_elem_n_sides, "_current_elem_n_sides incorrect");
  traceAssert(_intersected_extrema.atExtrema(), "Should be at extrema");
  debugRay("Called findExternalBoundarySide() on side ", _intersected_side);
  debugRay("  _intersected_side = ", _intersected_side);
  debugRay("  _intersected_extrema", _intersected_extrema);
 
  const auto & direction = (*_current_ray)->direction();
  const auto at_edge = _intersected_extrema.atEdge();
 
  // First, look for other sides on _current_elem that touch the intersected vertex/edge
  // that are on the boundary and are outgoing
  for (unsigned short s = 0; s < _current_elem_n_sides; ++s)
    if (!_current_elem->neighbor_ptr(s) && s != _intersected_side &&
        _current_elem->is_node_on_side(_intersected_extrema.first, s) &&
        (!at_edge || _current_elem->is_node_on_side(_intersected_extrema.second, s)) &&
        !_study.sideIsIncoming(_current_elem, s, direction, _tid))
    {
      debugRay("  Side ", s, " is a boundary side and the Ray exits");
      boundary_side = s;
      boundary_extrema = _intersected_extrema;
      boundary_elem = _current_elem;
      return;
    }
 
  // No luck in our element, so see if any neighbors at this vertex/edge are
  // on the boundary and are outgoing
  const auto & neighbors = getNeighbors(_current_elem, _intersected_extrema, _intersection_point);
  debugRay("Checking current element failed, now checking neighbors");
 
  debugRay("Found ", neighbors.size(), " candidate neighbors (including self)");
  for (const auto & neighbor_info : neighbors)
  {
    // This will be false when we have an edge neighbor that isn't a neighbor at
    // _intersection_point
    if (!neighbor_info._valid)
      continue;
 
    const Elem * neighbor = neighbor_info._elem;
    // We've already checked ourself
    if (neighbor == _current_elem)
      continue;
 
    debugRay("Checking neighbor ", neighbor->id());
 
    // Loop through the sides we touch and look for one that this Ray exits and is on the boundary
    for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
      if (neighbor_info._side_normals[i] * direction > TRACE_TOLERANCE &&
          !neighbor->neighbor_ptr(neighbor_info._sides[i]))
      {
        withinExtremaOnSide(
            neighbor, _intersection_point, neighbor_info._sides[i], _dim, boundary_extrema);
        traceAssert(boundary_extrema.atExtrema(), "Should be at extrema");
        boundary_side = neighbor_info._sides[i];
        boundary_elem = neighbor;
        return;
      }
  }
}

Referenced by trace().

getEdgeNeighbors() [1/2]

const std::vector< NeighborInfo > & TraceRay::getEdgeNeighbors	(	const Elem *	elem,
		const std::pair< const Node *, const Node * > &	vertices,
		const Point &	point
	)

Get the neighbors at an edge.

Parameters

elem	An elem that contains the edge
vertices	A pair of Nodes that are the vertices that contain the edge fully
point	The point on the edge where neighbors are desired

Definition at line 1900 of file TraceRay.C.

{
  traceAssert(elem, "Invalid elem");
  traceAssert(vertices.first, "Must be valid");
  traceAssert(vertices.second, "Must be valid");
 
  debugRay("Called getEdgeNeighbors() with:");
  debugRay("  elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
  debugRay("  vertices.first = ", vertices.first->id(), " at ", (Point)*vertices.first);
  debugRay("  vertices.second = ", vertices.second->id(), " at ", (Point)*vertices.second);
  debugRay("  point = ", point);
 
  traceAssert(elem->get_node_index(vertices.first) != libMesh::invalid_uint,
              "Doesn't contain vertex");
  traceAssert(elem->get_node_index(vertices.second) != libMesh::invalid_uint,
              "Doesn't contain vertex");
  traceAssert(isWithinSegment(
                  (Point)*vertices.first, (Point)*vertices.second, point, LOOSE_TRACE_TOLERANCE),
              "Point not within edge");
 
  ++_results[EDGE_NEIGHBOR_LOOKUPS];
 
  const auto ordered_vertices = vertices.first->id() < vertices.second->id()
                                    ? vertices
                                    : std::make_pair(vertices.second, vertices.first);
 
  // Look for the entry and build if necessary
  std::pair<bool, std::vector<NeighborInfo>> * entry;
  auto search = _edge_neighbors.find(ordered_vertices);
  if (search != _edge_neighbors.end())
    entry = &search->second;
  else
  {
    debugRay("Building edge neighbors");
    ++_results[EDGE_NEIGHBOR_BUILDS];
    entry = &_edge_neighbors
                 .emplace(ordered_vertices, std::make_pair(true, std::vector<NeighborInfo>()))
                 .first->second;
    findEdgeNeighbors(elem,
                      ordered_vertices.first,
                      ordered_vertices.second,
                      _neighbor_set,
                      _neighbor_untested_set,
                      _neighbor_next_untested_set,
                      _neighbor_active_neighbor_children,
                      entry->second);
 
    bool all_same_edge = true;
    for (auto & neighbor_info : entry->second)
    {
      traceAssert(neighbor_info._lower_bound <= neighbor_info._upper_bound,
                  "Bound order incorrect");
 
      // Fill the side normals
      for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
        neighbor_info._side_normals[i] =
            _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
 
      // See if the bounds are the same as the target edge
      if (neighbor_info._lower_bound != 0 || neighbor_info._upper_bound != 1)
        all_same_edge = false;
    }
    entry->first = all_same_edge;
  }
 
  // Means that all neighbors are not on the exact same edge, so we must
  // validate/invalidate based on if the neighbor's edge contains our point
  if (!entry->first)
  {
    const auto edge_length =
        ((Point)*ordered_vertices.first - (Point)*ordered_vertices.second).norm();
    const auto point_location = ((Point)*ordered_vertices.first - point).norm() / edge_length;
    for (auto & info : entry->second)
      info._valid = (info._lower_bound - TRACE_TOLERANCE) < point_location &&
                    point_location < (info._upper_bound + TRACE_TOLERANCE);
  }
 
  return entry->second;
}

Referenced by getEdgeNeighbors(), getNeighbors(), and trace().

getEdgeNeighbors() [2/2]

const std::vector< NeighborInfo > & TraceRay::getEdgeNeighbors	(	const Elem *	elem,
		const std::pair< unsigned short, unsigned short > &	vertices,
		const Point &	point
	)

Get the neighbors at an edge.

Parameters

elem	An elem that contains the edge
vertices	A pair of local vertex IDs of the vertices that contain the edge fully
point	The point on the edge where neighbors are desired

Definition at line 1983 of file TraceRay.C.

{
  debugRay("Called getEdgeNeighbors(), local index version with:");
  debugRay("  vertices.first = ", vertices.first);
  debugRay("  vertices.second = ", vertices.second);
  traceAssert(vertices.first < elem->n_vertices(),
              "Invalid vertex with ray " + std::to_string((*_current_ray)->id()));
  traceAssert(vertices.second < elem->n_vertices(), "Invalid vertex");
 
  return getEdgeNeighbors(
      elem, std::make_pair(elem->node_ptr(vertices.first), elem->node_ptr(vertices.second)), point);
}

getNeighbors()

const std::vector< NeighborInfo > & TraceRay::getNeighbors	(	const Elem *	elem,
		const ElemExtrema &	extrema,
		const Point &	point
	)

Get the point/vertex/edge neighbors depending on extrema.

Definition at line 1999 of file TraceRay.C.

{
  if (!extrema.atExtrema())
    return getPointNeighbors(elem, point);
  if (extrema.atVertex())
    return getVertexNeighbors(elem, extrema.vertex());
  return getEdgeNeighbors(elem, extrema.edgeVertices(), point);
}

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), findExternalBoundarySide(), and trace().

getPointNeighbors()

const std::vector< NeighborInfo > & TraceRay::getPointNeighbors	(	const Elem *	elem,
		const Point &	point
	)

Get the point neighbors.

Parameters

elem	The elem
point	The point

Definition at line 2009 of file TraceRay.C.

{
  traceAssert(elem, "Invalid elem");
 
  debugRay("Called getPointNeighbors()");
  debugRay(" elem = ", elem->id());
  debugRay(" point = ", point);
 
  ++_results[POINT_NEIGHBOR_BUILDS];
  _point_neighbor_helper.clear();
 
  findPointNeighbors(elem,
                     point,
                     _neighbor_set,
                     _neighbor_untested_set,
                     _neighbor_next_untested_set,
                     _neighbor_active_neighbor_children,
                     _point_neighbor_helper);
 
  // Fill the side normals
  for (auto & neighbor_info : _point_neighbor_helper)
    for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
      neighbor_info._side_normals[i] =
          _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
 
  return _point_neighbor_helper;
}

Referenced by getNeighbors(), and trace().

getVertexNeighbors() [1/2]

const std::vector< NeighborInfo > & TraceRay::getVertexNeighbors	(	const Elem *	elem,
		const Node *	vertex
	)

Gets the neighbors at a vertex.

Parameters

elem	An elem that contains the vertex
vertex	The Node that is the vertex

Definition at line 1848 of file TraceRay.C.

{
  traceAssert(elem, "Elem must be valid");
  traceAssert(vertex, "Vertex must be valid");
 
  debugRay("Called getVertexNeighbors() with:");
  debugRay("  elem->id() = ", elem->id(), " with centroid ", elem->vertex_average());
  debugRay("  vertex->id() = ", vertex->id(), ", at ", (Point)*vertex);
 
  traceAssert(elem->get_node_index(vertex) != libMesh::invalid_uint, "Doesn't contain node");
  traceAssert(elem->is_vertex(elem->get_node_index(vertex)), "Node is not a vertex");
 
  ++_results[VERTEX_NEIGHBOR_LOOKUPS];
 
  // Return the entry if we have it cached
  auto search = _vertex_neighbors.find(vertex);
  if (search != _vertex_neighbors.end())
    return search->second;
 
  ++_results[VERTEX_NEIGHBOR_BUILDS];
 
  // Make a new entry
  debugRay("Building vertex neighbors");
  std::vector<NeighborInfo> & entry =
      _vertex_neighbors.emplace(vertex, std::vector<NeighborInfo>()).first->second;
 
  findNodeNeighbors(elem,
                    vertex,
                    _neighbor_set,
                    _neighbor_untested_set,
                    _neighbor_next_untested_set,
                    _neighbor_active_neighbor_children,
                    entry);
 
  // Fill the side normals
  for (auto & neighbor_info : entry)
    for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
      neighbor_info._side_normals[i] =
          _study.getSideNormal(neighbor_info._elem, neighbor_info._sides[i], _tid);
 
  return entry;
}

Referenced by getNeighbors(), getVertexNeighbors(), and trace().

getVertexNeighbors() [2/2]

const std::vector< NeighborInfo > & TraceRay::getVertexNeighbors	(	const Elem *	elem,
		const unsigned short	vertex
	)

Gets the neighbors at a vertex.

Parameters

elem	An elem that contains the vertex
vertex	The local ID of the vertex on elem

Definition at line 1892 of file TraceRay.C.

{
  traceAssert(vertex < elem->n_vertices(), "Invalid vertex");
 
  return getVertexNeighbors(elem, elem->node_ptr(vertex));
}

meshChanged()

void TraceRay::meshChanged

(

)

Called on mesh change.

Definition at line 74 of file TraceRay.C.

{
  // Invalidate the vertex and edge neighbor caches
  _vertex_neighbors.clear();
  _edge_neighbors.clear();
}

moveThroughNeighbor()

TraceRay::ExitsElemResult TraceRay::moveThroughNeighbor ( const NeighborInfo &  neighbor_info, unsigned short &  incoming_side, Point &  intersection_point, unsigned short &  intersected_side, ElemExtrema &  intersected_extrema, Real &  intersection_distance  )

private

Sees if a Ray can move through a neighbor (vertex/edge/point)

Parameters

neighbor_info	The NeighborInfo for the neighbor
incoming_side	To be filled with the incoming side on the neighbor, if any
intersection_point	To be filled with the intersection point on the neighbor, if any
intersected_side	To be filled with the intersected side on the neighbor, if any
intersected_extrema	To be filled with the intersected extrema on the neighbor, if any
intersection_distance	To be filled with the intersection distance, if any

Returns

The trace result

Definition at line 697 of file TraceRay.C.

{
  if (!neighbor_info._valid)
    return NO_EXIT;
 
  const Elem * neighbor = neighbor_info._elem;
  debugRay("Checking neighbor ", neighbor->id(), " with centroid ", neighbor->vertex_average());
 
  // Find an entrant side (if any)
  incoming_side = RayTracingCommon::invalid_side;
  for (MooseIndex(neighbor_info._sides.size()) i = 0; i < neighbor_info._sides.size(); ++i)
    if (neighbor_info._side_normals[i] * (*_current_ray)->direction() < LOOSE_TRACE_TOLERANCE)
    {
      incoming_side = neighbor_info._sides[i];
      break;
    }
 
  // No entrant sides on this neighbor were found
  if (incoming_side == RayTracingCommon::invalid_side)
    return NO_EXIT;
 
  intersection_point = RayTracingCommon::invalid_point;
  intersected_side = RayTracingCommon::invalid_side;
  intersected_extrema.invalidate();
  intersection_distance = RayTracingCommon::invalid_distance;
 
  // See if there is a way through the neighbor element
  debugRay("Called exitsElem() from moveThroughNeighbor()");
  const auto exit_result =
      exitsElem(neighbor,
                neighbor->type(),
                incoming_side,
                intersection_point,
                intersected_side,
                intersected_extrema,
                intersection_distance,
                _backface_culling ? _study.getElemNormals(neighbor, _tid) : nullptr);
  debugRay("Done with exitsElem() from moveThroughNeighbor()");
 
  return exit_result;
}

Referenced by moveThroughNeighbors().

moveThroughNeighbors()

TraceRay::ExitsElemResult TraceRay::moveThroughNeighbors ( const std::vector< NeighborInfo > &  neighbors, const Elem *  last_elem, const Elem *&  best_elem, unsigned short &  best_elem_incoming_side  )

private

Moves the Ray though neighbors (vertex/edge/point)

Parameters

neighbors	The neighbors to try to move through
last_elem	The last element that was moved through - this one will be tried last
best_elem	The resulting element for which an intersection was found
best_elem_incoming_side	The incoming side on the resulting element for which

Returns

The trace result

Definition at line 584 of file TraceRay.C.

{
  ++_results[MOVED_THROUGH_NEIGHBORS];
 
  debugRay("Called moveThroughNeighbors() with ", neighbors.size(), " neighbors, and:");
  debugRay("  last_elem->id() = ", last_elem ? last_elem->id() : DofObject::invalid_id);
  debugRay("  _incoming_point = ", _incoming_point);
 
  traceAssert(!best_elem, "Best elem should be null");
  traceAssert(best_elem_incoming_side == RayTracingCommon::invalid_side,
              "Best elem side should be invalid");
  traceAssert(_intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  traceAssert(_intersected_side == RayTracingCommon::invalid_side, "Side should be invalid");
  traceAssert(_intersected_extrema.isInvalid(), "Extrema should be invalid");
  traceAssert(_intersection_distance == RayTracingCommon::invalid_distance,
              "Distance should be invalid");
 
  // Marker for the longest ray segment distance we've found in a neighbor. Start with a quite
  // small number because at this point we're desperate and will take anything!
  Real longest_distance = 1.0e-12;
  // Temporaries for the intersection checks
  unsigned short current_incoming_side;
  Point current_intersection_point;
  unsigned short current_intersected_side;
  ElemExtrema current_intersected_extrema;
  Real current_intersection_distance;
  TraceRay::ExitsElemResult best_exit_result = NO_EXIT;
 
  // The NeighborInfo for the last_elem
  // (to store so we can try it later if we fail for everyone else)
  const NeighborInfo * last_elem_info = nullptr;
 
  for (const NeighborInfo & neighbor_info : neighbors)
  {
    // If we're at the element we want to do last, skip it and store the info so
    // we can get to it later in case we fail at finding anything
    if (neighbor_info._elem == last_elem)
    {
      debugRay("Skipping last elem ", last_elem->id());
      last_elem_info = &neighbor_info;
      continue;
    }
 
    const auto exit_result = moveThroughNeighbor(neighbor_info,
                                                 current_incoming_side,
                                                 current_intersection_point,
                                                 current_intersected_side,
                                                 current_intersected_extrema,
                                                 current_intersection_distance);
 
    // Found one! See if it's the best way
    if (exit_result != NO_EXIT)
    {
      debugRay("Ray can exit through neighbor ", neighbor_info._elem->id());
 
      if (current_intersection_distance > longest_distance)
      {
        best_elem = neighbor_info._elem;
        best_elem_incoming_side = current_incoming_side;
        _intersection_point = current_intersection_point;
        _intersected_side = current_intersected_side;
        _intersected_extrema = current_intersected_extrema;
        _intersection_distance = current_intersection_distance;
        longest_distance = current_intersection_distance;
        best_exit_result = exit_result;
      }
    }
  }
 
  // Didn't find someone to exit, so try the last_elem (if any)
  if (!best_elem && last_elem_info)
  {
    const auto exit_result = moveThroughNeighbor(*last_elem_info,
                                                 current_incoming_side,
                                                 current_intersection_point,
                                                 current_intersected_side,
                                                 current_intersected_extrema,
                                                 current_intersection_distance);
 
    if (exit_result != NO_EXIT && current_intersection_distance > longest_distance)
    {
      debugRay("Ray can exit through last_elem ", last_elem->id());
 
      best_elem = last_elem;
      best_elem_incoming_side = current_incoming_side;
      _intersection_point = current_intersection_point;
      _intersected_side = current_intersected_side;
      _intersected_extrema = current_intersected_extrema;
      _intersection_distance = current_intersection_distance;
      best_exit_result = exit_result;
    }
  }
 
  debugRay("moveThroughNeighbors() best result:");
  debugRay("  best_elem = ", best_elem ? best_elem->id() : DofObject::invalid_id);
  debugRay("  best_elem_incoming_side = ", best_elem_incoming_side);
  debugRay("  _intersection_point = ", _intersection_point);
  debugRay("  _intersected_side = ", _intersected_side);
  debugRay("  _intersected_extrema = ", _intersected_extrema);
  debugRay("  _intersection_distance = ", _intersection_distance);
  if (best_elem)
  {
    debugRay("moveThroughNeighbors() next neighbor elem info:");
    debugRay(best_elem->get_info());
  }
 
  return best_exit_result;
}

Referenced by trace().

onBoundary()

void TraceRay::onBoundary ( const std::shared_ptr< Ray > &  ray, const bool  external  )

private

Called when a Ray hits a boundary.

Definition at line 2108 of file TraceRay.C.

{
  traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point),
              "Ray currentPoint() not set before onBoundary()");
 
  // Get the RayBCs on bnd_elems
  _study.getRayBCs(_on_boundary_ray_bcs, _boundary_elems, _tid, ray->id());
 
  // Store this information temprorarily because we are going to change it as we
  // apply each boundary condition
  const auto old_current_elem = _current_elem;
  const auto old_intersected_side = _intersected_side;
  const auto old_intersected_extrema = _intersected_extrema;
  const auto old_subdomain_id = _current_subdomain_id;
 
  // For each RayBC we found, apply it on the boundaries that we need
  for (RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
  {
    // First, find the boundaries this RayBC is valid on that are also in _boundary_elems.
    // We do this ahead of time so that we can pass in to RayBC::apply if the same
    // boundary condition is being appled multiple times on different boundaries.
    // This is useful in situations like reflection where multiple reflections are
    // necessary at a corner to perfectly reflect.
    _on_boundary_apply_index.clear();
    for (MooseIndex(_boundary_elems.size()) bnd_elems_i = 0; bnd_elems_i < _boundary_elems.size();
         ++bnd_elems_i)
      if (rbc->hasBoundary(_boundary_elems[bnd_elems_i].bnd_id))
        _on_boundary_apply_index.push_back(bnd_elems_i);
 
    traceAssert(!_on_boundary_apply_index.empty(), "Must not be empty");
 
    // Apply the RayBC on each of the relevant boundary elements
    for (const auto bnd_elems_i : _on_boundary_apply_index)
    {
      auto & bnd_elem = _boundary_elems[bnd_elems_i];
 
      debugRay("Calling ",
               rbc->type(),
               "::onBoundary for \"",
               rbc->name(),
               "\" on elem ",
               bnd_elem.elem->id(),
               " and side ",
               bnd_elem.side,
               " for bnd_id ",
               bnd_elem.bnd_id);
 
      _current_elem = bnd_elem.elem;
      _current_bnd_id = bnd_elem.bnd_id;
      _intersected_side = bnd_elem.side;
      _intersected_extrema = bnd_elem.extrema;
      _current_subdomain_id = _current_elem->subdomain_id();
 
      rbc->onBoundary(_on_boundary_apply_index.size());
      postRayTracingObject(ray, rbc);
    }
  }
 
  // Set this info back now that we're done applying BCs
  _current_elem = old_current_elem;
  _intersected_side = old_intersected_side;
  _intersected_extrema = old_intersected_extrema;
  _current_bnd_id = BoundaryInfo::invalid_id;
  _current_subdomain_id = old_subdomain_id;
 
  // When on an external boundary, the Ray must have been changed or killed.
  // Otherwise, we don't know what to do with it now! If this didn't happen,
  // output a detailed error message.
  if (external && !ray->trajectoryChanged() && ray->shouldContinue())
  {
    std::stringstream oss;
    oss << "Don't know what to do with a Ray after it hit an external\n";
    oss << "boundary at point " << _intersection_point << "!\n\n";
    oss << "When hitting an external RayBC, a Ray must either:\n";
    oss << "  Be killed by a RayBC\n";
    oss << "  Have its trajectory changed by the RayBC\n";
    oss << "by at least one of the executed RayBCs.\n\n";
    oss << "You need to either:\n";
    oss << "  Kill/change the Ray sooner with RayKernels, internal RayBCs, or a max distance\n";
    oss << "  Kill/change the Ray on the boundary with a RayBC\n\n";
    if (!_on_boundary_ray_bcs.empty())
    {
      oss << "RayBCs executed that did not kill or change the Ray:\n";
      for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
        for (const auto & bnd_elem : _boundary_elems)
          if (rbc->hasBoundary(bnd_elem.bnd_id))
            oss << "  " << rbc->typeAndName() << " on boundary " << bnd_elem.bnd_id << " ("
                << _mesh.getBoundaryName(bnd_elem.bnd_id) << ")\n";
      oss << "\n";
    }
    bool output_header = false;
    for (std::size_t i = 0; i < _boundary_elems.size(); ++i)
    {
      const auto bnd_id = _boundary_elems[i].bnd_id;
      bool found = false;
      for (const RayBoundaryConditionBase * rbc : _on_boundary_ray_bcs)
        if (rbc->hasBoundary(bnd_id))
        {
          found = true;
          break;
        }
 
      if (!found)
      {
        if (!output_header)
        {
          oss << "Boundaries that did not have any RayBCs:\n";
          output_header = true;
        }
        oss << "  " << bnd_id << " (" << _mesh.getBoundaryName(bnd_id) << ")\n";
      }
    }
 
    failTrace(oss.str(), _study.tolerateFailure(), __LINE__);
  }
}

Referenced by applyOnExternalBoundary(), and applyOnInternalBoundary().

onCompleteTrace()

void TraceRay::onCompleteTrace ( const std::shared_ptr< Ray > &  ray )

private

Called when a Ray is finished tracing (whenever !ray->shouldContinue())

Parameters

ray	The ray that is finished tracing

Definition at line 1653 of file TraceRay.C.

{
  for (RayKernelBase * rk : _study.currentRayKernels(_tid))
    rk->postTrace();
 
  debugRay("Called onCompleteTrace()\n", (*_current_ray)->getInfo());
  if (_intersection_distance > 0)
    possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
  possiblySaveDebugRayMesh();
 
  if (_current_cached_trace)
  {
    _current_cached_trace->_last = true;
 
    if (_intersection_distance > 0)
    {
      _current_cached_trace->addPoint(ray->currentPoint());
      if (_study.dataOnCacheTraces())
        _current_cached_trace->lastPoint()._data = ray->data();
      if (_study.auxDataOnCacheTraces())
        _current_cached_trace->lastPoint()._aux_data = ray->auxData();
    }
 
    mooseAssert(ray->stationary() == _current_cached_trace->stationary(), "Stationary mismatch");
  }
}

Referenced by trace().

onContinueTrace()

void TraceRay::onContinueTrace ( const std::shared_ptr< Ray > &  ray )

private

Called when a Ray is continuing to trace after segment.

Parameters

ray

The ray

Definition at line 1681 of file TraceRay.C.

{
  traceAssert(ray->shouldContinue(), "Ray must continue");
 
  if (_current_cached_trace && _study.segmentsOnCacheTraces() && _intersection_distance > 0)
  {
    _current_cached_trace->addPoint(ray->currentPoint());
    if (_study.dataOnCacheTraces())
      _current_cached_trace->lastPoint()._data = ray->data();
    if (_study.auxDataOnCacheTraces())
      _current_cached_trace->lastPoint()._aux_data = ray->auxData();
  }
}

Referenced by trace().

onSegment()

void TraceRay::onSegment ( const std::shared_ptr< Ray > &  ray )

private

Called on a single segment traced by a Ray.

Definition at line 2051 of file TraceRay.C.

{
  traceAssert((*_current_ray)->currentElem() == _current_elem, "Ray currentElem() incorrect");
  traceAssert((*_current_ray)->currentPoint() == _intersection_point,
              "Ray currentPoint() incorrect");
  traceAssert((*_current_ray)->currentIncomingSide() == _incoming_side,
              "Ray currentIncomingSide() incorrect");
#ifndef NDEBUG
  if (_study.verifyTraceIntersections())
  {
    if (_current_elem->has_affine_map())
      traceAssert(_current_elem->contains_point(_incoming_point),
                  "_current_elem does not contain incoming point");
 
    if (_intersected_side != RayTracingCommon::invalid_side &&
        !_study.sideIsNonPlanar(_current_elem, _intersected_side))
    {
      traceAssert(_elem_side_builder(*_current_elem, _intersected_side)
                      .close_to_point(_intersection_point, LOOSE_TRACE_TOLERANCE),
                  "Intersected point is not on intersected side");
      traceAssert(!_study.sideIsIncoming(
                      _current_elem, _intersected_side, (*_current_ray)->direction(), _tid),
                  "Intersected side is not outgoing");
    }
    if (_incoming_side != RayTracingCommon::invalid_side &&
        !_study.sideIsNonPlanar(_current_elem, _incoming_side))
    {
      traceAssert(_elem_side_builder(*_current_elem, _incoming_side)
                      .close_to_point(_incoming_point, LOOSE_TRACE_TOLERANCE),
                  "Incoming point is not on incoming side");
      if (ray->intersections() != 0 && ray->maxDistance() != 0)
        traceAssert(_study.sideIsIncoming(
                        _current_elem, _incoming_side, (*_current_ray)->direction(), _tid),
                    "Incoming side is not incoming");
    }
  }
#endif
  traceAssert(MooseUtils::absoluteFuzzyEqual(_intersection_distance,
                                             (_intersection_point - _incoming_point).norm()),
              "_intersection_distance is incorrect");
  traceAssert(_current_subdomain_id == _current_elem->subdomain_id(), "Subdomain incorrect");
  traceAssert(MooseUtils::absoluteFuzzyEqual((_incoming_point - _intersection_point).norm(),
                                             _intersection_distance),
              "Invalid intersection distance");
 
  _study.reinitSegment(
      _current_elem, _incoming_point, _intersection_point, _intersection_distance, _tid);
 
  const auto & rks = _study.currentRayKernels(_tid);
  for (auto & rk : rks)
  {
    rk->onSegment();
    postRayTracingObject(ray, rk);
  }
}

Referenced by trace().

onSubdomainChanged()

void TraceRay::onSubdomainChanged ( const std::shared_ptr< Ray > &  ray, const bool  same_ray  )

private

Called when the subdomain changes.

Parameters

ray	The current Ray
same_ray	Whether or not this is being called on the same Ray as previously

Definition at line 1750 of file TraceRay.C.

{
  debugRay("Calling onSubdomainChanged() on subdomain ", _current_elem->subdomain_id());
  debugRay("  _current_subdomain_id = ", _current_subdomain_id);
 
  _current_subdomain_id = _current_elem->subdomain_id();
  _current_subdomain_hmax = _study.subdomainHmax(_current_subdomain_id);
  _current_elem_type = _current_elem->type();
  _current_elem_n_sides = _current_elem->n_sides();
 
  if (_has_ray_kernels)
  {
    auto & current_ray_kernels = _study.currentRayKernels(_tid);
 
    // If we're still tracing the same Ray, keep track of our old RayKernels
    // so that we don't call preTrace() on them again
    if (same_ray)
      _old_ray_kernels.insert(current_ray_kernels.begin(), current_ray_kernels.end());
    // If we're not tracing the same Ray, we need to call preTrace() on everything
    else
      _old_ray_kernels.clear();
 
    // Call segmentSubdomainSetup to get new kernels etc
    _study.segmentSubdomainSetup(_current_subdomain_id, _tid, ray->id());
 
    // Call preTrace() on all of the RayKernels that need it
    for (RayKernelBase * rk : current_ray_kernels)
      // Haven't called preTrace() for this Ray on this RayKernel yet
      if (!_old_ray_kernels.count(rk))
        rk->preTrace();
  }
}

Referenced by trace().

onTrajectoryChanged()

void TraceRay::onTrajectoryChanged ( const std::shared_ptr< Ray > &  ray )

private

Called when a Ray's trajectory changes.

Parameters

ray

The ray

Definition at line 1721 of file TraceRay.C.

{
#ifndef NDEBUG
  if (_study.verifyTraceIntersections() &&
      (_intersected_extrema.atExtrema()
           ? !_current_elem->close_to_point(ray->currentPoint(), LOOSE_TRACE_TOLERANCE)
           : !_current_elem->contains_point(ray->currentPoint())))
    failTrace("Elem does not contain point after trajectory change",
              /* warning = */ false,
              __LINE__);
#endif
 
  traceAssert(ray->shouldContinue(), "Ray should continue when trajectory is being changed");
 
  ray->setTrajectoryChanged(false);
  ray->addTrajectoryChange();
 
  if (_current_cached_trace && !_study.segmentsOnCacheTraces())
  {
    if (_intersection_distance > 0)
      _current_cached_trace->addPoint(ray->currentPoint());
    if (_study.dataOnCacheTraces())
      _current_cached_trace->lastPoint()._data = ray->data();
    if (_study.auxDataOnCacheTraces())
      _current_cached_trace->lastPoint()._aux_data = ray->auxData();
  }
}

Referenced by trace().

possiblyAddToBoundaryElems()

void TraceRay::possiblyAddToBoundaryElems ( const Elem *  elem, const unsigned short  side, const std::vector< BoundaryID > &  bnd_ids, const ElemExtrema &  extrema  )

private

Helper for possibly storing boundary information in _boundary_elems, which is storage for boundary elements (elem, side, bnd_id) that need to have RayBCs applied to them.

For each BoundaryID in bnd_ids, a ConstBndElement will be added if said BoundaryID does not already exist in _boundary_elems.

Parameters

elem	The element to possibly add
side	The side to possibly add
bnd_ids	The BoundaryIDs to possibly add
extrema	The intersected elem extrema (vertex/edge), if any

Definition at line 916 of file TraceRay.C.

{
  if (!_study.sideIsNonPlanar(elem, side))
    traceAssert(extrema.isValid(elem, _intersection_point), "Extrema not correct");
 
  for (const auto bnd_id : bnd_ids)
  {
    bool found = false;
    for (const auto & bnd_elem : _boundary_elems)
      if (bnd_elem.bnd_id == bnd_id)
      {
        found = true;
        break;
      }
 
    if (!found)
    {
      debugRay("  Need to apply boundary on elem ",
               elem->id(),
               " and side ",
               side,
               " for bnd_id ",
               bnd_id);
 
      _boundary_elems.emplace_back(elem, side, bnd_id, extrema);
    }
  }
}

Referenced by applyOnExternalBoundary(), and applyOnInternalBoundary().

postRayTracingObject()

void TraceRay::postRayTracingObject ( const std::shared_ptr< Ray > &  ray, const RayTracingObject *  rto  )

private

Called after executing a RayTracingObject (RayBCs and RayKernels)

Verifies the configuration of ray->shouldContinue() and ray->trajectoryChanged()

Definition at line 2234 of file TraceRay.C.

{
  if (!ray->shouldContinue())
  {
    if (_should_continue)
      _should_continue = false;
  }
  else if (!_should_continue)
    failTrace(rto->typeAndName() +
                  " set a Ray to continue that was previously set to not continue.\n\n" +
                  "Once a Ray has been set to not continue, its continue status cannot change.",
              /* warning = */ false,
              __LINE__);
 
  if (!_should_continue && ray->trajectoryChanged())
    failTrace(rto->typeAndName() +
                  " changed the trajectory of a Ray that was set to not continue,\n" +
                  "or set a Ray whose trajectory was changed to not continue.",
              /* warning = */ false,
              __LINE__);
}

Referenced by onBoundary(), and onSegment().

preExecute()

void TraceRay::preExecute

(

)

Should be called immediately before calling any traces.

Definition at line 60 of file TraceRay.C.

{
  _current_subdomain_id = Elem::invalid_subdomain_id;
  _current_elem_type = libMesh::INVALID_ELEM;
 
  // Zero out all results
  for (auto & val : _results)
    val = 0;
 
  _has_ray_kernels = _study.hasRayKernels(_tid);
  _is_rectangular_domain = _study.isRectangularDomain();
}

results()

std::vector< unsigned long long int > TraceRay::results ( ) const

inline

Get the various results reported by this object, indexed by TraceRayResult.

Definition at line 93 of file TraceRay.h.

{ return _results; }

setBackfaceCulling()

void TraceRay::setBackfaceCulling ( const bool  backface_culling )

inline

Enable or disable the use of element normals for backface culling through the getElemNormals() method.

By default, the _study.getElemNormals() method is not implemented. This is because in order to be able to obtain the element normals often enough during a trace, one needs to implement a method to cache the normals. Otherwise, it is not worthwhile to use the normals during the trace because the time spent computing the normals on the fly will be greater than the time gain by using the normals in the tracing process.

Definition at line 105 of file TraceRay.h.

{ _backface_culling = backface_culling; }

Referenced by BackfaceCullingStudyTest::BackfaceCullingStudyTest().

storeExitsElemResult()

void TraceRay::storeExitsElemResult ( const ExitsElemResult  result )

private

Stores the result given by an intersection in _results as necessary.

Definition at line 2038 of file TraceRay.C.

{
  if (result == HIT_FACE)
    ++_results[FACE_HITS];
  else if (result == HIT_VERTEX)
    ++_results[VERTEX_HITS];
  else if (result == HIT_EDGE)
    ++_results[EDGE_HITS];
  else
    mooseError("Should not call storeExitsElemResult() with result ", result);
}

Referenced by trace().

subdomainHmax()

Real TraceRay::subdomainHmax ( const Elem *  elem ) const

private

Get the approximate subdomain hmax for an element.

Uses the cached value of _current_subdomain_hmax if the element is in the cached subdomain.

Definition at line 2226 of file TraceRay.C.

{
  const auto subdomain_id = elem->subdomain_id();
  return subdomain_id == _current_subdomain_id ? _current_subdomain_hmax
                                               : _study.subdomainHmax(subdomain_id);
}

Referenced by exitsElem().

trace()

void TraceRay::trace

(

const std::shared_ptr< Ray > &

ray

)

Traces a ray.

Definition at line 1018 of file TraceRay.C.

{
  mooseAssert(_study.currentlyPropagating(), "Should only use while propagating rays");
 
  _current_ray = &ray;
  _current_elem = ray->currentElem();
  _last_elem = nullptr;
  _incoming_point = ray->currentPoint();
  _incoming_side = ray->currentIncomingSide();
  _should_continue = true;
 
  _study.preTrace(_tid, ray);
 
  traceAssert(_current_elem, "Current element is not set");
  traceAssert(_current_elem->active(), "Current element is not active");
  traceAssert(!ray->invalidCurrentPoint(), "Current point is invalid");
  traceAssert(ray->shouldContinue(), "Ray should not continue");
  if (_study.verifyRays() && !ray->invalidCurrentIncomingSide() && ray->maxDistance() > 0 &&
      !_study.sideIsNonPlanar(_current_elem, _incoming_side) &&
      !_study.sideIsIncoming(_current_elem, _incoming_side, ray->direction(), _tid))
    failTrace("Ray incoming side is not incoming", /* warning = */ false, __LINE__);
 
#ifdef DEBUG_RAY_IF
  if (DEBUG_RAY_IF)
    libMesh::err << "\n\n";
#endif
  debugRay("At top of trace for Ray");
  debugRay("Top of trace loop Ray info\n", ray->getInfo());
  debugRay("Top of trace loop starting elem info\n", ray->currentElem()->get_info());
 
  // Invalidate this up front because it's copied immediately into _last_intersected_extrema
  _intersected_extrema.invalidate();
 
#ifdef DEBUG_RAY_MESH_IF
  _debug_mesh = nullptr;
  _debug_node_count = 0;
  if (DEBUG_RAY_MESH_IF)
  {
    _debug_mesh = new Mesh(_debug_comm, _dim);
    _debug_mesh->skip_partitioning(true);
  }
#endif
 
  // Caching trace along the way: init for this Ray
  if (_study.shouldCacheTrace(ray))
  {
    debugRay("Trying to init threaded cached trace");
 
    _current_cached_trace = &_study.initThreadedCachedTrace(ray, _tid);
 
    // Add starting data
    if (_study.dataOnCacheTraces())
      _current_cached_trace->lastPoint()._data = ray->data();
    if (_study.auxDataOnCacheTraces())
      _current_cached_trace->lastPoint()._aux_data = ray->auxData();
  }
  else
    _current_cached_trace = nullptr;
 
  // Need to call subdomain setup
  if (_current_elem->subdomain_id() != _current_subdomain_id || _study.rayDependentSubdomainSetup())
    onSubdomainChanged(ray, /* same_ray = */ false);
  // If we didn't change subdomains, we still need to call this on the RayKernels
  else
    for (RayKernelBase * rk : _study.currentRayKernels(_tid))
      rk->preTrace();
 
  // Ray tracing loop through each segment for a single Ray on this processor
  do
  {
#ifdef DEBUG_RAY_IF
    if (DEBUG_RAY_IF)
      libMesh::err << "\n\n";
#endif
    debugRay("At top of ray tracing loop");
    debugRay("  ray->id() = ", ray->id());
    debugRay("  _incoming_point = ", _incoming_point);
    debugRay("  _incoming_side = ", _incoming_side);
    debugRay("  _current_elem->id() = ", _current_elem->id());
    debugRay("  _current_elem->subdomain_id() = ", _current_elem->subdomain_id());
    debugRay("  _current_subdomain_id = ", _current_subdomain_id);
    debugRay("  _current_elem_type = ", Utility::enum_to_string(_current_elem_type));
    debugRay("Top of ray tracing loop Ray info\n", ray->getInfo());
    debugRay("Top of ray tracing loop current elem info\n", _current_elem->get_info());
 
    traceAssert(_current_ray == &ray, "Current ray mismatch");
 
    // Copy over in case we need to use it
    _last_intersected_extrema = _intersected_extrema;
    // Invalidate all intersections as we're tracing again
    _exits_elem = false;
    _intersection_point = RayTracingCommon::invalid_point;
    _intersected_side = RayTracingCommon::invalid_side;
    _intersected_extrema.invalidate();
    _intersection_distance = RayTracingCommon::invalid_distance;
 
    // Stationary ray
    if (ray->stationary())
    {
      mooseAssert(ray->invalidDirection(), "Should have an invalid direction");
      _exits_elem = true;
      _intersection_point = _incoming_point;
      _intersected_extrema = _last_intersected_extrema;
      _intersection_distance = 0;
      ++_results[ENDED_STATIONARY];
    }
    // If we haven't hit a vertex or an edge, do the normal exit algorithm first.
    // In the case of a Ray that previously moved through point neighbors due to
    // being at a vertex/edge, this will be true because we do not communicate the
    // vertex/edge intersection. The previous processor already set us up for
    // the best intersection because it already computed it and chose to
    // send it our way as such
    else if (!_last_intersected_extrema.atExtrema())
    {
      traceAssert(_current_elem->processor_id() == _pid, "Trace elem not on processor");
      debugRay("Didn't hit vertex or edge: doing normal exits elem check");
 
      if (_backface_culling)
        _current_normals = _study.getElemNormals(_current_elem, _tid);
 
      const auto exits_elem_result = exitsElem(_current_elem,
                                               _current_elem_type,
                                               _incoming_side,
                                               _intersection_point,
                                               _intersected_side,
                                               _intersected_extrema,
                                               _intersection_distance,
                                               _current_normals);
 
      if (exits_elem_result != NO_EXIT)
      {
        storeExitsElemResult(exits_elem_result);
        _exits_elem = true;
        ray->setCurrentPoint(_intersection_point);
      }
    }
    else
    {
      debugRay("Will not do normal exits elem check because at a vertex/edge");
    }
 
    // At this point, we either did a regular exit check and it failed, or we hit a vertex/edge
    // on the previous intersection and didn't bother to do a regular exit check on
    // _current_elem
    if (!_exits_elem)
    {
      debugRay("Moving through neighbors");
 
      // The element we want moveThroughNeighbors() to try last. That is - if all others fail,
      // this is the last resort. If we have a last intersection that tells us we're going
      // through a vertex or edge, we probably won't go back to that elem. If we don't, it means
      // that we failed the exitsElem() check on _current_elem and probably won't return to it
      // either.
      const Elem * move_through_neighbors_last = nullptr;
 
      // Get the neighbors
      const std::vector<NeighborInfo> * neighbors = nullptr;
      // If we have previous vertex/edge information, use it
      if (_last_intersected_extrema.atExtrema())
      {
        traceAssert(_last_elem, "Should be valid");
        move_through_neighbors_last = _last_elem;
        neighbors = &getNeighbors(_last_elem, _last_intersected_extrema, _incoming_point);
      }
      // Without previous vertex/edge information, let's try to find some. We could just do
      // a pure point neighbor check at the current point, but we'd rather be able to cache
      // this information so that someone else can use it - this requires the more unique
      // identifier of which vertex/edge we are at.
      else
      {
        debugRay("  Searching for vertex/edge hit with incoming side ",
                 _incoming_side,
                 " on elem ",
                 _current_elem->id(),
                 " at ",
                 _incoming_point);
 
        move_through_neighbors_last = _current_elem;
 
        // If we have side info: check for vertices on said side, otherwise, check everywhere
        const auto at_v = _incoming_side != RayTracingCommon::invalid_side
                              ? atVertexOnSide(_current_elem, _incoming_point, _incoming_side)
                              : atVertex(_current_elem, _incoming_point);
 
        if (at_v != RayTracingCommon::invalid_vertex)
          neighbors = &getVertexNeighbors(_current_elem, at_v);
        // If still nothing and in 3D, check if we're on an edge
        // TODO: Handle 2D with a similar check for if we're on a side
        else if (_dim == 3)
        {
          ElemExtrema extrema;
          if (_incoming_side != RayTracingCommon::invalid_side)
            withinEdgeOnSide(_current_elem, _incoming_point, _incoming_side, extrema);
          else
            withinEdge(_current_elem, _incoming_point, extrema);
 
          if (extrema.atEdge())
            neighbors = &getEdgeNeighbors(_current_elem, extrema, _incoming_point);
        }
 
        // If we still haven't found anything - let's try a last-ditch effort
        if (!neighbors || neighbors->empty())
          neighbors = &getPointNeighbors(_current_elem, _incoming_point);
 
        // Couldn't find anything
        if (neighbors->empty())
        {
          failTrace("Could not find neighbors to move through", _study.tolerateFailure(), __LINE__);
          return;
        }
      }
 
      // Move through a neighbor
      const Elem * best_neighbor = nullptr;
      auto best_neighbor_side = RayTracingCommon::invalid_side;
      const auto exits_elem_result = moveThroughNeighbors(
          *neighbors, move_through_neighbors_last, best_neighbor, best_neighbor_side);
 
      // If we didn't find anything... we're out of luck
      if (exits_elem_result == NO_EXIT)
      {
        failTrace("Could not find intersection after trying to move through point neighbors",
                  _study.tolerateFailure(),
                  __LINE__);
        return;
      }
 
      // At this point, we've successfully made it through an element and also started
      // through another element, so set that
      _exits_elem = true;
      _last_elem = _current_elem;
      _current_elem = best_neighbor;
      _incoming_side = best_neighbor_side;
      ray->setCurrentElem(best_neighbor);
      ray->setCurrentIncomingSide(best_neighbor_side);
      ray->setCurrentPoint(_intersection_point);
 
      // Don't own this element - return exits trace for this Ray on this proc
      if (best_neighbor->processor_id() != _pid)
      {
        // We've already computed the next intersection but said intersection is
        // on an elem on another processor. Therefore, the next proc will re-trace
        // this Ray on the perfect elem that we've picked (best_neighbor)
        ray->setCurrentPoint(_incoming_point);
        _intersection_distance = RayTracingCommon::invalid_distance;
 
        continueTraceOffProcessor(ray);
        return;
      }
 
      // Subdomain changed
      if (_current_elem->subdomain_id() != _current_subdomain_id)
        onSubdomainChanged(ray, /* same_ray = */ true);
 
      // Do own this element - tally the result as we're the ones tracing it
      storeExitsElemResult(exits_elem_result);
    }
 
    debugRay("Done with trace");
    debugRay("  _exits_elem: ", _exits_elem);
    debugRay("  _intersection_point: ", _intersection_point);
    debugRay("  _intersected_side: ", _intersected_side);
    debugRay("  _intersected_side centroid: ",
             _intersected_side == RayTracingCommon::invalid_side
                 ? RayTracingCommon::invalid_point
                 : _current_elem->side_ptr(_intersected_side)->vertex_average());
    debugRay("  _intersected_extrema = ", _intersected_extrema);
    debugRay("  _intersection_distance = ", _intersection_distance);
 
    // Increment intersections
    if (_intersection_distance > 0)
    {
      debugRay("Incrementing ray intersections by 1 to ", ray->intersections() + 1);
      ray->addIntersection();
      _results[INTERSECTIONS]++;
    }
 
    // Increment distance
    ray->addDistance(_intersection_distance);
    debugRay("Incremented ray distance by ", _intersection_distance);
 
    // The effective max distance that the Ray should travel - minimum of the two
    const auto max_distance = std::min(ray->maxDistance(), _study.rayMaxDistance());
    debugRay("Max distance checks");
    debugRay("  _study.rayMaxDistance() = ", _study.rayMaxDistance());
    debugRay("  ray->maxDistance() = ", ray->maxDistance());
    debugRay("  max_distance (effective) = ", max_distance);
 
    // At the maximum distance - distinguish between this and past the maximum
    // distance because in this case we're close enough to the intersection point
    // that we can keep it, its intersected side, and intersected extrema
    if (MooseUtils::absoluteFuzzyEqual(ray->distance(), max_distance))
    {
      debugRay("At max distance");
 
      ray->setShouldContinue(false);
      _should_continue = false;
    }
    // Past the max distance - need to remove the additional distance we traveled,
    // change the point and invalidate the intersection data (moves the Ray back)
    else if (ray->distance() > max_distance)
    {
      debugRay("Past max distance");
 
      // The distance past the max distance we have traveled
      const auto difference = ray->distance() - max_distance;
      traceAssert(difference > 0, "Negative distance change after past_max_distance");
 
      debugRay("Removing distance ", difference);
      ray->addDistance(-difference);
      debugRay("  New ray->distance() = ", ray->distance());
 
      _intersection_point -= ray->direction() * difference;
      _intersection_distance -= difference;
      _intersected_side = RayTracingCommon::invalid_side;
      _intersected_extrema.invalidate();
      ray->setCurrentPoint(_intersection_point);
 
      ray->setShouldContinue(false);
      _should_continue = false;
 
      traceAssert(_intersection_distance >= 0, "Negative _intersection_distance");
#ifndef NDEBUG
      if (_study.verifyTraceIntersections() && !_current_elem->contains_point(_intersection_point))
        failTrace("Does not contain point after past max distance",
                  /* warning = */ false,
                  __LINE__);
#endif
    }
 
    if (!_study.currentRayKernels(_tid).empty())
    {
      debugRay("Calling onSegment() with");
      debugRay("  current_elem->id() = ", _current_elem->id());
      debugRay("  _incoming_point = ", _incoming_point);
      debugRay("  _incoming_side = ", _incoming_side);
      debugRay("  _intersection_point = ", _intersection_point);
      debugRay("  _intersected_side = ", _intersected_side);
      debugRay("  _intersected_extrema = ", _intersected_extrema);
      debugRay("  _intersection_distance = ", _intersection_distance);
      onSegment(ray);
 
      _study.postOnSegment(_tid, ray);
 
      // RayKernel killed a Ray or we're at the end
      traceAssert(_should_continue == ray->shouldContinue(), "Should be the same");
      if (!_should_continue)
      {
        debugRay("RayKernel killed the ray or past max distance");
        traceAssert(!ray->trajectoryChanged(),
                    "RayKernels should not change trajectories of Rays at end");
 
        onCompleteTrace(ray);
        return;
      }
 
      // RayKernel moved a Ray
      if (ray->trajectoryChanged())
      {
        debugRay("RayKernel changed the Ray's trajectory");
        debugRay("  new direction = ", ray->direction());
        debugRay("  old incoming point = ", _incoming_point);
        debugRay("  new incoming point = ", ray->currentPoint());
        possiblyAddDebugRayMeshPoint(_incoming_point, ray->currentPoint());
 
        _incoming_point = ray->currentPoint();
        _incoming_side = RayTracingCommon::invalid_side;
        _intersected_extrema.invalidate();
        ray->setCurrentIncomingSide(_incoming_side);
 
        const auto new_intersection_distance = (ray->currentPoint() - _incoming_point).norm();
        ray->addDistance(-_intersection_distance + new_intersection_distance);
        _intersection_distance = new_intersection_distance;
 
        onTrajectoryChanged(ray);
        onContinueTrace(ray);
        continue;
      }
      else
        possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
    }
    else
    {
      _study.postOnSegment(_tid, ray);
 
      if (!_should_continue)
      {
        debugRay("Killing due to at end without RayKernels");
        traceAssert(!ray->shouldContinue(), "Ray shouldn't continue");
 
        onCompleteTrace(ray);
        return;
      }
    }
 
    // If at a vertex/on an edge and not on the boundary, we may actually be on the boundary,
    // just not on a boundary side. Check for that here. If the domain is rectangular we will
    // do a quick check against the bounding box as if we're not on the boundary of the
    // bounding box for a rectangular problem, we can skip this.
    // TODO: see how much the tolerance for the bounding box check can be tightened
    if (_dim > 1 && _intersected_extrema.atExtrema() &&
        _current_elem->neighbor_ptr(_intersected_side) &&
        (!_is_rectangular_domain ||
         onBoundingBoxBoundary(_study.boundingBox(),
                               _intersection_point,
                               _dim,
                               LOOSE_TRACE_TOLERANCE * _study.domainMaxLength())))
    {
      auto boundary_side = RayTracingCommon::invalid_side;
      ElemExtrema boundary_extrema;
      const Elem * boundary_elem = nullptr;
 
      findExternalBoundarySide(boundary_side, boundary_extrema, boundary_elem);
 
      if (boundary_elem)
      {
        // At this point, the new incoming side is very difficult to find
        // and may require some re-tracing backwards. Let's not do that -
        // we don't need it to continue the trace. This is reason why
        // _incoming_side is not available in RayBCs.
        _last_elem = _current_elem;
        _current_elem = boundary_elem;
        _intersected_side = boundary_side;
        _intersected_extrema = boundary_extrema;
        ray->setCurrentElem(_current_elem);
 
        debugRay("Found a neighbor boundary side with:");
        debugRay("  _current_elem->id() = ", _current_elem->id());
        debugRay("  _intersected_side = ", _intersected_side);
        debugRay("  _intersected_extrema = ", _intersected_extrema);
      }
    }
 
    // The next element
    const Elem * neighbor = nullptr;
 
    // Set up where to go next
    _incoming_point = _intersection_point;
    debugRay("Set _incoming point to ", _incoming_point);
 
    debugRay("Looking for neighbor on side ", _intersected_side, " of elem ", _current_elem->id());
 
    // Get the neighbor on that side
    // If the mesh has active elements of the same level, just grab the neighbor directly. In
    // this case, we can avoid a call to active() on the neighbor in getActiveNeighbor(), which
    // can be expensive depending on caching.
    neighbor = _study.hasSameLevelActiveElems()
                   ? _current_elem->neighbor_ptr(_intersected_side)
                   : getActiveNeighbor(_current_elem, _intersected_side, _intersection_point);
 
    // Found one - not on the boundary so set the next info
    if (neighbor)
    {
      traceAssert(neighbor->active(), "Inactive neighbor");
      traceAssert(_current_elem->subdomain_id() == _current_subdomain_id,
                  "_current_subdomain_id invalid");
 
      // If the mesh has active elements of the same level, don't use
      // neighbor->which_neighbor_am_i(), which calls active() and an n_sides() virtual call
      if (_study.hasSameLevelActiveElems())
      {
        // If the subdomain hasn't changed, we can guarantee n_sides is the same as
        // _current_elem. Prefer this because neighbor->n_sides() is an expensive virtual. Even
        // better, subdomain_id() isn't a virtual!
        const unsigned short n_sides = neighbor->subdomain_id() == _current_subdomain_id
                                           ? _current_elem_n_sides
                                           : neighbor->n_sides();
        traceAssert(n_sides == neighbor->n_sides(), "n_sides incorrect");
 
        for (_incoming_side = 0; _incoming_side < n_sides; ++_incoming_side)
          if (neighbor->neighbor_ptr(_incoming_side) == _current_elem)
            break;
      }
      else
        _incoming_side = neighbor->which_neighbor_am_i(_current_elem);
 
      _last_elem = _current_elem;
      _current_elem = neighbor;
      ray->setCurrentElem(neighbor);
      ray->setCurrentIncomingSide(_incoming_side);
 
      debugRay("Next elem: ", neighbor->id(), " with centroid ", neighbor->vertex_average());
      debugRay("Next _incoming_side: ",
               _incoming_side,
               " with centroid ",
               neighbor->side_ptr(_incoming_side)->vertex_average());
      traceAssert(_last_elem->subdomain_id() == _current_subdomain_id,
                  "_current_subdomain_id invalid");
 
      // Whether or not the subdomain changed
      const bool subdomain_changed = neighbor->subdomain_id() != _current_subdomain_id;
 
      // See if we hit any internal sides leaving this element and apply
      // We require that all internal RayBCs be on internal sidesets that have different
      // subdomain ids on each side. Therefore, only check if at vertex/edge or if
      // the subdomain changes
      if (_study.hasInternalSidesets() && (subdomain_changed || _intersected_extrema.atExtrema()))
      {
        applyOnInternalBoundary(ray);
 
        // Internal RayBC killed a Ray
        if (!_should_continue)
        {
          traceAssert(!ray->shouldContinue(), "Should be the same");
          debugRay("Internal RayBC killed the ray");
 
          onCompleteTrace(ray);
          return;
        }
 
        // Internal RayBC changed the Ray
        if (ray->trajectoryChanged())
        {
          debugRay("Internal RayBC changed the trajectory:");
          debugRay("  new direction = ", ray->direction());
 
          // Is this side still incoming?
          const auto normal = _study.getSideNormal(_current_elem, _incoming_side, _tid);
          const auto dot = normal * ray->direction();
          debugRay("Dot product with new direction and side = ", dot);
          if (dot > -TRACE_TOLERANCE)
          {
            _incoming_side = _last_elem->which_neighbor_am_i(_current_elem);
            _current_elem = _last_elem;
            neighbor = _last_elem;
            ray->setCurrentElem(_current_elem);
            ray->setCurrentIncomingSide(_incoming_side);
            debugRay("  Dot > 0 (Ray turned around): Setting _current_elem = ",
                     _current_elem->id(),
                     " and _incoming_side = ",
                     _incoming_side);
          }
 
          onTrajectoryChanged(ray);
        }
 
        traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point, TRACE_TOLERANCE),
                    "Internal RayBC changed the Ray point");
      }
 
      // Neighbor is off processor
      // If we hit at a vertex/edge, we will continue and let the move through neighbor exit
      // figure out who to send to
      if (neighbor->processor_id() != _pid)
      {
        if (_intersected_extrema.atExtrema())
        {
          debugRay("Neighbor is off processor but continuing to move through neighbors");
        }
        else
        {
          continueTraceOffProcessor(ray);
          return;
        }
      }
 
      // Neighbor is on processor, call subdomain setup if needed
      if (subdomain_changed)
        onSubdomainChanged(ray, /* same_ray = */ true);
    }
    // No neighbor found: on the boundary
    else
    {
      debugRay("No neighbor found - on the boundary");
 
      // Apply boundary conditions
      applyOnExternalBoundary(ray);
 
      if (_current_elem == ray->currentElem())
        traceAssert(ray->currentPoint().absolute_fuzzy_equals(_intersection_point, TRACE_TOLERANCE),
                    "RayBC changed the Ray point");
 
      // Quit tracing if the Ray was killed by a BC
      if (!_should_continue)
      {
        traceAssert(!ray->shouldContinue(), "Should be the same");
        debugRay("Exiting due to death by BC");
 
        onCompleteTrace(ray);
        return;
      }
      // RayBC changed the direction of the Ray
      if (ray->trajectoryChanged())
      {
        possiblyAddDebugRayMeshPoint(_incoming_point, _intersection_point);
 
        _last_elem = _current_elem;
        // Direction changed
        if (_current_elem == ray->currentElem())
        {
          debugRay("RayBC reflected the ray");
          debugRay("  new direction = ", ray->direction());
          traceAssert(ray->direction() *
                              _study.getSideNormal(_current_elem, _intersected_side, _tid) <
                          TRACE_TOLERANCE,
                      "Reflected ray is not incoming");
 
          _incoming_point = _intersection_point;
          _incoming_side = _intersected_side;
          ray->setCurrentPoint(_incoming_point);
          ray->setCurrentIncomingSide(_incoming_side);
        }
        // Position changed (PeriodicRayBC)
        else
        {
          debugRay("RayBC moved the ray");
          debugRay("  new point = ", ray->currentPoint());
          debugRay("  new pid = ", ray->currentElem()->processor_id());
          debugRay("  new elem id = ", ray->currentElem()->id());
          debugRay("  new side = ", ray->currentIncomingSide());
 
          _current_elem = ray->currentElem();
          _incoming_point = ray->currentPoint();
          _incoming_side = ray->currentIncomingSide();
          _intersected_extrema.invalidate();
 
          if (_current_elem->processor_id() != _pid)
          {
            continueTraceOffProcessor(ray);
            return;
          }
        }
        onTrajectoryChanged(ray);
      }
    }
 
    onContinueTrace(ray);
 
  } while (true);
 
  // If a trace made its way down here and didn't return... it failed
  failTrace("Could not find an intersection", _study.tolerateFailure(), __LINE__);
}

Member Data Documentation

_backface_culling

bool TraceRay::_backface_culling

private

Whether or not to use element normals for backface culling.

Definition at line 449 of file TraceRay.h.

Referenced by exitsElem(), moveThroughNeighbor(), setBackfaceCulling(), and trace().

_boundary_elems

std::vector<TraceRayBndElement> TraceRay::_boundary_elems

private

Boundary elements that need RayBCs to be applied.

Definition at line 497 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), onBoundary(), and possiblyAddToBoundaryElems().

_boundary_ids

std::vector<BoundaryID> TraceRay::_boundary_ids

private

Reusable vector for calling _boundary_info.boundary_ids()

Definition at line 495 of file TraceRay.h.

Referenced by applyOnExternalBoundary().

_boundary_info

const BoundaryInfo& TraceRay::_boundary_info

private

The BoundaryInfo for the mesh.

Definition at line 442 of file TraceRay.h.

Referenced by applyOnExternalBoundary().

_current_bnd_id

BoundaryID TraceRay::_current_bnd_id

private

The current BoundaryID (used when calling RayBoundaryConditionBase::onBoundary())

Definition at line 486 of file TraceRay.h.

Referenced by currentBoundaryID(), and onBoundary().

_current_cached_trace

TraceData* TraceRay::_current_cached_trace

private

The TraceData for the current cached trace (if any)

Definition at line 452 of file TraceRay.h.

Referenced by continueTraceOffProcessor(), onCompleteTrace(), onContinueTrace(), onTrajectoryChanged(), and trace().

_current_elem

const Elem* TraceRay::_current_elem

private

The element the current Ray is being traced in.

Definition at line 457 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), continueTraceOffProcessor(), currentElem(), failTraceMessage(), findExternalBoundarySide(), onBoundary(), onSegment(), onSubdomainChanged(), onTrajectoryChanged(), and trace().

_current_elem_n_sides

unsigned short TraceRay::_current_elem_n_sides

private

The number of sides on the current elem, used to avoid elem->n_sides() virtual calls.

Definition at line 467 of file TraceRay.h.

Referenced by failTraceMessage(), findExternalBoundarySide(), onSubdomainChanged(), and trace().

_current_elem_type

libMesh::ElemType TraceRay::_current_elem_type

private

The current elem type (constant on subdomain), used to avoid elem->type() calls.

Definition at line 465 of file TraceRay.h.

Referenced by failTraceMessage(), onSubdomainChanged(), preExecute(), and trace().

_current_normals

const Point* TraceRay::_current_normals

private

The normals for the current element for backface culling (pointer to the first normal - optional)

Definition at line 492 of file TraceRay.h.

Referenced by trace().

_current_ray

const std::shared_ptr<Ray>* TraceRay::_current_ray

private

The current ray being traced.

Definition at line 455 of file TraceRay.h.

Referenced by currentRay(), and trace().

_current_subdomain_hmax

Real TraceRay::_current_subdomain_hmax

private

The current subdomain hmax.

Definition at line 463 of file TraceRay.h.

Referenced by onSubdomainChanged(), and subdomainHmax().

_current_subdomain_id

SubdomainID TraceRay::_current_subdomain_id

private

The current SubdomainID.

Definition at line 461 of file TraceRay.h.

Referenced by currentSubdomainID(), failTraceMessage(), onBoundary(), onSegment(), onSubdomainChanged(), preExecute(), subdomainHmax(), and trace().

_debug_comm

Parallel::Communicator TraceRay::_debug_comm

private

Definition at line 534 of file TraceRay.h.

Referenced by trace().

_debug_mesh

Mesh* TraceRay::_debug_mesh

private

Definition at line 533 of file TraceRay.h.

Referenced by trace().

_debug_node_count

unsigned int TraceRay::_debug_node_count

private

Definition at line 535 of file TraceRay.h.

Referenced by trace().

_dim

const unsigned int TraceRay::_dim

private

The mesh dimension.

Definition at line 440 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), findExternalBoundarySide(), and trace().

_edge_neighbors

std::unordered_map<std::pair<const Node *, const Node *>, std::pair<bool, std::vector<NeighborInfo> > > TraceRay::_edge_neighbors

private

The cached edge neighbors.

Definition at line 504 of file TraceRay.h.

Referenced by getEdgeNeighbors(), and meshChanged().

_elem_side_builder

libMesh::ElemSideBuilder TraceRay::_elem_side_builder

private

Helper for building element sides without excessive allocation.

Definition at line 517 of file TraceRay.h.

Referenced by exitsElem(), and onSegment().

_exits_elem

bool TraceRay::_exits_elem

private

Whether or not the current trace exits an element.

Definition at line 489 of file TraceRay.h.

Referenced by failTraceMessage(), and trace().

_has_ray_kernels

bool TraceRay::_has_ray_kernels

private

Whether or not the RayTracingStudy has any RayKernels.

Definition at line 525 of file TraceRay.h.

Referenced by onSubdomainChanged(), and preExecute().

_incoming_point

Point TraceRay::_incoming_point

private

The incoming point of the current Ray.

Definition at line 469 of file TraceRay.h.

Referenced by applyOnInternalBoundary(), continueTraceOffProcessor(), currentIncomingPoint(), exitsElem(), failTraceMessage(), moveThroughNeighbors(), onCompleteTrace(), onSegment(), and trace().

_incoming_side

unsigned short TraceRay::_incoming_side

private

The incoming side of the current Ray.

Definition at line 471 of file TraceRay.h.

Referenced by applyOnInternalBoundary(), continueTraceOffProcessor(), currentIncomingSide(), exitsElem(), failTraceMessage(), onSegment(), and trace().

_intersected_extrema

ElemExtrema TraceRay::_intersected_extrema

private

The work point for the intersected vertex/edge vertices of the current Ray, if any.

Definition at line 480 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), currentIntersectedExtrema(), failTraceMessage(), findExternalBoundarySide(), moveThroughNeighbors(), onBoundary(), onTrajectoryChanged(), and trace().

_intersected_side

unsigned short TraceRay::_intersected_side

private

The work point for the intersected side of the current Ray.

Definition at line 478 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), currentIntersectedSide(), failTraceMessage(), findExternalBoundarySide(), moveThroughNeighbors(), onBoundary(), onSegment(), and trace().

_intersection_distance

Real TraceRay::_intersection_distance

private

The work point for the intersection distance of the current Ray.

Definition at line 484 of file TraceRay.h.

Referenced by continueTraceOffProcessor(), currentIntersectionDistance(), moveThroughNeighbors(), onCompleteTrace(), onContinueTrace(), onSegment(), onTrajectoryChanged(), and trace().

_intersection_point

Point TraceRay::_intersection_point

private

The work point for the intersection of the current Ray.

Definition at line 476 of file TraceRay.h.

Referenced by applyOnExternalBoundary(), applyOnInternalBoundary(), continueTraceOffProcessor(), currentIntersectionPoint(), failTraceMessage(), findExternalBoundarySide(), moveThroughNeighbors(), onBoundary(), onCompleteTrace(), onSegment(), possiblyAddToBoundaryElems(), and trace().

_is_rectangular_domain

bool TraceRay::_is_rectangular_domain

private

Whether or not the domain is rectangular (defined perfectly by its bounding box)

Definition at line 527 of file TraceRay.h.

Referenced by preExecute(), and trace().

_last_elem

const Elem* TraceRay::_last_elem

private

The last element the current Ray was traced in.

Definition at line 459 of file TraceRay.h.

Referenced by applyOnInternalBoundary(), and trace().

_last_intersected_extrema

ElemExtrema TraceRay::_last_intersected_extrema

private

The intersected vertex/edge vertices for the previous intersection, if any.

Definition at line 482 of file TraceRay.h.

Referenced by trace().

_mesh

MooseMesh& TraceRay::_mesh

private

The MooseMesh.

Definition at line 438 of file TraceRay.h.

Referenced by onBoundary().

_neighbor_active_neighbor_children

std::vector<const Elem *> TraceRay::_neighbor_active_neighbor_children

private

Definition at line 511 of file TraceRay.h.

Referenced by getEdgeNeighbors(), getPointNeighbors(), and getVertexNeighbors().

_neighbor_next_untested_set

MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> TraceRay::_neighbor_next_untested_set

private

Definition at line 510 of file TraceRay.h.

Referenced by getEdgeNeighbors(), getPointNeighbors(), and getVertexNeighbors().

_neighbor_set

MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> TraceRay::_neighbor_set

private

Definition at line 508 of file TraceRay.h.

Referenced by getEdgeNeighbors(), getPointNeighbors(), and getVertexNeighbors().

_neighbor_untested_set

MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> TraceRay::_neighbor_untested_set

private

Definition at line 509 of file TraceRay.h.

Referenced by getEdgeNeighbors(), getPointNeighbors(), and getVertexNeighbors().

_old_ray_kernels

std::set<RayKernelBase *> TraceRay::_old_ray_kernels

private

Helper for avoiding calling preTrace() on the same RayKernel multiple times.

Definition at line 530 of file TraceRay.h.

Referenced by onSubdomainChanged().

_on_boundary_apply_index

std::vector<std::size_t> TraceRay::_on_boundary_apply_index

private

Reusable for which boundary elements to apply for a specific RayBC in onBoundary()

Definition at line 522 of file TraceRay.h.

Referenced by onBoundary().

_on_boundary_ray_bcs

std::vector<RayBoundaryConditionBase *> TraceRay::_on_boundary_ray_bcs

private

Reusable for getting the RayBCs in onBoundary()

Definition at line 520 of file TraceRay.h.

Referenced by onBoundary().

_pid

const processor_id_type TraceRay::_pid

private

The processor id.

Definition at line 444 of file TraceRay.h.

Referenced by continueTraceOffProcessor(), failTraceMessage(), and trace().

_point_neighbor_helper

std::vector<NeighborInfo> TraceRay::_point_neighbor_helper

private

Reusable for building neighbors.

Definition at line 507 of file TraceRay.h.

Referenced by getPointNeighbors().

_results

std::vector<unsigned long long int> TraceRay::_results

private

Results over all of the local traces, indexed by TraceRayResult.

Definition at line 514 of file TraceRay.h.

Referenced by exitsElem(), failTrace(), getEdgeNeighbors(), getPointNeighbors(), getVertexNeighbors(), moveThroughNeighbors(), preExecute(), results(), storeExitsElemResult(), and trace().

_should_continue

bool TraceRay::_should_continue

private

Whether or not the current Ray should continue.

Definition at line 473 of file TraceRay.h.

Referenced by failTrace(), postRayTracingObject(), and trace().

_study

RayTracingStudy& TraceRay::_study

private

The RayTracingStudy.

Definition at line 436 of file TraceRay.h.

Referenced by applyOnInternalBoundary(), continueTraceOffProcessor(), exitsElem(), exitsElem(), findExternalBoundarySide(), getEdgeNeighbors(), getPointNeighbors(), getVertexNeighbors(), moveThroughNeighbor(), onBoundary(), onCompleteTrace(), onContinueTrace(), onSegment(), onSubdomainChanged(), onTrajectoryChanged(), possiblyAddToBoundaryElems(), preExecute(), subdomainHmax(), and trace().

_tid

const THREAD_ID TraceRay::_tid

private

The thread id.

Definition at line 446 of file TraceRay.h.

Referenced by exitsElem(), failTraceMessage(), findExternalBoundarySide(), getEdgeNeighbors(), getPointNeighbors(), getVertexNeighbors(), moveThroughNeighbor(), onBoundary(), onCompleteTrace(), onSegment(), onSubdomainChanged(), preExecute(), and trace().

_vertex_neighbors

std::unordered_map<const Node *, std::vector<NeighborInfo> > TraceRay::_vertex_neighbors

private

The cached vertex neighbors.

Definition at line 500 of file TraceRay.h.

Referenced by getVertexNeighbors(), and meshChanged().

---

The documentation for this class was generated from the following files:

• ray_tracing/include/raytracing/TraceRay.h
• ray_tracing/src/raytracing/TraceRay.C