AdaptiveRayContainmentCheck Class Reference

Ray-casting point-in-solid engine over a closed surface mesh represented as a collection of SurfaceElement wrappers. More...

#include <AdaptiveRayContainmentCheck.h>

Public Member Functions
	AdaptiveRayContainmentCheck (const std::vector< std::unique_ptr< SurfaceElement > > &bd_elements, const std::vector< Point > &centroids, const SurfaceGeometry::RayDirectionOptions &ray_options, const Real eps_on_surface=libMesh::TOLERANCE, const int leaf_max_size=10, const FileName &obb_file_name="", const FileName &ray_file_name="", const libMesh::Parallel::Communicator *comm=nullptr)
SurfaceGeometry::SurfaceSide	sideness (const Point &p) const
	Main function: Determine if a point is inside the geometry.
const Point &	rayDirection () const
	The resolved ray direction actually used for shooting: the (normalized) user direction for a user-selected ray, or the PCA-selected direction for an auto ray.

Private Member Functions
bool	rayIntersectGeometry (const Point &ray_start, const Point &ray_end, const SurfaceElement *elem) const
	Ray-element intersection (e.g., ray-line for 2D, ray-triangle for 3D)
bool	isOutsideBoundingBox (const Point &query_point) const
	Check if point is outside global bounding box.
bool	isOutsideRayBBox (const Point &orig, const Point &dir, const Ball &ball) const
	Check if element center is outside ray bounding box.
bool	isOutsideBoundingRegion (const Point &orig, const Point &dir, const Ball &ball) const
	Check if element center is outside ray bounding circle/sphere.
BoundingBox	computeGlobalBoundingBox ()
	Compute the global bounding box of all boundary elements.
void	preparePCASVD ()
void	initializeRayDirection ()
	Finalizes the ray direction and the matching bounding box.
std::optional< SurfaceGeometry::SurfaceSide >	sidenessFromRayPair (const Point &p, const std::array< Point, 2 > &ray_starts) const
	Determine sideness from a pair of opposite rays.
bool	isOnSurface (const Point &p) const
	True if p lies on the surface (within _eps_on_surface), i.e.
int	countCrossings (const Point &ray_start, const Point &ray_end, const bool use_primary_direction=true) const
	Count how many times the segment from ray_start to ray_end crosses the surface.
template<typename CrossingTest >
int	countFilteredCrossings (const Point &ray_start, const Point &ray_end, const bool use_primary_direction, CrossingTest is_crossing) const
	Shared traversal for the 2D and 3D crossing counts: walk the candidate elements (KD-tree candidates for the primary direction, all elements otherwise), apply the ray-bbox and bounding-region pruning, and count those for which is_crossing(surface) is true.
int	countCrossings2D (const Point &ray_start, const Point &ray_end, const bool use_primary_direction) const
	2D crossing count using a half-open side-based crossing rule: an edge is counted when its two endpoints lie on strictly opposite sides of the ray line (a vertex or collinear edge exactly on the line is thus counted consistently, with no tolerance) and the crossing is on the ray_start side of the query point.
Point	rayStartOutsideAABB (const Point &point, const Point &unit_direction, const bool inverted) const
	Ray start strictly outside the global AABB along unit_direction, for any direction.
Point	rayStartOutsideOBB (const Point &point, const Point &ray_direction, const unsigned int obb_axis, const bool inverted=false) const
	Computes the starting point of an OBB-based ray (auto/PCA policy) for a given query point.
Point	projectPointOntoPlane (const Point &point_to_project, const Point &plane_point, const Point &plane_normal) const
	Orthogonally project point_to_project onto the plane defined by plane_point and unit normal plane_normal.
void	buildObbKdtreeAndMaxProjectedDiagonal (const Real expand_box_length)
	Constructs an oriented bounding box (OBB) using the results of PCA and the KD-tree.
std::vector< unsigned int >	collectCandidateElementIDs (const Point &query_point) const
	Use the kd-tree to collect candidate element IDs to check intersections.

Private Attributes
Point	_centroid_nodal_points
	The centroid of the boundary elements' node points (prepare inside this class).
std::vector< Point >	_projected_centroids
	Projected centroids of the elements in the boundary mesh.
std::unique_ptr< KDTree >	_kd_tree
	The KDTree is constructed using the projected centroids of the elements in the boundary mesh.
Real	_max_projected_diag_length
	The maximum diagonal length of the projected bounding box from the boundary elements.
const std::vector< std::unique_ptr< SurfaceElement > > &	_bd_elements
	pass into the constructor for the surface elements
const std::vector< Point > &	_centroids
	pass into the constructor for the surface element centroids
int	_dim = -1
	The dimension of the embedding mesh.
std::size_t	_num_elements = 0
	The number of elements in the boundary mesh.
Point	_ray_direction
	Ray shooting direction.
Real	_eps_on_surface
	Epsilon value for checking if a point is on the surface of the geometry.
int	_leaf_max_size
	Configures KDTree leaf node size for performance tuning.
FileName	_obb_file_name
	The file name for the OBB.
FileName	_ray_file_name
	The file name for the ray.
const libMesh::Parallel::Communicator *	_comm = nullptr
	Communicator used only for writing the debug OBB/ray mesh files.
OrientedBoundingBox	_obb_bounds
	The oriented bounding box (OBB).
BoundingBox	_bounds
	The bounding box AABB.
bool	_auto_ray_direction = false
	Whether the ray direction is auto-selected via PCA (true) or user-selected (false).
