Add SOR, cylinder and cone to in-memory shaping.

LLNL · gunney1 · Jan 14, 2025 · Jul 19, 2024 · Jul 31, 2024 · Jul 31, 2024
commit 63358e9505c1a5a2c40ece10ec093df0f975a6b3
diff --git a/src/axom/klee/Geometry.cpp b/src/axom/klee/Geometry.cpp
@@ -45,16 +45,36 @@ Geometry::Geometry(const TransformableGeometryProperties &startProperties,
 }
 
 Geometry::Geometry(const TransformableGeometryProperties &startProperties,
-                   const axom::primal::Sphere<double, 3>& sphere,
+                   const Sphere3D& sphere,
                    axom::IndexType levelOfRefinement,
                    std::shared_ptr<GeometryOperator const> operator_)
   : m_startProperties(startProperties)
   , m_format("sphere3D")
   , m_path()
   , m_meshGroup(nullptr)
   , m_topology()
+  , m_levelOfRefinement(levelOfRefinement)
   , m_sphere(sphere)
+  , m_operator(std::move(operator_))
+{
+}
+
+Geometry::Geometry(const TransformableGeometryProperties &startProperties,
+                   const axom::Array<double, 2>& discreteFunction,
+                   const Point3D& vorBase,
+                   const Vector3D& vorDirection,
+                   axom::IndexType levelOfRefinement,
+                   std::shared_ptr<GeometryOperator const> operator_)
+  : m_startProperties(startProperties)
+  , m_format("vor3D")
+  , m_path()
+  , m_meshGroup(nullptr)
+  , m_topology()
   , m_levelOfRefinement(levelOfRefinement)
+  , m_sphere()
+  , m_discreteFunction(discreteFunction)
+  , m_vorBase(vorBase)
+  , m_vorDirection(vorDirection)
   , m_operator(std::move(operator_))
 {
 }

