SOR implementation for MeshClipper

src/axom/quest/CMakeLists.txt

@@ -150,11 +150,22 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE)
         if(RAJA_FOUND)
             list(APPEND quest_headers MeshClipper.hpp
                                       MeshClipperStrategy.hpp
+                                      detail/clipping/Plane3DClipper.hpp
                                       detail/clipping/TetClipper.hpp
+                                      detail/clipping/HexClipper.hpp
+                                      detail/clipping/SphereClipper.hpp
+                                      detail/clipping/FSorClipper.hpp
+                                      detail/clipping/SorClipper.hpp
                                       detail/clipping/MeshClipperImpl.hpp)
             list(APPEND quest_sources MeshClipper.cpp
+                                      MeshClipperStrategy.cpp
+                                      detail/clipping/Plane3DClipper.cpp
                                       detail/clipping/TetClipper.cpp
-                                      MeshClipperStrategy.cpp)
+                                      detail/clipping/HexClipper.cpp
+                                      detail/clipping/SphereClipper.cpp
+                                      detail/clipping/FSorClipper.cpp
+                                      detail/clipping/SorClipper.cpp
+                                      detail/clipping/MeshClipperImpl.hpp)
         endif()
     endif()
 

src/axom/quest/IntersectionShaper.hpp

@@ -825,6 +825,7 @@ class IntersectionShaper : public Shaper
 
       // Generate the Octahedra
       // (octahedra m_octs will be on device)
+      m_octs = axom::Array<OctahedronType>(0, 0, axom::execution_space<ExecSpace>::allocatorID());
       const bool disc_status =
         axom::quest::discretize<ExecSpace>(polyline, polyline_size, m_level, m_octs, m_octcount);
 

src/axom/quest/detail/Discretize_detail.hpp

@@ -7,6 +7,7 @@
 #define AXOM_QUEST_DISCRETIZE_DETAIL_
 
 #include "axom/primal/constants.hpp"
+#include "math.h"
 
 namespace
 {
@@ -152,12 +153,11 @@ int discrSeg(const Point2D &a, const Point2D &b, int levels, axom::ArrayView<Oct
   int hostAllocID = axom::execution_space<axom::SEQ_EXEC>::allocatorID();
 
   // Assert input assumptions
-  SLIC_ASSERT(b[0] - a[0] >= 0);
   SLIC_ASSERT(a[1] >= 0);
   SLIC_ASSERT(b[1] >= 0);
 
   // Deal with degenerate segments
-  if(b[0] - a[0] < axom::primal::PRIMAL_TINY)
+  if(fabs(b[0] - a[0]) < axom::primal::PRIMAL_TINY)
   {
     return 0;
   }
@@ -256,7 +256,7 @@ namespace quest
  * less than the polyline's length).  That is exponential growth.  Use
  * appropriate caution.
  *
- * This routine initializes an Array pointed to by \a out.
+ * This routine resizes and populates an Array pointed to by \a out.
  */
 template <typename ExecSpace>
 bool discretize(const axom::ArrayView<Point2D> &polyline,
@@ -265,19 +265,18 @@ bool discretize(const axom::ArrayView<Point2D> &polyline,
                 axom::Array<OctType> &out,
                 int &octcount)
 {
-  int allocId = axom::execution_space<ExecSpace>::allocatorID();
+  SLIC_ERROR_IF(!axom::execution_space<ExecSpace>::usesAllocId(out.getAllocatorID()),
+                axom::fmt::format("Execution space {} cannot access allocator id {}",
+                                  axom::execution_space<ExecSpace>::name(),
+                                  out.getAllocatorID()));
+
   // Check for invalid input.  If any segment is invalid, exit returning false.
   bool stillValid = true;
   int segmentcount = pointcount - 1;
   for(int seg = 0; seg < segmentcount && stillValid; ++seg)
   {
     const Point2D &a = polyline[seg];
     const Point2D &b = polyline[seg + 1];
-    // invalid if a.x > b.x
-    if(a[0] > b[0])
-    {
-      stillValid = false;
-    }
     if(a[1] < 0 || b[1] < 0)
     {
       stillValid = false;
@@ -293,7 +292,8 @@ bool discretize(const axom::ArrayView<Point2D> &polyline,
   // That was the octahedron count for one segment.  Multiply by the number
   // of segments we will compute.
   int totaloctcount = segoctcount * segmentcount;
-  out = axom::Array<OctType>(totaloctcount, totaloctcount, allocId);
+  out.empty();
+  out.resize(axom::ArrayOptions::Uninitialized(), totaloctcount);
   axom::ArrayView<OctType> out_view = out.view();
   octcount = 0;
 
@@ -303,6 +303,8 @@ bool discretize(const axom::ArrayView<Point2D> &polyline,
       discrSeg<ExecSpace>(polyline[seg], polyline[seg + 1], levels, out_view, octcount);
     octcount += segment_prism_count;
   }
+  // octcount may be < totaloctcount if there are degenerate segments.
+  out.resize(octcount);
 
   // TODO check for errors in each segment's computation
   return true;