bool	_build_obb = false
	When the ray direction is auto-selected (PCA) we build an Oriented Bounding Box (OBB); a user-selected ray uses only the global Axis-Aligned Bounding Box (AABB) instead.
Point	_plane_origin = Point(0.0, 0.0, 0.0)
	The origin of the plane used to ensure that every projected point is correctly aligned and lies on the same plane.
Point	_max_variance_vector
	max variance vector
Point	_second_variance_vector
	second max variance vector
Point	_min_variance_vector
	min variance vector (only used for 3D)

Detailed Description

Ray-casting point-in-solid engine over a closed surface mesh represented as a collection of SurfaceElement wrappers.

The shooting direction follows one of two explicit policies (RayDirectionMode):

• AUTO_PCA: the engine auto-selects a robust direction via PCA, builds an oriented bounding box (OBB), and may use a parity-tie fallback along other variance axes.
• USER_SPECIFIED: the caller's direction is used exactly. There is no PCA, no OBB, and no fallback; the engine uses a global axis-aligned bounding box (AABB) and, on a genuinely ambiguous (grazing/tangent) query, errors rather than silently changing direction. The direction is validated only to be finite, non-zero, and in-plane for a 2D surface.

Definition at line 36 of file AdaptiveRayContainmentCheck.h.

Constructor & Destructor Documentation

AdaptiveRayContainmentCheck()

AdaptiveRayContainmentCheck::AdaptiveRayContainmentCheck	(	const std::vector< std::unique_ptr< SurfaceElement > > &	bd_elements,
		const std::vector< Point > &	centroids,
		const SurfaceGeometry::RayDirectionOptions &	ray_options,
		const Real	eps_on_surface = libMesh::TOLERANCE,
		const int	leaf_max_size = 10,
		const FileName &	obb_file_name = "",
		const FileName &	ray_file_name = "",
		const libMesh::Parallel::Communicator *	comm = nullptr
	)

Definition at line 19 of file AdaptiveRayContainmentCheck.C.

  : _bd_elements(bd_elements),
    _centroids(centroids),
    _ray_direction(ray_options.direction),
    _eps_on_surface(eps_on_surface),
    _leaf_max_size(leaf_max_size),
    _obb_file_name(obb_file_name),
    _ray_file_name(ray_file_name),
    _comm(comm),
    _plane_origin(Point(0.0, 0.0, 0.0))
{
  // Runtime (not mooseAssert) check: this precondition depends on caller-supplied input and
  // guards the _bd_elements[0] dereference below. mooseAssert is compiled out in opt builds,
  // which would leave that dereference as out-of-bounds undefined behavior on an empty set.
  if (_bd_elements.empty())
    mooseError(
        "AdaptiveRayContainmentCheck: boundary elements must not be empty or uninitialized.");
  _num_elements = _bd_elements.size();
  _dim = _bd_elements[0]->expectedEmbeddingMeshDim();
 
  _auto_ray_direction = (ray_options.mode == SurfaceGeometry::RayDirectionMode::AUTO_PCA);
 
  if (!_auto_ray_direction)
  {
    // USER_SPECIFIED policy: use the direction exactly. Validate only that it is usable, then
    // normalize. Degeneracy (a ray grazing a vertex/edge or tangent to the surface) is the
    // user's responsibility; the engine never auto-corrects or switches the direction.
    for (const auto i : make_range(3u))
      if (!std::isfinite(_ray_direction(i)))
        mooseError(
            "AdaptiveRayContainmentCheck: a user-selected ray_direction must be finite; got ",
            _ray_direction,
            ".");
    if (MooseUtils::absoluteFuzzyEqual(_ray_direction.norm(), 0.0))
      mooseError("AdaptiveRayContainmentCheck: a user-selected ray_direction must be non-zero.");
    if (_dim == 2 && !MooseUtils::absoluteFuzzyEqual(_ray_direction(2), 0.0))
      mooseError("AdaptiveRayContainmentCheck: a user-selected ray_direction for a 2D surface must "
                 "lie in the mesh plane (its z component must be zero); got ",
                 _ray_direction,
                 ".");
    _ray_direction = _ray_direction.unit();
  }
 
  // PCA is only used to auto-select the direction; a user-selected ray never runs it.
  if (_auto_ray_direction)
    preparePCASVD();
  buildObbKdtreeAndMaxProjectedDiagonal(
      1e-2 /*safe protect: expanded box length in each direction and both sides*/);
}

Member Function Documentation

buildObbKdtreeAndMaxProjectedDiagonal()

void AdaptiveRayContainmentCheck::buildObbKdtreeAndMaxProjectedDiagonal ( const Real  expand_box_length )

private

Constructs an oriented bounding box (OBB) using the results of PCA and the KD-tree.

During this process, it also finds the maximum projected diagonal length.

Definition at line 506 of file AdaptiveRayContainmentCheck.C.