diff --git a/src/axom/klee/Geometry.hpp b/src/axom/klee/Geometry.hpp
@@ -56,8 +56,12 @@ inline bool operator!=(const TransformableGeometryProperties &lhs,
 class Geometry
 {
 public:
+  using Point3D = axom::primal::Point<double, 3>;
+  using Vector3D = axom::primal::Vector<double, 3>;
+  using Sphere3D = axom::primal::Sphere<double, 3>;
+
   /**
-   * Create a new Geometry object based on a file representation.
+   * Create a Geometry object based on a file representation.
    *
    * \param startProperties the transformable properties before any
    * operators are applied
@@ -71,7 +75,7 @@ class Geometry
            std::shared_ptr<GeometryOperator const> operator_);
 
   /**
-   * Create a new Geometry object based on a blueprint tetrahedral mesh.
+   * Create a Geometry object based on a blueprint tetrahedral mesh.
    *
    * \param startProperties the transformable properties before any
    * operators are applied
@@ -88,7 +92,7 @@ class Geometry
            std::shared_ptr<GeometryOperator const> operator_);
 
   /**
-   * Create a new sphere Geometry object.
+   * Create a sphere Geometry object.
    *
    * \param startProperties the transformable properties before any
    * operators are applied
@@ -104,16 +108,42 @@ class Geometry
            axom::IndexType levelOfRefinement,
            std::shared_ptr<GeometryOperator const> operator_);
 
+  /**
+   * Create a volume-of-revolution (VOR) Geometry object.
+   *
+   * \param startProperties the transformable properties before any
+   * operators are applied
+   * \param discreteFunction Discrete function describing the surface
+   *        of revolution.
+   * \param vorBase Coordinates of the base of the VOR.
+   * \param vorDirection VOR axis, in the direction of increasing z.
+   * \param levelOfRefinement Number of refinement levels to use for
+   *        discretizing the sphere.
+   * \param operator_ a possibly null operator to apply to the geometry.
+   *
+   * The \c discreteFunction should be an Nx2 array, interpreted as
+   * (z,r) pairs, where z is the axial distance and r is the radius.
+   * The \c vorBase coordinates corresponds to z=0.
+   * \c vorAxis should point in the direction of increasing z.
+   *
+   * \internal TODO: Is this the simplex requirement overly restrictive?
+   */
+  Geometry(const TransformableGeometryProperties &startProperties,
+           const axom::Array<double, 2>& discreteFunction,
+           const Point3D& vorBase,
+           const Vector3D& vorDirection,
+           axom::IndexType levelOfRefinement,
+           std::shared_ptr<GeometryOperator const> operator_);
+
   /**
    * Get the format in which the geometry is specified.
    *
    * Values are:
    * - "c2c" = C2C file
    * - "proe" = ProE file
    * - "memory-blueprint" = Blueprint tetrahedral mesh in memory
-   * - "memory-xy" = discretized 2D, non-negative function as an
-   *   array of points in memory
    * - "sphere3D" = 3D sphere, as \c primal::Sphere<double,3>
+   * - "vor3D" = 3D volume of revolution.
    * - "cone3D" = 3D cone, as \c primal::Cone<double,3>
    *   "cylinder3D" = 3D cylinder, as \c primal::Cylinder<double,3>
    *
@@ -202,6 +232,14 @@ class Geometry
     return m_sphere;
   }
 
+  /**
+   @brief Get the discrete function.
+  */
+  axom::ArrayView<const double, 2> getDiscreteFunction() const
+  {
+    return m_discreteFunction.view();
+  }
+
 private:
   TransformableGeometryProperties m_startProperties;
 
@@ -217,24 +255,20 @@ class Geometry
   //!@brief Topology of the blueprint simplex mesh, if it's in memory.
   std::string m_topology;
 
-  //!@brief The analytical sphere, if used.
-  axom::primal::Sphere<double, 3> m_sphere;
-
-#if 0
-  //!@brief The analytical cylinder, if used.
-  axom::primal::Sphere<double, 3> m_cylinder;
-
-  //!@brief The analytical cone, if used.
-  axom::primal::Sphere<double, 3> m_cone;
-
-  //!@brief The discrete 2D curve, if used.
-  axom::Array<primal::Point<double, 2> m_discreteFunction;
-#endif
-
   //!@brief Level of refinement for discretizing analytical shapes
   // and surfaces of revolutions.
   axom::IndexType m_levelOfRefinement = 0;
 
+  //!@brief The analytical sphere, if used.
+  Sphere3D m_sphere;
+
+  /*!
+    @brief The discrete 2D function, as an Nx2 array, if used.
+  */
+  axom::Array<double, 2> m_discreteFunction;
+  Point3D m_vorBase;
+  Vector3D m_vorDirection;
+
   std::shared_ptr<const GeometryOperator> m_operator;
 };
 

diff --git a/src/axom/quest/DiscreteShape.cpp b/src/axom/quest/DiscreteShape.cpp
@@ -103,12 +103,13 @@ std::shared_ptr<mint::Mesh> DiscreteShape::createMeshRepresentation()
 {
   if (m_meshRep) { return m_meshRep; }
 
+  using TetMesh = axom::mint::UnstructuredMesh<axom::mint::Topology::SINGLE_SHAPE>;
+
   const axom::klee::Geometry& geometry = m_shape.getGeometry();
   const auto& geometryFormat = geometry.getFormat();
 
   if (geometryFormat == "memory-blueprint")
   {
-    // TODO: For readability, move most of this into a function.
     // Put the in-memory geometry in m_meshRep.
     axom::sidre::Group* inputGroup = geometry.getBlueprintMesh();
     axom::sidre::Group* rootGroup = inputGroup->getDataStore()->getRoot();
@@ -171,11 +172,10 @@ std::shared_ptr<mint::Mesh> DiscreteShape::createMeshRepresentation()
         }
       });
 
