Merge pull request #1111 from LLNL/feature/han12/shaping_proe

Add Pro/E input to intersection shaping driver

RELEASE-NOTES.md

@@ -62,6 +62,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
   striding layouts.
 - Adds an `ArrayView::subspan()` overload for multi-dimensional subspans
 - Adds an `axom::utilities::insertionSort()` method.
+- Quest: Adds Pro/E tetrahedral meshes as input to the `IntersectionShaper`
 
 ### Changed
 - Fixed bug in `mint::mesh::UnstructuredMesh` constructors, affecting capacity.

src/axom/primal/geometry/Octahedron.hpp

@@ -72,10 +72,7 @@ class Octahedron
   using PointType = Point<T, NDIMS>;
   using VectorType = Vector<T, NDIMS>;
 
-  enum
-  {
-    NUM_OCT_VERTS = 6
-  };
+  static constexpr int NUM_VERTS = 6;
 
 public:
   /*!
@@ -121,7 +118,7 @@ class Octahedron
   AXOM_HOST_DEVICE
   explicit Octahedron(const PointType* pts)
   {
-    for(int i = 0; i < NUM_OCT_VERTS; i++)
+    for(int i = 0; i < NUM_VERTS; i++)
     {
       m_points[i] = pts[i];
     }
@@ -135,9 +132,9 @@ class Octahedron
   AXOM_HOST_DEVICE
   explicit Octahedron(const axom::ArrayView<PointType> pts)
   {
-    SLIC_ASSERT(pts.size() == NUM_OCT_VERTS);
+    SLIC_ASSERT(pts.size() == NUM_VERTS);
 
-    for(int i = 0; i < NUM_OCT_VERTS; i++)
+    for(int i = 0; i < NUM_VERTS; i++)
     {
       m_points[i] = pts[i];
     }
@@ -151,7 +148,7 @@ class Octahedron
   AXOM_HOST_DEVICE
   explicit Octahedron(std::initializer_list<PointType> pts)
   {
-    SLIC_ASSERT(pts.size() == NUM_OCT_VERTS);
+    SLIC_ASSERT(pts.size() == NUM_VERTS);
 
     int i = 0;
     for(const auto& pt : pts)
@@ -168,7 +165,7 @@ class Octahedron
    */
   AXOM_HOST_DEVICE PointType& operator[](int idx)
   {
-    SLIC_ASSERT(idx >= 0 && idx < NUM_OCT_VERTS);
+    SLIC_ASSERT(idx >= 0 && idx < NUM_VERTS);
     return m_points[idx];
   }
 
@@ -179,7 +176,7 @@ class Octahedron
    */
   AXOM_HOST_DEVICE const PointType& operator[](int idx) const
   {
-    SLIC_ASSERT(idx >= 0 && idx < NUM_OCT_VERTS);
+    SLIC_ASSERT(idx >= 0 && idx < NUM_VERTS);
     return m_points[idx];
   }
 