{
  if (!_centroids.empty())
    mooseAssert(_centroids.size() >= 3, "Need at least three points.");
 
  // (a) Prepare KD-tree data (optional) and track PCA-space extents
 
  // Initialize the ray direction if not set
  initializeRayDirection();
 
  // Global min / max along the three PCA axes
  Real u_min = std::numeric_limits<Real>::max();
  Real u_max = std::numeric_limits<Real>::lowest();
  Real v_min = u_min, v_max = u_max;
  Real w_min = u_min, w_max = u_max;
 
  _projected_centroids.resize(_num_elements);
 
  _max_projected_diag_length = 0.0;
 
  for (const auto i : make_range(_num_elements))
  {
    // Per-element KD-tree data
    {
      const Point & pt =
          (!_centroids.empty() ? _centroids[i] : _bd_elements[i]->elem().vertex_average());
 
      _projected_centroids[i] = projectPointOntoPlane(pt, _plane_origin, _ray_direction);
 
      if (_dim == 2) // flatten Z in 2-D mode
        _projected_centroids[i](2) = 0.0;
 
      _max_projected_diag_length =
          std::max(_max_projected_diag_length,
                   _bd_elements[i]->getProjectedBoundingBoxDiagonal(_ray_direction));
    }
 
    if (_build_obb)
    {
      // Update PCA-space bounding box
      const Elem & e = _bd_elements[i]->elem();
 
      for (const auto j : make_range(e.n_nodes()))
      {
        const Point d = *(e.node_ptr(j)) - _centroid_nodal_points;
        const Real u = d * _max_variance_vector;
        const Real v = d * _second_variance_vector;
        const Real w = (_dim == 3) ? d * _min_variance_vector : 0.0;
 
        u_min = std::min(u_min, u);
        u_max = std::max(u_max, u);
        v_min = std::min(v_min, v);
        v_max = std::max(v_max, v);
        if (_dim == 3)
        {
          w_min = std::min(w_min, w);
          w_max = std::max(w_max, w);
        }
      }
    }
  }
 
  if (_build_obb)
  {
    // (b) Build the oriented bounding box (OBB)
    const Point min_corner =
        _centroid_nodal_points + u_min * _max_variance_vector + v_min * _second_variance_vector +
        ((_dim == 3) ? w_min * _min_variance_vector : Point()) -
        expand_box_length * _max_variance_vector - expand_box_length * _second_variance_vector -
        ((_dim == 3) ? expand_box_length * _min_variance_vector : Point());
 
    std::vector<std::pair<Point, Point>> axis_pairs{
        {min_corner,
         min_corner + (u_max - u_min) * _max_variance_vector +
             2 * expand_box_length *
                 _max_variance_vector /*2 because we subtract 1 in min_corner*/},
 
        {min_corner,
         min_corner + (v_max - v_min) * _second_variance_vector +
             2 * expand_box_length * _second_variance_vector}};
 
    if (_dim == 3)
      axis_pairs.emplace_back(min_corner,
                              min_corner + (w_max - w_min) * _min_variance_vector +
                                  2 * expand_box_length * _min_variance_vector);
 
    _obb_bounds = OrientedBoundingBox(axis_pairs);
 
    if (_obb_file_name != "")
    {
      if (!_comm)
        mooseError("A communicator is required to write the OBB mesh file '", _obb_file_name, "'.");
      std::filesystem::path obb_path(_obb_file_name.c_str());
      _obb_bounds.writeMesh(obb_path, *_comm);
    }
 
    if (_ray_file_name != "")
    {
      if (!_comm)
        mooseError("A communicator is required to write the ray mesh file '", _ray_file_name, "'.");
      std::filesystem::path ray_path(_ray_file_name.c_str());
      _obb_bounds.writeRayAlongShortestAxis(ray_path, *_comm);
    }
  }
 
  // (c) Finalise KD-tree
  _kd_tree = std::make_unique<KDTree>(_projected_centroids, _leaf_max_size);
}

Referenced by AdaptiveRayContainmentCheck().

collectCandidateElementIDs()

std::vector< unsigned int > AdaptiveRayContainmentCheck::collectCandidateElementIDs ( const Point &  query_point ) const

private

Use the kd-tree to collect candidate element IDs to check intersections.

Definition at line 616 of file AdaptiveRayContainmentCheck.C.

{
  std::vector<unsigned int> elem_ids;
 
  // KD-tree radius search in projected PCA space
  Point proj = projectPointOntoPlane(query_point, _plane_origin, _ray_direction);
  if (_dim == 2)
    proj(2) = 0.0; // flatten Z for 2-D
 
  std::vector<nanoflann::ResultItem<std::size_t, Real>> matches;
  _kd_tree->radiusSearch(proj, _max_projected_diag_length, matches);
 
  elem_ids.reserve(matches.size());
  for (const auto & m : matches)
    elem_ids.push_back(static_cast<unsigned int>(m.first));
 
  return elem_ids;
}

Referenced by countFilteredCrossings(), and isOnSurface().

computeGlobalBoundingBox()

BoundingBox AdaptiveRayContainmentCheck::computeGlobalBoundingBox ( )

private

Compute the global bounding box of all boundary elements.

Definition at line 310 of file AdaptiveRayContainmentCheck.C.

{
  const auto & first_elem = _bd_elements[0]->elem();
  BoundingBox bbox = first_elem.loose_bounding_box();
 
  for (const auto & bd_elem : _bd_elements)
    bbox.union_with(bd_elem->elem().loose_bounding_box());
 
  const Real eps = _eps_on_surface;
  Point min_pt = bbox.min();
  Point max_pt = bbox.max();
 
  for (const auto d : make_range(3u))
  {
    min_pt(d) -= eps;
    max_pt(d) += eps;
  }
 
  return BoundingBox(min_pt, max_pt);
}

Referenced by initializeRayDirection().

countCrossings()

int AdaptiveRayContainmentCheck::countCrossings ( const Point &  ray_start, const Point &  ray_end, const bool  use_primary_direction = true  ) const

private

Count how many times the segment from ray_start to ray_end crosses the surface.

