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 #include "gtest/gtest.h"
7
8 #include "axom/config.hpp"
9
10 #include "axom/core/Types.hpp"
11 #include "axom/core/execution/for_all.hpp" // for for_all, execution policies
12 #include "axom/core/memory_management.hpp"
13
14 #include "axom/primal/geometry/Point.hpp"
15 #include "axom/primal/geometry/BoundingBox.hpp"
16 #include "axom/primal/geometry/Triangle.hpp"
17 #include "axom/primal/geometry/Plane.hpp"
18
19 #include "axom/primal/operators/clip.hpp"
20
21 #include <limits>
22
23 namespace Primal3D
24 {
25 using PointType = axom::primal::Point<double, 3>;
26 using VectorType = axom::primal::Vector<double, 3>;
27 using BoundingBoxType = axom::primal::BoundingBox<double, 3>;
28 using TriangleType = axom::primal::Triangle<double, 3>;
29 using TetrahedronType = axom::primal::Tetrahedron<double, 3>;
30 using OctahedronType = axom::primal::Octahedron<double, 3>;
31 using PolyhedronType = axom::primal::Polyhedron<double, 3>;
32 using PolygonType = axom::primal::Polygon<double, 3>;
33 using PlaneType = axom::primal::Plane<double, 3>;
34 using PolyhedronType = axom::primal::Polyhedron<double, 3>;
35 } // namespace Primal3D
36
TEST(primal_clip,simple_clip)37 TEST(primal_clip, simple_clip)
38 {
39 using namespace Primal3D;
40 BoundingBoxType bbox;
41 bbox.addPoint(PointType::zero());
42 bbox.addPoint(PointType::ones());
43
44 PointType points[] = {
45 PointType::make_point(2, 2, 2),
46 PointType::make_point(2, 2, 4),
47 PointType::make_point(2, 4, 2),
48 PointType::make_point(-100, -100, 0.5),
49 PointType::make_point(-100, 100, 0.5),
50 PointType::make_point(100, 0, 0.5),
51 PointType::make_point(0.25, 0.25, 0.5),
52 PointType::make_point(0.75, 0.25, 0.5),
53 PointType::make_point(0.66, 0.5, 0.5),
54 PointType::make_point(1.5, 0.5, 0.5),
55 };
56
57 {
58 TriangleType tri(points[0], points[1], points[2]);
59
60 PolygonType poly = axom::primal::clip(tri, bbox);
61 EXPECT_EQ(0, poly.numVertices());
62 }
63
64 {
65 TriangleType tri(points[3], points[4], points[5]);
66
67 PolygonType poly = axom::primal::clip(tri, bbox);
68 EXPECT_EQ(4, poly.numVertices());
69
70 SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox
71 << " is polygon" << poly);
72 }
73
74 {
75 TriangleType tri(points[3], points[4], points[5]);
76
77 PolygonType poly = axom::primal::clip(tri, bbox);
78 EXPECT_EQ(4, poly.numVertices());
79
80 EXPECT_EQ(PointType(.5), poly.centroid());
81
82 SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox
83 << " is polygon" << poly);
84 }
85
86 {
87 TriangleType tri(points[6], points[7], points[9]);
88
89 PolygonType poly = axom::primal::clip(tri, bbox);
90 EXPECT_EQ(4, poly.numVertices());
91
92 SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox
93 << " is polygon" << poly);
94 }
95 }
96
TEST(primal_clip,unit_simplex)97 TEST(primal_clip, unit_simplex)
98 {
99 using namespace Primal3D;
100 double delta = 1e-5;
101
102 // Test the "unit simplex", and a jittered version
103 PointType points[] = {PointType::make_point(1, 0, 0),
104 PointType::make_point(0, 1, 0),
105 PointType::make_point(0, 0, 1),
106 PointType::make_point(1 + delta, delta, delta),
107 PointType::make_point(delta, 1 + delta, delta),
108 PointType::make_point(delta, delta, 1 + delta)};
109
110 BoundingBoxType bbox;
111 bbox.addPoint(PointType::zero());
112 bbox.addPoint(PointType(.75));
113
114 // intersection of this triangle and cube is a hexagon
115 {
116 TriangleType tri(points[0], points[1], points[2]);
117
118 PolygonType poly = axom::primal::clip(tri, bbox);
119 EXPECT_EQ(6, poly.numVertices());
120
121 SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox
122 << " is polygon" << poly);
123 }
124 {
125 TriangleType tri(points[3], points[4], points[5]);
126
127 PolygonType poly = axom::primal::clip(tri, bbox);
128 EXPECT_EQ(6, poly.numVertices());
129
130 SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox
131 << " is polygon" << poly);
132 }
133 }
134
TEST(primal_clip,boundingBoxOptimization)135 TEST(primal_clip, boundingBoxOptimization)
136 {
137 using namespace Primal3D;
138
139 SLIC_INFO("Checking correctness of optimization for skipping clipping "
140 << " of planes that the triangle's bounding box doesn't cover");
141
142 const double VAL1 = 3.;
143 const double VAL2 = 2.;
144
145 BoundingBoxType bbox;
146 bbox.addPoint(PointType(-1.));
147 bbox.addPoint(PointType(1.));
148
149 PointType midpoint = PointType::zero();
150
151 PointType points[] = {
152 PointType::make_point(VAL1, VAL2, 0),
153 PointType::make_point(-VAL1, VAL2, 0),
154 PointType::make_point(VAL1, -VAL2, 0),
155 PointType::make_point(-VAL1, -VAL2, 0),
156
157 PointType::make_point(VAL1, 0, VAL2),
158 PointType::make_point(-VAL1, 0, VAL2),
159 PointType::make_point(VAL1, 0, -VAL2),
160 PointType::make_point(-VAL1, 0, -VAL2),
161
162 PointType::make_point(0, VAL2, VAL1),
163 PointType::make_point(0, VAL2, -VAL1),
164 PointType::make_point(0, -VAL2, VAL1),
165 PointType::make_point(0, -VAL2, -VAL1),
166
167 PointType::make_point(0, VAL1, VAL2),
168 PointType::make_point(0, -VAL1, VAL2),
169 PointType::make_point(0, VAL1, -VAL2),
170 PointType::make_point(0, -VAL1, -VAL2),
171 };
172
173 for(int i = 0; i < 16; i += 2)
174 {
175 TriangleType tri(midpoint, points[i], points[i + 1]);
176 PolygonType poly = axom::primal::clip(tri, bbox);
177 SLIC_INFO(poly);
178 EXPECT_EQ(5, poly.numVertices());
179 }
180 }
181
TEST(primal_clip,experimentalData)182 TEST(primal_clip, experimentalData)
183 {
184 using namespace Primal3D;
185
186 const double EPS = 1e-8;
187
188 // Triangle 248 from sphere mesh
189 TriangleType tri(PointType::make_point(0.405431, 3.91921, 3.07821),
190 PointType::make_point(1.06511, 3.96325, 2.85626),
191 PointType::make_point(0.656002, 4.32465, 2.42221));
192
193 // Block index {grid pt: (19,29,24); level: 5} from InOutOctree
194 BoundingBoxType box12(PointType::make_point(0.937594, 4.06291, 2.50025),
195 PointType::make_point(1.25012, 4.37544, 2.81278));
196
197 PolygonType poly = axom::primal::clip(tri, box12);
198 EXPECT_EQ(3, poly.numVertices());
199
200 SLIC_INFO("Intersection of triangle "
201 << tri << " \n\t and bounding box " << box12 << " \n\t is polygon"
202 << poly << " with centroid " << poly.centroid());
203
204 // Check that the polygon vertices are on the triangle
205 for(int i = 0; i < poly.numVertices(); ++i)
206 {
207 PointType bary = tri.physToBarycentric(poly[i]);
208 PointType reconstructed = tri.baryToPhysical(bary);
209
210 SLIC_INFO("Testing clipped polygon point "
211 << poly[i] << "-- w/ barycentric coords " << bary
212 << "\n\t-- reconstructed point is: " << reconstructed << "...\n");
213
214 double barySum = bary[0] + bary[1] + bary[2];
215 EXPECT_NEAR(1., barySum, EPS);
216
217 for(int dim = 0; dim < 3; ++dim)
218 {
219 EXPECT_GE(bary[dim], -EPS);
220 EXPECT_NEAR(poly[i][dim], reconstructed[dim], EPS);
221 }
222
223 EXPECT_TRUE(box12.contains(poly[i]));
224 }
225
226 // Check that the polygon centroid is on the triangle
227 {
228 PointType centroid = poly.centroid();
229 PointType bary = tri.physToBarycentric(centroid);
230 PointType reconstructed = tri.baryToPhysical(bary);
231
232 SLIC_INFO("Testing clipped polygon centroid "
233 << centroid << "-- w/ barycentric coords " << bary
234 << "\n\t-- reconstructed point is: " << reconstructed << "...\n");
235
236 double barySum = bary[0] + bary[1] + bary[2];
237 EXPECT_NEAR(1., barySum, EPS);
238
239 for(int dim = 0; dim < 3; ++dim)
240 {
241 EXPECT_GE(bary[dim], -EPS);
242 EXPECT_NEAR(centroid[dim], reconstructed[dim], EPS);
243 }
244 EXPECT_TRUE(box12.contains(centroid));
245 }
246 }
247
248 template <typename ExecPolicy>
unit_check_poly_clip()249 void unit_check_poly_clip()
250 {
251 using namespace Primal3D;
252
253 constexpr double EPS = 1.e-24;
254
255 // Create a square in 3D space as our "polyhedron", and a plane splitting it
256 // between the middle. The polyhedron after a call to poly_clip_vertices
257 // should have two new vertices where the plane crosses the 0-1 and 2-3 edges:
258 // |
259 // 3 ---|--- 2 3 -- 2' -- 2
260 // | | | | |
261 // | p-> | --> | |
262 // | | | | |
263 // 0 ---|--- 1 0 -- 1' -- 1
264 // |
265 // In addition, vertices 0 and 3 should be marked as clipped.
266
267 PolyhedronType square;
268 square.addVertex({0.0, 0.0, 0.0});
269 square.addVertex({1.0, 0.0, 0.0});
270 square.addVertex({1.0, 1.0, 0.0});
271 square.addVertex({0.0, 1.0, 0.0});
272
273 square.addNeighbors(0, {1, 3});
274 square.addNeighbors(1, {0, 2});
275 square.addNeighbors(2, {1, 3});
276 square.addNeighbors(3, {0, 2});
277
278 PlaneType plane(VectorType {1, 0, 0}, PointType {0.5, 0.0, 0.0});
279
280 const int current_allocator = axom::getDefaultAllocatorID();
281 axom::setDefaultAllocator(axom::execution_space<ExecPolicy>::allocatorID());
282
283 PolyhedronType* out_square = axom::allocate<PolyhedronType>(1);
284 out_square[0] = square;
285
286 unsigned int* out_clipped = axom::allocate<unsigned int>(1);
287 out_clipped[0] = 0;
288
289 axom::for_all<ExecPolicy>(
290 0,
291 1,
292 AXOM_LAMBDA(int /* idx */) {
293 axom::primal::detail::poly_clip_vertices(out_square[0],
294 plane,
295 EPS,
296 out_clipped[0]);
297 });
298
299 const PolyhedronType& clippedSquare = out_square[0];
300
301 EXPECT_EQ(clippedSquare.numVertices(), 6);
302 // Check that existing vertices were not modified
303 EXPECT_EQ(clippedSquare[0], square[0]);
304 EXPECT_EQ(clippedSquare[1], square[1]);
305 EXPECT_EQ(clippedSquare[2], square[2]);
306 EXPECT_EQ(clippedSquare[3], square[3]);
307
308 int lo_idx = -1, hi_idx = -1;
309 for(int i = 4; i < 6; i++)
310 {
311 if(clippedSquare[i] == PointType {0.5, 0.0, 0.0})
312 {
313 lo_idx = i;
314 }
315 else if(clippedSquare[i] == PointType {0.5, 1.0, 0.0})
316 {
317 hi_idx = i;
318 }
319 }
320
321 // Check that plane-crossing vertices were generated
322 EXPECT_NE(lo_idx, -1);
323 EXPECT_NE(hi_idx, -1);
324
325 // Check that vertices outside the plane were marked as clipped
326 EXPECT_EQ(out_clipped[0], (1 | (1 << 3)));
327
328 // Generate sets of expected neighbors
329 std::vector<std::set<int>> expectedNbrs(6);
330 expectedNbrs[0] = {3, lo_idx};
331 expectedNbrs[1] = {2, lo_idx};
332 expectedNbrs[2] = {1, hi_idx};
333 expectedNbrs[3] = {0, hi_idx};
334 expectedNbrs[lo_idx] = {0, 1};
335 expectedNbrs[hi_idx] = {2, 3};
336
337 // Check connectivity
338 for(int iv = 0; iv < 6; iv++)
339 {
340 // Each vertex should only have two neighbors
341 EXPECT_EQ(clippedSquare.getNumNeighbors(iv), 2);
342 std::set<int> resultNbrs {clippedSquare.getNeighbors(iv)[0],
343 clippedSquare.getNeighbors(iv)[1]};
344 // Check that neighbors match expected connectivity
345 EXPECT_EQ(expectedNbrs[iv], resultNbrs);
346 }
347
348 axom::deallocate(out_square);
349 axom::deallocate(out_clipped);
350 axom::setDefaultAllocator(current_allocator);
351 }
352
353 #if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE)
354 template <typename ExecPolicy>
check_oct_tet_clip(double EPS)355 void check_oct_tet_clip(double EPS)
356 {
357 using namespace Primal3D;
358
359 umpire::ResourceManager& rm = umpire::ResourceManager::getInstance();
360
361 // Save current/default allocator
362 const int current_allocator = axom::getDefaultAllocatorID();
363
364 // Determine new allocator (for CUDA policy, set to Unified)
365 umpire::Allocator allocator =
366 rm.getAllocator(axom::execution_space<ExecPolicy>::allocatorID());
367
368 // Set new default to device
369 axom::setDefaultAllocator(allocator.getId());
370
371 // Allocate memory for shapes
372 TetrahedronType* tet = axom::allocate<TetrahedronType>(1);
373 OctahedronType* oct = axom::allocate<OctahedronType>(1);
374
375 PolyhedronType* res = (axom::execution_space<ExecPolicy>::onDevice()
376 ? axom::allocate<PolyhedronType>(
377 1,
378 rm.getAllocator(umpire::resource::Unified).getId())
379 : axom::allocate<PolyhedronType>(1));
380
381 tet[0] = TetrahedronType(PointType({1, 0, 0}),
382 PointType({1, 1, 0}),
383 PointType({0, 1, 0}),
384 PointType({1, 0, 1}));
385
386 oct[0] = OctahedronType(PointType({1, 0, 0}),
387 PointType({1, 1, 0}),
388 PointType({0, 1, 0}),
389 PointType({0, 1, 1}),
390 PointType({0, 0, 1}),
391 PointType({1, 0, 1}));
392
393 axom::for_all<ExecPolicy>(
394 1,
395 AXOM_LAMBDA(int i) { res[i] = axom::primal::clip(oct[i], tet[i]); });
396
397 EXPECT_NEAR(0.1666, res[0].volume(), EPS);
398
399 axom::deallocate(tet);
400 axom::deallocate(oct);
401 axom::deallocate(res);
402
403 axom::setDefaultAllocator(current_allocator);
404 }
405
TEST(primal_clip,unit_poly_clip_vertices_sequential)406 TEST(primal_clip, unit_poly_clip_vertices_sequential)
407 {
408 unit_check_poly_clip<axom::SEQ_EXEC>();
409 }
410
TEST(primal_clip,clip_oct_tet_sequential)411 TEST(primal_clip, clip_oct_tet_sequential)
412 {
413 const double EPS = 1e-4;
414 check_oct_tet_clip<axom::SEQ_EXEC>(EPS);
415 }
416
417 #ifdef AXOM_USE_OPENMP
TEST(primal_clip,unit_poly_clip_vertices_omp)418 TEST(primal_clip, unit_poly_clip_vertices_omp)
419 {
420 unit_check_poly_clip<axom::OMP_EXEC>();
421 }
422
TEST(primal_clip,clip_oct_tet_omp)423 TEST(primal_clip, clip_oct_tet_omp)
424 {
425 const double EPS = 1e-4;
426 check_oct_tet_clip<axom::OMP_EXEC>(EPS);
427 }
428 #endif /* AXOM_USE_OPENMP */
429
430 #if defined(AXOM_USE_CUDA)
AXOM_CUDA_TEST(primal_clip,unit_poly_clip_vertices_gpu)431 AXOM_CUDA_TEST(primal_clip, unit_poly_clip_vertices_gpu)
432 {
433 unit_check_poly_clip<axom::CUDA_EXEC<256>>();
434 }
435
AXOM_CUDA_TEST(primal_clip,clip_oct_tet_cuda)436 AXOM_CUDA_TEST(primal_clip, clip_oct_tet_cuda)
437 {
438 const double EPS = 1e-4;
439 check_oct_tet_clip<axom::CUDA_EXEC<256>>(EPS);
440 }
441 #endif /* AXOM_USE_CUDA */
442 #endif /* AXOM_USE_RAJA && AXOM_USE_UMPIRE */
443
444 // Tetrahedron does not clip octahedron.
TEST(primal_clip,oct_tet_clip_nonintersect)445 TEST(primal_clip, oct_tet_clip_nonintersect)
446 {
447 using namespace Primal3D;
448
449 TetrahedronType tet(PointType({-1, -1, -1}),
450 PointType({-1, 0, 0}),
451 PointType({-1, -1, 0}),
452 PointType({0, 0, 0}));
453 OctahedronType oct(PointType({1, 0, 0}),
454 PointType({1, 1, 0}),
455 PointType({0, 1, 0}),
456 PointType({0, 1, 1}),
457 PointType({0, 0, 1}),
458 PointType({1, 0, 1}));
459
460 PolyhedronType poly = axom::primal::clip(oct, tet);
461 EXPECT_EQ(0.0, poly.volume());
462 }
463
464 // Tetrahedron is encapsulated by the octahedron
TEST(primal_clip,oct_tet_clip_encapsulate)465 TEST(primal_clip, oct_tet_clip_encapsulate)
466 {
467 using namespace Primal3D;
468 const double EPS = 1e-4;
469
470 TetrahedronType tet(PointType({1, 0, 0}),
471 PointType({1, 1, 0}),
472 PointType({0, 1, 0}),
473 PointType({1, 0, 1}));
474 OctahedronType oct(PointType({1, 0, 0}),
475 PointType({1, 1, 0}),
476 PointType({0, 1, 0}),
477 PointType({0, 1, 1}),
478 PointType({0, 0, 1}),
479 PointType({1, 0, 1}));
480
481 PolyhedronType poly = axom::primal::clip(oct, tet);
482
483 // Expected result should be 0.666 / 4 = 0.1666, volume of tet.
484 EXPECT_NEAR(0.1666, poly.volume(), EPS);
485 }
486
487 // Octahedron is encapsulated inside the tetrahedron
TEST(primal_clip,oct_tet_clip_encapsulate_inv)488 TEST(primal_clip, oct_tet_clip_encapsulate_inv)
489 {
490 using namespace Primal3D;
491 const double EPS = 1e-4;
492
493 TetrahedronType tet(PointType({0, 0, 0}),
494 PointType({0, 2, 0}),
495 PointType({0, 0, 2}),
496 PointType({2, 0, 0}));
497 OctahedronType oct(PointType({1, 0, 0}),
498 PointType({1, 1, 0}),
499 PointType({0, 1, 0}),
500 PointType({0, 1, 1}),
501 PointType({0, 0, 1}),
502 PointType({1, 0, 1}));
503
504 PolyhedronType poly = axom::primal::clip(oct, tet);
505
506 // Expected result should be 0.6666, volume of oct.
507 EXPECT_NEAR(0.6666, poly.volume(), EPS);
508 }
509
510 // Half of the octahedron is clipped by the tetrahedron
TEST(primal_clip,oct_tet_clip_half)511 TEST(primal_clip, oct_tet_clip_half)
512 {
513 using namespace Primal3D;
514 const double EPS = 1e-4;
515
516 TetrahedronType tet(PointType({0.5, 0.5, 2}),
517 PointType({2, -1, 0}),
518 PointType({-1, -1, 0}),
519 PointType({-1, 2, 0}));
520 OctahedronType oct(PointType({1, 0, 0}),
521 PointType({1, 1, 0}),
522 PointType({0, 1, 0}),
523 PointType({0, 1, 1}),
524 PointType({0, 0, 1}),
525 PointType({1, 0, 1}));
526
527 PolyhedronType poly = axom::primal::clip(oct, tet);
528
529 // Expected result should be 0.3333, half the volume of oct.
530 EXPECT_NEAR(0.3333, poly.volume(), EPS);
531 }
532
533 // Octahedron is adjacent to tetrahedron
TEST(primal_clip,oct_tet_clip_adjacent)534 TEST(primal_clip, oct_tet_clip_adjacent)
535 {
536 using namespace Primal3D;
537
538 TetrahedronType tet(PointType({0, -1, 0}),
539 PointType({0, 0, 1}),
540 PointType({1, 0, 1}),
541 PointType({1, 0, 0}));
542 OctahedronType oct(PointType({1, 0, 0}),
543 PointType({1, 1, 0}),
544 PointType({0, 1, 0}),
545 PointType({0, 1, 1}),
546 PointType({0, 0, 1}),
547 PointType({1, 0, 1}));
548
549 PolyhedronType poly = axom::primal::clip(oct, tet);
550
551 EXPECT_EQ(0.0, poly.volume());
552 }
553
554 // Tetrahedron clips octahedron at a single vertex
TEST(primal_clip,oct_tet_clip_point)555 TEST(primal_clip, oct_tet_clip_point)
556 {
557 using namespace Primal3D;
558
559 TetrahedronType tet(PointType({-1, -1, 0}),
560 PointType({-0.5, 0.5, 0}),
561 PointType({0, 0, 2}),
562 PointType({0.5, -0.5, 0}));
563 OctahedronType oct(PointType({1, 0, 0}),
564 PointType({1, 1, 0}),
565 PointType({0, 1, 0}),
566 PointType({0, 1, 1}),
567 PointType({0, 0, 1}),
568 PointType({1, 0, 1}));
569
570 PolyhedronType poly = axom::primal::clip(oct, tet);
571
572 EXPECT_EQ(0.0, poly.volume());
573 }
574
TEST(primal_clip,oct_tet_clip_special_case_1)575 TEST(primal_clip, oct_tet_clip_special_case_1)
576 {
577 using namespace Primal3D;
578 const double EPS = 1e-4;
579
580 TetrahedronType tet(PointType({0.5, 0.5, 0.5}),
581 PointType({1, 1, 0}),
582 PointType({1, 0, 0}),
583 PointType({0.5, 0.5, 0}));
584
585 OctahedronType oct(PointType({0.5, 0.853553, 0.146447}),
586 PointType({0.853553, 0.853553, 0.5}),
587 PointType({0.853553, 0.5, 0.146447}),
588 PointType({1, 0.5, 0.5}),
589 PointType({0.5, 0.5, 0}),
590 PointType({0.5, 1, 0.5}));
591
592 // NOTE: Order of vertices 1,2 and 4,5 are flipped due
593 // to winding being opposite of what volume() expects
594 PolyhedronType octPoly;
595 octPoly.addVertex(PointType({0.5, 0.853553, 0.146447}));
596 octPoly.addVertex(PointType({0.853553, 0.5, 0.146447}));
597 octPoly.addVertex(PointType({0.853553, 0.853553, 0.5}));
598 octPoly.addVertex(PointType({1, 0.5, 0.5}));
599 octPoly.addVertex(PointType({0.5, 1, 0.5}));
600 octPoly.addVertex(PointType({0.5, 0.5, 0}));
601
602 octPoly.addNeighbors(0, {1, 5, 4, 2});
603 octPoly.addNeighbors(1, {0, 2, 3, 5});
604 octPoly.addNeighbors(2, {0, 4, 3, 1});
605 octPoly.addNeighbors(3, {1, 2, 4, 5});
606 octPoly.addNeighbors(4, {0, 5, 3, 2});
607 octPoly.addNeighbors(5, {0, 1, 3, 4});
608
609 EXPECT_NEAR(0.0251, octPoly.volume(), EPS);
610
611 PolyhedronType poly = axom::primal::clip(oct, tet);
612
613 EXPECT_NEAR(0.0041, poly.volume(), EPS);
614 }
615
TEST(primal_clip,oct_tet_clip_special_case_2)616 TEST(primal_clip, oct_tet_clip_special_case_2)
617 {
618 using namespace Primal3D;
619 const double EPS = 1e-4;
620
621 TetrahedronType tet(PointType({0.5, 0.5, 0.5}),
622 PointType({0, 1, 0}),
623 PointType({1, 1, 0}),
624 PointType({0.5, 0.5, 0}));
625
626 OctahedronType oct(PointType({0.5, 0.853553, 0.146447}),
627 PointType({0.853553, 0.853553, 0.5}),
628 PointType({0.853553, 0.5, 0.146447}),
629 PointType({1, 0.5, 0.5}),
630 PointType({0.5, 0.5, 0}),
631 PointType({0.5, 1, 0.5}));
632
633 // NOTE: Order of vertices 1,2 and 4,5 are flipped due
634 // to winding being opposite of what volume() expects
635 PolyhedronType octPoly;
636 octPoly.addVertex(PointType({0.5, 0.853553, 0.146447}));
637 octPoly.addVertex(PointType({0.853553, 0.5, 0.146447}));
638 octPoly.addVertex(PointType({0.853553, 0.853553, 0.5}));
639 octPoly.addVertex(PointType({1, 0.5, 0.5}));
640 octPoly.addVertex(PointType({0.5, 1, 0.5}));
641 octPoly.addVertex(PointType({0.5, 0.5, 0}));
642
643 octPoly.addNeighbors(0, {1, 5, 4, 2});
644 octPoly.addNeighbors(1, {0, 2, 3, 5});
645 octPoly.addNeighbors(2, {0, 4, 3, 1});
646 octPoly.addNeighbors(3, {1, 2, 4, 5});
647 octPoly.addNeighbors(4, {0, 5, 3, 2});
648 octPoly.addNeighbors(5, {0, 1, 3, 4});
649
650 EXPECT_NEAR(0.0251, octPoly.volume(), EPS);
651
652 PolyhedronType poly = axom::primal::clip(oct, tet);
653
654 EXPECT_NEAR(0.0041, poly.volume(), EPS);
655 }
656
657 //------------------------------------------------------------------------------
658 #include "axom/slic/core/SimpleLogger.hpp"
659 using axom::slic::SimpleLogger;
660
main(int argc,char * argv[])661 int main(int argc, char* argv[])
662 {
663 ::testing::InitGoogleTest(&argc, argv);
664
665 SimpleLogger logger; // create & initialize test logger,
666
667 int result = RUN_ALL_TESTS();
668
669 return result;
670 }
671