@@ -192,14 +189,14 @@ class Octahedron
   bool equals(const Octahedron& other, double eps = 1.e-24) const
   {
     // Two octs are equal if each vertex is closer than eps to a vertex of the other.
-    int matched[NUM_OCT_VERTS];
-    for(int i = 0; i < NUM_OCT_VERTS; ++i)
+    int matched[NUM_VERTS];
+    for(int i = 0; i < NUM_VERTS; ++i)
     {
       matched[i] = 0;
     }
-    for(int ourvert = 0; ourvert < NUM_OCT_VERTS; ++ourvert)
+    for(int ourvert = 0; ourvert < NUM_VERTS; ++ourvert)
     {
-      for(int theirvert = 0; theirvert < NUM_OCT_VERTS; ++theirvert)
+      for(int theirvert = 0; theirvert < NUM_VERTS; ++theirvert)
       {
         if(!matched[theirvert] &&
            squared_distance(m_points[ourvert], other[theirvert]) < eps)
@@ -209,34 +206,12 @@ class Octahedron
       }
     }
     int matchedcount = 0;
-    for(int i = 0; i < NUM_OCT_VERTS; ++i)
+    for(int i = 0; i < NUM_VERTS; ++i)
     {
       matchedcount += matched[i];
     }
 
-    return (matchedcount == NUM_OCT_VERTS);
-  }
-
-  /*!
-   * \brief Determines whether the octahedron has duplicate vertices.
-   * \return True if there are duplicate vertices; False otherwise.
-   */
-  AXOM_HOST_DEVICE
-  bool has_duplicate_vertices() const
-  {
-    for(int i = 0; i < 6; i++)
-    {
-      for(int j = i + 1; j < NUM_OCT_VERTS; j++)
-      {
-        // operator= for Point does not want to play nice...
-        if(m_points[i][0] == m_points[j][0] &&
-           m_points[i][1] == m_points[j][1] && m_points[i][2] == m_points[j][2])
-        {
-          return true;
-        }
-      }
-    }
-    return false;
+    return (matchedcount == NUM_VERTS);
   }
 
   /*!
@@ -253,7 +228,7 @@ class Octahedron
   }
 
 private:
-  PointType m_points[NUM_OCT_VERTS];
+  PointType m_points[NUM_VERTS];
 };
 
 //------------------------------------------------------------------------------

src/axom/primal/geometry/Polyhedron.hpp

@@ -511,7 +511,7 @@ class Polyhedron
   {
     double retVol = 0.0;
 
-    if(!isValid())
+    if(!isValid() || hasDuplicateVertices())
     {
       return retVol;
     }
@@ -607,6 +607,31 @@ class Polyhedron
   AXOM_HOST_DEVICE
   bool isValid() const { return m_num_vertices >= 4; }
 
+  /*!
+   * \brief Check if any vertices are duplicate, within a tolerance
+   *
+   * \param [in] eps The tolerance
+   *
+   * \return True, if the polyhedron has duplicate vertices, False otherwise
+   */
+  AXOM_HOST_DEVICE
+  bool hasDuplicateVertices(double eps = 1.e-10) const
+  {
+    for(int i = 0; i < m_num_vertices; i++)
+    {
+      for(int j = i + 1; j < m_num_vertices; j++)
+      {
+        if(axom::utilities::isNearlyEqual(m_vertices[i][0], m_vertices[j][0], eps) &&
+           axom::utilities::isNearlyEqual(m_vertices[i][1], m_vertices[j][1], eps) &&
+           axom::utilities::isNearlyEqual(m_vertices[i][2], m_vertices[j][2], eps))
+        {
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+
   /*!
    * \brief Check if vertex neighbors information is available
    *

src/axom/primal/geometry/Tetrahedron.hpp

@@ -38,7 +38,7 @@ class Tetrahedron
   using VectorType = Vector<T, NDIMS>;
   using SphereType = Sphere<T, NDIMS>;
 
-  static constexpr int NUM_TET_VERTS = 4;
+  static constexpr int NUM_VERTS = 4;
 
 public:
   /// \brief Default constructor. Creates a degenerate tetrahedron.
@@ -76,7 +76,7 @@ class Tetrahedron
   AXOM_HOST_DEVICE
   explicit Tetrahedron(const PointType* pts)
   {
-    for(int i = 0; i < NUM_TET_VERTS; i++)
+    for(int i = 0; i < NUM_VERTS; i++)
     {
       m_points[i] = pts[i];
     }
@@ -90,9 +90,9 @@ class Tetrahedron
   AXOM_HOST_DEVICE
   explicit Tetrahedron(const axom::ArrayView<PointType> pts)
   {
-    SLIC_ASSERT(pts.size() == NUM_TET_VERTS);
+    SLIC_ASSERT(pts.size() == NUM_VERTS);
 
-    for(int i = 0; i < NUM_TET_VERTS; i++)
+    for(int i = 0; i < NUM_VERTS; i++)
     {
       m_points[i] = pts[i];
     }
@@ -106,7 +106,7 @@ class Tetrahedron
   AXOM_HOST_DEVICE
   explicit Tetrahedron(std::initializer_list<PointType> pts)
   {
-    SLIC_ASSERT(pts.size() == NUM_TET_VERTS);
+    SLIC_ASSERT(pts.size() == NUM_VERTS);
 
     int i = 0;
     for(const auto& pt : pts)
@@ -123,7 +123,7 @@ class Tetrahedron
    */
   AXOM_HOST_DEVICE PointType& operator[](int idx)
   {
-    SLIC_ASSERT(idx >= 0 && idx < NUM_TET_VERTS);
+    SLIC_ASSERT(idx >= 0 && idx < NUM_VERTS);
     return m_points[idx];
   }
 
@@ -134,7 +134,7 @@ class Tetrahedron
    */
   AXOM_HOST_DEVICE const PointType& operator[](int idx) const
   {
-    SLIC_ASSERT(idx >= 0 && idx < NUM_TET_VERTS);
+    SLIC_ASSERT(idx >= 0 && idx < NUM_VERTS);
     return m_points[idx];
   }
 
@@ -220,7 +220,7 @@ class Tetrahedron
       "Barycentric coordinates must sum to (near) one.");
 
     PointType res;
-    for(int i = 0; i < NUM_TET_VERTS; ++i)
+    for(int i = 0; i < NUM_VERTS; ++i)
     {
       res.array() += bary[i] * m_points[i].array();
     }

src/axom/primal/operators/compute_bounding_box.hpp

@@ -172,6 +172,23 @@ AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
   return res;
 }
 
+/*!
+ * \brief Creates a bounding box around a Tetrahedron
+ *
+ * \param [in] tet The Tetrahedron
+ */
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
+  const Tetrahedron<T, NDIMS> &tet)
+{
+  BoundingBox<T, NDIMS> res(tet[0]);
+  for(int i = 1; i < 4; i++)
+  {
+    res.addPoint(tet[i]);
+  }
+  return res;
+}
+
 }  // namespace primal
 }  // namespace axom
 

