https://mooseframework.inl.gov
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
TriangleManifold Class Reference

Utility for querying point containment against a closed triangulated surface mesh. More...

#include <TriangleManifold.h>

Public Member Functions

 TriangleManifold (MeshBase &mesh, const Real surface_tolerance)
 Build a manifold classifier from a prepared surface mesh. More...
 
bool contains (const Point &point) const
 
const libMesh::BoundingBoxboundingBox () const
 
std::size_t numTriangles () const
 

Private Types

enum  RayIntersection { RayIntersection::Miss, RayIntersection::Hit, RayIntersection::Ambiguous }
 Result of intersecting the positive x-direction ray with a triangle. More...
 

Private Member Functions

void finalize ()
 Complete post-parse validation and acceleration-structure setup. More...
 
void buildCandidateGrid ()
 Build the yz-plane lookup grid used to accelerate +x ray queries. More...
 
bool pointInsideBoundingBox (const Point &point) const
 Cheap global bounding-box rejection for containment queries. More...
 
bool pointOnSurface (const Point &point) const
 Detect whether a query point lies on or extremely near the manifold surface. More...
 
RayIntersection rayIntersectsTriangle (const Point &point, const libMesh::Elem &tri) const
 Intersect a positive x-direction ray with a single triangle. More...
 
bool containsBySolidAngle (const Point &point) const
 Robust fallback containment query based on accumulated solid angle. More...
 
std::vector< dof_id_typerayCandidates (const Point &point) const
 Get the subset of triangles whose yz extents may intersect the query ray. More...
 

Private Attributes

std::size_t _num_y_cells = 1
 Number of yz-grid bins in the y direction. More...
 
std::size_t _num_z_cells = 1
 Number of yz-grid bins in the z direction. More...
 
Real _y_min = 0.0
 Minimum global y coordinate used to map query points into yz-grid bins. More...
 
Real _z_min = 0.0
 Minimum global z coordinate used to map query points into yz-grid bins. More...
 
Real _y_cell_size = 1.0
 Width of one yz-grid cell in the y direction. More...
 
Real _z_cell_size = 1.0
 Width of one yz-grid cell in the z direction. More...
 
std::unordered_map< dof_id_type, std::vector< dof_id_type > > _ray_grid
 Lookup from packed yz-grid cell index to triangles that could intersect the +x query ray. More...
 
MeshBase & _mesh
 
const Real _surface_tolerance
 Absolute tolerance used throughout validation and geometric classification. More...
 
const libMesh::BoundingBox _bounding_box
 Global bounding box of the transformed manifold. More...
 
const std::unique_ptr< libMesh::PointLocatorBase_point_locator
 Pre-built point locator for fast proximity-to-surface detection. More...
 

Detailed Description

Utility for querying point containment against a closed triangulated surface mesh.

The class takes ownership of a reference to a prepared 2D surface mesh, validates that it forms a closed 2-manifold (every triangle has exactly three neighbors, all elements are Tri3, and the mesh is consistently oriented), builds a lightweight yz-plane acceleration grid, and exposes contains() for point-in-solid queries.

Containment is resolved in three stages:

  1. Cheap global bounding-box rejection.
  2. Near-surface detection via a pre-built point locator (points within surface_tolerance of the mesh surface are treated as inside).
  3. Odd/even parity counting on a fixed +x ray, with automatic fallback to a solid-angle accumulation test when the ray grazes a triangle edge or vertex.

The referenced mesh must outlive this object. Any geometric transforms (scale, rotation, translation) should be applied to the mesh before constructing a TriangleManifold.

Definition at line 44 of file TriangleManifold.h.

Member Enumeration Documentation

◆ RayIntersection

Result of intersecting the positive x-direction ray with a triangle.

Ambiguous is returned when the hit is too close to an edge, vertex, or the ray origin. In that case we abandon parity counting and fall back to the more expensive but robust solid-angle test.

Enumerator
Miss 
Hit 
Ambiguous 

Definition at line 81 of file TriangleManifold.h.

82  {
83  Miss,
84  Hit,
85  Ambiguous
86  };

Constructor & Destructor Documentation

◆ TriangleManifold()

TriangleManifold::TriangleManifold ( MeshBase &  mesh,
const Real  surface_tolerance 
)

