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 QUEST_TEST_UTILITIES_HPP_
7 #define QUEST_TEST_UTILITIES_HPP_
8
9 #include "axom/core.hpp"
10 #include "axom/slic.hpp"
11 #include "axom/primal.hpp"
12 #include "axom/mint.hpp"
13
14 #include <math.h>
15
16 /*!
17 * \file quest_test_utilities.hpp
18 *
19 * This file contains several utility functions to facilitate in the development
20 * of unit tests in quest, e.g., generating a random points and simple
21 * unstructured mesh types, such as, the surface mesh of a sphere with a
22 * given center and radius and an octahedron.
23 */
24
25 namespace axom
26 {
27 namespace quest
28 {
29 namespace utilities
30 {
31 namespace primal = axom::primal;
32 namespace mint = axom::mint;
33
34 constexpr double DEG_TO_RAD = M_PI / 180.;
35
36 /*!
37 * \brief Gets a surface mesh instance for the sphere.
38 *
39 * \param [in,out] mesh pointer to the mesh instance, where to store the mesh
40 * \param [in] SPHERE_CENTER the center of the sphere
41 * \param [in] RADIUS the radius of the sphere
42 * \param [in] THETA_RES the theta resolution for generating the sphere
43 * \param [in] PHI_RES the phi resolution for generating the sphere
44 *
45 * \pre mesh != nullptr
46 * \pre mesh->getDimension() == 3
47 * \pre mesh->getMeshType() == mint::UNSTRUCTURED_MESH
48 * \pre mesh->hasMixedCellTypes() == false
49 * \pre mesh->getCellType() == mint::TRIANGLE
50 */
getSphereSurfaceMesh(mint::UnstructuredMesh<mint::SINGLE_SHAPE> * mesh,const double * SPHERE_CENTER,double RADIUS,int THETA_RES,int PHI_RES)51 void getSphereSurfaceMesh(mint::UnstructuredMesh<mint::SINGLE_SHAPE>* mesh,
52 const double* SPHERE_CENTER,
53 double RADIUS,
54 int THETA_RES,
55 int PHI_RES)
56 {
57 SLIC_ASSERT(mesh != nullptr);
58 SLIC_ASSERT(mesh->getDimension() == 3);
59 SLIC_ASSERT(mesh->getMeshType() == mint::UNSTRUCTURED_MESH);
60 SLIC_ASSERT(mesh->hasMixedCellTypes() == false);
61 SLIC_ASSERT(mesh->getCellType() == mint::TRIANGLE);
62
63 const IndexType totalNodes = THETA_RES * PHI_RES;
64 const IndexType totalCells = (2 * THETA_RES) + (THETA_RES * PHI_RES);
65 mesh->reserve(totalNodes, totalCells);
66
67 const double theta_start = 0;
68 const double theta_end = 360 * DEG_TO_RAD;
69 const double phi_start = 0;
70 const double phi_end = 180 * DEG_TO_RAD;
71
72 double x[3];
73 double n[3];
74 axom::IndexType c[3];
75
76 // North pole point
77 x[0] = SPHERE_CENTER[0];
78 x[1] = SPHERE_CENTER[1];
79 x[2] = SPHERE_CENTER[2] + RADIUS;
80 mesh->appendNode(x[0], x[1], x[2]);
81
82 // South pole point
83 x[2] = SPHERE_CENTER[2] - RADIUS;
84 mesh->appendNode(x[0], x[1], x[2]);
85
86 // Calculate spacing
87 const double dphi = (phi_end - phi_start) / (static_cast<double>(PHI_RES - 1));
88 const double dtheta =
89 (theta_end - theta_start) / (static_cast<double>(THETA_RES - 1));
90
91 // Generate points
92 for(int i = 0; i < THETA_RES; ++i)
93 {
94 const double theta = theta_start + i * dtheta;
95
96 for(int j = 0; j < PHI_RES - 2; ++j)
97 {
98 const double phi = phi_start + j * dphi;
99 const double radius = RADIUS * sin(phi);
100
101 n[0] = radius * cos(theta);
102 n[1] = radius * sin(theta);
103 n[2] = RADIUS * cos(phi);
104 x[0] = n[0] + SPHERE_CENTER[0];
105 x[1] = n[1] + SPHERE_CENTER[1];
106 x[2] = n[2] + SPHERE_CENTER[2];
107
108 mesh->appendNode(x[0], x[1], x[2]);
109
110 } // END for all j
111
112 } // END for all i
113
114 const int phiResolution = PHI_RES - 2; // taking in to account two pole points.
115 int stride = phiResolution * THETA_RES;
116
117 // Generate mesh connectivity around north pole
118 for(int i = 0; i < THETA_RES; ++i)
119 {
120 c[2] = phiResolution * i + /* number of poles */ 2;
121 c[1] = (phiResolution * (i + 1) % stride) + /* number of poles */ 2;
122 c[0] = 0;
123 mesh->appendCell(c);
124 } // END for
125
126 // Generate mesh connectivity around south pole
127 int offset = PHI_RES - 1;
128 for(int i = 0; i < THETA_RES; ++i)
129 {
130 c[2] = phiResolution * i + offset;
131 c[1] = (phiResolution * (i + 1) % stride) + offset;
132 c[0] = 1;
133 mesh->appendCell(c);
134 }
135
136 // Generate mesh connectivity in between poles
137 for(int i = 0; i < THETA_RES; ++i)
138 {
139 for(int j = 0; j < PHI_RES - 3; ++j)
140 {
141 c[0] = phiResolution * i + j + 2;
142 c[1] = c[0] + 1;
143 c[2] = ((phiResolution * (i + 1) + j) % stride) + 3;
144
145 mesh->appendCell(c);
146
147 c[1] = c[2];
148 c[2] = c[1] - 1;
149 mesh->appendCell(c);
150 } // END for all j
151 } // END for all i
152 }
153
154 /*!
155 * \brief Simple utility to generate a Point whose entries
156 * are random values in the range [beg, end]
157 */
158 template <int DIM>
randomSpacePt(double beg,double end)159 primal::Point<double, DIM> randomSpacePt(double beg, double end)
160 {
161 primal::Point<double, DIM> pt;
162 for(int i = 0; i < DIM; ++i) pt[i] = axom::utilities::random_real(beg, end);
163
164 return pt;
165 }
166
167 /*!
168 * \brief Simple utility to find the centroid of two points
169 */
170 template <int DIM>
getCentroid(const primal::Point<double,DIM> & pt0,const primal::Point<double,DIM> & pt1)171 primal::Point<double, DIM> getCentroid(const primal::Point<double, DIM>& pt0,
172 const primal::Point<double, DIM>& pt1)
173 {
174 return primal::Point<double, DIM>((pt0.array() + pt1.array()) / 2.);
175 }
176
177 /*!
178 * \brief Simple utility to find the centroid of three points
179 */
180 template <int DIM>
getCentroid(const primal::Point<double,DIM> & pt0,const primal::Point<double,DIM> & pt1,const primal::Point<double,DIM> & pt2)181 primal::Point<double, DIM> getCentroid(const primal::Point<double, DIM>& pt0,
182 const primal::Point<double, DIM>& pt1,
183 const primal::Point<double, DIM>& pt2)
184 {
185 return primal::Point<double, DIM>((pt0.array() + pt1.array() + pt2.array()) /
186 3.);
187 }
188
189 /*!
190 * \brief Simple utility to find the centroid of four points
191 */
192 template <int DIM>
getCentroid(const primal::Point<double,DIM> & pt0,const primal::Point<double,DIM> & pt1,const primal::Point<double,DIM> & pt2,const primal::Point<double,DIM> & pt3)193 primal::Point<double, DIM> getCentroid(const primal::Point<double, DIM>& pt0,
194 const primal::Point<double, DIM>& pt1,
195 const primal::Point<double, DIM>& pt2,
196 const primal::Point<double, DIM>& pt3)
197 {
198 return primal::Point<double, DIM>(
199 (pt0.array() + pt1.array() + pt2.array() + pt3.array()) / 4.);
200 }
201
202 /*!
203 * \brief Utility function to generate a triangle mesh of an octahedron
204 * Vertices of the octahedron are at +-i, +-j and +-k.
205 * \note The caller must delete the mesh
206 */
make_octahedron_mesh()207 axom::mint::Mesh* make_octahedron_mesh()
208 {
209 using VertexIndex = axom::IndexType;
210 using SpacePt = primal::Point<double, 3>;
211 using SpaceTriangle = primal::Triangle<double, 3>;
212
213 enum
214 {
215 POS_X,
216 NEG_X,
217 POS_Y,
218 NEG_Y,
219 POS_Z,
220 NEG_Z
221 };
222
223 // The six vertices of the octahedron
224 const int NUM_VERTS = 6;
225 SpacePt verts[NUM_VERTS] = {SpacePt::make_point(1., 0., 0.),
226 SpacePt::make_point(-1., 0., 0.),
227 SpacePt::make_point(0, 1., 0.),
228 SpacePt::make_point(0, -1., 0.),
229 SpacePt::make_point(0, 0, 1.),
230 SpacePt::make_point(0, 0, -1.)};
231
232 // The eight triangles of the octahedron
233 // Explicit representation of triangle-vertex incidence relation
234 // Note: We are orienting the triangles with normals pointing outside
235 const int NUM_TRIS = 8;
236 const int VERTS_PER_TRI = 3;
237 VertexIndex tvRelation[NUM_TRIS * VERTS_PER_TRI] = {
238 POS_Z, POS_X, POS_Y, POS_Z, POS_Y, NEG_X, POS_Z, NEG_X,
239 NEG_Y, POS_Z, NEG_Y, POS_X, NEG_Z, POS_Y, POS_X, NEG_Z,
240 NEG_X, POS_Y, NEG_Z, NEG_Y, NEG_X, NEG_Z, POS_X, NEG_Y};
241
242 // First, confirm that all triangle normals point away from the origin
243 for(int i = 0; i < NUM_TRIS; ++i)
244 {
245 int baseIndex = i * VERTS_PER_TRI;
246 SpaceTriangle tri(verts[tvRelation[baseIndex + 0]],
247 verts[tvRelation[baseIndex + 1]],
248 verts[tvRelation[baseIndex + 2]]);
249
250 SLIC_ASSERT(primal::ON_NEGATIVE_SIDE == primal::orientation(SpacePt(), tri));
251 }
252
253 // Now create an unstructured triangle mesh from the two arrays
254 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
255 UMesh* triMesh = new UMesh(3, mint::TRIANGLE);
256
257 // insert verts
258 for(int i = 0; i < NUM_VERTS; ++i)
259 {
260 triMesh->appendNode(verts[i][0], verts[i][1], verts[i][2]);
261 }
262
263 // insert triangles
264 for(int i = 0; i < NUM_TRIS; ++i)
265 {
266 triMesh->appendCell(&tvRelation[i * VERTS_PER_TRI]);
267 }
268
269 SLIC_ASSERT(NUM_VERTS == triMesh->getNumberOfNodes());
270 SLIC_ASSERT(NUM_TRIS == triMesh->getNumberOfCells());
271
272 return triMesh;
273 }
274
275 /*!
276 * \brief Utility function to generate the triangle mesh of a tetrahedron
277 *
278 * Vertices are close to, but not exactly: (0, 0, 20),
279 * (-18.21, 4.88, -6.66), (4.88, -18.21, -6.66), (13.33, 13.33, -6.66)
280 *
281 * \note The caller must delete the mesh
282 */
make_tetrahedron_mesh()283 mint::Mesh* make_tetrahedron_mesh()
284 {
285 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
286
287 UMesh* surface_mesh = new UMesh(3, mint::TRIANGLE);
288 surface_mesh->appendNode(-0.000003, -0.000003, 19.999999);
289 surface_mesh->appendNode(-18.213671, 4.880339, -6.666668);
290 surface_mesh->appendNode(4.880339, -18.213671, -6.666668);
291 surface_mesh->appendNode(13.333334, 13.333334, -6.666663);
292 axom::IndexType cell[3];
293 cell[0] = 0;
294 cell[1] = 1;
295 cell[2] = 2;
296 surface_mesh->appendCell(cell);
297 cell[0] = 0;
298 cell[1] = 3;
299 cell[2] = 1;
300 surface_mesh->appendCell(cell);
301 cell[0] = 0;
302 cell[1] = 2;
303 cell[2] = 3;
304 surface_mesh->appendCell(cell);
305 cell[0] = 1;
306 cell[1] = 3;
307 cell[2] = 2;
308 surface_mesh->appendCell(cell);
309
310 return surface_mesh;
311 }
312
313 /*!
314 * \brief Generate the triangle mesh of a cracked tetrahedron
315 *
316 * Vertices are those of the tetrahedron generated by
317 * make_tetrahedron_mesh(), except that one corner of one face is
318 * displaced out away from its two neighbors. This mesh is not
319 * watertight and has no self-intersections.
320 *
321 * \note The caller must delete the mesh
322 */
make_crackedtet_mesh()323 mint::Mesh* make_crackedtet_mesh()
324 {
325 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
326
327 UMesh* surface_mesh = new UMesh(3, mint::TRIANGLE);
328 surface_mesh->appendNode(-0.000003, -0.000003, 19.999999);
329 surface_mesh->appendNode(-18.213671, 4.880339, -6.666668);
330 surface_mesh->appendNode(4.880339, -18.213671, -6.666668);
331 surface_mesh->appendNode(13.333334, 13.333334, -6.666663);
332 surface_mesh->appendNode(-0.200003, -0.100003, 18.999999);
333 axom::IndexType cell[3];
334 cell[0] = 4;
335 cell[1] = 1;
336 cell[2] = 2;
337 surface_mesh->appendCell(cell);
338 cell[0] = 0;
339 cell[1] = 3;
340 cell[2] = 1;
341 surface_mesh->appendCell(cell);
342 cell[0] = 0;
343 cell[1] = 2;
344 cell[2] = 3;
345 surface_mesh->appendCell(cell);
346 cell[0] = 1;
347 cell[1] = 3;
348 cell[2] = 2;
349 surface_mesh->appendCell(cell);
350
351 return surface_mesh;
352 }
353
354 /*!
355 * \brief Generate the triangle mesh of a caved-in tetrahedron
356 *
357 * Vertices are those of the tetrahedron generated by
358 * make_tetrahedron_mesh(), except that one corner of one face is
359 * displaced into its two neighbors. This mesh is not watertight
360 * and has two self-intersections.
361 *
362 * \note The caller must delete the mesh
363 */
make_cavedtet_mesh()364 mint::Mesh* make_cavedtet_mesh()
365 {
366 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
367
368 UMesh* surface_mesh = new UMesh(3, mint::TRIANGLE);
369 surface_mesh->appendNode(2.00003, 1.00003, 18.999999);
370 surface_mesh->appendNode(-18.213671, 4.880339, -6.666668);
371 surface_mesh->appendNode(4.880339, -18.213671, -6.666668);
372 surface_mesh->appendNode(-0.000003, -0.000003, 19.999999);
373 surface_mesh->appendNode(13.333334, 13.333334, -6.666663);
374 axom::IndexType cell[3];
375 cell[0] = 0;
376 cell[1] = 1;
377 cell[2] = 2;
378 surface_mesh->appendCell(cell);
379 cell[0] = 3;
380 cell[1] = 4;
381 cell[2] = 1;
382 surface_mesh->appendCell(cell);
383 cell[0] = 3;
384 cell[1] = 2;
385 cell[2] = 4;
386 surface_mesh->appendCell(cell);
387 cell[0] = 1;
388 cell[1] = 4;
389 cell[2] = 2;
390 surface_mesh->appendCell(cell);
391
392 return surface_mesh;
393 }
394
395 /*!
396 * \brief Generate the triangle mesh of a caved-in tetrahedron
397 * with added degenerate triangles
398 *
399 * Vertices are those of the tetrahedron generated by
400 * make_tetrahedron_mesh(), except that one corner of one face is
401 * displaced into its two neighbors. This mesh is not watertight
402 * has two self-intersections, and has two degenerate triangles.
403 *
404 * \note The caller must delete the mesh
405 */
make_degen_cavedtet_mesh()406 mint::Mesh* make_degen_cavedtet_mesh()
407 {
408 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
409
410 UMesh* surface_mesh = new UMesh(3, mint::TRIANGLE);
411 surface_mesh->appendNode(2.00003, 1.00003, 18.999999);
412 surface_mesh->appendNode(-18.213671, 4.880339, -6.666668);
413 surface_mesh->appendNode(4.880339, -18.213671, -6.666668);
414 surface_mesh->appendNode(-0.000003, -0.000003, 19.999999);
415 surface_mesh->appendNode(13.333334, 13.333334, -6.666663);
416 axom::IndexType cell[3];
417 cell[0] = 0;
418 cell[1] = 1;
419 cell[2] = 2;
420 surface_mesh->appendCell(cell);
421 cell[0] = 3;
422 cell[1] = 4;
423 cell[2] = 1;
424 surface_mesh->appendCell(cell);
425 cell[0] = 3;
426 cell[1] = 2;
427 cell[2] = 4;
428 surface_mesh->appendCell(cell);
429 cell[0] = 1;
430 cell[1] = 4;
431 cell[2] = 2;
432 surface_mesh->appendCell(cell);
433 cell[0] = 0;
434 cell[1] = 0;
435 cell[2] = 0;
436 surface_mesh->appendCell(cell);
437 cell[0] = 3;
438 cell[1] = 4;
439 cell[2] = 3;
440 surface_mesh->appendCell(cell);
441
442 return surface_mesh;
443 }
444
445 /*!
446 * \brief Utility function to generate a mesh consisting of the boundary segments of a circle
447 *
448 * \note The caller must delete the mesh
449 */
make_circle_mesh_2d(double radius,int num_segments)450 mint::Mesh* make_circle_mesh_2d(double radius, int num_segments)
451 {
452 constexpr int DIM = 2;
453 using UMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
454
455 UMesh* circle_mesh = new UMesh(DIM, mint::SEGMENT);
456
457 // Generate vertices on the boundary of a circle in CCW order
458 for(int i = 0; i < num_segments; ++i)
459 {
460 double theta = 2. * M_PI * i / static_cast<double>(num_segments);
461 circle_mesh->appendNode(radius * cos(theta), radius * sin(theta));
462
463 int next_idx = (i + 1) == num_segments ? 0 : (i + 1);
464 axom::IndexType cell[2] = {i, next_idx};
465 circle_mesh->appendCell(cell);
466 }
467
468 return circle_mesh;
469 }
470
471 } // end namespace utilities
472 } // end namespace quest
473 } // end namespace axom
474
475 #endif // QUEST_TEST_UTILITIES_HPP_
476