-    // TODO: Set this up with sidre::Group to have data on device.
-    axom::mint::UnstructuredMesh<axom::mint::Topology::SINGLE_SHAPE>* tetMesh = nullptr;
+    TetMesh* tetMesh = nullptr;
     if (m_sidreGroup != nullptr)
     {
-      tetMesh = new axom::mint::UnstructuredMesh<axom::mint::Topology::SINGLE_SHAPE>(
+      tetMesh = new TetMesh(
         3,
         axom::mint::CellType::TET,
         m_sidreGroup,
@@ -184,7 +184,7 @@ std::shared_ptr<mint::Mesh> DiscreteShape::createMeshRepresentation()
     }
     else
     {
-      tetMesh = new axom::mint::UnstructuredMesh<axom::mint::Topology::SINGLE_SHAPE>(
+      tetMesh = new TetMesh(
         3,
         axom::mint::CellType::TET,
         nodeCount,
@@ -196,6 +196,60 @@ std::shared_ptr<mint::Mesh> DiscreteShape::createMeshRepresentation()
 
     applyTransforms();
   }
+  if (geometryFormat == "vor3D")
+  {
+    // Construct the tet m_meshRep from the volume-of-revolution.
+    auto& vorGeom = m_shape.getGeometry();
+    const auto& discreteFcn = vorGeom.getDiscreteFunction();
+
+    // Generate the Octahedra
+    axom::Array<OctType> octs;
+    int octCount = 0;
+    axom::ArrayView<Point2D> polyline((Point2D*)discreteFcn.data(),
+                                      discreteFcn.shape()[0]);
+    const bool good = axom::quest::discretize<axom::SEQ_EXEC>(
+      polyline,
+      int(polyline.size()),
+      m_shape.getGeometry().getLevelOfRefinement(),
+      octs,
+      octCount);
+    SLIC_ASSERT(good);
+
+    // Dump discretized octs as a tet mesh
+    axom::mint::Mesh* mesh;
+    axom::quest::mesh_from_discretized_polyline(
+      octs.view(),
+      octCount,
+      polyline.size() - 1,
+      mesh);
+
+    if (m_sidreGroup)
+    {
+      auto* tetMesh = static_cast<TetMesh*>(mesh);
+      // Don't directly use tetMesh, because we want the data in Sidre.
+      axom::mint::UnstructuredMesh<axom::mint::SINGLE_SHAPE>*
+        siMesh = new axom::mint::UnstructuredMesh<axom::mint::SINGLE_SHAPE>(
+          tetMesh->getDimension(),
+          tetMesh->getCellType(),
+          m_sidreGroup,
+          tetMesh->getTopologyName(),
+          tetMesh->getCoordsetName(),
+          tetMesh->getNumberOfNodes(),
+          tetMesh->getNumberOfCells());
+      siMesh->appendNodes( tetMesh->getCoordinateArray(0)
+                         , tetMesh->getCoordinateArray(1)
+                         , tetMesh->getCoordinateArray(2)
+                         , tetMesh->getNumberOfNodes() );
+      siMesh->appendCells( tetMesh->getCellNodesArray()
+                         , tetMesh->getNumberOfCells() );
+      m_meshRep.reset(siMesh);
+    } else {
+      m_meshRep.reset(mesh);
+    }
+
+    // Transform the coordinates of the linearized mesh.
+    applyTransforms();
+  }
 
   if (m_meshRep)
   {

diff --git a/src/axom/quest/Discretize.hpp b/src/axom/quest/Discretize.hpp
@@ -69,7 +69,7 @@ bool discretize(const SphereType& s,
  * This routine initializes an Array, \a out.
  */
 template <typename ExecSpace>
-bool discretize(const axom::Array<Point2D>& polyline,
+bool discretize(const axom::ArrayView<Point2D>& polyline,
                 int len,
                 int levels,
                 axom::Array<OctType>& out,

diff --git a/src/axom/quest/IntersectionShaper.hpp b/src/axom/quest/IntersectionShaper.hpp
@@ -317,6 +317,7 @@ class IntersectionShaper : public Shaper
   using Point3D = primal::Point<double, 3>;
   using TetrahedronType = primal::Tetrahedron<double, 3>;
   using SegmentMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
+  using TetMesh = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
 
   using RuntimePolicy = axom::runtime_policy::Policy;
 
@@ -610,21 +611,15 @@ class IntersectionShaper : public Shaper
 
     std::string shapeFormat = shape.getGeometry().getFormat();
 
+    // C2C mesh is not discretized into tets, but all others are.
     if(shapeFormat == "c2c")
     {
       prepareC2CCells<ExecSpace>();
     }
-
-    else if(shapeFormat == "proe")
-    {
-      prepareTetCells<ExecSpace>();
-    }
-
-    else if(shapeFormat == "memory-blueprint" || shapeFormat == "sphere3D")
+    else if(surfaceMeshIsTet())
     {
       prepareTetCells<ExecSpace>();
     }
-
     else
     {
       SLIC_ERROR(
@@ -1492,8 +1487,8 @@ class IntersectionShaper : public Shaper
     AXOM_ANNOTATE_SCOPE("runShapeQuery");
     const std::string shapeFormat = shape.getGeometry().getFormat();
 
-    // Testing separate workflow for Pro/E
-    if(shapeFormat == "proe" || shapeFormat == "memory-blueprint" || shapeFormat == "sphere3D")
+    // C2C mesh is not discretized into tets, but all others are.
+    if(surfaceMeshIsTet())
     {
       switch(m_execPolicy)
       {
@@ -2029,6 +2024,14 @@ class IntersectionShaper : public Shaper
     return volFrac;
   }
 
+  bool surfaceMeshIsTet() const {
+    bool isTet = m_surfaceMesh != nullptr &&
+      m_surfaceMesh->getDimension() == 3 &&
+      !m_surfaceMesh->hasMixedCellTypes() &&
+      m_surfaceMesh->getCellType() == mint::TET;
+    return isTet;
+  }
+
 private:
   RuntimePolicy m_execPolicy {RuntimePolicy::seq};
   int m_level {DEFAULT_CIRCLE_REFINEMENT_LEVEL};

diff --git a/src/axom/quest/Shaper.cpp b/src/axom/quest/Shaper.cpp
@@ -93,7 +93,7 @@ bool Shaper::isValidFormat(const std::string& format) const
 {
   return (format == "stl" || format == "proe" || format == "c2c" ||
           format == "memory-blueprint" || format == "sphere3D" ||
-          format == "none");
+          format == "vor3D" || format == "none");
 }
 
 void Shaper::loadShape(const klee::Shape& shape)

diff --git a/src/axom/quest/detail/Discretize_detail.hpp b/src/axom/quest/detail/Discretize_detail.hpp
@@ -264,7 +264,7 @@ namespace quest
  * This routine initializes an Array pointed to by \a out.
  */
 template <typename ExecSpace>
-bool discretize(const axom::Array<Point2D> &polyline,
+bool discretize(const axom::ArrayView<Point2D> &polyline,
                 int pointcount,
                 int levels,
                 axom::Array<OctType> &out,

diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt
@@ -252,15 +252,16 @@ if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI
     if(AXOM_ENABLE_TESTS AND AXOM_DATA_DIR)
         # 3D shaping tests.  No 2D shaping in this feature, at least for now.
         set(_nranks 1)
-        set(_testshapes "tet" "sphere")
+        set(_testshapes "tetmesh" "sphere" "cyl" "cone" "vor")
         foreach(_testshape ${_testshapes})
           set(_testname "quest_shape_in_memory_${_testshape}")
           axom_add_test(
             NAME ${_testname}
             COMMAND quest_shape_in_memory_ex
+                       --testShape ${_testshape}
                        --method intersection --refinements 2
                        --scale 0.75 0.75 0.75
-                       inline_mesh --min -1 -1 -1 --max 5 5 5 --res 20 20 20 -d 3
+                       inline_mesh --min -2 -2 -2 --max 4 4 4 --res 20 20 20 -d 3
             NUM_MPI_TASKS ${_nranks})
         endforeach()
     endif()