1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level COPYRIGHT file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 #include "axom/config.hpp"
7 #include "axom/core.hpp"
8 #include "axom/slic.hpp"
9 #include "axom/mint.hpp"
10 #include "axom/primal.hpp"
11 #include "axom/quest/Discretize.hpp"
12
13 using SphereType = axom::primal::Sphere<double, 3>;
14 using OctType = axom::primal::Octahedron<double, 3>;
15 using Point3D = axom::primal::Point<double, 3>;
16 using Point2D = axom::primal::Point<double, 2>;
17 using NAType = axom::primal::NumericArray<double, 3>;
18
19 // gtest includes
20 #include "gtest/gtest.h"
21
22 #include <vector>
23
24 //------------------------------------------------------------------------------
check_generation(OctType * & standard,OctType * & test,int generation,int offset,int count)25 bool check_generation(OctType*& standard,
26 OctType*& test,
27 int generation,
28 int offset,
29 int count)
30 {
31 std::vector<int> matched(count, 0);
32
33 for(int i = 0; i < count; ++i)
34 {
35 int has_matched = -1;
36 for(int j = 0; j < count; ++j)
37 {
38 if(!matched[j] && standard[offset + i].equals(test[offset + j]))
39 {
40 matched[j] = 1;
41 has_matched = offset + j;
42 }
43 }
44 if(has_matched < 0)
45 {
46 std::cout << "Gen " << generation << " standard oct " << (offset + i)
47 << " didn't match: " << standard[offset + i] << std::endl;
48 }
49 }
50
51 int matchcount = 0;
52 for(int i = 0; i < count; ++i)
53 {
54 matchcount += matched[i];
55 }
56
57 bool matches = (matchcount == count);
58
59 if(!matches)
60 {
61 // a little more reporting here
62 std::cout << "Generation " << generation << " had " << (count - matchcount)
63 << " test octahedra not matched to standard octahedra:"
64 << std::endl;
65 for(int i = 0; i < count; ++i)
66 {
67 if(!matched[i])
68 {
69 std::cout << "Test oct " << (offset + i)
70 << " not matched:" << test[offset + i] << std::endl;
71 }
72 }
73 }
74
75 return matches;
76 }
77
78 /* Return a handwritten list of the octahedra discretizing a one-segment polyline.
79 *
80 * The routine allocates and returns an array of octahedrons pointed to by out.
81 * The caller must free that array.
82 */
discretized_segment(Point2D a,Point2D b,OctType * & out)83 void discretized_segment(Point2D a, Point2D b, OctType*& out)
84 {
85 // We're going to return three generations in the out-vector:
86 // one in the first generation, three in the second (covering each of the
87 // side-walls of the first generation), eight in the third (covering the
88 // two exposed side-walls in each of the second-gen octs).
89 //
90 // These octahedra are actually triangular prisms (if the polyline
91 // is parallel to the X-axis) or truncated tetrahedra. They have two parallel
92 // triangle end-caps and three quadrilateral side-walls. Because we store
93 // them in Octahedron objects, those side-walls are actually two triangles
94 // that are coplanar.
95 constexpr int ZEROTH_GEN_COUNT = 1;
96 constexpr int FIRST_GEN_COUNT = 3;
97 constexpr int SECOND_GEN_COUNT = 6;
98 out = axom::allocate<OctType>(ZEROTH_GEN_COUNT + FIRST_GEN_COUNT +
99 SECOND_GEN_COUNT);
100
101 // The first generation puts a triangle in the end-discs of the truncated
102 // cones with vertices at 12, 4, and 8 o'clock.
103 // The second generation puts three triangles in the end-discs;
104 // the first has vertices at 12, 2, and 4 o'clock,
105 // the second at 4, 6, and 8 o'clock,
106 // the third at 8, 10, and 12 o'clock.
107 // The third generation puts six more triangles in the ends;
108 // first at 12, 1, 2 o'clock,
109 // second at 2, 3, 4 o'clock,
110 // third at 4, 5, 6, o'clock,
111 // fourth at 6, 7, 8 o'clock,
112 // fifth at 8, 9, 10 o'clock,
113 // sixth at 10, 11, 12 o'clock.
114 // These are constructed with 1, 1/2, and sqrt(3)/2 in discs perpendicular
115 // to the X-axis at a[0] and b[0]. We're looking out the X-axis so a
116 // clockwise spin corresponds to the right-hand spin in the function
117 // from_segment() in Discretize.cpp.
118
119 const double SQ_3_2 = sqrt(3.) / 2.;
120
121 // clang-format off
122 Point3D Ia {a[0], -0.5*a[1], SQ_3_2*a[1]};
123 Point3D IIa {a[0], -SQ_3_2*a[1], 0.5*a[1]};
124 Point3D IIIa {a[0], -1*a[1], 0*a[1]};
125 Point3D IVa {a[0], -SQ_3_2*a[1], -0.5*a[1]};
126 Point3D Va {a[0], -0.5*a[1], -SQ_3_2*a[1]};
127 Point3D VIa {a[0], 0*a[1], -1*a[1]};
128 Point3D VIIa {a[0], 0.5*a[1], -SQ_3_2*a[1]};
129 Point3D VIIIa{a[0], SQ_3_2*a[1], -0.5*a[1]};
130 Point3D IXa {a[0], 1*a[1], 0*a[1]};
131 Point3D Xa {a[0], SQ_3_2*a[1], 0.5*a[1]};
132 Point3D XIa {a[0], 0.5*a[1], SQ_3_2*a[1]};
133 Point3D XIIa {a[0], 0*a[1], 1*a[1]};
134
135 Point3D Ib {b[0], -0.5*b[1], SQ_3_2*b[1]};
136 Point3D IIb {b[0], -SQ_3_2*b[1], 0.5*b[1]};
137 Point3D IIIb {b[0], -1*b[1], 0*b[1]};
138 Point3D IVb {b[0], -SQ_3_2*b[1], -0.5*b[1]};
139 Point3D Vb {b[0], -0.5*b[1], -SQ_3_2*b[1]};
140 Point3D VIb {b[0], 0*b[1], -1*b[1]};
141 Point3D VIIb {b[0], 0.5*b[1], -SQ_3_2*b[1]};
142 Point3D VIIIb{b[0], SQ_3_2*b[1], -0.5*b[1]};
143 Point3D IXb {b[0], 1*b[1], 0*b[1]};
144 Point3D Xb {b[0], SQ_3_2*b[1], 0.5*b[1]};
145 Point3D XIb {b[0], 0.5*b[1], SQ_3_2*b[1]};
146 Point3D XIIb {b[0], 0*b[1], 1*b[1]};
147 // clang-format on
148
149 // First generation: a single prism
150 out[0] = OctType(XIIa, IVb, IVa, VIIIb, VIIIa, XIIb);
151
152 // Second generation: a prism on the outside of each side wall.
153 out[1] = OctType(IVb, VIa, VIb, VIIIa, VIIIb, IVa);
154 out[2] = OctType(XIIb, IIa, IIb, IVa, IVb, XIIa);
155 out[3] = OctType(VIIIb, Xa, Xb, XIIa, XIIb, VIIIa);
156
157 // Third generation: a prism on the outside of each exterior side wall.
158 out[4] = OctType(XIIa, Ib, Ia, IIb, IIa, XIIb);
159 out[5] = OctType(IIa, IIIb, IIIa, IVb, IVa, IIb);
160 out[6] = OctType(IVa, Vb, Va, VIb, VIa, IVb);
161 out[7] = OctType(VIa, VIIb, VIIa, VIIIb, VIIIa, VIb);
162 out[8] = OctType(VIIIa, IXb, IXa, Xb, Xa, VIIIb);
163 out[9] = OctType(Xa, XIb, XIa, XIIb, XIIa, Xb);
164 }
165
166 template <typename ExecPolicy>
degenerate_segment_test(const char * label,Point2D * polyline,int len,bool expsuccess)167 void degenerate_segment_test(const char* label,
168 Point2D* polyline,
169 int len,
170 bool expsuccess)
171 {
172 constexpr int gens = 2;
173 int octcount = 0;
174
175 SCOPED_TRACE(label);
176
177 OctType* generated = nullptr;
178 if(expsuccess)
179 {
180 EXPECT_TRUE(
181 axom::quest::discretize<ExecPolicy>(polyline, len, gens, generated, octcount));
182 }
183 else
184 {
185 EXPECT_FALSE(
186 axom::quest::discretize<ExecPolicy>(polyline, len, gens, generated, octcount));
187 }
188 EXPECT_EQ(0, octcount);
189
190 axom::deallocate(generated);
191 }
192
193 //------------------------------------------------------------------------------
194 template <typename ExecPolicy>
run_degen_segment_tests()195 void run_degen_segment_tests()
196 {
197 // Test each of the three generations.
198 // We don't know what order they'll be in, but we do know how many octahedra
199 // will be in each generation.
200
201 Point2D* polyline = axom::allocate<Point2D>(2);
202
203 polyline[0] = {0., 0.};
204 polyline[1] = {0., 0.};
205 degenerate_segment_test<ExecPolicy>("a.x == b.x, a.y == b.y", polyline, 2, true);
206
207 polyline[0] = {1., 0.};
208 polyline[1] = {1., 1.};
209 degenerate_segment_test<ExecPolicy>("a.x == b.x, a.y != b.y", polyline, 2, true);
210
211 polyline[0] = {1., -0.1};
212 polyline[1] = {1.5, 1.};
213 degenerate_segment_test<ExecPolicy>("a.y < 0", polyline, 2, false);
214
215 polyline[0] = {1., -0.1};
216 polyline[1] = {1.5, 1.};
217 degenerate_segment_test<ExecPolicy>("a.y < 0", polyline, 2, false);
218
219 polyline[0] = {1., 1.};
220 polyline[1] = {1.5, -.1};
221 degenerate_segment_test<ExecPolicy>("b.y < 0", polyline, 2, false);
222
223 polyline[0] = {.5, 1.};
224 polyline[1] = {0., 1.};
225 degenerate_segment_test<ExecPolicy>("a.x > b.x", polyline, 2, false);
226
227 axom::deallocate(polyline);
228 }
229
230 template <typename ExecPolicy>
segment_test(const char * label,Point2D * polyline,int len)231 void segment_test(const char* label, Point2D* polyline, int len)
232 {
233 SCOPED_TRACE(label);
234
235 // Test each of the three generations.
236 constexpr int generations = 2;
237
238 // We don't know what order they'll be in, but we do know how many octahedra
239 // will be in each generation.
240 constexpr int ZEROTH_GEN_COUNT = 1;
241 constexpr int FIRST_GEN_COUNT = 3;
242 constexpr int SECOND_GEN_COUNT = 6;
243
244 int generation = 0;
245
246 // The discretized_segment() routine produces a list of 10 hand-calculated
247 // octahedra (three generations) that discretize the surface of revolution
248 // (SoR) produced by revolving a one-segment polyline around the positive
249 // X-axis.
250 OctType* handcut = nullptr;
251 discretized_segment(polyline[0], polyline[1], handcut);
252
253 OctType* generated = nullptr;
254 int octcount = 0;
255 axom::quest::discretize<ExecPolicy>(polyline,
256 len,
257 generations,
258 generated,
259 octcount);
260
261 EXPECT_TRUE(
262 check_generation(handcut, generated, generation, 0, ZEROTH_GEN_COUNT));
263 generation += 1;
264 EXPECT_TRUE(check_generation(handcut,
265 generated,
266 generation,
267 ZEROTH_GEN_COUNT,
268 FIRST_GEN_COUNT));
269 generation += 1;
270 EXPECT_TRUE(check_generation(handcut,
271 generated,
272 generation,
273 ZEROTH_GEN_COUNT + FIRST_GEN_COUNT,
274 SECOND_GEN_COUNT));
275 axom::deallocate(generated);
276 axom::deallocate(handcut);
277 }
278
279 //------------------------------------------------------------------------------
280 template <typename ExecPolicy>
run_single_segment_tests()281 void run_single_segment_tests()
282 {
283 Point2D* polyline = axom::allocate<Point2D>(2);
284
285 polyline[0] = Point2D {0.5, 0.};
286 polyline[1] = Point2D {1.8, 0.8};
287 segment_test<ExecPolicy>("Cone (pointing left)", polyline, 2);
288
289 polyline[0] = Point2D {1.0, 1.0};
290 polyline[1] = Point2D {1.8, 0.};
291 segment_test<ExecPolicy>("Cone (pointing right)", polyline, 2);
292
293 polyline[0] = Point2D {1.0, 1.0};
294 polyline[1] = Point2D {1.8, 0.8};
295 segment_test<ExecPolicy>("Truncated cone", polyline, 2);
296
297 polyline[0] = Point2D {1., 1.};
298 polyline[1] = Point2D {3., 1.};
299 segment_test<ExecPolicy>("Cylinder", polyline, 2);
300
301 polyline[0] = Point2D {-.4, 1.2};
302 polyline[1] = Point2D {1.2, 1.};
303 segment_test<ExecPolicy>("a.x < 0, b.x > 0", polyline, 2);
304
305 axom::deallocate(polyline);
306 }
307
308 template <typename ExecPolicy>
multi_segment_test(const char * label,Point2D * polyline,int len)309 void multi_segment_test(const char* label, Point2D* polyline, int len)
310 {
311 SCOPED_TRACE(label);
312
313 // Test each of the three generations.
314 constexpr int generations = 2;
315
316 // We don't know what order they'll be in, but we do know how many octahedra
317 // will be in each generation.
318 constexpr int ZEROTH_GEN_COUNT = 1;
319 constexpr int FIRST_GEN_COUNT = 3;
320 constexpr int SECOND_GEN_COUNT = 6;
321 constexpr int TOTAL_COUNT =
322 ZEROTH_GEN_COUNT + FIRST_GEN_COUNT + SECOND_GEN_COUNT;
323
324 int generation = 0;
325
326 OctType* generated = nullptr;
327 int octcount = 0;
328 axom::quest::discretize<ExecPolicy>(polyline,
329 len,
330 generations,
331 generated,
332 octcount);
333
334 int segcount = len - 1;
335 OctType* handcut = new OctType[segcount * TOTAL_COUNT];
336
337 for(int segidx = 0; segidx < segcount; ++segidx)
338 {
339 std::stringstream seglabel;
340 seglabel << "Segment " << segidx;
341 SCOPED_TRACE(seglabel.str());
342
343 // The discretized_segment() routine produces a list of 10 hand-calculated
344 // octahedra (three generations) that discretize the surface of revolution
345 // (SoR) produced by revolving a one-segment polyline around the positive
346 // X-axis.
347 OctType* octpointer = &handcut[segidx * TOTAL_COUNT];
348 discretized_segment(polyline[segidx], polyline[segidx + 1], octpointer);
349
350 EXPECT_TRUE(check_generation(octpointer,
351 generated,
352 generation,
353 segidx * TOTAL_COUNT + 0,
354 ZEROTH_GEN_COUNT));
355 generation += 1;
356 EXPECT_TRUE(check_generation(octpointer,
357 generated,
358 generation,
359 segidx * TOTAL_COUNT + ZEROTH_GEN_COUNT,
360 FIRST_GEN_COUNT));
361 generation += 1;
362 EXPECT_TRUE(check_generation(
363 octpointer,
364 generated,
365 generation,
366 segidx * TOTAL_COUNT + ZEROTH_GEN_COUNT + FIRST_GEN_COUNT,
367 SECOND_GEN_COUNT));
368 }
369
370 axom::deallocate(handcut);
371 axom::deallocate(generated);
372 }
373
374 //------------------------------------------------------------------------------
375 template <typename ExecPolicy>
run_multi_segment_tests()376 void run_multi_segment_tests()
377 {
378 constexpr int pointcount = 5;
379 Point2D* polyline = axom::allocate<Point2D>(pointcount);
380
381 polyline[0] = Point2D {1.0, 0.5};
382 polyline[1] = Point2D {1.6, 0.3};
383 polyline[2] = Point2D {1.7, 2.0};
384 polyline[3] = Point2D {3.1, 1.5};
385 polyline[4] = Point2D {4.0, 2.0};
386 segment_test<ExecPolicy>("multi-segment", polyline, pointcount);
387
388 axom::deallocate(polyline);
389 }
390
391 //------------------------------------------------------------------------------
392 enum ReflectDimension
393 {
394 X = 0,
395 Y,
396 Z
397 };
398
399 /* Reflect an octahedron in a given dimension.
400 */
reflect(ReflectDimension d,OctType o)401 OctType reflect(ReflectDimension d, OctType o)
402 {
403 OctType out(o);
404 for(int i = 0; i < OctType::NUM_OCT_VERTS; ++i)
405 {
406 Point3D& pt = out[i];
407 pt[d] *= -1;
408 }
409 return out;
410 }
411
412 /* Return a handwritten list of the octahedra discretizing the unit sphere.
413 *
414 * The routine allocates and returns an array of octahedrons pointed to by out.
415 * The caller must free that array.
416 */
discretized_sphere(OctType * & out)417 void discretized_sphere(OctType*& out)
418 {
419 // We're going to return three generations in the out-vector:
420 // one in the first generation, eight in the second (covering each
421 // of the faces of the first generation), 32 in the third (covering
422 // the four exposed faces in each of the second-gen octs).
423 constexpr int ZEROTH_GEN_COUNT = 1;
424 constexpr int FIRST_GEN_COUNT = 8;
425 constexpr int SECOND_GEN_COUNT = 32;
426 out = axom::allocate<OctType>(ZEROTH_GEN_COUNT + FIRST_GEN_COUNT +
427 SECOND_GEN_COUNT);
428
429 // First generation: one octahedron, with vertices on the unit vectors.
430 NAType ihat({1., 0., 0.});
431 NAType jhat({0., 1., 0.});
432 NAType khat({0., 0., 1.});
433
434 OctType center(Point3D(ihat),
435 Point3D(jhat),
436 Point3D(khat),
437 Point3D(-1 * ihat),
438 Point3D(-1 * jhat),
439 Point3D(-1 * khat));
440 out[0] = center;
441
442 // Second generation: eight octs, one for each face of the unit oct.
443 // Point ij is halfway between (1, 0, 0) and (0, 1, 0),
444 // point jk is halfway between (0, 1, 0) and (0, 0, 1),
445 // point ki is halfway between (0, 0, 1) and (1, 0, 0).
446 Point3D ij(NAType({M_SQRT1_2, M_SQRT1_2, 0.}));
447 Point3D jk(NAType({0., M_SQRT1_2, M_SQRT1_2}));
448 Point3D ki(NAType({M_SQRT1_2, 0., M_SQRT1_2}));
449
450 OctType second_gen(Point3D(ihat), Point3D(jhat), Point3D(khat), jk, ki, ij);
451 out[1] = second_gen;
452 out[2] = reflect(X, second_gen);
453 out[3] = reflect(Y, second_gen);
454 out[4] = reflect(Z, second_gen);
455 out[5] = reflect(Y, reflect(X, second_gen));
456 out[6] = reflect(Z, reflect(Y, second_gen));
457 out[7] = reflect(X, reflect(Z, second_gen));
458 out[8] = reflect(X, reflect(Y, reflect(Z, second_gen)));
459
460 // Third generation: 32 new octs, one for each exposed face of the previous
461 // generation.
462 double SQRT1_6 = 1. / sqrt(6.);
463 // There are three interior points, derived from ij, jk, and ki.
464 // Point a is halfway between ij and ki, at (1/sqrt(6))(2, 1, 1).
465 Point3D a(SQRT1_6 * NAType({2, 1, 1}));
466 // Point b is halfway between ij and jk, at (1/sqrt(6))(1, 2, 1).
467 Point3D b(SQRT1_6 * NAType({1, 2, 1}));
468 // Point c is halfway between jk and ki, at (1/sqrt(6))(1, 1, 2).
469 Point3D c(SQRT1_6 * NAType({1, 1, 2}));
470
471 // There are six edge points, derived from the original corner points and
472 // ij, jk, and ki.
473 // Point d is halfway between ihat and ij, at
474 // (1/sqrt(4 + 2 sqrt(2)))(1+sqrt(2), 1, 0)
475 double FACTOR_3G = 1. / sqrt(2. * M_SQRT2 + 4);
476 Point3D d(FACTOR_3G * NAType({1. + M_SQRT2, 1., 0}));
477 // Point e splits jhat and ij, at (1/sqrt(2 sqrt(2) + 2))(1, 1+sqrt(2), 0)
478 Point3D e(FACTOR_3G * NAType({1., 1. + M_SQRT2, 0}));
479 // Point f splits jhat and jk, at (1/sqrt(2 sqrt(2) + 2))(0, 1+sqrt(2), 1)
480 Point3D f(FACTOR_3G * NAType({0, 1. + M_SQRT2, 1.}));
481 // Point g splits khat and jk, at (1/sqrt(2 sqrt(2) + 2))(0, 1, 1+sqrt(2))
482 Point3D g(FACTOR_3G * NAType({0, 1., 1. + M_SQRT2}));
483 // Point m splits khat and ki, at (1/sqrt(2 sqrt(2) + 2))(1, 0, 1+sqrt(2))
484 Point3D m(FACTOR_3G * NAType({1., 0, 1. + M_SQRT2}));
485 // Point n splits ihat and ki, at (1/sqrt(2 sqrt(2) + 2))(1+sqrt(2), 0, 1)
486 Point3D n(FACTOR_3G * NAType({1. + M_SQRT2, 0, 1.}));
487
488 int offset = ZEROTH_GEN_COUNT + FIRST_GEN_COUNT;
489 // Here's the first octant of third-generation octahedra (octant 0).
490 int octant = 0;
491 // First, the interior oct
492 out[offset + 0 + octant * 4] = OctType(ij, jk, ki, c, a, b);
493 // The one next to ihat
494 out[offset + 1 + octant * 4] = OctType(Point3D(ihat), ij, ki, a, n, d);
495 // The one next to jhat
496 out[offset + 2 + octant * 4] = OctType(Point3D(jhat), jk, ij, b, e, f);
497 // The one next to khat
498 out[offset + 3 + octant * 4] = OctType(Point3D(khat), ki, jk, c, g, m);
499
500 // Now we get to transform these into all seven remaining octants.
501 // Reflect in X
502 octant = 1;
503 ReflectDimension rd0 = X;
504 out[offset + 0 + octant * 4] = reflect(rd0, out[offset + 0]);
505 out[offset + 1 + octant * 4] = reflect(rd0, out[offset + 1]);
506 out[offset + 2 + octant * 4] = reflect(rd0, out[offset + 2]);
507 out[offset + 3 + octant * 4] = reflect(rd0, out[offset + 3]);
508 // Reflect in Y
509 octant = 2;
510 rd0 = Y;
511 out[offset + 0 + octant * 4] = reflect(rd0, out[offset + 0]);
512 out[offset + 1 + octant * 4] = reflect(rd0, out[offset + 1]);
513 out[offset + 2 + octant * 4] = reflect(rd0, out[offset + 2]);
514 out[offset + 3 + octant * 4] = reflect(rd0, out[offset + 3]);
515 // Reflect in Z
516 octant = 3;
517 rd0 = Z;
518 out[offset + 0 + octant * 4] = reflect(rd0, out[offset + 0]);
519 out[offset + 1 + octant * 4] = reflect(rd0, out[offset + 1]);
520 out[offset + 2 + octant * 4] = reflect(rd0, out[offset + 2]);
521 out[offset + 3 + octant * 4] = reflect(rd0, out[offset + 3]);
522 // Reflect in X, then Y
523 octant = 4;
524 rd0 = X;
525 ReflectDimension rd1 = Y;
526 out[offset + 0 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 0]));
527 out[offset + 1 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 1]));
528 out[offset + 2 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 2]));
529 out[offset + 3 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 3]));
530 // Reflect in Y, then Z
531 octant = 5;
532 rd0 = Y;
533 rd1 = Z;
534 out[offset + 0 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 0]));
535 out[offset + 1 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 1]));
536 out[offset + 2 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 2]));
537 out[offset + 3 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 3]));
538 // Reflect in Z, then X
539 octant = 6;
540 rd0 = Z;
541 rd1 = X;
542 out[offset + 0 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 0]));
543 out[offset + 1 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 1]));
544 out[offset + 2 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 2]));
545 out[offset + 3 + octant * 4] = reflect(rd1, reflect(rd0, out[offset + 3]));
546 // And the grand finale: reflect in Z, then Y, then X
547 octant = 7;
548 rd0 = Z;
549 rd1 = Y;
550 ReflectDimension rd2 = X;
551 out[offset + 0 + octant * 4] =
552 reflect(rd2, reflect(rd1, reflect(rd0, out[offset + 0])));
553 out[offset + 1 + octant * 4] =
554 reflect(rd2, reflect(rd1, reflect(rd0, out[offset + 1])));
555 out[offset + 2 + octant * 4] =
556 reflect(rd2, reflect(rd1, reflect(rd0, out[offset + 2])));
557 out[offset + 3 + octant * 4] =
558 reflect(rd2, reflect(rd1, reflect(rd0, out[offset + 3])));
559 }
560
561 //------------------------------------------------------------------------------
TEST(quest_discretize,sphere_test)562 TEST(quest_discretize, sphere_test)
563 {
564 // The discretize() routine chops up a given sphere into the specified
565 // number of generations of octahedra. Here, we'll discretize the unit
566 // sphere in three generations, to match the hand-calculated octs from
567 // discretized_sphere().
568 SphereType sph; // Unit sphere at the origin
569 constexpr int generations = 2;
570 OctType* generated = nullptr;
571 int octcount = 0;
572 axom::quest::discretize(sph, generations, generated, octcount);
573
574 // The discretized_sphere() routine produces a list of 41 hand-calculated
575 // octahedra (three generations) that discretize the unit sphere.
576 OctType* handcut = nullptr;
577 discretized_sphere(handcut);
578
579 // Test each of the three generations.
580 // We don't know what order they'll be in, but we do know how many octahedra
581 // will be in each generation.
582 constexpr int ZEROTH_GEN_COUNT = 1;
583 constexpr int FIRST_GEN_COUNT = 8;
584 constexpr int SECOND_GEN_COUNT = 32;
585 int generation = 0;
586
587 EXPECT_TRUE(
588 check_generation(handcut, generated, generation, 0, ZEROTH_GEN_COUNT));
589 generation += 1;
590 EXPECT_TRUE(check_generation(handcut,
591 generated,
592 generation,
593 ZEROTH_GEN_COUNT,
594 FIRST_GEN_COUNT));
595 generation += 1;
596 EXPECT_TRUE(check_generation(handcut,
597 generated,
598 generation,
599 ZEROTH_GEN_COUNT + FIRST_GEN_COUNT,
600 SECOND_GEN_COUNT));
601
602 axom::deallocate(generated);
603 axom::deallocate(handcut);
604 }
605
606 //------------------------------------------------------------------------------
TEST(quest_discretize,degenerate_sphere_test)607 TEST(quest_discretize, degenerate_sphere_test)
608 {
609 constexpr int generations = 2;
610 OctType* generated = nullptr;
611 int octcount = 0;
612
613 {
614 SCOPED_TRACE("Negative sphere radius");
615 SphereType sph(-.2); // BAD sphere at the origin with negative radius
616 EXPECT_FALSE(axom::quest::discretize(sph, generations, generated, octcount));
617 EXPECT_EQ(0, octcount);
618
619 axom::deallocate(generated);
620 }
621
622 {
623 SCOPED_TRACE("Zero sphere radius");
624 SphereType sph(0.); // Degenerate sphere at the origin with zero radius
625 EXPECT_TRUE(axom::quest::discretize(sph, generations, generated, octcount));
626 EXPECT_EQ(0, octcount);
627
628 axom::deallocate(generated);
629 }
630 }
631
632 //------------------------------------------------------------------------------
TEST(quest_discretize,degenerate_segment_test)633 TEST(quest_discretize, degenerate_segment_test)
634 {
635 SLIC_INFO("Discretizing sequentially");
636 {
637 SCOPED_TRACE("sequential execution");
638 run_degen_segment_tests<axom::SEQ_EXEC>();
639 }
640
641 #if defined(AXOM_USE_OPENMP) && defined(AXOM_USE_RAJA)
642 SLIC_INFO("Discretizing with OpenMP");
643 {
644 SCOPED_TRACE("OpenMP execution");
645 run_degen_segment_tests<axom::OMP_EXEC>();
646 }
647 #endif
648
649 #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_RAJA) && \
650 defined(AXOM_USE_UMPIRE) && defined(__CUDACC__)
651 SLIC_INFO("Discretizing with CUDA");
652 {
653 SCOPED_TRACE("32-wide CUDA execution");
654 run_degen_segment_tests<axom::CUDA_EXEC<32>>();
655 }
656
657 {
658 SCOPED_TRACE("64-wide CUDA execution");
659 run_degen_segment_tests<axom::CUDA_EXEC<64>>();
660 }
661
662 {
663 SCOPED_TRACE("128-wide CUDA execution");
664 run_degen_segment_tests<axom::CUDA_EXEC<128>>();
665 }
666
667 {
668 SCOPED_TRACE("256-wide CUDA execution");
669 run_degen_segment_tests<axom::CUDA_EXEC<256>>();
670 }
671
672 {
673 SCOPED_TRACE("512-wide CUDA execution");
674 run_degen_segment_tests<axom::CUDA_EXEC<512>>();
675 }
676 #endif
677 }
678
679 //------------------------------------------------------------------------------
TEST(quest_discretize,segment_test)680 TEST(quest_discretize, segment_test)
681 {
682 SLIC_INFO("Discretizing sequentially");
683 {
684 SCOPED_TRACE("sequential execution");
685 run_single_segment_tests<axom::SEQ_EXEC>();
686 }
687
688 #if defined(AXOM_USE_OPENMP) && defined(AXOM_USE_RAJA)
689 SLIC_INFO("Discretizing with OpenMP");
690 {
691 SCOPED_TRACE("OpenMP execution");
692 run_single_segment_tests<axom::OMP_EXEC>();
693 }
694 #endif
695
696 #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_RAJA) && \
697 defined(AXOM_USE_UMPIRE) && defined(__CUDACC__)
698 SLIC_INFO("Discretizing with CUDA");
699 {
700 SCOPED_TRACE("32-wide CUDA execution");
701 run_single_segment_tests<axom::CUDA_EXEC<32>>();
702 }
703
704 {
705 SCOPED_TRACE("64-wide CUDA execution");
706 run_single_segment_tests<axom::CUDA_EXEC<64>>();
707 }
708
709 {
710 SCOPED_TRACE("128-wide CUDA execution");
711 run_single_segment_tests<axom::CUDA_EXEC<128>>();
712 }
713
714 {
715 SCOPED_TRACE("256-wide CUDA execution");
716 run_single_segment_tests<axom::CUDA_EXEC<256>>();
717 }
718
719 {
720 SCOPED_TRACE("512-wide CUDA execution");
721 run_single_segment_tests<axom::CUDA_EXEC<512>>();
722 }
723 #endif
724 }
725
726 //------------------------------------------------------------------------------
TEST(quest_discretize,multi_segment_test)727 TEST(quest_discretize, multi_segment_test)
728 {
729 SLIC_INFO("Discretizing sequentially");
730 {
731 SCOPED_TRACE("sequential execution");
732 run_multi_segment_tests<axom::SEQ_EXEC>();
733 }
734
735 #if defined(AXOM_USE_OPENMP) && defined(AXOM_USE_RAJA)
736 SLIC_INFO("Discretizing with OpenMP");
737 {
738 SCOPED_TRACE("OpenMP execution");
739 run_multi_segment_tests<axom::OMP_EXEC>();
740 }
741 #endif
742
743 #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_RAJA) && \
744 defined(AXOM_USE_UMPIRE) && defined(__CUDACC__)
745 SLIC_INFO("Discretizing with CUDA");
746 {
747 SCOPED_TRACE("32-wide CUDA execution");
748 run_multi_segment_tests<axom::CUDA_EXEC<32>>();
749 }
750
751 {
752 SCOPED_TRACE("64-wide CUDA execution");
753 run_multi_segment_tests<axom::CUDA_EXEC<64>>();
754 }
755
756 {
757 SCOPED_TRACE("128-wide CUDA execution");
758 run_multi_segment_tests<axom::CUDA_EXEC<128>>();
759 }
760
761 {
762 SCOPED_TRACE("256-wide CUDA execution");
763 run_multi_segment_tests<axom::CUDA_EXEC<256>>();
764 }
765
766 {
767 SCOPED_TRACE("512-wide CUDA execution");
768 run_multi_segment_tests<axom::CUDA_EXEC<512>>();
769 }
770 #endif
771 }
772
773 //------------------------------------------------------------------------------
TEST(quest_discretize,to_tet_mesh)774 TEST(quest_discretize, to_tet_mesh)
775 {
776 constexpr int pointcount = 5;
777 constexpr int segcount = pointcount - 1;
778 constexpr int generations = 2;
779 Point2D* polyline = axom::allocate<Point2D>(pointcount);
780
781 polyline[0] = Point2D {1.0, 0.5};
782 polyline[1] = Point2D {1.6, 0.3};
783 polyline[2] = Point2D {1.7, 2.0};
784 polyline[3] = Point2D {3.1, 1.5};
785 polyline[4] = Point2D {4.0, 2.0};
786
787 OctType* generated = nullptr;
788 int octcount = 0;
789 axom::quest::discretize<axom::SEQ_EXEC>(polyline,
790 pointcount,
791 generations,
792 generated,
793 octcount);
794
795 axom::mint::Mesh* mesh;
796 axom::quest::mesh_from_discretized_polyline(generated, octcount, segcount, mesh);
797 axom::mint::write_vtk(mesh, "tet_mesh.vtk");
798
799 delete mesh;
800 axom::deallocate(polyline);
801 }
802
803 //------------------------------------------------------------------------------
804 //------------------------------------------------------------------------------
main(int argc,char * argv[])805 int main(int argc, char* argv[])
806 {
807 int result = 0;
808
809 ::testing::InitGoogleTest(&argc, argv);
810
811 axom::slic::SimpleLogger logger;
812 result = RUN_ALL_TESTS();
813
814 return result;
815 }
816