1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 #ifndef AXOM_QUEST_MESH_TESTER_HPP_
7 #define AXOM_QUEST_MESH_TESTER_HPP_
8
9 // Axom includes
10 #include "axom/config.hpp"
11 #include "axom/core.hpp"
12 #include "axom/primal.hpp"
13 #include "axom/spin.hpp"
14 #include "axom/mint.hpp"
15
16 // C++ includes
17 #include <cmath>
18 #include <algorithm>
19 #include <vector>
20 #include <unordered_map>
21 #include <functional> // for std::hash
22
23 // BVH includes
24 #if defined(AXOM_USE_RAJA)
25 #include "RAJA/RAJA.hpp"
26 #include "axom/spin/BVH.hpp"
27 #include "axom/mint/execution/internal/structured_exec.hpp"
28 #endif
29
30 // Umpire
31 #if defined(AXOM_USE_UMPIRE)
32 #include "umpire/strategy/QuickPool.hpp"
33 #endif
34
35 // C/C++ includes
36 #include <utility>
37 #include <vector>
38
39 /*!
40 * \file MeshTester.hpp
41 * \brief Defines functions to test Quest meshes for common defects.
42 */
43
44 namespace axom
45 {
46 namespace quest
47 {
48 namespace detail
49 {
50 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
51 using Triangle3 = primal::Triangle<double, 3>;
52 using SpatialBoundingBox = primal::BoundingBox<double, 3>;
53 using UniformGrid3 = spin::UniformGrid<int, 3>;
54 using Point3 = primal::Point<double, 3>;
55 } // namespace detail
56
57 /*! Enumeration indicating mesh watertightness */
58 enum class WatertightStatus : signed char
59 {
60 WATERTIGHT = 0, ///< Each edge in a surface mesh is incident in two cells
61 NOT_WATERTIGHT, ///< Each edge is incident in one or two cells
62 CHECK_FAILED ///< Calculation failed (possibly a non-manifold mesh)
63 };
64
getMeshTriangle(axom::IndexType i,detail::UMesh * surface_mesh)65 inline detail::Triangle3 getMeshTriangle(axom::IndexType i,
66 detail::UMesh* surface_mesh)
67 {
68 SLIC_ASSERT(surface_mesh->getNumberOfCellNodes(i) == 3);
69
70 detail::Triangle3 tri;
71
72 const axom::IndexType* triCell = surface_mesh->getCellNodeIDs(i);
73
74 const double* x = surface_mesh->getCoordinateArray(mint::X_COORDINATE);
75 const double* y = surface_mesh->getCoordinateArray(mint::Y_COORDINATE);
76 const double* z = surface_mesh->getCoordinateArray(mint::Z_COORDINATE);
77
78 for(int n = 0; n < 3; ++n)
79 {
80 const axom::IndexType nodeIdx = triCell[n];
81 tri[n][0] = x[nodeIdx];
82 tri[n][1] = y[nodeIdx];
83 tri[n][2] = z[nodeIdx];
84 }
85
86 return tri;
87 }
88
89 /// \name Mesh test and repair
90 /// @{
91
92 /*!
93 * \brief Find self-intersections and degenerate triangles in a surface mesh
94 * utilizing a Bounding Volume Hierarchy.
95 *
96 * \param [in] surface_mesh A triangle mesh in three dimensions
97 * \param [out] intersection Pairs of indices of intersecting mesh triangles
98 * \param [out] degenerateIndices indices of degenerate mesh triangles
99 * \param [in] intersectionThreshold Tolerance threshold for triangle
100 * intersection tests (default: 1E-8)
101 * After running this function over a surface mesh, intersection will be filled
102 * with pairs of indices of intersecting triangles and degenerateIndices will
103 * be filled with the indices of the degenerate triangles in the mesh.
104 * Triangles that share vertex pairs (adjacent triangles in a watertight
105 * surface mesh) are not reported as intersecting. Degenerate triangles
106 * are not reported as intersecting other triangles.
107 *
108 */
109 #if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE)
110 template <typename ExecSpace, typename FloatType>
findTriMeshIntersectionsBVH(mint::UnstructuredMesh<mint::SINGLE_SHAPE> * surface_mesh,std::vector<std::pair<int,int>> & intersections,std::vector<int> & degenerateIndices,double intersectionThreshold=1E-8)111 void findTriMeshIntersectionsBVH(
112 mint::UnstructuredMesh<mint::SINGLE_SHAPE>* surface_mesh,
113 std::vector<std::pair<int, int>>& intersections,
114 std::vector<int>& degenerateIndices,
115 double intersectionThreshold = 1E-8)
116 {
117 AXOM_PERF_MARK_FUNCTION("findTriMeshIntersectionsBVH");
118
119 SLIC_INFO("Running BVH intersection algorithm "
120 << " in execution Space: "
121 << axom::execution_space<ExecSpace>::name());
122
123 constexpr size_t POOL_SIZE = (1024 * 1024 * 1024) + 1;
124 umpire::ResourceManager& rm = umpire::ResourceManager::getInstance();
125
126 // Use device pool for CUDA policy, host pool otherwise
127 umpire::Allocator allocator =
128 (axom::execution_space<ExecSpace>::onDevice()
129 ? rm.getAllocator(umpire::resource::Device)
130 : rm.getAllocator(axom::execution_space<ExecSpace>::allocatorID()));
131
132 umpire::Allocator pool_allocator =
133 rm.makeAllocator<umpire::strategy::QuickPool>(allocator.getName() + "_POOL",
134 allocator,
135 POOL_SIZE);
136
137 const int poolID = pool_allocator.getId();
138
139 // Get allocator
140 const int current_allocator = axom::getDefaultAllocatorID();
141 axom::setDefaultAllocator(poolID);
142
143 constexpr int NDIMS = 3;
144 const int ncells = surface_mesh->getNumberOfCells();
145
146 using BoxType = typename primal::BoundingBox<FloatType, NDIMS>;
147
148 int* ZERO =
149 axom::allocate<int>(1, getUmpireResourceAllocatorID(umpire::resource::Host));
150 ZERO[0] = 0;
151
152 detail::Triangle3* tris = axom::allocate<detail::Triangle3>(ncells);
153
154 // Marks each cell/triangle as degenerate (1) or not (0)
155 int* degenerate = axom::allocate<int>(ncells);
156
157 // Each access-aligned bounding box represented by 2 (x,y,z) points
158 BoxType* aabbs = axom::allocate<BoxType>(ncells);
159
160 // Initialize the bounding box for each Triangle and marks
161 // if the Triangle is degenerate.
162 AXOM_PERF_MARK_SECTION(
163 "init_tri_bb",
164 mint::for_all_cells<ExecSpace, mint::xargs::coords>(
165 surface_mesh,
166 AXOM_LAMBDA(IndexType cellIdx,
167 numerics::Matrix<double> & coords,
168 const IndexType* AXOM_UNUSED_PARAM(nodeIds)) {
169 detail::Triangle3 tri;
170
171 for(IndexType inode = 0; inode < 3; ++inode)
172 {
173 const double* node = coords.getColumn(inode);
174 tri[inode][0] = node[mint::X_COORDINATE];
175 tri[inode][1] = node[mint::Y_COORDINATE];
176 tri[inode][2] = node[mint::Z_COORDINATE];
177 } // END for all cells nodes
178
179 degenerate[cellIdx] = (tri.degenerate() ? 1 : 0);
180
181 tris[cellIdx] = tri;
182
183 aabbs[cellIdx] = compute_bounding_box(tri);
184 }););
185
186 // Copy degenerate data back to host
187 int* host_degenerate =
188 axom::allocate<int>(ncells,
189 getUmpireResourceAllocatorID(umpire::resource::Host));
190 axom::copy(host_degenerate, degenerate, ncells * sizeof(int));
191
192 // Return degenerateIndices
193 for(int i = 0; i < ncells; i++)
194 {
195 if(host_degenerate[i] == 1)
196 {
197 degenerateIndices.push_back(i);
198 }
199 }
200
201 // Construct BVH
202 axom::spin::BVH<NDIMS, ExecSpace, FloatType> bvh;
203 bvh.setAllocatorID(poolID);
204 bvh.initialize(aabbs, ncells);
205
206 // Run find algorithm
207 IndexType* offsets = axom::allocate<IndexType>(ncells);
208 IndexType* counts = axom::allocate<IndexType>(ncells);
209 IndexType* candidates = nullptr;
210 bvh.findBoundingBoxes(offsets,
211 counts,
212 candidates,
213 ncells,
214 reinterpret_cast<BoxType*>(aabbs));
215
216 // Get the total number of candidates
217 using REDUCE_POL = typename axom::execution_space<ExecSpace>::reduce_policy;
218
219 RAJA::ReduceSum<REDUCE_POL, int> totalCandidates(0);
220 for_all<ExecSpace>(
221 ncells,
222 AXOM_LAMBDA(IndexType i) { totalCandidates += counts[i]; });
223
224 //Deallocate no longer needed variables
225 axom::deallocate(aabbs);
226
227 axom::deallocate(degenerate);
228 axom::deallocate(host_degenerate);
229
230 int* intersection_pairs = axom::allocate<int>(totalCandidates.get() * 2);
231
232 using ATOMIC_POL = typename axom::execution_space<ExecSpace>::atomic_policy;
233 int* counter = axom::allocate<int>(1);
234 axom::copy(counter, ZERO, sizeof(int));
235
236 // Initialize triangle indices and valid candidates
237 IndexType* indices = axom::allocate<IndexType>(totalCandidates.get());
238 IndexType* validCandidates = axom::allocate<IndexType>(totalCandidates.get());
239 int* numValidCandidates = axom::allocate<int>(1);
240 axom::copy(numValidCandidates, ZERO, sizeof(int));
241
242 AXOM_PERF_MARK_SECTION(
243 "init_candidates",
244 for_all<ExecSpace>(
245 ncells,
246 AXOM_LAMBDA(IndexType i) {
247 for(int j = 0; j < counts[i]; j++)
248 {
249 if(i < candidates[offsets[i] + j])
250 {
251 auto idx = RAJA::atomicAdd<ATOMIC_POL>(numValidCandidates, 1);
252 indices[idx] = i;
253 validCandidates[idx] = candidates[offsets[i] + j];
254 }
255 }
256 }););
257
258 // Copy numValidCandidates back to host
259 int* host_numValidCandidates =
260 axom::allocate<int>(1, getUmpireResourceAllocatorID(umpire::resource::Host));
261 axom::copy(host_numValidCandidates, numValidCandidates, sizeof(int));
262
263 AXOM_PERF_MARK_SECTION("find_tri_pairs",
264 for_all<ExecSpace>(
265 *host_numValidCandidates,
266 AXOM_LAMBDA(IndexType i) {
267 int index = indices[i];
268 int candidate = validCandidates[i];
269 if(primal::intersect(tris[index],
270 tris[candidate],
271 false,
272 intersectionThreshold))
273 {
274 auto idx = RAJA::atomicAdd<ATOMIC_POL>(counter, 2);
275 intersection_pairs[idx] = index;
276 intersection_pairs[idx + 1] = candidate;
277 }
278 }););
279
280 // Copy intersection pairs and counter back to host
281 int* host_counter =
282 axom::allocate<int>(1, getUmpireResourceAllocatorID(umpire::resource::Host));
283 axom::copy(host_counter, counter, sizeof(int));
284
285 int* host_intersection_pairs =
286 axom::allocate<int>(totalCandidates.get() * 2,
287 getUmpireResourceAllocatorID(umpire::resource::Host));
288 axom::copy(host_intersection_pairs,
289 intersection_pairs,
290 totalCandidates.get() * 2 * sizeof(int));
291
292 // Initialize pairs of clashes
293 for(int i = 0; i < host_counter[0]; i += 2)
294 {
295 intersections.push_back(std::make_pair(host_intersection_pairs[i],
296 host_intersection_pairs[i + 1]));
297 }
298
299 // Deallocate
300 axom::deallocate(tris);
301
302 axom::deallocate(offsets);
303 axom::deallocate(counts);
304 axom::deallocate(candidates);
305 axom::deallocate(indices);
306 axom::deallocate(validCandidates);
307 axom::deallocate(host_numValidCandidates);
308
309 axom::deallocate(intersection_pairs);
310 axom::deallocate(host_intersection_pairs);
311 axom::deallocate(counter);
312 axom::deallocate(host_counter);
313
314 axom::setDefaultAllocator(current_allocator);
315 }
316 #endif
317
318 /*!
319 * \brief Find self-intersections and degenerate triangles in a surface mesh
320 * utilizing a Uniform Grid.
321 *
322 * \param [in] surface_mesh A triangle mesh in three dimensions
323 * \param [out] intersection Pairs of indices of intersecting mesh triangles
324 * \param [out] degenerateIndices indices of degenerate mesh triangles
325 * \param [in] spatialIndexResolution The grid resolution for the index
326 * structure (default: 0)
327 * \param [in] intersectionThreshold Tolerance threshold for triangle
328 * intersection tests (default: 1E-8)
329 *
330 * After running this function over a surface mesh, intersection will be filled
331 * with pairs of indices of intersecting triangles and degenerateIndices will
332 * be filled with the indices of the degenerate triangles in the mesh.
333 * Triangles that share vertex pairs (adjacent triangles in a watertight
334 * surface mesh) are not reported as intersecting. Degenerate triangles
335 * are not reported as intersecting other triangles.
336 *
337 * This function uses a quest::UniformGrid spatial index. Input
338 * spatialIndexResolution specifies the bin size for the UniformGrid. The
339 * default value of 0 causes this routine to calculate a heuristic bin size
340 * based on the cube root of the number of cells in the mesh.
341 */
342 void findTriMeshIntersections(mint::UnstructuredMesh<mint::SINGLE_SHAPE>* surface_mesh,
343 std::vector<std::pair<int, int>>& intersections,
344 std::vector<int>& degenerateIndices,
345 int spatialIndexResolution = 0,
346 double intersectionThreshold = 1E-8);
347
348 /*!
349 * \brief Check a surface mesh for holes using its face relation.
350 *
351 * \param [in] surface_mesh A surface mesh in three dimensions
352 *
353 * \returns status If the mesh is watertight, is not watertight, or
354 * if an error occurred (possibly due to non-manifold mesh).
355 *
356 * \note This method marks the cells on the boundary by creating a new
357 * cell-centered field variable, called "boundary", on the given input mesh.
358 *
359 * \note This function computes the mesh's cell-face and face-vertex relations.
360 * For large meshes, this can take a long time. The relations are used to
361 * check for holes, and remain cached with the mesh after this function
362 * finishes.
363 */
364 WatertightStatus isSurfaceMeshWatertight(
365 mint::UnstructuredMesh<mint::SINGLE_SHAPE>* surface_mesh);
366
367 /*!
368 * \brief Mesh repair function to weld vertices that are closer than \a eps
369 *
370 * \param [in,out] surface_mesh A pointer to a pointer to a triangle mesh
371 * \param [in] eps Distance threshold for welding vertices (using the max norm)
372 *
373 * \pre \a eps must be greater than zero
374 * \pre \a surface_mesh is a pointer to a pointer to a non-null triangle mesh.
375 * \post The triangles of \a surface_mesh are reindexed using the welded
376 * vertices and degenerate triangles are removed. The mesh can still contain
377 * vertices that are not referenced by any triangles.
378 *
379 * This utility function repairs an input triangle mesh (embedded in three
380 * dimensional space) by 'welding' vertices that are closer than \a eps.
381 * The vertices are quantized to an integer lattice with spacing \a eps
382 * and vertices that fall into the same cell on this lattice are identified.
383 * All identified vertices are given the coordinates of the first such vertex
384 * and all incident triangles use the same index for this vertex.
385 *
386 * The input mesh can be a "soup of triangles", where the vertices
387 * of adjacent triangles have distinct indices. After running this function,
388 * vertices that are closer than \a eps are welded, and their incident
389 * triangles use the new vertex indices. Thus, the output is an
390 * "indexed triangle mesh".
391 *
392 * This function also removes degenerate triangles from the mesh. These
393 * are triangles without three distinct vertices after the welding.
394 *
395 * \note This function is destructive. It modifies the input triangle
396 * mesh in place.
397 * \note The distance metric in this function uses the "max" norm (l_inf).
398 */
399 void weldTriMeshVertices(mint::UnstructuredMesh<mint::SINGLE_SHAPE>** surface_mesh,
400 double eps);
401
402 /// @}
403
404 } // end namespace quest
405 } // end namespace axom
406
407 #endif // AXOM_QUEST_MESH_TESTER_HPP_
408