src/axom/primal/operators/split.hpp

@@ -54,11 +54,11 @@ void split(const Octahedron<Tp, NDIMS>& oct,
 
   // Step 1: Find the centroid
   NumArray c;  // ctor fills with 0
-  for(int i = 0; i < Oct::NUM_OCT_VERTS; ++i)
+  for(int i = 0; i < Oct::NUM_VERTS; ++i)
   {
     c += oct[i].array();
   }
-  c = c / (double)Oct::NUM_OCT_VERTS;
+  c = c / static_cast<double>(Oct::NUM_VERTS);
   typename Oct::PointType C(c);
 
   // Step 2: Now store the new tets.  The documentation for the Octahedron class

src/axom/primal/tests/primal_clip.cpp

@@ -1260,6 +1260,36 @@ TEST(primal_clip, tet_tet_clip_split)
   EXPECT_NEAR(0.3333, tet_volumes, EPS);
 }
 
+TEST(primal_clip, tet_tet_clip_special_case_1)
+{
+  using namespace Primal3D;
+  constexpr double EPS = 1e-10;
+  constexpr bool CHECK_SIGN = true;
+
+  // Tets do not intersect, but share a face
+  TetrahedronType tet1(PointType {0.5, 0.5, -0.125},
+                       PointType {0, 0, -0.25},
+                       PointType {1, 0, -0.25},
+                       PointType {0.5, 0, -0.125});
+
+  TetrahedronType tet2(PointType {0.1875, 0.0625, -0.234375},
+                       PointType {0.125, 0.125, -0.25},
+                       PointType {0.125, 0, -0.25},
+                       PointType {0.125, 0.0625, -0.234375});
+
+  PolyhedronType poly = axom::primal::clip(tet1, tet2, EPS, CHECK_SIGN);
+
+  EXPECT_NEAR(0.00, poly.volume(), EPS);
+  EXPECT_NEAR(
+    0.00,
+    axom::primal::intersection_volume<double>(tet2, tet1, EPS, CHECK_SIGN),
+    EPS);
+  EXPECT_NEAR(
+    0.00,
+    axom::primal::intersection_volume<double>(tet1, tet2, EPS, CHECK_SIGN),
+    EPS);
+}
+
 //------------------------------------------------------------------------------
 int main(int argc, char* argv[])
 {

src/axom/primal/tests/primal_tetrahedron.cpp

@@ -388,7 +388,7 @@ TEST_F(TetrahedronTest, tetrahedron_roundtrip_bary_to_physical)
   using CoordType = TetrahedronTest::CoordType;
   using QPoint = TetrahedronTest::QPoint;
   using QTet = TetrahedronTest::QTet;
-  using RPoint = primal::Point<CoordType, QTet::NUM_TET_VERTS>;
+  using RPoint = primal::Point<CoordType, QTet::NUM_VERTS>;
 
   // Test tets
   std::vector<QTet> tets = {this->getTet(0),

src/axom/quest/Discretize.cpp

@@ -110,7 +110,10 @@ OctType new_inscribed_oct(const SphereType& sphere, OctType& o, int s, int t, in
  * This routine allocates an array pointed to by \a out.  The caller is responsible
  * to free the array.
  */
-bool discretize(const SphereType& sphere, int levels, OctType*& out, int& octcount)
+bool discretize(const SphereType& sphere,
+                int levels,
+                axom::Array<OctType>& out,
+                int& octcount)
 {
   // Check input.  Negative radius: return false.
   if(sphere.getRadius() < 0)
@@ -126,7 +129,7 @@ bool discretize(const SphereType& sphere, int levels, OctType*& out, int& octcou
 
   octcount = count_sphere_octahedra(levels);
 
-  out = axom::allocate<OctType>(octcount);
+  out = axom::Array<OctType>(octcount, octcount);
 
   // index points to an octahedron of the last generation.  We'll generate
   // new octahedra based on out[index].
@@ -215,12 +218,12 @@ double octPolyVolume(const OctType& o)
 }  // namespace
 
 //------------------------------------------------------------------------------
-int mesh_from_discretized_polyline(const OctType* octs,
+int mesh_from_discretized_polyline(axom::ArrayView<OctType>& octs,
                                    int octcount,
                                    int segcount,
                                    mint::Mesh*& mesh)
 {
-  SLIC_ASSERT(octs != nullptr);
+  SLIC_ASSERT(octs.data() != nullptr);
 
   const int tetcount = 8 * octcount;
   const int vertcount = 4 * tetcount;