TraceRayTools.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

// MOOSE includes
#include "StaticallyAllocatedSet.h"
#include "Conversion.h"
#include "MooseTypes.h"
#include "libMeshReducedNamespace.h"

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

// libMesh includes
#include "libmesh/face.h"
#include "libmesh/cell_hex.h"
#include "libmesh/cell_prism.h"
#include "libmesh/cell_pyramid.h"
#include "libmesh/cell_tet.h"
#include "libmesh/remote_elem.h"

namespace libMesh
{
class Mesh;
class Edge;
class Cell;
}

namespace TraceRayTools
{

extern const std::set<int> TRACEABLE_ELEMTYPES;
extern const std::set<int> ADAPTIVITY_TRACEABLE_ELEMTYPES;

const Real TRACE_TOLERANCE = 1e-8;
const Real LOOSE_TRACE_TOLERANCE = 1e-5;

const int INVALID_INT = std::numeric_limits<int>::max();

enum SegmentVertices
{
  SEGMENT_VERTEX_0 = 0,
  SEGMENT_VERTEX_1 = 1,
  SEGMENT_VERTEX_NONE = 7
};

bool lineLineIntersect2D(const Point & start,
                         const Point & direction,
                         const Real length,
                         const Point & v0,
                         const Point & v1,
                         Point & intersection_point,
                         Real & intersection_distance,
                         SegmentVertices & segment_vertex
#ifdef DEBUG_RAY_INTERSECTIONS
                         ,
                         const bool debug
#endif
);

bool isWithinSegment(const Point & segment1,
                     const Point & segment2,
                     const Point & point,
                     const Real tolerance = TRACE_TOLERANCE);

bool isWithinSegment(const Point & segment1,
                     const Point & segment2,
                     const Real segment_length,
                     const Point & point,
                     const Real tolerance = TRACE_TOLERANCE);

void findPointNeighbors(
    const Elem * const elem,
    const Point & point,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
    std::vector<const Elem *> active_neighbor_children,
    std::vector<NeighborInfo> & info);

void findNodeNeighbors(
    const Elem * const elem,
    const Node * const node,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
    std::vector<const Elem *> active_neighbor_children,
    std::vector<NeighborInfo> & info);

void findEdgeNeighbors(
    const Elem * const elem,
    const Node * const node1,
    const Node * const node2,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
    std::vector<const Elem *> active_neighbor_children,
    std::vector<NeighborInfo> & info);

template <typename KeepFunctor>
void
findNeighbors(
    const Elem * const elem,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & neighbor_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & untested_set,
    MooseUtils::StaticallyAllocatedSet<const Elem *, MAX_POINT_NEIGHBORS> & next_untested_set,
    std::vector<const Elem *> active_neighbor_children,
    KeepFunctor & keep_functor)
{
  mooseAssert(elem->active(), "Inactive element");

  neighbor_set.clear();
  untested_set.clear();
  next_untested_set.clear();

  untested_set.insert(elem);

  while (!untested_set.empty())
  {
    // Loop over all the elements in the patch that haven't already been tested
    for (const Elem * const test_elem : untested_set)
      for (const auto * const current_neighbor : test_elem->neighbor_ptr_range())
        if (current_neighbor && current_neighbor != remote_elem &&
            current_neighbor != elem) // we have a real neighbor on elem side
        {
          if (current_neighbor->active()) // ... if it is active
          {
            if (!neighbor_set.contains(current_neighbor) && keep_functor(current_neighbor))
            {
              next_untested_set.insert(current_neighbor);
              neighbor_set.insert(current_neighbor);
            }
          }
#ifdef LIBMESH_ENABLE_AMR
          else // ... the neighbor is *not* active,
          {    // ... so add *all* neighboring
               // active children that touch p
            current_neighbor->active_family_tree_by_neighbor(active_neighbor_children, test_elem);

            for (const Elem * current_child : active_neighbor_children)
              if (!neighbor_set.contains(current_child) && keep_functor(current_child) &&
                  current_child != elem)
              {
                next_untested_set.insert(current_child);
                neighbor_set.insert(current_child);
              }
          }
#endif // #ifdef LIBMESH_ENABLE_AMR
        }
    untested_set.swap(next_untested_set);
    next_untested_set.clear();
  }
}

template <typename T>
bool
findEdgeNeighborsWithinEdgeInternal(const Elem * const candidate,
                                    const Elem * const elem,
                                    const Node * const vertex1,
                                    const Node * const vertex2,
                                    const Real edge_length,
                                    std::vector<NeighborInfo> & info)
{
  // If we have the type already: use it so we can avoid all of the virtual
  // calls that we would be making on Elem
  mooseAssert(dynamic_cast<const T *>(candidate), "Incorrect elem type");
  const T * const T_candidate = static_cast<const T *>(candidate);

  // Local index of our two vertices of interest (if any)
  auto v1 = T_candidate->get_node_index(vertex1);
  auto v2 = T_candidate->get_node_index(vertex2);
  // Whether or not we have said vertices
  const bool has_v1 = v1 != invalid_uint;
  const bool has_v2 = v2 != invalid_uint;

  // If the candidate has both vertices, it shares the complete edge
  if (has_v1 && has_v2)
  {
    // Add the sides that contain said edge
    for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
      if (T_candidate->is_node_on_edge(v1, e) && T_candidate->is_node_on_edge(v2, e))
      {
        std::vector<unsigned short> sides(2);
        sides[0] = T::edge_sides_map[e][0];
        sides[1] = T::edge_sides_map[e][1];
        info.emplace_back(candidate, std::move(sides), 0, 1);
        return true;
      }

    mooseError("Failed to find a side that the vertices are on");
  }

  const auto n_vertices = T_candidate->n_vertices();

  // If we only have one of the vertices, we can still be contained within the edge if we have
  // another vertex that is within the edge
  if (has_v1 || has_v2)
  {
    // Local index of the vertex that the candidate and the target edge have in common
    const auto common_v = has_v1 ? v1 : v2;

    // See if another vertex that isn't the common node is contained
    MooseIndex(n_vertices) other_v;
    for (other_v = 0; other_v < n_vertices; ++other_v)
      if (other_v != common_v &&
          isWithinSegment(*vertex1, *vertex2, edge_length, T_candidate->point(other_v)))
        break;

    // If we have the common vertex and another vertex within the target edge, use the sides
    // that contain both of those vertices
    if (other_v != n_vertices)
    {
      for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
        if (T_candidate->is_node_on_edge(common_v, e) && T_candidate->is_node_on_edge(other_v, e))
        {
          std::vector<unsigned short> sides(2);
          sides[0] = T::edge_sides_map[e][0];
          sides[1] = T::edge_sides_map[e][1];

          info.emplace_back(
              candidate,
              std::move(sides),
              has_v1 ? 0 : (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length,
              has_v1 ? (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length : 1);
          return true;
        }
      mooseError("Failed to find a side that the vertices are on");
    }
    else if (T_candidate->level() < elem->level())
    {
      for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
        if (T_candidate->is_node_on_edge(common_v, e) &&
            isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
                            T_candidate->point(T::edge_nodes_map[e][1]),
                            has_v1 ? *vertex2 : *vertex1))
        {
          other_v = T::edge_nodes_map[e][0] == common_v ? T::edge_nodes_map[e][1]
                                                        : T::edge_nodes_map[e][0];
          std::vector<unsigned short> sides(2);
          sides[0] = T::edge_sides_map[e][0];
          sides[1] = T::edge_sides_map[e][1];

          info.emplace_back(
              candidate,
              std::move(sides),
              has_v1 ? 0 : (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length,
              has_v1 ? (T_candidate->point(other_v) - (Point)*vertex1).norm() / edge_length : 1);
          return true;
        }

      return false;
    }
  }
  // When the candidate level is less, one of our vertices could be in a candidate's face and the
  // other could be within one of the candidate's edges
  else if (T_candidate->level() < elem->level())
  {
    auto v1_edge = RayTracingCommon::invalid_edge;
    auto v2_edge = RayTracingCommon::invalid_edge;

    // See if any of the edges contain one of the vertices
    for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
    {
      if (v1_edge == RayTracingCommon::invalid_edge &&
          isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
                          T_candidate->point(T::edge_nodes_map[e][1]),
                          *vertex1))
        v1_edge = e;
      if (v2_edge == RayTracingCommon::invalid_edge &&
          isWithinSegment(T_candidate->point(T::edge_nodes_map[e][0]),
                          T_candidate->point(T::edge_nodes_map[e][1]),
                          *vertex2))
        v2_edge = e;
    }

    const auto v1_within = v1_edge != RayTracingCommon::invalid_edge;
    const auto v2_within = v2_edge != RayTracingCommon::invalid_edge;

    if (v1_within && v2_within)
    {
      mooseAssert(v1_edge == v2_edge, "Vertices should be contained in same edge");

      std::vector<unsigned short> sides(2);
      sides[0] = T::edge_sides_map[v1_edge][0];
      sides[1] = T::edge_sides_map[v1_edge][1];
      info.emplace_back(candidate, std::move(sides), 0, 1);
      return true;
    }
    else if (v1_within || v2_within)
    {
      const auto in_edge = v1_within ? v1_edge : v2_edge;
      const Point & check_point = v1_within ? *vertex2 : *vertex1;

      std::vector<unsigned short> sides(1);
      Real lower_bound = 0;
      Real upper_bound = 1;

      if (candidate->side_ptr(T::edge_sides_map[in_edge][0])->contains_point(check_point))
        sides[0] = T::edge_sides_map[in_edge][0];
      else if (candidate->side_ptr(T::edge_sides_map[in_edge][1])->contains_point(check_point))
        sides[0] = T::edge_sides_map[in_edge][1];
      else
      {
        const Real point_bound = v1_within ? 0 : 1;
        lower_bound = point_bound;
        upper_bound = point_bound;
        sides[0] = T::edge_sides_map[in_edge][0];
        sides.push_back(T::edge_sides_map[in_edge][1]);
      }

      info.emplace_back(candidate, std::move(sides), lower_bound, upper_bound);
      return true;
    }

    return false;
  }
  // If we have neither of the vertices, and the candidate's edge can't fully contain [vertex1,
  // vertex2], check for other vertices that are contained within the edge
  else
  {
    for (v1 = 0; v1 < n_vertices; ++v1)
      if (isWithinSegment(*vertex1, *vertex2, edge_length, candidate->point(v1)))
      {
        for (v2 = v1 + 1; v2 < n_vertices; ++v2)
          if (isWithinSegment(*vertex1, *vertex2, edge_length, candidate->point(v2)))
            break;
        break;
      }

    // No vertices are contained within the edge
    if (v1 == n_vertices)
      return false;

    // Only one vertex contained: add the sides associated with that vertex
    if (v2 >= n_vertices)
    {
      std::vector<unsigned short> sides;
      sides.reserve(2);

      for (MooseIndex(T::num_sides) s = 0; s < T::num_sides; ++s)
        if (T_candidate->is_node_on_side(v1, s))
          sides.push_back(s);

      if (sides.empty())
        mooseError("Failed to find sides that the vertex is on");

      // Bounds are the same because we only touch at a point
      const auto bound = (T_candidate->point(v1) - (Point)*vertex1).norm() / edge_length;
      info.emplace_back(candidate, std::move(sides), bound, bound);
      return true;
    }
    // Two vertices contained: find the sides that contain both vertices
    else
    {
      for (MooseIndex(T::num_edges) e = 0; e < T::num_edges; ++e)
        if (T_candidate->is_node_on_edge(v1, e) && T_candidate->is_node_on_edge(v2, e))
        {
          auto lower_bound = (T_candidate->point(v1) - (Point)*vertex1).norm() / edge_length;
          auto upper_bound = (T_candidate->point(v2) - (Point)*vertex1).norm() / edge_length;
          if (lower_bound > upper_bound)
          {
            const auto temp = lower_bound;
            lower_bound = upper_bound;
            upper_bound = temp;
          }

          std::vector<unsigned short> sides(2);
          sides[0] = T::edge_sides_map[e][0];
          sides[1] = T::edge_sides_map[e][1];
          info.emplace_back(candidate, std::move(sides), lower_bound, upper_bound);
          return true;
        }

      mooseError("Failed to find sides that the vertices are on");
    }
  }

  return false;
}

const Elem *
childContainingPointOnSide(const Elem * elem, const Point & point, const unsigned short side);

const Elem * getActiveNeighbor(const Elem * elem, const unsigned short side, const Point & point);

bool intersectTriangle(const Point & start,
                       const Point & direction,
                       const Elem * const elem,
                       const unsigned short v0,
                       const unsigned short v1,
                       const unsigned short v2,
                       Real & intersection_distance,
                       ElemExtrema & intersected_extrema,
                       const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                       ,
                       const bool debug
#endif
);

bool intersectQuad(const Point & start,
                   const Point & direction,
                   const Elem * const elem,
                   const unsigned short v00,
                   const unsigned short v10,
                   const unsigned short v11,
                   const unsigned short v01,
                   Real & intersection_distance,
                   ElemExtrema & intersected_extrema,
                   const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                   ,
                   const bool debug
#endif
);

template <typename T>
typename std::enable_if<std::is_base_of<libMesh::Face, T>::value, bool>::type
sideIntersectedByLine(const Elem * elem,
                      const Point & start_point,
                      const Point & direction,
                      const unsigned short side,
                      const Real max_length,
                      Point & intersection_point,
                      Real & intersection_distance,
                      ElemExtrema & intersected_extrema,
                      const Real /* hmax */
#ifdef DEBUG_RAY_INTERSECTIONS
                      ,
                      const bool debug
#endif
)
{
  mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");

  SegmentVertices segment_vertex = SEGMENT_VERTEX_NONE;

  const bool intersected = lineLineIntersect2D(start_point,
                                               direction,
                                               max_length,
                                               elem->point(T::side_nodes_map[side][0]),
                                               elem->point(T::side_nodes_map[side][1]),
                                               intersection_point,
                                               intersection_distance,
                                               segment_vertex
#ifdef DEBUG_RAY_INTERSECTIONS
                                               ,
                                               debug
#endif
  );

  if (segment_vertex != SEGMENT_VERTEX_NONE)
  {
    intersected_extrema.setVertex(T::side_nodes_map[side][segment_vertex]);
    mooseAssert(
        intersected_extrema.vertexPoint(elem).absolute_fuzzy_equals(intersection_point,
                                                                    3 * TRACE_TOLERANCE),
        "Doesn't intersect vertex at: " + Moose::stringify(intersected_extrema.vertexPoint(elem)) +
            " tentative intersection point: " + Moose::stringify(intersection_point));
  }

  return intersected;
}

template <typename T>
typename std::enable_if<std::is_base_of<Hex, T>::value, bool>::type
sideIntersectedByLine(const Elem * elem,
                      const Point & start_point,
                      const Point & direction,
                      const unsigned short side,
                      const Real,
                      Point & intersection_point,
                      Real & intersection_distance,
                      ElemExtrema & intersected_extrema,
                      const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                      ,
                      const bool debug
#endif
)
{
  mooseAssert(elem->first_order_equivalent_type(elem->type()) == HEX8, "Not a Hex");
  mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");

  const bool intersected = intersectQuad(start_point,
                                         direction,
                                         elem,
                                         T::side_nodes_map[side][3],
                                         T::side_nodes_map[side][2],
                                         T::side_nodes_map[side][1],
                                         T::side_nodes_map[side][0],
                                         intersection_distance,
                                         intersected_extrema,
                                         hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                         ,
                                         debug
#endif
  );

  if (intersected)
    intersection_point = start_point + intersection_distance * direction;

  return intersected;
}

template <typename T>
typename std::enable_if<std::is_base_of<Tet, T>::value, bool>::type
sideIntersectedByLine(const Elem * elem,
                      const Point & start_point,
                      const Point & direction,
                      const unsigned short side,
                      const Real /* max_length */,
                      Point & intersection_point,
                      Real & intersection_distance,
                      ElemExtrema & intersected_extrema,
                      const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                      ,
                      const bool debug
#endif
)
{
  mooseAssert(elem->first_order_equivalent_type(elem->type()) == TET4, "Not a Tet");
  mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");

  const bool intersected = intersectTriangle(start_point,
                                             direction,
                                             elem,
                                             T::side_nodes_map[side][2],
                                             T::side_nodes_map[side][1],
                                             T::side_nodes_map[side][0],
                                             intersection_distance,
                                             intersected_extrema,
                                             hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                             ,
                                             debug
#endif
  );

  if (intersected)
    intersection_point = start_point + intersection_distance * direction;

  return intersected;
}

template <typename T>
typename std::enable_if<std::is_base_of<libMesh::Pyramid, T>::value, bool>::type
sideIntersectedByLine(const Elem * elem,
                      const Point & start_point,
                      const Point & direction,
                      const unsigned short side,
                      const Real /* max_length */,
                      Point & intersection_point,
                      Real & intersection_distance,
                      ElemExtrema & intersected_extrema,
                      const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                      ,
                      const bool debug
#endif
)
{
  mooseAssert(elem->first_order_equivalent_type(elem->type()) == PYRAMID5, "Not a Pyramid");
  mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");

  const bool intersected = side < 4 ? intersectTriangle(start_point,
                                                        direction,
                                                        elem,
                                                        T::side_nodes_map[side][2],
                                                        T::side_nodes_map[side][1],
                                                        T::side_nodes_map[side][0],
                                                        intersection_distance,
                                                        intersected_extrema,
                                                        hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                                        ,
                                                        debug
#endif
                                                        )
                                    : intersectQuad(start_point,
                                                    direction,
                                                    elem,
                                                    T::side_nodes_map[side][3],
                                                    T::side_nodes_map[side][2],
                                                    T::side_nodes_map[side][1],
                                                    T::side_nodes_map[side][0],
                                                    intersection_distance,
                                                    intersected_extrema,
                                                    hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                                    ,
                                                    debug
#endif
                                      );

  if (intersected)
    intersection_point = start_point + intersection_distance * direction;

  return intersected;
}

template <typename T>
typename std::enable_if<std::is_base_of<libMesh::Prism, T>::value, bool>::type
sideIntersectedByLine(const Elem * elem,
                      const Point & start_point,
                      const Point & direction,
                      const unsigned short side,
                      const Real /* max_length */,
                      Point & intersection_point,
                      Real & intersection_distance,
                      ElemExtrema & intersected_extrema,
                      const Real hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                      ,
                      const bool debug
#endif
)
{
  mooseAssert(elem->first_order_equivalent_type(elem->type()) == PRISM6, "Not a Prism");
  mooseAssert(intersection_point == RayTracingCommon::invalid_point, "Point should be invalid");
  mooseAssert(intersected_extrema.isInvalid(), "Should be invalid");

  const bool intersected = (side == 0 || side == 4) ? intersectTriangle(start_point,
                                                                        direction,
                                                                        elem,
                                                                        T::side_nodes_map[side][2],
                                                                        T::side_nodes_map[side][1],
                                                                        T::side_nodes_map[side][0],
                                                                        intersection_distance,
                                                                        intersected_extrema,
                                                                        hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                                                        ,
                                                                        debug
#endif
                                                                        )
                                                    : intersectQuad(start_point,
                                                                    direction,
                                                                    elem,
                                                                    T::side_nodes_map[side][3],
                                                                    T::side_nodes_map[side][2],
                                                                    T::side_nodes_map[side][1],
                                                                    T::side_nodes_map[side][0],
                                                                    intersection_distance,
                                                                    intersected_extrema,
                                                                    hmax
#ifdef DEBUG_RAY_INTERSECTIONS
                                                                    ,
                                                                    debug
#endif
                                                      );

  if (intersected)
    intersection_point = start_point + intersection_distance * direction;

  return intersected;
}

bool isTraceableElem(const Elem * elem);
bool isAdaptivityTraceableElem(const Elem * elem);

unsigned short atVertex(const Elem * elem, const Point & point);

unsigned short atVertexOnSide(const Elem * elem, const Point & point, const unsigned short side);

template <typename T>
inline typename std::enable_if<!std::is_base_of<libMesh::Pyramid, T>::value &&
                                   !std::is_base_of<libMesh::Prism, T>::value,
                               unsigned short>::type
nodesPerSide(const unsigned short)
{
  return T::nodes_per_side;
}

template <typename T>
inline typename std::enable_if<std::is_base_of<Pyramid, T>::value, unsigned short>::type
nodesPerSide(const unsigned short side)
{
  return T::nodes_per_side - (side != 4);
}

template <typename T>
inline typename std::enable_if<std::is_base_of<libMesh::Prism, T>::value, unsigned short>::type
nodesPerSide(const unsigned short side)
{
  return T::nodes_per_side - (side == 0 || side == 4);
}

template <typename T>
typename std::enable_if<!std::is_base_of<Edge, T>::value, unsigned short>::type
atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side);

template <typename T>
typename std::enable_if<std::is_base_of<Edge, T>::value, unsigned short>::type
atVertexOnSideTempl(const Elem * elem, const Point & point, const unsigned short side);

template <typename T>
bool withinEdgeTempl(const Elem * elem,
                     const Point & point,
                     ElemExtrema & extrema,
                     const Real tolerance = TRACE_TOLERANCE);

bool withinEdge(const Elem * elem,
                const Point & point,
                ElemExtrema & extrema,
                const Real tolerance = TRACE_TOLERANCE);

bool withinEdgeOnSide(const Elem * const elem,
                      const Point & point,
                      const unsigned short side,
                      ElemExtrema & extrema);

bool withinExtremaOnSide(const Elem * const elem,
                         const Point & point,
                         const unsigned short side,
                         const unsigned int dim,
                         ElemExtrema & extrema);

template <typename T>
typename std::enable_if<std::is_base_of<libMesh::Cell, T>::value, bool>::type withinEdgeOnSideTempl(
    const Elem * const elem, const Point & point, const unsigned short side, ElemExtrema & extrema);

template <typename T>
typename std::enable_if<!std::is_base_of<libMesh::Cell, T>::value, bool>::type
withinEdgeOnSideTempl(const Elem * const, const Point &, const unsigned short, ElemExtrema &)
{
  mooseError("Should not call withinEdgeOnSideTempl() with a non-Cell derived Elem");
}

bool onBoundingBoxBoundary(const BoundingBox & bbox,
                           const Point & point,
                           const unsigned int dim,
                           const Real tolerance);

}