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