The primary direction uses KD-tree candidates; PCA fallback directions scan all elements because the KD-tree projection is aligned only with the primary direction. Callers must rule out the on-surface case (via isOnSurface) before interpreting the count.

Definition at line 196 of file AdaptiveRayContainmentCheck.C.

{
  // The 2D count uses a half-open side-based crossing rule, which counts a ray passing exactly
  // through a boundary vertex or along a collinear edge correctly without any tolerance. The 3D
  // path is unchanged.
  if (_dim == 2)
    return countCrossings2D(ray_start, ray_end, use_primary_direction);
 
  const auto ray_hits_surface = [this, &ray_start, &ray_end](const SurfaceElement * surface)
  { return rayIntersectGeometry(ray_start, ray_end, surface); };
 
  return countFilteredCrossings(ray_start, ray_end, use_primary_direction, ray_hits_surface);
}

Referenced by sideness(), and sidenessFromRayPair().

countCrossings2D()

int AdaptiveRayContainmentCheck::countCrossings2D ( const Point &  ray_start, const Point &  ray_end, const bool  use_primary_direction  ) const

private

2D crossing count using a half-open side-based crossing rule: an edge is counted when its two endpoints lie on strictly opposite sides of the ray line (a vertex or collinear edge exactly on the line is thus counted consistently, with no tolerance) and the crossing is on the ray_start side of the query point.

See countCrossings() for the parameters.

Definition at line 213 of file AdaptiveRayContainmentCheck.C.

{
  // The ray goes from ray_start (placed outside the geometry) to the query point p = ray_end, so
  // the crossings on the segment are exactly those of the half-line from p toward ray_start.
  const Point & p = ray_end;
  const Point dir = ray_end - ray_start; // start -> p
 
  const auto edge_crosses_ray = [&p, &dir](const SurfaceElement * surface)
  {
    const Elem & e = surface->elem();
    const Point a = e.point(0) - p;
    const Point b = e.point(1) - p;
 
    // Signed perpendicular position of each edge endpoint relative to the ray line. The strict ">"
    // on both endpoints is the half-open ("one end closed, one end open") crossing convention: an
    // endpoint exactly on the ray line is assigned one fixed side, so a vertex shared by two edges
    // is counted by exactly one of them, an edge collinear with the ray (both sides zero) is
    // skipped, and a tangential touch cancels to an even count. No tolerance is needed.
    const Real side_a = dir.cross(a)(2);
    const Real side_b = dir.cross(b)(2);
    if ((side_a > 0.0) == (side_b > 0.0))
      return false;
 
    // The edge crosses the ray line; keep it only if the crossing lies on the ray_start side of p
    // (on the segment [ray_start, p]), i.e. opposite the +dir direction that points from start to
    // p.
    const Real t = side_a / (side_a - side_b);
    const Point crossing = a + t * (b - a); // crossing point relative to p
    return dir * crossing < 0.0;
  };
 
  return countFilteredCrossings(ray_start, ray_end, use_primary_direction, edge_crosses_ray);
}

Referenced by countCrossings().

countFilteredCrossings()

template<typename CrossingTest >

int AdaptiveRayContainmentCheck::countFilteredCrossings ( const Point &  ray_start, const Point &  ray_end, const bool  use_primary_direction, CrossingTest  is_crossing  ) const

private

Shared traversal for the 2D and 3D crossing counts: walk the candidate elements (KD-tree candidates for the primary direction, all elements otherwise), apply the ray-bbox and bounding-region pruning, and count those for which is_crossing(surface) is true.

The per-dimension crossing test is supplied by the caller so the pruning loop lives in exactly one place. Defined in the .C: only the two in-file callers instantiate it.

Definition at line 167 of file AdaptiveRayContainmentCheck.C.

{
  const auto candidate_ids =
      use_primary_direction ? collectCandidateElementIDs(ray_end) : std::vector<unsigned int>{};
  const auto num_candidates = use_primary_direction ? candidate_ids.size() : _num_elements;
  const Point segment_direction = ray_end - ray_start;
 
  int count = 0;
  for (const auto candidate : make_range(num_candidates))
  {
    const auto elem_id = use_primary_direction ? candidate_ids[candidate] : candidate;
    const auto & surface = _bd_elements[elem_id].get();
    const auto ball = surface->computeBoundingBall();
 
    if (isOutsideRayBBox(ray_start, segment_direction, ball))
      continue;
    if (isOutsideBoundingRegion(ray_start, segment_direction, ball))
      continue;
 
    if (is_crossing(surface))
      ++count;
  }
  return count;
}

Referenced by countCrossings(), and countCrossings2D().

initializeRayDirection()

void AdaptiveRayContainmentCheck::initializeRayDirection ( )

private

Finalizes the ray direction and the matching bounding box.

Finalize the ray direction (auto -> PCA, user -> as given) and its bounding box.

For an auto ray it sets _ray_direction to the PCA-selected direction and requests the OBB; for a user-selected ray it retains the user's direction and builds the global AABB.

Definition at line 376 of file AdaptiveRayContainmentCheck.C.

{
  if (_auto_ray_direction)
  {
    // Auto ray: adopt the PCA-selected direction and use the oriented bounding box.
    _ray_direction = (_dim == 3) ? _min_variance_vector : _second_variance_vector;
    _build_obb = true;
  }
  else
    // User-selected axis-aligned ray: keep the user's (normalized) direction and use a
    // global axis-aligned bounding box.
    _bounds = computeGlobalBoundingBox();
}

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

