Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "TriangleManifold.h"
11 :
12 : #include "GeometryUtils.h"
13 : #include "MooseError.h"
14 :
15 : #include "libmesh/mesh_tools.h"
16 : #include "libmesh/unstructured_mesh.h"
17 :
18 : #include <algorithm>
19 : #include <cmath>
20 : #include <cstdint>
21 : #include <cstring>
22 :
23 : namespace TriangleManifoldUtils
24 : {
25 : std::uint64_t
26 57294 : packCell(const std::size_t y_index, const std::size_t z_index)
27 : {
28 : // The acceleration structure is indexed in yz because all rays travel in the +x direction.
29 57294 : return (static_cast<std::uint64_t>(y_index) << 32) | static_cast<std::uint64_t>(z_index);
30 : }
31 :
32 : std::size_t
33 126332 : cellIndex(const Real value, const Real min_value, const Real cell_size, const std::size_t num_cells)
34 : {
35 : // Clamp to the valid range so near-boundary roundoff does not produce an invalid cell id.
36 126332 : if (num_cells <= 1)
37 0 : return 0;
38 :
39 126332 : const auto index = static_cast<long long>(std::floor((value - min_value) / cell_size));
40 126332 : if (index < 0)
41 0 : return 0;
42 126332 : if (index >= static_cast<long long>(num_cells))
43 3316 : return num_cells - 1;
44 123016 : return static_cast<std::size_t>(index);
45 : }
46 :
47 : class SurfaceChecker : public libMesh::MeshTetInterface
48 : {
49 : public:
50 95 : explicit SurfaceChecker(UnstructuredMesh & mesh) : libMesh::MeshTetInterface(mesh) {}
51 0 : void triangulate() override { mooseError("SurfaceChecker is not meant for triangulation."); }
52 : std::string improveAndValidate();
53 : };
54 : }
55 :
56 95 : TriangleManifold::TriangleManifold(MeshBase & mesh, const Real surface_tolerance)
57 95 : : _mesh(mesh),
58 95 : _surface_tolerance(surface_tolerance),
59 95 : _bounding_box(MeshTools::create_bounding_box(_mesh)),
60 190 : _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 95 : finalize();
68 :
69 81 : _point_locator->set_close_to_point_tol(_surface_tolerance);
70 97 : }
71 :
72 : bool
73 21434 : TriangleManifold::contains(const Point & point) const
74 : {
75 : // Most points are rejected here before we perform any per-triangle work.
76 21434 : if (!pointInsideBoundingBox(point))
77 3888 : return false;
78 :
79 : // By design, near-surface points count as inside for downstream subdomain tagging.
80 17546 : if (pointOnSurface(point))
81 4 : return true;
82 :
83 : // Candidate filtering keeps the parity test from touching every triangle in large surfaces.
84 17542 : const auto candidates = rayCandidates(point);
85 17542 : unsigned int num_hits = 0;
86 :
87 731582 : for (const auto triangle_index : candidates)
88 : {
89 714320 : const auto tri = _mesh.elem_ptr(triangle_index);
90 : // Triangles wholly behind the ray origin cannot contribute to the parity count.
91 714320 : if (tri->loose_bounding_box().max()(0) < point(0) - _surface_tolerance)
92 334238 : continue;
93 :
94 380082 : switch (rayIntersectsTriangle(point, *tri))
95 : {
96 371740 : case RayIntersection::Miss:
97 371740 : break;
98 8062 : case RayIntersection::Hit:
99 : // Standard odd/even counting on a closed surface.
100 8062 : ++num_hits;
101 8062 : break;
102 280 : case RayIntersection::Ambiguous:
103 : // Edge and vertex grazing hits are exactly where parity counting becomes brittle.
104 280 : return containsBySolidAngle(point);
105 : }
106 : }
107 :
108 17262 : return num_hits % 2;
109 17542 : }
110 :
111 : void
112 95 : TriangleManifold::finalize()
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 95 : auto umesh = dynamic_cast<UnstructuredMesh *>(&_mesh);
117 95 : TriangleManifoldUtils::SurfaceChecker checker(*umesh);
118 95 : auto msg = checker.improveAndValidate();
119 95 : if (!msg.empty())
120 14 : mooseError(
121 : "The inputted surface mesh cannot be treated as manifold for the following reasons:\n",
122 : msg);
123 :
124 81 : buildCandidateGrid();
125 97 : }
126 :
127 : void
128 81 : TriangleManifold::buildCandidateGrid()
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 =
133 81 : std::max(boundingBox().max()(1) - boundingBox().min()(1), _surface_tolerance);
134 : const auto extent_z =
135 81 : std::max(boundingBox().max()(2) - boundingBox().min()(2), _surface_tolerance);
136 :
137 : // Choose a roughly square yz grid scaled by aspect ratio so candidate lists stay short.
138 : const auto target_cells =
139 81 : std::max<dof_id_type>(1, static_cast<dof_id_type>(std::sqrt(numTriangles())));
140 81 : const auto aspect = std::sqrt(extent_y / extent_z);
141 :
142 162 : _num_y_cells = std::max<dof_id_type>(
143 81 : 1, static_cast<dof_id_type>(std::lround(std::sqrt(target_cells) * aspect)));
144 162 : _num_z_cells = std::max<std::size_t>(
145 81 : 1, static_cast<std::size_t>(std::ceil(static_cast<Real>(target_cells) / _num_y_cells)));
146 :
147 81 : _y_min = boundingBox().min()(1);
148 81 : _z_min = boundingBox().min()(2);
149 81 : _y_cell_size = extent_y / _num_y_cells;
150 81 : _z_cell_size = extent_z / _num_z_cells;
151 :
152 22893 : for (const auto elem : _mesh.active_element_ptr_range())
153 : {
154 22812 : const auto triangle_index = elem->id();
155 22812 : const auto bbox = elem->loose_bounding_box();
156 :
157 : const auto y_start =
158 22812 : TriangleManifoldUtils::cellIndex(bbox.min()(1), _y_min, _y_cell_size, _num_y_cells);
159 : const auto y_stop =
160 22812 : TriangleManifoldUtils::cellIndex(bbox.max()(1), _y_min, _y_cell_size, _num_y_cells);
161 : const auto z_start =
162 22812 : TriangleManifoldUtils::cellIndex(bbox.min()(2), _z_min, _z_cell_size, _num_z_cells);
163 : const auto z_stop =
164 22812 : TriangleManifoldUtils::cellIndex(bbox.max()(2), _z_min, _z_cell_size, _num_z_cells);
165 :
166 : // Insert the triangle into every grid cell touched by its yz projection.
167 56752 : for (const auto iy : make_range(y_start, y_stop + 1))
168 73692 : for (const auto iz : make_range(z_start, z_stop + 1))
169 39752 : _ray_grid[TriangleManifoldUtils::packCell(iy, iz)].push_back(triangle_index);
170 81 : }
171 81 : }
172 :
173 : bool
174 21434 : TriangleManifold::pointInsideBoundingBox(const Point & point) const
175 : {
176 : // Inflate the global box by the tolerance so near-surface points are not rejected too early.
177 21434 : return point(0) >= boundingBox().min()(0) - _surface_tolerance &&
178 20792 : point(0) <= boundingBox().max()(0) + _surface_tolerance &&
179 20310 : point(1) >= boundingBox().min()(1) - _surface_tolerance &&
180 20002 : point(1) <= boundingBox().max()(1) + _surface_tolerance &&
181 61766 : point(2) >= boundingBox().min()(2) - _surface_tolerance &&
182 40974 : point(2) <= boundingBox().max()(2) + _surface_tolerance;
183 : }
184 :
185 : bool
186 17546 : TriangleManifold::pointOnSurface(const Point & point) const
187 : {
188 17546 : return (*_point_locator)(point) != nullptr;
189 : }
190 :
191 : TriangleManifold::RayIntersection
192 380082 : TriangleManifold::rayIntersectsTriangle(const Point & point, const libMesh::Elem & tri) const
193 : {
194 : // This is a fixed-direction Moller-Trumbore style ray/triangle intersection test.
195 380082 : static const Point direction(1.0, 0.0, 0.0);
196 :
197 380082 : const Point edge1 = tri.node_ref(1) - tri.node_ref(0);
198 380082 : const Point edge2 = tri.node_ref(2) - tri.node_ref(0);
199 380082 : const Point h = direction.cross(edge2);
200 380082 : const Real determinant = edge1 * h;
201 380082 : const Real characteristic_length = std::max(edge1.norm(), edge2.norm());
202 :
203 380082 : if (std::abs(determinant) <= _surface_tolerance * characteristic_length)
204 : // Nearly parallel triangles are ignored because they do not provide a stable parity event.
205 175192 : return RayIntersection::Miss;
206 :
207 204890 : const Real inv_determinant = 1.0 / determinant;
208 204890 : const Point s = point - tri.node_ref(0);
209 204890 : const Real u = inv_determinant * (s * h);
210 204890 : if (u < 0.0 || u > 1.0)
211 154856 : return RayIntersection::Miss;
212 :
213 50034 : const Point q = s.cross(edge1);
214 50034 : const Real v = inv_determinant * (direction * q);
215 50034 : if (v < 0.0 || u + v > 1.0)
216 41338 : return RayIntersection::Miss;
217 :
218 8696 : const Real t = inv_determinant * (edge2 * q);
219 8696 : if (t < 0.0)
220 : // Intersections behind the ray origin do not contribute to +x parity counting.
221 354 : return RayIntersection::Miss;
222 8342 : if (t <= _surface_tolerance)
223 : // Hits too close to the ray origin are treated as ambiguous boundary situations.
224 0 : return RayIntersection::Ambiguous;
225 :
226 8342 : const Point intersection = point + t * direction;
227 8342 : const auto tolerance_sq = _surface_tolerance * _surface_tolerance;
228 8342 : if (geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(0), tri.node_ref(1)) <=
229 8064 : tolerance_sq ||
230 8064 : geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(1), tri.node_ref(2)) <=
231 16406 : tolerance_sq ||
232 8064 : 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.
235 280 : return RayIntersection::Ambiguous;
236 :
237 8062 : return RayIntersection::Hit;
238 : }
239 :
240 : bool
241 280 : TriangleManifold::containsBySolidAngle(const Point & point) const
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 280 : Real total_angle = 0.0;
246 3640 : for (const auto elem : _mesh.active_element_ptr_range())
247 3360 : total_angle +=
248 3640 : geom_utils::solidAngle(point, elem->node_ref(0), elem->node_ref(1), elem->node_ref(2));
249 :
250 280 : return std::abs(total_angle) > 2.0 * libMesh::pi;
251 : }
252 :
253 : std::vector<dof_id_type>
254 17542 : TriangleManifold::rayCandidates(const Point & point) const
255 : {
256 : // If the query is outside the manifold's yz extent, the fixed +x ray cannot hit anything.
257 17542 : if (point(1) < boundingBox().min()(1) - _surface_tolerance ||
258 17542 : point(1) > boundingBox().max()(1) + _surface_tolerance ||
259 52626 : point(2) < boundingBox().min()(2) - _surface_tolerance ||
260 17542 : point(2) > boundingBox().max()(2) + _surface_tolerance)
261 0 : return {};
262 :
263 : const auto y_index =
264 17542 : TriangleManifoldUtils::cellIndex(point(1), _y_min, _y_cell_size, _num_y_cells);
265 : const auto z_index =
266 17542 : TriangleManifoldUtils::cellIndex(point(2), _z_min, _z_cell_size, _num_z_cells);
267 17542 : const auto it = _ray_grid.find(TriangleManifoldUtils::packCell(y_index, z_index));
268 17542 : if (it == _ray_grid.end())
269 5888 : return {};
270 :
271 : // Return by value to keep the helper simple and avoid exposing internal storage.
272 11654 : return it->second;
273 : }
274 :
275 : std::string
276 95 : TriangleManifoldUtils::SurfaceChecker::improveAndValidate()
277 : {
278 : using libMesh::MeshTetInterface;
279 95 : auto result = improve_hull_integrity();
280 :
281 95 : if (result.empty())
282 162 : return "";
283 :
284 14 : std::ostringstream err_msg;
285 14 : if (result.count(NON_TRI3))
286 5 : err_msg << "- At least one non-Tri3 element was found.\n" << std::endl;
287 14 : if (result.count(MISSING_NEIGHBOR))
288 9 : err_msg << "- At least one triangle without three neighbors was found.\n" << std::endl;
289 14 : if (result.count(EMPTY_MESH))
290 2 : err_msg << "- The surface mesh was empty\n" << std::endl;
291 14 : if (result.count(MISSING_BACKLINK))
292 3 : err_msg << "- At least one triangle neighbor without a return neighbor link was found.\n";
293 14 : if (result.count(BAD_NEIGHBOR_NODES))
294 3 : err_msg << "- At least one triangle neighbor without expected node links was found.\n";
295 14 : if (result.count(NON_ORIENTED))
296 3 : err_msg << "- At least one triangle neighbor with an inconsistent orientation was found.\n";
297 14 : if (result.count(BAD_NEIGHBOR_LINKS))
298 : err_msg << "- At least one triangle neighbor with inconsistent node and neighbor links was "
299 0 : "found.\n";
300 14 : if (result.count(DEGENERATE_ELEMENT))
301 : err_msg << "- At least one input triangle is degenerate, with near-zero area relative to the "
302 0 : "manifold.\n";
303 14 : if (result.count(DEGENERATE_MESH))
304 4 : err_msg << "- Mesh is degenerate, with zero thickness in at least one direction.\n";
305 14 : return err_msg.str();
306 95 : }
|