Build a manifold classifier from a prepared surface mesh.

Parameters
meshSerialized 2D surface mesh that defines the closed manifold. Must outlive this object. Any desired transforms should already be applied to the mesh.
surface_toleranceAbsolute tolerance used for manifold validation and near-surface classification. Choose this relative to the mesh length scale and expected coordinate noise from the export pipeline.

Definition at line 56 of file TriangleManifold.C.

57  : _mesh(mesh),
58  _surface_tolerance(surface_tolerance),
59  _bounding_box(MeshTools::create_bounding_box(_mesh)),
60  _point_locator(_mesh.sub_point_locator())
61 {
62  mooseAssert(_surface_tolerance > 0.0, "surface_tolerance must be strictly positive.");
63  mooseAssert(mesh.is_serial(), "Input manifold mesh must be serialized.");
64  mooseAssert(mesh.mesh_dimension() == 2, "Manifold mesh must be a surface.");
65 
66  // Finish topology validation and acceleration-structure setup before queries are allowed.
67  finalize();
68 
69  _point_locator->set_close_to_point_tol(_surface_tolerance);
70 }
const std::unique_ptr< libMesh::PointLocatorBase > _point_locator
Pre-built point locator for fast proximity-to-surface detection.
MeshBase & mesh
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
const libMesh::BoundingBox _bounding_box
Global bounding box of the transformed manifold.
void finalize()
Complete post-parse validation and acceleration-structure setup.

Member Function Documentation

◆ boundingBox()

const libMesh::BoundingBox& TriangleManifold::boundingBox ( ) const
inline
Returns
The manifold bounding box.

Definition at line 66 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), pointInsideBoundingBox(), and rayCandidates().

66 { return _bounding_box; };
const libMesh::BoundingBox _bounding_box
Global bounding box of the transformed manifold.

◆ buildCandidateGrid()

void TriangleManifold::buildCandidateGrid ( )
private

Build the yz-plane lookup grid used to accelerate +x ray queries.

Definition at line 128 of file TriangleManifold.C.

Referenced by finalize().

129 {
130  // Since every containment ray travels in +x, only y and z are needed to choose candidate
131  // triangles for the parity test.
132  const auto extent_y =
134  const auto extent_z =
136 
137  // Choose a roughly square yz grid scaled by aspect ratio so candidate lists stay short.
138  const auto target_cells =
139  std::max<dof_id_type>(1, static_cast<dof_id_type>(std::sqrt(numTriangles())));
140  const auto aspect = std::sqrt(extent_y / extent_z);
141 
142  _num_y_cells = std::max<dof_id_type>(
143  1, static_cast<dof_id_type>(std::lround(std::sqrt(target_cells) * aspect)));
144  _num_z_cells = std::max<std::size_t>(
145  1, static_cast<std::size_t>(std::ceil(static_cast<Real>(target_cells) / _num_y_cells)));
146 
147  _y_min = boundingBox().min()(1);
148  _z_min = boundingBox().min()(2);
149  _y_cell_size = extent_y / _num_y_cells;
150  _z_cell_size = extent_z / _num_z_cells;
151 
152  for (const auto elem : _mesh.active_element_ptr_range())
153  {
154  const auto triangle_index = elem->id();
155  const auto bbox = elem->loose_bounding_box();
156 
157  const auto y_start =
159  const auto y_stop =
161  const auto z_start =
163  const auto z_stop =
165 
166  // Insert the triangle into every grid cell touched by its yz projection.
167  for (const auto iy : make_range(y_start, y_stop + 1))
168  for (const auto iz : make_range(z_start, z_stop + 1))
169  _ray_grid[TriangleManifoldUtils::packCell(iy, iz)].push_back(triangle_index);
170  }
171 }
Real _z_min
Minimum global z coordinate used to map query points into yz-grid bins.
Real _z_cell_size
Width of one yz-grid cell in the z direction.
Real _y_cell_size
Width of one yz-grid cell in the y direction.
Real _y_min
Minimum global y coordinate used to map query points into yz-grid bins.
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
auto max(const L &left, const R &right)
std::size_t _num_z_cells
Number of yz-grid bins in the z direction.
const Point & min() const
std::uint64_t packCell(const std::size_t y_index, const std::size_t z_index)
std::unordered_map< dof_id_type, std::vector< dof_id_type > > _ray_grid
Lookup from packed yz-grid cell index to triangles that could intersect the +x query ray...
std::size_t cellIndex(const Real value, const Real min_value, const Real cell_size, const std::size_t num_cells)
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
const libMesh::BoundingBox & boundingBox() const
std::size_t numTriangles() const
IntRange< T > make_range(T beg, T end)
auto min(const L &left, const R &right)
uint8_t dof_id_type
std::size_t _num_y_cells
Number of yz-grid bins in the y direction.