isOnSurface()

bool AdaptiveRayContainmentCheck::isOnSurface ( const Point &  p ) const

private

True if p lies on the surface (within _eps_on_surface), i.e.

some candidate element contains it. This is a property of p alone, independent of any ray direction.

Definition at line 157 of file AdaptiveRayContainmentCheck.C.

{
  for (const auto elem_id : collectCandidateElementIDs(p))
    if (_bd_elements[elem_id].get()->elem().contains_point(p, _eps_on_surface))
      return true;
  return false;
}

Referenced by sideness().

isOutsideBoundingBox()

bool AdaptiveRayContainmentCheck::isOutsideBoundingBox ( const Point &  query_point ) const

private

Check if point is outside global bounding box.

Definition at line 259 of file AdaptiveRayContainmentCheck.C.

{
  return (_build_obb) ? !_obb_bounds.contains(query_point, _eps_on_surface)
                      : !_bounds.contains_point(query_point);
}

Referenced by sideness().

isOutsideBoundingRegion()

bool AdaptiveRayContainmentCheck::isOutsideBoundingRegion ( const Point &  orig, const Point &  dir, const Ball &  ball  ) const

private

Check if element center is outside ray bounding circle/sphere.

Definition at line 290 of file AdaptiveRayContainmentCheck.C.

{
  const auto & center = ball.center();
  const auto radius = ball.radius();
 
  const auto w = center - orig;
 
  Real b = (w * dir) / (dir * dir);
  Point Pb = orig + b * dir;
 
  Real distance_squared = 0.0;
  for (const auto i : make_range(_dim))
    distance_squared += Utility::pow<2>(Pb(i) - center(i));
 
  return (distance_squared > radius * radius);
}

Referenced by countFilteredCrossings().

isOutsideRayBBox()

bool AdaptiveRayContainmentCheck::isOutsideRayBBox ( const Point &  orig, const Point &  dir, const Ball &  ball  ) const

private

Check if element center is outside ray bounding box.

Definition at line 266 of file AdaptiveRayContainmentCheck.C.

{
  Point lb, ub;
  const auto & center = ball.center();
  const Real radius = ball.radius();
 
  for (const auto i : make_range(_dim))
  {
    lb(i) = std::min(orig(i), orig(i) + dir(i)) - radius;
    ub(i) = std::max(orig(i), orig(i) + dir(i)) + radius;
  }
 
  for (const auto i : make_range(_dim))
  {
    if (center(i) < lb(i) || center(i) > ub(i))
      return true;
  }
 
  return false;
}

Referenced by countFilteredCrossings().

preparePCASVD()

void AdaptiveRayContainmentCheck::preparePCASVD ( )

private

Definition at line 430 of file AdaptiveRayContainmentCheck.C.

{
  Point centroid_sum;
 
  std::vector<Point> nodal_points;
  for (const auto & elem : _bd_elements)
  {
    const auto & e = elem->elem();
    for (const auto i : make_range(e.n_nodes()))
    {
      const Node * node = e.node_ptr(i);
      mooseAssert(node, "Node pointer is null!");
      nodal_points.push_back(*node);
      centroid_sum += *node;
    }
  }
 
  const unsigned int N = nodal_points.size();
  mooseAssert(N >= 3, "At least 3 points required");
 
  // (a) Compute the centroid
  _centroid_nodal_points = centroid_sum / static_cast<Real>(N);
 
  // (b) Build the mean-centered matrix X (N x 3)
  DenseMatrix<Real> X(N, 3);
  for (const auto i : make_range(N))
  {
    const Point d = nodal_points[i] - _centroid_nodal_points;
    X(i, 0) = d(0);
    X(i, 1) = d(1);
    X(i, 2) = d(2);
  }
 
  // (c) Perform SVD: X = U * sigma * V^T
  DenseVector<Real> sigma;
  DenseMatrix<Real> U, VT;
  X.svd(sigma, U, VT); // VT is 3x3, each row is a principal direction
 
  // (d) Extract principal directions
  _max_variance_vector = Point(VT(0, 0), VT(0, 1), VT(0, 2));    // max variance
  _second_variance_vector = Point(VT(1, 0), VT(1, 1), VT(1, 2)); // second largest variance
  _min_variance_vector = Point(VT(2, 0), VT(2, 1), VT(2, 2));    // min variance
 
  // (e) normalize them to be safe (unit() returns a normalized copy)
  _max_variance_vector = _max_variance_vector.unit();
  _second_variance_vector = _second_variance_vector.unit();
  _min_variance_vector = _min_variance_vector.unit();
 
  // (f) Canonicalize the sign of each principal direction. SVD singular vectors
  // are only defined up to sign, and LAPACK can return opposite signs on
  // different platforms or versions. Because the auto-selected ray direction is
  // one of these vectors, an unstable sign makes the in-out classification of
  // borderline elements non-reproducible across platforms. Fix a deterministic
  // convention: make the largest-magnitude component positive (ties broken by
  // the lowest index).
  auto canonicalize_sign = [](Point & v)
  {
    unsigned int i_max = 0;
    for (const auto i : make_range(1, 3))
      if (std::abs(v(i)) > std::abs(v(i_max)))
        i_max = i;
    if (v(i_max) < 0.0)
      v *= -1.0;
  };
  canonicalize_sign(_max_variance_vector);
  canonicalize_sign(_second_variance_vector);
  canonicalize_sign(_min_variance_vector);
 
  mooseAssert(
      MooseUtils::absoluteFuzzyEqual(_max_variance_vector * _second_variance_vector, 0.0) &&
          MooseUtils::absoluteFuzzyEqual(_max_variance_vector * _min_variance_vector, 0.0) &&
          MooseUtils::absoluteFuzzyEqual(_second_variance_vector * _min_variance_vector, 0.0),
      "Principal directions are not orthogonal.");
}

