SOR implementation for MeshClipper

RELEASE-NOTES.md

@@ -48,6 +48,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 
 
 ###  Changed
+- Version of `quest::discretize` that approximates a surface-of-revolution from a polyline
+  now respects the allocator ID of the output `Array`.  It no longer resets the ID to the
+  execution space default.
 - Evaluation methods for line integrals in `axom::primal` have been generalized, and 
   `evaluate_scalar_line_integral` has been renamed to `evaluate_line_integral`.
 - Treatment of materials on strided-structured Blueprint meshes has changed in `axom::mir`.

src/axom/primal/geometry/CoordinateTransformer.hpp

@@ -8,6 +8,7 @@
 
 #include "axom/core/numerics/Matrix.hpp"
 #include "axom/core/utilities/Utilities.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/core/numerics/floating_point_limits.hpp"
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/primal/geometry/Vector.hpp"
@@ -157,7 +158,7 @@ class CoordinateTransformer
    *
    * Validity can be checked with isValid().
    */
-  AXOM_HOST_DEVICE void setInvalid() { m_P[0][0] = std::numeric_limits<T>::quiet_NaN(); }
+  AXOM_HOST_DEVICE void setInvalid() { m_P[0][0] = axom::numeric_limits<T>::quiet_NaN(); }
 
   //! @brief Whether transformer is valid.
   AXOM_HOST_DEVICE bool isValid() { return !std::isnan(m_P[0][0]); }
@@ -377,7 +378,7 @@ class CoordinateTransformer
    * M = [ P v ]
    *     [ 0 1 ]
    *
-   * Store m_P[0][0] = std::numeric_limits<T>::quiet_NaN() to set this
+   * Store m_P[0][0] = axom::numeric_limits<T>::quiet_NaN() to set this
    * CoordinateTransformer as invalid.  See \a setInvalid()
    */
   Matrx m_P;
@@ -411,7 +412,7 @@ class CoordinateTransformer
     }
     else
     {
-      m[0][0] = std::numeric_limits<T>::quiet_NaN();
+      m[0][0] = axom::numeric_limits<T>::quiet_NaN();
     }
   }
 

src/axom/quest/CMakeLists.txt

@@ -154,13 +154,17 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE)
                                       detail/clipping/TetClipper.hpp
                                       detail/clipping/HexClipper.hpp
                                       detail/clipping/SphereClipper.hpp
+                                      detail/clipping/MonotonicZSORClipper.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
                                       detail/clipping/HexClipper.cpp
                                       detail/clipping/SphereClipper.cpp
+                                      detail/clipping/MonotonicZSORClipper.cpp
+                                      detail/clipping/SORClipper.cpp
                                       detail/clipping/MeshClipperImpl.hpp)
         endif()
     endif()

src/axom/quest/IntersectionShaper.hpp

@@ -824,7 +824,8 @@ class IntersectionShaper : public Shaper
         axom::fmt::format("{:-^80}", axom::fmt::format(" Refinement level set to {} ", m_level)));
 
       // Generate the Octahedra
-      // (octahedra m_octs will be on device)
+      // (Set m_octs's allocator id to where we want its data to live.)
+      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.clear();
+  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;

src/axom/quest/detail/clipping/MeshClipperImpl.hpp

@@ -109,12 +109,19 @@ class MeshClipperImpl : public MeshClipper::Impl
       });
   }
 
-  //! @brief Make a list of indices where labels have value LABEL_ON.
+  /*!
+   * @brief Make a list of indices where labels have value LABEL_ON,
+   * stored in the same allocator id as the labels.
+   */
   void collectOnIndices(const axom::ArrayView<LabelType>& labels,
                         axom::Array<axom::IndexType>& onIndices) override
   {
     if(labels.empty())
     {
+      if(onIndices.getAllocatorID() != labels.getAllocatorID())
+      {
+        onIndices = axom::Array<IndexType>(0, 0, labels.getAllocatorID());
+      }
       return;
     };