◆ contains()

bool TriangleManifold::contains ( const Point &  point) const
Returns
True if the point is inside the manifold or lies on the surface.

Definition at line 73 of file TriangleManifold.C.

Referenced by ManifoldSubdomainGenerator::generate().

74 {
75  // Most points are rejected here before we perform any per-triangle work.
77  return false;
78 
79  // By design, near-surface points count as inside for downstream subdomain tagging.
80  if (pointOnSurface(point))
81  return true;
82 
83  // Candidate filtering keeps the parity test from touching every triangle in large surfaces.
84  const auto candidates = rayCandidates(point);
85  unsigned int num_hits = 0;
86 
87  for (const auto triangle_index : candidates)
88  {
89  const auto tri = _mesh.elem_ptr(triangle_index);
90  // Triangles wholly behind the ray origin cannot contribute to the parity count.
91  if (tri->loose_bounding_box().max()(0) < point(0) - _surface_tolerance)
92  continue;
93 
94  switch (rayIntersectsTriangle(point, *tri))
95  {
97  break;
99  // Standard odd/even counting on a closed surface.
100  ++num_hits;
101  break;
103  // Edge and vertex grazing hits are exactly where parity counting becomes brittle.
104  return containsBySolidAngle(point);
105  }
106  }
107 
108  return num_hits % 2;
109 }
bool containsBySolidAngle(const Point &point) const
Robust fallback containment query based on accumulated solid angle.
bool pointOnSurface(const Point &point) const
Detect whether a query point lies on or extremely near the manifold surface.
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
bool pointInsideBoundingBox(const Point &point) const
Cheap global bounding-box rejection for containment queries.
std::vector< dof_id_type > rayCandidates(const Point &point) const
Get the subset of triangles whose yz extents may intersect the query ray.
RayIntersection rayIntersectsTriangle(const Point &point, const libMesh::Elem &tri) const
Intersect a positive x-direction ray with a single triangle.

◆ containsBySolidAngle()

bool TriangleManifold::containsBySolidAngle ( const Point &  point) const
private

Robust fallback containment query based on accumulated solid angle.

Definition at line 241 of file TriangleManifold.C.

Referenced by contains().