Referenced by AdaptiveRayContainmentCheck().

projectPointOntoPlane()

Point AdaptiveRayContainmentCheck::projectPointOntoPlane ( const Point &  point_to_project, const Point &  plane_point, const Point &  plane_normal  ) const

private

Orthogonally project point_to_project onto the plane defined by plane_point and unit normal plane_normal.

plane_normal is assumed to be a unit vector.

Definition at line 419 of file AdaptiveRayContainmentCheck.C.

{
  // Delegate to libMesh::Plane::closest_point, which returns the orthogonal
  // projection of the point onto the plane. `plane_normal` is assumed to be a
  // unit vector (closest_point does not normalize it).
  return libMesh::Plane(plane_point, plane_normal).closest_point(point_to_project);
}

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), collectCandidateElementIDs(), and rayStartOutsideOBB().

rayDirection()

const Point & AdaptiveRayContainmentCheck::rayDirection ( ) const

inline

The resolved ray direction actually used for shooting: the (normalized) user direction for a user-selected ray, or the PCA-selected direction for an auto ray.

Definition at line 53 of file AdaptiveRayContainmentCheck.h.

{ return _ray_direction; }

rayIntersectGeometry()

bool AdaptiveRayContainmentCheck::rayIntersectGeometry ( const Point &  ray_start, const Point &  ray_end, const SurfaceElement *  elem  ) const

private

Ray-element intersection (e.g., ray-line for 2D, ray-triangle for 3D)

Definition at line 250 of file AdaptiveRayContainmentCheck.C.

{
  LineSegment ray_segment(ray_start, ray_end);
  return elem->intersect(ray_segment);
}

Referenced by countCrossings().

rayStartOutsideAABB()

Point AdaptiveRayContainmentCheck::rayStartOutsideAABB ( const Point &  point, const Point &  unit_direction, const bool  inverted  ) const

private

Ray start strictly outside the global AABB along unit_direction, for any direction.

The 8 AABB corners are projected onto the direction to find the box extent, and the start is placed just past the far side (a scale-aware padding), moving only the distance needed. inverted shoots from the opposite side. unit_direction must be normalized.

Definition at line 391 of file AdaptiveRayContainmentCheck.C.

{
  // Project the 8 AABB corners onto the direction to find the box's extent along it. The ray
  // start is then placed just past the far side, so it is provably outside the box while moving
  // only the distance needed (a full-diagonal displacement would make unnecessarily long rays).
  const Point & lo = _bounds.min();
  const Point & hi = _bounds.max();
 
  Real min_projection = std::numeric_limits<Real>::max();
  Real max_projection = std::numeric_limits<Real>::lowest();
  for (const auto c : make_range(8u))
  {
    const Point corner(
        (c & 1u) ? hi(0) : lo(0), (c & 2u) ? hi(1) : lo(1), (c & 4u) ? hi(2) : lo(2));
    const Real projection = corner * unit_direction;
    min_projection = std::min(min_projection, projection);
    max_projection = std::max(max_projection, projection);
  }
 
  // Scale-aware padding so the start never lands exactly on the box boundary.
  const Real padding = _eps_on_surface + 1e-2 * (max_projection - min_projection);
  const Real target = inverted ? max_projection + padding : min_projection - padding;
  return point + (target - point * unit_direction) * unit_direction;
}

Referenced by sideness().

rayStartOutsideOBB()

Point AdaptiveRayContainmentCheck::rayStartOutsideOBB ( const Point &  point, const Point &  ray_direction, const unsigned int  obb_axis, const bool  inverted = false  ) const

private

Computes the starting point of an OBB-based ray (auto/PCA policy) for a given query point.

The point is projected onto the OBB face normal to ray_direction, selected by obb_axis, and displaced along that same direction so the ray originates outside the geometry.

Definition at line 332 of file AdaptiveRayContainmentCheck.C.

{
  mooseAssert(_build_obb,
              "AdaptiveRayContainmentCheck::rayStartOutsideOBB: OBB-based ray start is only used "
              "by the auto (PCA) policy.");
  mooseAssert(obb_axis < static_cast<unsigned int>(_dim),
              "AdaptiveRayContainmentCheck::rayStartOutsideOBB: invalid OBB axis index.");
 
  const Real axis_length = _obb_bounds.getAxisLength(obb_axis);
  const Real half_axis_length = axis_length / 2.0;
 
  // The projection below lands the point exactly on the OBB face perpendicular to the ray
  // (ray_direction == the OBB axis direction), so any positive outward step is provably
  // outside the box. Use a scale-aware padding, consistent with rayStartOutsideAABB: an
  // absolute floor (also guards a near-zero axis) plus a small fraction of the box extent
  // to stay clear of floating-point noise at large coordinate scales.
  const Real padding = _eps_on_surface + 1e-2 * axis_length;
 
  Point projection_plane_corner;
  Real direction_multiplier;
 
  if (_obb_bounds.getProjectedLength(point, obb_axis) < half_axis_length)
  {
    projection_plane_corner =
        inverted ? _obb_bounds.getMaximalCorner() : _obb_bounds.getMinimalCorner();
    direction_multiplier = inverted ? -1.0 : 1.0;
  }
  else
  {
    projection_plane_corner =
        inverted ? _obb_bounds.getMinimalCorner() : _obb_bounds.getMaximalCorner();
    direction_multiplier = inverted ? 1.0 : -1.0;
  }
 
  const Point projected_point =
      projectPointOntoPlane(point, projection_plane_corner, ray_direction);
  return projected_point - padding * direction_multiplier * ray_direction;
}

