TraceRay.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

// Local includes
#include "DebugRay.h"
#include "ElemExtrema.h"
#include "NeighborInfo.h"
#include "RayTracingCommon.h"
#include "TraceRayBndElement.h"

// MOOSE Includes
#include "MooseMesh.h"
#include "MooseHashing.h"
#include "MooseTypes.h"
#include "StaticallyAllocatedSet.h"

// libMesh includes
#include "libmesh/enum_elem_type.h"
#include "libmesh/elem_side_builder.h"

// Forward declarations
class Ray;
class RayBoundaryConditionBase;
class RayKernelBase;
class RayTracingStudy;
class RayTracingObject;
struct TraceData;
namespace libMesh
{
class BoundaryInfo;
class Edge;
class Mesh;
}

class TraceRay
{
public:
  TraceRay(RayTracingStudy & study, const THREAD_ID tid);

  virtual ~TraceRay() {}

  enum TraceRayResult
  {
    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
  };

  void preExecute();

  void meshChanged();

  void trace(const std::shared_ptr<Ray> & ray);

  std::vector<unsigned long long int> results() const { return _results; }

  void setBackfaceCulling(const bool backface_culling) { _backface_culling = backface_culling; }

  const std::shared_ptr<Ray> * const & currentRay() const { return _current_ray; }
  const Elem * const & currentElem() const { return _current_elem; }
  const Point & currentIntersectionPoint() const { return _intersection_point; }
  const Point & currentIncomingPoint() const { return _incoming_point; }
  const unsigned short & currentIncomingSide() const { return _incoming_side; }
  const unsigned short & currentIntersectedSide() const { return _intersected_side; }
  const ElemExtrema & currentIntersectedExtrema() const { return _intersected_extrema; }
  const Real & currentIntersectionDistance() const { return _intersection_distance; }
  const BoundaryID & currentBoundaryID() const { return _current_bnd_id; }
  const SubdomainID & currentSubdomainID() const { return _current_subdomain_id; }

  const std::vector<NeighborInfo> & getVertexNeighbors(const Elem * elem, const Node * vertex);
  const std::vector<NeighborInfo> & getVertexNeighbors(const Elem * elem,
                                                       const unsigned short vertex);
  const std::vector<NeighborInfo> &
  getEdgeNeighbors(const Elem * elem,
                   const std::pair<const Node *, const Node *> & vertices,
                   const Point & point);
  const std::vector<NeighborInfo> &
  getEdgeNeighbors(const Elem * elem,
                   const std::pair<unsigned short, unsigned short> & vertices,
                   const Point & point);
  const std::vector<NeighborInfo> & getPointNeighbors(const Elem * elem, const Point & point);
  const std::vector<NeighborInfo> &
  getNeighbors(const Elem * elem, const ElemExtrema & extrema, const Point & point);

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

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

  void onCompleteTrace(const std::shared_ptr<Ray> & ray);
  void onContinueTrace(const std::shared_ptr<Ray> & /* ray */);
  void onTrajectoryChanged(const std::shared_ptr<Ray> & ray);
  void onSubdomainChanged(const std::shared_ptr<Ray> & ray, const bool same_ray);

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

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

  enum ExitsElemResult
  {
    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
  };

  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);

  template <typename T, typename FirstOrderT>
  typename 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);

  template <typename T, typename FirstOrderT>
  typename 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);

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

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

  void applyOnExternalBoundary(const std::shared_ptr<Ray> & ray);
  void applyOnInternalBoundary(const std::shared_ptr<Ray> & ray);
  void possiblyAddToBoundaryElems(const Elem * elem,
                                  const unsigned short side,
                                  const std::vector<BoundaryID> & bnd_ids,
                                  const ElemExtrema & extrema);
  void continueTraceOffProcessor(const std::shared_ptr<Ray> & ray);

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

  void storeExitsElemResult(const ExitsElemResult result);

  Real subdomainHmax(const Elem * elem) const;

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

  RayTracingStudy & _study;
  MooseMesh & _mesh;
  const unsigned int _dim;
  const BoundaryInfo & _boundary_info;
  const processor_id_type _pid;
  const THREAD_ID _tid;

  bool _backface_culling;

  TraceData * _current_cached_trace;

  const std::shared_ptr<Ray> * _current_ray;
  const Elem * _current_elem;
  const Elem * _last_elem;
  SubdomainID _current_subdomain_id;
  Real _current_subdomain_hmax;
  libMesh::ElemType _current_elem_type;
  unsigned short _current_elem_n_sides;
  Point _incoming_point;
  unsigned short _incoming_side;
  bool _should_continue;

  Point _intersection_point;
  unsigned short _intersected_side;
  ElemExtrema _intersected_extrema;
  ElemExtrema _last_intersected_extrema;
  Real _intersection_distance;
  BoundaryID _current_bnd_id;

  bool _exits_elem;

  const Point * _current_normals;

  std::vector<BoundaryID> _boundary_ids;
  std::vector<TraceRayBndElement> _boundary_elems;

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

  std::vector<NeighborInfo> _point_neighbor_helper;
  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;

  libMesh::ElemSideBuilder _elem_side_builder;

  std::vector<RayBoundaryConditionBase *> _on_boundary_ray_bcs;
  std::vector<std::size_t> _on_boundary_apply_index;

  bool _has_ray_kernels;
  bool _is_rectangular_domain;

  std::set<RayKernelBase *> _old_ray_kernels;

#ifdef DEBUG_RAY_MESH_IF
  Mesh * _debug_mesh;
  Parallel::Communicator _debug_comm;
  unsigned int _debug_node_count;
#endif
};

// Helper for a mooseAssert() with useful trace information
#ifdef NDEBUG
#define traceAssert(asserted, msg) ((void)0)
#else
#define traceAssert(asserted, msg)                                                                 \
  do                                                                                               \
  {                                                                                                \
    if (!(asserted))                                                                               \
      mooseAssert(asserted, failTraceMessage(msg, __LINE__));                                      \
  } while (0)
#endif