242 {
243  // For a closed oriented mesh, the total solid angle is approximately +-4\pi inside and 0
244  // outside. Taking parity over components handles nested shells cleanly.
245  Real total_angle = 0.0;
246  for (const auto elem : _mesh.active_element_ptr_range())
247  total_angle +=
248  geom_utils::solidAngle(point, elem->node_ref(0), elem->node_ref(1), elem->node_ref(2));
249 
250  return std::abs(total_angle) > 2.0 * libMesh::pi;
251 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:50
Real solidAngle(const Point &point, const Point &v0, const Point &v1, const Point &v2)
Compute the signed solid angle subtended by one oriented triangle at the query point.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Real pi

◆ finalize()

void TriangleManifold::finalize ( )
private

Complete post-parse validation and acceleration-structure setup.

Definition at line 112 of file TriangleManifold.C.

Referenced by TriangleManifold().

113 {
114  // Validate that the mesh can be used as a manifold
115  // We use the same logic as determining if the mesh can the tetrahedralized.
116  auto umesh = dynamic_cast<UnstructuredMesh *>(&_mesh);
118  auto msg = checker.improveAndValidate();
119  if (!msg.empty())
120  mooseError(
121  "The inputted surface mesh cannot be treated as manifold for the following reasons:\n",
122  msg);
123 
125 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:311
void buildCandidateGrid()
Build the yz-plane lookup grid used to accelerate +x ray queries.

◆ numTriangles()

std::size_t TriangleManifold::numTriangles ( ) const
inline
Returns
The number of triangles in the loaded manifold.

Definition at line 71 of file TriangleManifold.h.

Referenced by buildCandidateGrid().

71 { return _mesh.n_active_elem(); }

◆ pointInsideBoundingBox()

bool TriangleManifold::pointInsideBoundingBox ( const Point &  point) const
private

Cheap global bounding-box rejection for containment queries.

Definition at line 174 of file TriangleManifold.C.

Referenced by contains().

175 {
176  // Inflate the global box by the tolerance so near-surface points are not rejected too early.
177  return point(0) >= boundingBox().min()(0) - _surface_tolerance &&
178  point(0) <= boundingBox().max()(0) + _surface_tolerance &&
179  point(1) >= boundingBox().min()(1) - _surface_tolerance &&
180  point(1) <= boundingBox().max()(1) + _surface_tolerance &&
181  point(2) >= boundingBox().min()(2) - _surface_tolerance &&
182  point(2) <= boundingBox().max()(2) + _surface_tolerance;
183 }
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
const Point & min() const
const libMesh::BoundingBox & boundingBox() const
const Point & max() const

◆ pointOnSurface()

bool TriangleManifold::pointOnSurface ( const Point &  point) const
private

Detect whether a query point lies on or extremely near the manifold surface.

Definition at line 186 of file TriangleManifold.C.

Referenced by contains().

187 {
188  return (*_point_locator)(point) != nullptr;
189 }
const std::unique_ptr< libMesh::PointLocatorBase > _point_locator
Pre-built point locator for fast proximity-to-surface detection.

◆ rayCandidates()

std::vector< dof_id_type > TriangleManifold::rayCandidates ( const Point &  point) const
private

Get the subset of triangles whose yz extents may intersect the query ray.

Definition at line 254 of file TriangleManifold.C.

Referenced by contains().

255 {
256  // If the query is outside the manifold's yz extent, the fixed +x ray cannot hit anything.
257  if (point(1) < boundingBox().min()(1) - _surface_tolerance ||
258  point(1) > boundingBox().max()(1) + _surface_tolerance ||
259  point(2) < boundingBox().min()(2) - _surface_tolerance ||
260  point(2) > boundingBox().max()(2) + _surface_tolerance)
261  return {};
262 
263  const auto y_index =
265  const auto z_index =
267  const auto it = _ray_grid.find(TriangleManifoldUtils::packCell(y_index, z_index));
268  if (it == _ray_grid.end())
269  return {};
270 
271  // Return by value to keep the helper simple and avoid exposing internal storage.
272  return it->second;
273 }
Real _z_min
Minimum global z coordinate used to map query points into yz-grid bins.
Real _z_cell_size
Width of one yz-grid cell in the z direction.
Real _y_cell_size
Width of one yz-grid cell in the y direction.
Real _y_min
Minimum global y coordinate used to map query points into yz-grid bins.
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
auto max(const L &left, const R &right)
std::size_t _num_z_cells
Number of yz-grid bins in the z direction.
std::uint64_t packCell(const std::size_t y_index, const std::size_t z_index)
std::unordered_map< dof_id_type, std::vector< dof_id_type > > _ray_grid
Lookup from packed yz-grid cell index to triangles that could intersect the +x query ray...
std::size_t cellIndex(const Real value, const Real min_value, const Real cell_size, const std::size_t num_cells)
const libMesh::BoundingBox & boundingBox() const
auto min(const L &left, const R &right)
std::size_t _num_y_cells
Number of yz-grid bins in the y direction.

◆ rayIntersectsTriangle()

TriangleManifold::RayIntersection TriangleManifold::rayIntersectsTriangle ( const Point &  point,
const libMesh::Elem tri 
) const
private

Intersect a positive x-direction ray with a single triangle.

Definition at line 192 of file TriangleManifold.C.

Referenced by contains().

193 {
194  // This is a fixed-direction Moller-Trumbore style ray/triangle intersection test.
195  static const Point direction(1.0, 0.0, 0.0);
196 
197  const Point edge1 = tri.node_ref(1) - tri.node_ref(0);
198  const Point edge2 = tri.node_ref(2) - tri.node_ref(0);
199  const Point h = direction.cross(edge2);
200  const Real determinant = edge1 * h;
201  const Real characteristic_length = std::max(edge1.norm(), edge2.norm());
202 
203  if (std::abs(determinant) <= _surface_tolerance * characteristic_length)
204  // Nearly parallel triangles are ignored because they do not provide a stable parity event.
205  return RayIntersection::Miss;
206 
207  const Real inv_determinant = 1.0 / determinant;
208  const Point s = point - tri.node_ref(0);
209  const Real u = inv_determinant * (s * h);
210  if (u < 0.0 || u > 1.0)
211  return RayIntersection::Miss;
212 
213  const Point q = s.cross(edge1);
214  const Real v = inv_determinant * (direction * q);
215  if (v < 0.0 || u + v > 1.0)
216  return RayIntersection::Miss;
217 
218  const Real t = inv_determinant * (edge2 * q);
219  if (t < 0.0)
220  // Intersections behind the ray origin do not contribute to +x parity counting.
221  return RayIntersection::Miss;
222  if (t <= _surface_tolerance)
223  // Hits too close to the ray origin are treated as ambiguous boundary situations.
225 
226  const Point intersection = point + t * direction;
227  const auto tolerance_sq = _surface_tolerance * _surface_tolerance;
228  if (geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(0), tri.node_ref(1)) <=
229  tolerance_sq ||
230  geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(1), tri.node_ref(2)) <=
231  tolerance_sq ||
232  geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(2), tri.node_ref(0)) <=
233  tolerance_sq)
234  // Edge and vertex hits are where odd/even parity counting is the least reliable.
236 
237  return RayIntersection::Hit;
238 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:50
const Real _surface_tolerance
Absolute tolerance used throughout validation and geometric classification.
auto max(const L &left, const R &right)
const Node & node_ref(const unsigned int i) const
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Real pointSegmentDistanceSq(const Point &point, const Point &a, const Point &b)
Compute the squared distance from a point to a 3-D line segment.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