Referenced by sideness().

sideness()

SurfaceGeometry::SurfaceSide AdaptiveRayContainmentCheck::sideness

(

const Point &

p

)

const

Main function: Determine if a point is inside the geometry.

Definition at line 78 of file AdaptiveRayContainmentCheck.C.

{
  if (isOutsideBoundingBox(p))
    return SurfaceGeometry::SurfaceSide::OUTSIDE;
 
  // Whether p sits on the surface is a property of p alone, so decide it once here rather than
  // re-testing it on every ray cast below.
  if (isOnSurface(p))
    return SurfaceGeometry::SurfaceSide::ON;
 
  const std::array<Point, 2> ray_starts =
      _auto_ray_direction
          ? std::array<Point, 2>{rayStartOutsideOBB(p, _ray_direction, _dim - 1, false),
                                 rayStartOutsideOBB(p, _ray_direction, _dim - 1, true)}
          : std::array<Point, 2>{rayStartOutsideAABB(p, _ray_direction, false),
                                 rayStartOutsideAABB(p, _ray_direction, true)};
 
  if (const auto side = sidenessFromRayPair(p, ray_starts))
    return *side;
 
  if (!_auto_ray_direction)
  {
    std::ostringstream oss;
    oss << p;
    mooseError("AdaptiveRayContainmentCheck: the user-selected ray_direction ",
               _ray_direction,
               " gives an ambiguous (grazing or tangent) intersection at point ",
               oss.str(),
               "; choose a different ray_direction or use the auto (pca_ray) method.");
  }
 
  // Defensive fallback. With the half-open crossing count, the two opposite primary rays always
  // agree in parity for a well-formed closed surface, so this path is not reached for valid input;
  // it is kept as a safety net should a future or degenerate case ever produce a parity
  // disagreement. Probe the remaining PCA variance directions: if any probe ray escapes without
  // crossing the surface, the point is outside. Otherwise the query is undecidable.
  for (const auto obb_axis : make_range(static_cast<unsigned int>(_dim - 1)))
  {
    const Point fallback_direction = _obb_bounds.getAxisDirection(obb_axis);
    const std::array<Point, 2> probe_starts = {
        rayStartOutsideOBB(p, fallback_direction, obb_axis, false),
        rayStartOutsideOBB(p, fallback_direction, obb_axis, true)};
 
    for (const auto & probe_start : probe_starts)
      if (countCrossings(probe_start, p, false) == 0)
        return SurfaceGeometry::SurfaceSide::OUTSIDE;
  }
 
  std::ostringstream oss;
  oss << p;
  mooseError("AdaptiveRayContainmentCheck: No decision could be made for point " + oss.str());
}

sidenessFromRayPair()

std::optional< SurfaceGeometry::SurfaceSide > AdaptiveRayContainmentCheck::sidenessFromRayPair ( const Point &  p, const std::array< Point, 2 > &  ray_starts  ) const

private

Determine sideness from a pair of opposite rays.

A conflicting parity returns std::nullopt so the caller can apply the selected ambiguity policy.

Definition at line 132 of file AdaptiveRayContainmentCheck.C.

{
  std::array<int, 2> counts = {0, 0};
 
  // Shoot the two (opposite) rays and count intersections with the boundary elements. p has
  // already been ruled on-surface by sideness(), so the count is always meaningful here.
  for (const auto i : make_range(2))
  {
    counts[i] = countCrossings(ray_starts[i], p);
 
    // A ray that never crosses the closed surface proves the point is outside.
    if (counts[i] == 0)
      return SurfaceGeometry::SurfaceSide::OUTSIDE;
  }
 
  // Consistent parity gives a definite decision; conflicting parity is undecided (nullopt) and
  // left to the caller's policy.
  if ((counts[0] % 2) == (counts[1] % 2))
    return (counts[0] % 2 == 1) ? SurfaceGeometry::SurfaceSide::INSIDE
                                : SurfaceGeometry::SurfaceSide::OUTSIDE;
  return std::nullopt;
}

Referenced by sideness().

Member Data Documentation

_auto_ray_direction

bool AdaptiveRayContainmentCheck::_auto_ray_direction = false

private

Whether the ray direction is auto-selected via PCA (true) or user-selected (false).

Set once in the constructor from the SurfaceGeometry::RayDirectionOptions mode. Only an auto ray runs PCA and may have its direction chosen automatically.

Definition at line 112 of file AdaptiveRayContainmentCheck.h.

Referenced by AdaptiveRayContainmentCheck(), initializeRayDirection(), and sideness().

_bd_elements

const std::vector<std::unique_ptr<SurfaceElement> >& AdaptiveRayContainmentCheck::_bd_elements

private

pass into the constructor for the surface elements

Definition at line 73 of file AdaptiveRayContainmentCheck.h.

Referenced by AdaptiveRayContainmentCheck(), buildObbKdtreeAndMaxProjectedDiagonal(), computeGlobalBoundingBox(), countFilteredCrossings(), isOnSurface(), and preparePCASVD().

_bounds

BoundingBox AdaptiveRayContainmentCheck::_bounds

private

The bounding box AABB.

Definition at line 107 of file AdaptiveRayContainmentCheck.h.

Referenced by initializeRayDirection(), isOutsideBoundingBox(), and rayStartOutsideAABB().

_build_obb

bool AdaptiveRayContainmentCheck::_build_obb = false

private

When the ray direction is auto-selected (PCA) we build an Oriented Bounding Box (OBB); a user-selected ray uses only the global Axis-Aligned Bounding Box (AABB) instead.

Definition at line 116 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), initializeRayDirection(), isOutsideBoundingBox(), and rayStartOutsideOBB().

_centroid_nodal_points

Point AdaptiveRayContainmentCheck::_centroid_nodal_points

private

The centroid of the boundary elements' node points (prepare inside this class).

Definition at line 58 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), and preparePCASVD().

_centroids

const std::vector<Point>& AdaptiveRayContainmentCheck::_centroids

private

pass into the constructor for the surface element centroids

Definition at line 76 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_comm

const libMesh::Parallel::Communicator* AdaptiveRayContainmentCheck::_comm = nullptr

private

Communicator used only for writing the debug OBB/ray mesh files.

Null when no debug output is requested.

Definition at line 101 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_dim

int AdaptiveRayContainmentCheck::_dim = -1

private

The dimension of the embedding mesh.

Definition at line 79 of file AdaptiveRayContainmentCheck.h.

Referenced by AdaptiveRayContainmentCheck(), buildObbKdtreeAndMaxProjectedDiagonal(), collectCandidateElementIDs(), countCrossings(), initializeRayDirection(), isOutsideBoundingRegion(), isOutsideRayBBox(), rayStartOutsideOBB(), and sideness().

_eps_on_surface

Real AdaptiveRayContainmentCheck::_eps_on_surface

private

Epsilon value for checking if a point is on the surface of the geometry.

Definition at line 88 of file AdaptiveRayContainmentCheck.h.

Referenced by computeGlobalBoundingBox(), isOnSurface(), isOutsideBoundingBox(), rayStartOutsideAABB(), and rayStartOutsideOBB().

_kd_tree

std::unique_ptr<KDTree> AdaptiveRayContainmentCheck::_kd_tree

private

The KDTree is constructed using the projected centroids of the elements in the boundary mesh.

(prepare inside this class)

Definition at line 67 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), and collectCandidateElementIDs().

_leaf_max_size

int AdaptiveRayContainmentCheck::_leaf_max_size

private

Configures KDTree leaf node size for performance tuning.

Definition at line 91 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_max_projected_diag_length

Real AdaptiveRayContainmentCheck::_max_projected_diag_length

private

The maximum diagonal length of the projected bounding box from the boundary elements.

Definition at line 70 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), and collectCandidateElementIDs().

_max_variance_vector

Point AdaptiveRayContainmentCheck::_max_variance_vector

private

max variance vector

Definition at line 123 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), and preparePCASVD().

_min_variance_vector

Point AdaptiveRayContainmentCheck::_min_variance_vector

private

min variance vector (only used for 3D)

Definition at line 129 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), initializeRayDirection(), and preparePCASVD().

_num_elements

std::size_t AdaptiveRayContainmentCheck::_num_elements = 0

private

The number of elements in the boundary mesh.

Definition at line 82 of file AdaptiveRayContainmentCheck.h.

Referenced by AdaptiveRayContainmentCheck(), buildObbKdtreeAndMaxProjectedDiagonal(), and countFilteredCrossings().

_obb_bounds

OrientedBoundingBox AdaptiveRayContainmentCheck::_obb_bounds

private

The oriented bounding box (OBB).

Definition at line 104 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), isOutsideBoundingBox(), rayStartOutsideOBB(), and sideness().

_obb_file_name

FileName AdaptiveRayContainmentCheck::_obb_file_name

private

The file name for the OBB.

Definition at line 94 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_plane_origin

Point AdaptiveRayContainmentCheck::_plane_origin = Point(0.0, 0.0, 0.0)

private

The origin of the plane used to ensure that every projected point is correctly aligned and lies on the same plane.

This is for fast k-d tree searching to select the candidate elements to check the intersection.

Definition at line 120 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), and collectCandidateElementIDs().

_projected_centroids

std::vector<Point> AdaptiveRayContainmentCheck::_projected_centroids

private

Projected centroids of the elements in the boundary mesh.

(prepare inside this class) This is used to construct the KDTree in this struct for fast searching.

Definition at line 63 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_ray_direction

Point AdaptiveRayContainmentCheck::_ray_direction

private

Ray shooting direction.

Definition at line 85 of file AdaptiveRayContainmentCheck.h.

Referenced by AdaptiveRayContainmentCheck(), buildObbKdtreeAndMaxProjectedDiagonal(), collectCandidateElementIDs(), initializeRayDirection(), rayDirection(), and sideness().

_ray_file_name

FileName AdaptiveRayContainmentCheck::_ray_file_name

private

The file name for the ray.

Definition at line 97 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal().

_second_variance_vector

Point AdaptiveRayContainmentCheck::_second_variance_vector

private

second max variance vector

Definition at line 126 of file AdaptiveRayContainmentCheck.h.

Referenced by buildObbKdtreeAndMaxProjectedDiagonal(), initializeRayDirection(), and preparePCASVD().

---

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

• include/utils/AdaptiveRayContainmentCheck.h
• src/utils/AdaptiveRayContainmentCheck.C