Member Data Documentation

◆ _bounding_box

const libMesh::BoundingBox TriangleManifold::_bounding_box
private

Global bounding box of the transformed manifold.

Definition at line 136 of file TriangleManifold.h.

Referenced by boundingBox().

◆ _mesh

MeshBase& TriangleManifold::_mesh
private

◆ _num_y_cells

std::size_t TriangleManifold::_num_y_cells = 1
private

Number of yz-grid bins in the y direction.

Definition at line 110 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _num_z_cells

std::size_t TriangleManifold::_num_z_cells = 1
private

Number of yz-grid bins in the z direction.

Definition at line 113 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _point_locator

const std::unique_ptr<libMesh::PointLocatorBase> TriangleManifold::_point_locator
private

Pre-built point locator for fast proximity-to-surface detection.

Definition at line 139 of file TriangleManifold.h.

Referenced by pointOnSurface(), and TriangleManifold().

◆ _ray_grid

std::unordered_map<dof_id_type, std::vector<dof_id_type> > TriangleManifold::_ray_grid
private

Lookup from packed yz-grid cell index to triangles that could intersect the +x query ray.

Definition at line 128 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _surface_tolerance

const Real TriangleManifold::_surface_tolerance
private

Absolute tolerance used throughout validation and geometric classification.

Definition at line 133 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), contains(), pointInsideBoundingBox(), rayCandidates(), rayIntersectsTriangle(), and TriangleManifold().

◆ _y_cell_size

Real TriangleManifold::_y_cell_size = 1.0
private

Width of one yz-grid cell in the y direction.

Definition at line 122 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _y_min

Real TriangleManifold::_y_min = 0.0
private

Minimum global y coordinate used to map query points into yz-grid bins.

Definition at line 116 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _z_cell_size

Real TriangleManifold::_z_cell_size = 1.0
private

Width of one yz-grid cell in the z direction.

Definition at line 125 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().

◆ _z_min

Real TriangleManifold::_z_min = 0.0
private

Minimum global z coordinate used to map query points into yz-grid bins.

Definition at line 119 of file TriangleManifold.h.

Referenced by buildCandidateGrid(), and rayCandidates().


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