diff --git a/src/axom/klee/Geometry.cpp b/src/axom/klee/Geometry.cpp index 35fe74940c..9f9171530d 100644 --- a/src/axom/klee/Geometry.cpp +++ b/src/axom/klee/Geometry.cpp @@ -6,6 +6,7 @@ #include "axom/klee/Geometry.hpp" #include "axom/klee/GeometryOperators.hpp" +#include "conduit_blueprint_mesh.hpp" #include @@ -13,22 +14,126 @@ namespace axom { namespace klee { -bool operator==(const TransformableGeometryProperties &lhs, - const TransformableGeometryProperties &rhs) +bool operator==(const TransformableGeometryProperties& lhs, + const TransformableGeometryProperties& rhs) { return lhs.dimensions == rhs.dimensions && lhs.units == rhs.units; } -Geometry::Geometry(const TransformableGeometryProperties &startProperties, +Geometry::Geometry(const TransformableGeometryProperties& startProperties, std::string format, std::string path, std::shared_ptr operator_) : m_startProperties(startProperties) , m_format(std::move(format)) , m_path(std::move(path)) + , m_levelOfRefinement(0) , m_operator(std::move(operator_)) { } +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const axom::sidre::Group* meshGroup, + const std::string& topology, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("blueprint-tets") + , m_path() + , m_meshGroup(meshGroup) + , m_topology(topology) + , m_levelOfRefinement(0) + , m_operator(std::move(operator_)) +{ +#ifdef AXOM_DEBUG + SLIC_ASSERT_MSG(isBlueprintTetMesh(m_meshGroup), + "Mesh provided to Geometry is not a valid blueprint " + "unstructured tetrahedral mesh."); +#endif +} + +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const axom::primal::Tetrahedron& tet, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("tet3D") + , m_path() + , m_meshGroup(nullptr) + , m_topology() + , m_tet(tet) + , m_levelOfRefinement(0) + , m_operator(std::move(operator_)) +{ } + +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const axom::primal::Hexahedron& hex, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("hex3D") + , m_path() + , m_meshGroup(nullptr) + , m_topology() + , m_hex(hex) + , m_levelOfRefinement(0) + , m_operator(std::move(operator_)) +{ } + +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const Sphere3D& sphere, + axom::IndexType levelOfRefinement, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("sphere3D") + , m_path() + , m_meshGroup(nullptr) + , m_topology() + , m_sphere(sphere) + , m_levelOfRefinement(levelOfRefinement) + , m_operator(std::move(operator_)) +{ } + +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const axom::Array& discreteFunction, + const Point3D& vorBase, + const Vector3D& vorDirection, + axom::IndexType levelOfRefinement, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("vor3D") + , m_path() + , m_meshGroup(nullptr) + , m_topology() + , m_sphere() + , m_discreteFunction(discreteFunction) + , m_vorBase(vorBase) + , m_vorDirection(vorDirection) + , m_levelOfRefinement(levelOfRefinement) + , m_operator(std::move(operator_)) +{ } + +Geometry::Geometry(const TransformableGeometryProperties& startProperties, + const axom::primal::Plane& plane, + std::shared_ptr operator_) + : m_startProperties(startProperties) + , m_format("plane3D") + , m_path() + , m_meshGroup(nullptr) + , m_topology() + , m_plane(plane) + , m_levelOfRefinement(0) + , m_operator(std::move(operator_)) +{ } + +bool Geometry::hasGeometry() const +{ + bool isInMemory = m_format == "blueprint-tets" || m_format == "sphere3D" || + m_format == "tet3D" || m_format == "hex3D" || m_format == "plane3D" || + m_format == "cone3D" || m_format == "cylinder3D"; + if(isInMemory) + { + return true; + } + return !m_path.empty(); +} + TransformableGeometryProperties Geometry::getEndProperties() const { if(m_operator) @@ -38,5 +143,69 @@ TransformableGeometryProperties Geometry::getEndProperties() const return m_startProperties; } +const axom::sidre::Group* Geometry::getBlueprintMesh() const +{ + SLIC_ASSERT_MSG( + m_meshGroup, + axom::fmt::format( + "The Geometry format '{}' is not specified " + "as a blueprint mesh and/or has not been converted into one.", + m_format)); + return m_meshGroup; +} + +const std::string& Geometry::getBlueprintTopology() const +{ + SLIC_ASSERT_MSG( + m_meshGroup, + axom::fmt::format( + "The Geometry format '{}' is not specified " + "as a blueprint mesh and/or has not been converted into one.", + m_format)); + return m_topology; +} + +bool Geometry::isBlueprintTetMesh(const axom::sidre::Group* meshGroup) const +{ + conduit::Node bpMesh; + meshGroup->createNativeLayout(bpMesh); + + conduit::Node info; + bool isValid = conduit::blueprint::mesh::verify(bpMesh, info); + // The above call to verify causes the crash, even though this + // function is never entered. + if(!isValid) + { + return false; + } + + const auto& topology = + bpMesh.fetch_existing(axom::fmt::format("topologies/{}", m_topology)); + + std::string coordsetName = topology.fetch_existing("coordset").as_string(); + const auto& coordSet = + bpMesh.fetch_existing(axom::fmt::format("coordsets/{}", coordsetName)); + + auto dim = conduit::blueprint::mesh::coordset::dims(coordSet); + if(dim != 3) + { + return false; + } + + auto topoType = topology.fetch_existing("type").as_string(); + if(topoType != "unstructured") + { + return false; + } + + auto shapeType = topology.fetch_existing("elements/shape").as_string(); + if(shapeType != "tet") + { + return false; + } + + return true; +} + } // namespace klee } // namespace axom diff --git a/src/axom/klee/Geometry.hpp b/src/axom/klee/Geometry.hpp index 2e8e63abe8..d794ee7132 100644 --- a/src/axom/klee/Geometry.hpp +++ b/src/axom/klee/Geometry.hpp @@ -11,6 +11,8 @@ #include "axom/klee/Dimensions.hpp" #include "axom/klee/Units.hpp" +#include "axom/primal.hpp" +#include "axom/sidre.hpp" namespace axom { @@ -54,8 +56,15 @@ inline bool operator!=(const TransformableGeometryProperties &lhs, class Geometry { public: + using Point3D = axom::primal::Point; + using Vector3D = axom::primal::Vector; + using Sphere3D = axom::primal::Sphere; + using Tet3D = axom::primal::Tetrahedron; + using Hex3D = axom::primal::Hexahedron; + using Plane3D = axom::primal::Plane; + /** - * Create a new Geometry object. + * Create a Geometry object based on a file representation. * * \param startProperties the transformable properties before any * operators are applied @@ -69,21 +78,164 @@ class Geometry std::shared_ptr operator_); /** - * Get the format in which the geometry is specified. + * Create a Geometry object based on a blueprint tetrahedral mesh. + * + * \param startProperties the transformable properties before any + * operators are applied + * \param meshGroup the geometry in blueprint format. + * The elements should be segments, triangles or tetrahedra. + * \param topology The blueprint topology to use. + * \param operator_ a possibly null operator to apply to the geometry. + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::sidre::Group *meshGroup, + const std::string &topology, + std::shared_ptr operator_); + + /** + * Create a tetrahedron Geometry object. + * + * \param startProperties the transformable properties before any + * operators are applied + * \param tet Tetrahedron + * \param operator_ a possibly null operator to apply to the geometry. + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::primal::Tetrahedron &tet, + std::shared_ptr operator_); + + /** + * Create a hexahedron Geometry object. + * + * \param startProperties the transformable properties before any + * operators are applied + * \param hex Hexahedron + * \param operator_ a possibly null operator to apply to the geometry. + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::primal::Hexahedron &hex, + std::shared_ptr operator_); + + /** + * Create a sphere Geometry object. + * + * \param startProperties the transformable properties before any + * operators are applied + * \param sphere Analytical sphere specifications + * \param levelOfRefinement Number of refinement levels to use for + * discretizing the sphere. + * \param operator_ a possibly null operator to apply to the geometry. + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::primal::Sphere &sphere, + axom::IndexType levelOfRefinement, + std::shared_ptr 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 VOR. + * \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. + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::Array &discreteFunction, + const Point3D &vorBase, + const Vector3D &vorDirection, + axom::IndexType levelOfRefinement, + std::shared_ptr operator_); + + /** + * Create a planar Geometry object. + * + * \param startProperties the transformable properties before any + * operators are applied + * \param tet Tetrahedron + * \param operator_ a possibly null operator to apply to the geometry. + * + * The space on the positive normal side of the plane is considered + * "inside the shape". + */ + Geometry(const TransformableGeometryProperties &startProperties, + const axom::primal::Plane &plane, + std::shared_ptr operator_); + + /** + * @brief Get the format in which the geometry was specified. + * + * The format is determined by the constructor used. + * Values are: + * - "c2c" = C2C file + * - "proe" = ProE file + * - "blueprint-tets" = Blueprint tetrahedral mesh in memory + * - "tet3D" = 3D tetrahedron (4 points) + * - "sphere3D" = 3D sphere, as \c primal::Sphere + * - "vor3D" = 3D volume of revolution. + * - "cone3D" = 3D cone, as \c primal::Cone + * - "cylinder3D" = 3D cylinder, as \c primal::Cylinder + * - "hex3D" = 3D hexahedron (8 points) + * - "plane3D" = 3D plane * * \return the format of the shape + * + * TODO: Depending on the specified geometry, some members are not + * used. It may clarify if we make each supported geometry a + * subclass. */ const std::string &getFormat() const { return m_format; } /** - * Get the path at which to find the specification of the geometry + * Get the path at which to find the specification of the geometry, + * for geometries stored in files. * * \return the path to the geometry file */ const std::string &getPath() const { return m_path; } - /// Predicate that returns true when the shape has an associated geometry - bool hasGeometry() const { return !m_path.empty(); } + /** + * @brief Return the blueprint mesh, for formats that are specified + * by a blueprint mesh or have been converted to a blueprint mesh. + */ + const axom::sidre::Group *getBlueprintMesh() const; + + /** + * @brief Return the blueprint mesh topology, for formats that are specified + * by a blueprint mesh or have been converted to a blueprint mesh. + */ + const std::string &getBlueprintTopology() const; + + /** + @brief Return the VOR axis direction. + */ + const Vector3D getVorDirection() const { return m_vorDirection; } + + /** + @brief Return the 3D coordinates of the VOR base. + */ + const Point3D getVorBaseCoords() const { return m_vorBase; } + + /*! @brief Predicate that returns true when the shape has an associated geometry + + A false means that this is set up to determine volume fractions without + computing on the geometry. + + TODO: We should just create a new format to represent getting + volume fractions without geometries. Or move this logic into + Shape, because it's confusing to have a Geometry that has no + geometry. + */ + bool hasGeometry() const; /** * Get a GeometryOperator to apply to this geometry. Can be null. @@ -113,11 +265,90 @@ class Geometry */ TransformableGeometryProperties getEndProperties() const; + /** + @brief Return the number of levels of refinement for discretization + of analytical curves. + + This number is unused for geometries that are specified in discrete + form. + */ + axom::IndexType getLevelOfRefinement() const { return m_levelOfRefinement; } + + /** + @brief Return the tet geometry, when the Geometry + represents a tetrahedron. + */ + const axom::primal::Tetrahedron &getTet() const { return m_tet; } + + /** + @brief Return the hex geometry, when the Geometry + represents a hexahedron. + */ + const axom::primal::Hexahedron &getHex() const { return m_hex; } + + /** + @brief Return the sphere geometry, when the Geometry + represents an alalytical sphere. + */ + const axom::primal::Sphere &getSphere() const { return m_sphere; } + + /** + @brief Return the plane geometry, when the Geometry + represents a plane. + */ + const axom::primal::Plane &getPlane() const { return m_plane; } + + /** + @brief Get the discrete function used in volumes of revolution. + */ + axom::ArrayView getDiscreteFunction() const + { + return m_discreteFunction.view(); + } + private: TransformableGeometryProperties m_startProperties; + + //!@brief Geometry format. std::string m_format; + + //!@brief Geometry file path, if it's file-based. std::string m_path; + + //!@brief Geometry blueprint simplex mesh, when/if it's in memory. + const axom::sidre::Group *m_meshGroup; + + //!@brief Topology of the blueprint simplex mesh, if it's in memory. + std::string m_topology; + + //!@brief The tetrahedron, if used. + Tet3D m_tet; + + //!@brief The hexahedron, if used. + Hex3D m_hex; + + //!@brief The plane, if used. + Plane3D m_plane; + + //!@brief The analytical sphere, if used. + Sphere3D m_sphere; + + //! @brief The discrete 2D function, as an Nx2 array, if used. + axom::Array m_discreteFunction; + + //!@brief The point corresponding to z=0 on the VOR axis. + Point3D m_vorBase; + + //!@brief VOR axis in the direction of increasing z. + Vector3D m_vorDirection; + + //!@brief Level of refinement for discretizing curved + // analytical shapes and surfaces of revolutions. + axom::IndexType m_levelOfRefinement = 0; + std::shared_ptr m_operator; + + bool isBlueprintTetMesh(const axom::sidre::Group *meshGroup) const; }; } // namespace klee diff --git a/src/axom/klee/GeometryOperatorsIO.hpp b/src/axom/klee/GeometryOperatorsIO.hpp index 04bf55f138..38eaaf7e43 100644 --- a/src/axom/klee/GeometryOperatorsIO.hpp +++ b/src/axom/klee/GeometryOperatorsIO.hpp @@ -91,7 +91,7 @@ class GeometryOperatorData private: Path m_path; std::vector m_singleOperatorData; -}; +}; // GeometryOperatorData /** * Data for a named operator. diff --git a/src/axom/mint/mesh/Mesh.hpp b/src/axom/mint/mesh/Mesh.hpp index f77e5e2279..d876329eb8 100644 --- a/src/axom/mint/mesh/Mesh.hpp +++ b/src/axom/mint/mesh/Mesh.hpp @@ -539,6 +539,12 @@ class Mesh */ inline sidre::Group* getSidreGroup() { return m_group; } + /*! + * \brief Return a const pointer to the sidre::Group associated with this Mesh + * instance or nullptr if none exists. + */ + inline const sidre::Group* getSidreGroup() const { return m_group; } + /*! * \brief Return the name of the topology associated with this Mesh instance, * the return value is undefined if the mesh is not in sidre. diff --git a/src/axom/primal/geometry/Hexahedron.hpp b/src/axom/primal/geometry/Hexahedron.hpp index c16a55886b..424a3f4c36 100644 --- a/src/axom/primal/geometry/Hexahedron.hpp +++ b/src/axom/primal/geometry/Hexahedron.hpp @@ -248,7 +248,7 @@ class Hexahedron * \note Assumes tets is pre-allocated */ AXOM_HOST_DEVICE - void triangulate(axom::StackArray& tets) + void triangulate(axom::StackArray& tets) const { // Hex center (hc) PointType hc = vertexMean(); diff --git a/src/axom/primal/geometry/NumericArray.hpp b/src/axom/primal/geometry/NumericArray.hpp index 3e40fdb0e4..0775826276 100644 --- a/src/axom/primal/geometry/NumericArray.hpp +++ b/src/axom/primal/geometry/NumericArray.hpp @@ -44,8 +44,8 @@ AXOM_HOST_DEVICE bool operator==(const NumericArray& lhs, * \return status true if lhs!=rhs, otherwise, false. */ template -bool operator!=(const NumericArray& lhs, - const NumericArray& rhs); +AXOM_HOST_DEVICE bool operator!=(const NumericArray& lhs, + const NumericArray& rhs); /*! * \brief Performs component-wise addition of two numeric arrays. @@ -73,7 +73,7 @@ AXOM_HOST_DEVICE NumericArray operator-(const NumericArray& lh * \result C resulting numeric array from unary negation. */ template -NumericArray operator-(const NumericArray& arr); +AXOM_HOST_DEVICE NumericArray operator-(const NumericArray& arr); /*! * \brief Scalar multiplication a numeric array; Scalar on rhs. @@ -82,7 +82,8 @@ NumericArray operator-(const NumericArray& arr); * \return C resutling numeric array, \f$ \ni: C_i = scalar*arr_i, \forall i\f$ */ template -NumericArray operator*(const NumericArray& arr, double scalar); +AXOM_HOST_DEVICE NumericArray operator*(const NumericArray& arr, + double scalar); /*! * \brief Scalar multiplication a numeric array; Scalar on lhs. @@ -91,7 +92,8 @@ NumericArray operator*(const NumericArray& arr, double scalar) * \return C resulting numeric array, \f$ \ni: C_i = scalar*arr_i, \forall i\f$ */ template -NumericArray operator*(double scalar, const NumericArray& arr); +AXOM_HOST_DEVICE NumericArray operator*(double scalar, + const NumericArray& arr); /*! * \brief Component-wise multiplication of NumericArrays @@ -111,8 +113,8 @@ AXOM_HOST_DEVICE NumericArray operator*(const NumericArray& lh * \pre \f$ rhs_i != 0.0, \forall i \f$ */ template -NumericArray operator/(const NumericArray& lhs, - const NumericArray& rhs); +AXOM_HOST_DEVICE NumericArray operator/(const NumericArray& lhs, + const NumericArray& rhs); /*! * \brief Scalar division of NumericArray; Scalar on rhs @@ -122,7 +124,8 @@ NumericArray operator/(const NumericArray& lhs, * \pre scalar != 0.0 */ template -NumericArray operator/(const NumericArray& arr, double scalar); +AXOM_HOST_DEVICE NumericArray operator/(const NumericArray& arr, + double scalar); /*! * \brief Coordinate-wise absolute value on the NumericArray @@ -131,7 +134,7 @@ NumericArray operator/(const NumericArray& arr, double scalar) * \return A NumericArray whose coordinates are the absolute value of arr */ template -NumericArray abs(const NumericArray& arr); +AXOM_HOST_DEVICE NumericArray abs(const NumericArray& arr); /*! * \brief Overloaded output operator for numeric arrays @@ -310,6 +313,7 @@ class NumericArray // NOLINT * \pre forall i, arr[i] != 0 * \return A reference to the NumericArray instance after cwise division. */ + AXOM_HOST_DEVICE NumericArray& operator/=(const NumericArray& arr); /*! @@ -535,7 +539,7 @@ AXOM_HOST_DEVICE inline NumericArray& NumericArray::operator*= //------------------------------------------------------------------------------ template -inline NumericArray& NumericArray::operator/=( +AXOM_HOST_DEVICE inline NumericArray& NumericArray::operator/=( const NumericArray& v) { for(int i = 0; i < SIZE; ++i) @@ -715,7 +719,8 @@ AXOM_HOST_DEVICE bool operator==(const NumericArray& lhs, //------------------------------------------------------------------------------ template -bool operator!=(const NumericArray& lhs, const NumericArray& rhs) +AXOM_HOST_DEVICE bool operator!=(const NumericArray& lhs, + const NumericArray& rhs) { return !(lhs == rhs); } @@ -730,8 +735,9 @@ std::ostream& operator<<(std::ostream& os, const NumericArray& arr) //------------------------------------------------------------------------------ template -inline NumericArray operator*(const NumericArray& arr, - double scalar) +AXOM_HOST_DEVICE inline NumericArray operator*( + const NumericArray& arr, + double scalar) { NumericArray result(arr); result *= scalar; @@ -740,8 +746,9 @@ inline NumericArray operator*(const NumericArray& arr, //------------------------------------------------------------------------------ template -inline NumericArray operator*(double scalar, - const NumericArray& arr) +AXOM_HOST_DEVICE inline NumericArray operator*( + double scalar, + const NumericArray& arr) { NumericArray result(arr); result *= scalar; @@ -772,8 +779,9 @@ AXOM_HOST_DEVICE inline NumericArray operator*( //------------------------------------------------------------------------------ template -inline NumericArray operator/(const NumericArray& lhs, - const NumericArray& rhs) +AXOM_HOST_DEVICE inline NumericArray operator/( + const NumericArray& lhs, + const NumericArray& rhs) { NumericArray result(lhs); result /= rhs; @@ -782,8 +790,9 @@ inline NumericArray operator/(const NumericArray& lhs, //------------------------------------------------------------------------------ template -inline NumericArray operator/(const NumericArray& arr, - double scalar) +AXOM_HOST_DEVICE inline NumericArray operator/( + const NumericArray& arr, + double scalar) { NumericArray result(arr); result /= scalar; @@ -803,7 +812,8 @@ AXOM_HOST_DEVICE inline NumericArray operator-( //------------------------------------------------------------------------------ template -inline NumericArray operator-(const NumericArray& arr) +AXOM_HOST_DEVICE inline NumericArray operator-( + const NumericArray& arr) { NumericArray result; result -= arr; @@ -812,7 +822,7 @@ inline NumericArray operator-(const NumericArray& arr) //------------------------------------------------------------------------------ template -inline NumericArray abs(const NumericArray& arr) +AXOM_HOST_DEVICE inline NumericArray abs(const NumericArray& arr) { NumericArray result(arr); diff --git a/src/axom/primal/operators/split.hpp b/src/axom/primal/operators/split.hpp index 0fea650b4f..2425df78a7 100644 --- a/src/axom/primal/operators/split.hpp +++ b/src/axom/primal/operators/split.hpp @@ -78,6 +78,58 @@ void split(const Octahedron& oct, out.push_back(Tet(oct[Q], oct[S], oct[U], C)); }; +/*! + * \brief Splits an Octahedron into eight Tetrahedrons + * + * \tparam Tp the coordinate type, such double or float + * \tparam NDIMS the number of spatial dimensions (must be 3). + * \param [in] oct The Octahedron to split + * \param [out] outPtr C-style array of 8 Tetrahedron objects; + * the fragments of oct are appended to out. + * + * \pre NDIMS == 3 + * + * The tets are produced by putting a vertex at the centroid of the oct + * and drawing an edge from each vertex to the centroid. + */ +template +AXOM_HOST_DEVICE void split(const Octahedron& oct, + Tetrahedron* outPtr) +{ + // Implemented for three dimensions + SLIC_ASSERT(NDIMS == 3); + + // Type aliases + using NumArray = NumericArray; + using Oct = Octahedron; + using Tet = Tetrahedron; + + // Step 1: Find the centroid + NumArray c; // ctor fills with 0 + for(int i = 0; i < Oct::NUM_VERTS; ++i) + { + c += oct[i].array(); + } + c = c / static_cast(Oct::NUM_VERTS); + typename Oct::PointType C(c); + + // Step 2: Now store the new tets. The documentation for the Octahedron class + // shows how the points are arranged to imply the faces. + + // clang-format off + enum OctVerts {P, Q, R, S, T, U}; + // clang-format on + + outPtr[0] = Tet(oct[P], oct[R], oct[Q], C); + outPtr[1] = Tet(oct[Q], oct[R], oct[S], C); + outPtr[2] = Tet(oct[R], oct[T], oct[S], C); + outPtr[3] = Tet(oct[S], oct[T], oct[U], C); + outPtr[4] = Tet(oct[T], oct[P], oct[U], C); + outPtr[5] = Tet(oct[U], oct[P], oct[Q], C); + outPtr[6] = Tet(oct[P], oct[T], oct[R], C); + outPtr[7] = Tet(oct[Q], oct[S], oct[U], C); +}; + } // namespace primal } // namespace axom diff --git a/src/axom/quest/CMakeLists.txt b/src/axom/quest/CMakeLists.txt index 079b8cbdd9..c9dd6f67da 100644 --- a/src/axom/quest/CMakeLists.txt +++ b/src/axom/quest/CMakeLists.txt @@ -134,8 +134,10 @@ if(MFEM_FOUND AND AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE AND AXOM_ENABLE_MFEM_SI list(APPEND quest_headers Shaper.hpp SamplingShaper.hpp IntersectionShaper.hpp + DiscreteShape.hpp detail/shaping/shaping_helpers.hpp) list(APPEND quest_sources Shaper.cpp + DiscreteShape.cpp detail/shaping/shaping_helpers.cpp) list(APPEND quest_depends_on klee) endif() diff --git a/src/axom/quest/DiscreteShape.cpp b/src/axom/quest/DiscreteShape.cpp new file mode 100644 index 0000000000..cc45a2915d --- /dev/null +++ b/src/axom/quest/DiscreteShape.cpp @@ -0,0 +1,805 @@ +// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "axom/quest/DiscreteShape.hpp" +#include "axom/quest/Discretize.hpp" +#include "axom/mint/mesh/UnstructuredMesh.hpp" +#include "axom/klee/GeometryOperators.hpp" +#include "axom/core/utilities/StringUtilities.hpp" +#include "axom/quest/interface/internal/QuestHelpers.hpp" + +#include +#include +#include + +namespace axom +{ +namespace quest +{ +namespace internal +{ +/*! + * \brief Implementation of a GeometryOperatorVisitor for processing klee shape operators + * + * This class extracts the matrix form of supported operators and marks the operator as unvalid otherwise + * To use, check the \a isValid() function after visiting and then call the \a getMatrix() function. + */ +class AffineMatrixVisitor : public klee::GeometryOperatorVisitor +{ +public: + AffineMatrixVisitor() : m_matrix(4, 4) { } + + void visit(const klee::Translation& translation) override + { + m_matrix = translation.toMatrix(); + m_isValid = true; + } + void visit(const klee::Rotation& rotation) override + { + m_matrix = rotation.toMatrix(); + m_isValid = true; + } + void visit(const klee::Scale& scale) override + { + m_matrix = scale.toMatrix(); + m_isValid = true; + } + void visit(const klee::UnitConverter& converter) override + { + m_matrix = converter.toMatrix(); + m_isValid = true; + } + + void visit(const klee::CompositeOperator&) override + { + SLIC_WARNING("CompositeOperator not supported for Shaper query"); + m_isValid = false; + } + void visit(const klee::SliceOperator&) override + { + SLIC_WARNING("SliceOperator not yet supported for Shaper query"); + m_isValid = false; + } + + const numerics::Matrix& getMatrix() const { return m_matrix; } + bool isValid() const { return m_isValid; } + +private: + bool m_isValid {false}; + numerics::Matrix m_matrix; +}; + +} // end namespace internal + +// These were needed for linking - but why? They are constexpr. +constexpr int DiscreteShape::DEFAULT_SAMPLES_PER_KNOT_SPAN; +constexpr double DiscreteShape::MINIMUM_PERCENT_ERROR; +constexpr double DiscreteShape::MAXIMUM_PERCENT_ERROR; +constexpr double DiscreteShape::DEFAULT_VERTEX_WELD_THRESHOLD; + +DiscreteShape::DiscreteShape(const axom::klee::Shape& shape, + axom::sidre::Group* parentGroup, + const std::string& prefixPath) + : m_shape(shape) + , m_sidreGroup(nullptr) + , m_refinementType(DiscreteShape::RefinementUniformSegments) + , m_percentError( + utilities::clampVal(0.0, MINIMUM_PERCENT_ERROR, MAXIMUM_PERCENT_ERROR)) +{ + setPrefixPath(prefixPath); + setParentGroup(parentGroup); + + if(parentGroup == nullptr && + m_shape.getGeometry().getFormat() == "blueprint-tets") + { + SLIC_ERROR( + "DiscreteShape: Support for Blueprint-mesh shape format currently " + "requires a non-null parentGroup in the constructor. This restriction " + "can be lifted with additional coding, so please file an issue with the " + "Axom team if you need this feature."); + } +} + +void DiscreteShape::clearInternalData() { m_meshRep.reset(); } + +std::shared_ptr DiscreteShape::createMeshRepresentation() +{ + if(m_meshRep) + { + return m_meshRep; + } + + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + const auto& geometryFormat = geometry.getFormat(); + + if(geometryFormat == "blueprint-tets") + { + createBlueprintTetsRepresentation(); + } + else if(geometryFormat == "tet3D") + { + createTetRepresentation(); + } + else if(geometryFormat == "hex3D") + { + createHexRepresentation(); + } + else if(geometryFormat == "plane3D") + { + createPlaneRepresentation(); + } + else if(geometryFormat == "sphere3D") + { + createSphereRepresentation(); + } + if(geometryFormat == "vor3D") + { + createVORRepresentation(); + } + + if(m_meshRep) + { + return m_meshRep; + } + + if(!m_shape.getGeometry().hasGeometry()) + { + // If shape has no geometry, there's nothing to discretize. + SLIC_DEBUG( + axom::fmt::format("Current shape '{}' of material '{}' has no geometry", + m_shape.getName(), + m_shape.getMaterial())); + return m_meshRep; + } + + // We handled all the non-file formats. The rest are file formats. + const std::string& file_format = geometryFormat; + + std::string shapePath = resolvePath(); + SLIC_INFO("Reading file: " << shapePath << "..."); + + // Initialize revolved volume. + m_revolvedVolume = 0.; + + if(utilities::string::endsWith(shapePath, ".stl")) + { + SLIC_ASSERT_MSG( + file_format == "stl", + axom::fmt::format(" '{}' format requires .stl file type", file_format)); + + axom::mint::Mesh* meshRep = nullptr; + quest::internal::read_stl_mesh(shapePath, meshRep, m_comm); + m_meshRep.reset(meshRep); + // Transform the coordinates of the linearized mesh. + applyTransforms(); + } + else if(utilities::string::endsWith(shapePath, ".proe")) + { + SLIC_ASSERT_MSG( + file_format == "proe", + axom::fmt::format(" '{}' format requires .proe file type", file_format)); + + axom::mint::Mesh* meshRep = nullptr; + quest::internal::read_pro_e_mesh(shapePath, meshRep, m_comm); + m_meshRep.reset(meshRep); + } +#ifdef AXOM_USE_C2C + else if(utilities::string::endsWith(shapePath, ".contour")) + { + SLIC_ASSERT_MSG( + file_format == "c2c", + axom::fmt::format(" '{}' format requires .contour file type", file_format)); + + // Get the transforms that are being applied to the mesh. Get them + // as a single concatenated matrix. + auto transform = getTransforms(); + + // Pass in the transform so any transformations can figure into + // computing the revolved volume. + axom::mint::Mesh* meshRep = nullptr; + if(m_refinementType == DiscreteShape::RefinementDynamic && + m_percentError > MINIMUM_PERCENT_ERROR) + { + quest::internal::read_c2c_mesh_non_uniform(shapePath, + transform, + m_percentError, + m_vertexWeldThreshold, + meshRep, + m_revolvedVolume, // output arg + m_comm); + } + else + { + quest::internal::read_c2c_mesh_uniform(shapePath, + transform, + m_samplesPerKnotSpan, + m_vertexWeldThreshold, + meshRep, + m_revolvedVolume, // output arg + m_comm); + } + m_meshRep.reset(meshRep); + + // Transform the coordinates of the linearized mesh. + applyTransforms(); + } +#endif + else + { + SLIC_ERROR( + axom::fmt::format("Unsupported filetype for this Axom configuration. " + "Provided file was '{}', with format '{}'", + shapePath, + file_format)); + } + + return m_meshRep; +} + +void DiscreteShape::createBlueprintTetsRepresentation() +{ + SLIC_ASSERT(m_sidreGroup != nullptr); + + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + + // Put the in-memory geometry in m_meshRep. + const axom::sidre::Group* inputGroup = geometry.getBlueprintMesh(); +#ifdef AXOM_USE_UMPIRE + int allocID = inputGroup->getDefaultAllocatorID(); +#else + int allocID = axom::execution_space::allocatorID(); +#endif + + std::string modName = inputGroup->getName() + "_modified"; + while(m_sidreGroup->hasGroup(modName)) + { + modName = modName + "-"; + } + + /* + We use a sidre::Group::deepCopyGroup as a convenient way to copy + the mesh. This requires a non-null m_sidreGroup. This + restriction can be lifted if we implement code to copy the + inputGroup directly into an mint::UnstructuredMesh without putting + it into another sidre::Group. We currently don't have that. + */ + axom::sidre::Group* modGroup = m_sidreGroup->createGroup(modName); + modGroup->deepCopyGroup(inputGroup, allocID); + + m_meshRep.reset( + axom::mint::getMesh(modGroup->getGroup(inputGroup->getName()), + m_shape.getGeometry().getBlueprintTopology())); + + // Transform the coordinates of the linearized mesh. + applyTransforms(); +} + +void DiscreteShape::createTetRepresentation() +{ + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + + const auto& tet = geometry.getTet(); + + const axom::IndexType tetCount = 1; + const axom::IndexType nodeCount = 4; + axom::Array nodeCoords(nodeCount, 3); + axom::Array connectivity(tetCount, 4); + + for(int iNode = 0; iNode < 4; ++iNode) + { + const auto& coords = tet[iNode]; + nodeCoords[iNode][0] = coords[0]; + nodeCoords[iNode][1] = coords[1]; + nodeCoords[iNode][2] = coords[2]; + connectivity[0][iNode] = iNode; + } + + TetMesh* tetMesh = nullptr; + if(m_sidreGroup != nullptr) + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + m_sidreGroup, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + else + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + tetMesh->appendNodes((double*)nodeCoords.data(), nodeCoords.shape()[0]); + tetMesh->appendCells(connectivity.data(), connectivity.shape()[0]); + m_meshRep.reset(tetMesh); + + applyTransforms(); +} + +void DiscreteShape::createHexRepresentation() +{ + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + const auto& hex = geometry.getHex(); + axom::StackArray tets; + hex.triangulate(tets); + + const axom::IndexType tetCount = HexType::NUM_TRIANGULATE; + const axom::IndexType nodeCount = HexType::NUM_TRIANGULATE * 4; + axom::Array nodeCoords(nodeCount, 3); + auto nodeCoordsView = nodeCoords.view(); + axom::Array connectivity(tetCount, 4); + auto connectivityView = connectivity.view(); + // NOTE: This is not much computation, so just run on host. + axom::for_all( + tetCount, + AXOM_LAMBDA(axom::IndexType iTet) { + const auto& tet = tets[iTet]; + for(int i = 0; i < 4; ++i) + { + axom::IndexType iNode = iTet * 4 + i; + const auto& coords = tet[i]; + nodeCoordsView[iNode][0] = coords[0]; + nodeCoordsView[iNode][1] = coords[1]; + nodeCoordsView[iNode][2] = coords[2]; + connectivityView[iTet][i] = iNode; + } + }); + + TetMesh* tetMesh = nullptr; + if(m_sidreGroup != nullptr) + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + m_sidreGroup, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + else + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + tetMesh->appendNodes((double*)nodeCoords.data(), nodeCoords.shape()[0]); + tetMesh->appendCells(connectivity.data(), connectivity.shape()[0]); + m_meshRep.reset(tetMesh); + + applyTransforms(); +} + +void DiscreteShape::createPlaneRepresentation() +{ + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + + const auto& plane = geometry.getPlane(); + // Generate a big bounding hex on the positive side of the plane. + /* + TODO: Use a specialized plane representation. Plane is very + simple, but representing it as huge tets has 2 problems. + 1. The shaper runs slow because the huge tets contains huge + numbers of cells, disbling attempts at divide-and-conquer + with the BVH. + + 2. When the tets are orders of magnitude larger then the + computational mesh, I've seen weird overlap computation + errors that affect certain resolutions and don't follow + any trends. I think it's due to round-off, but I've not + verified. + */ + axom::primal::Hexahedron boundingHex; + const double len = 1e2; // Big enough to contain anticipated mesh, + // but too big invites weird round-off errors. + // We should compute based on the mesh. + // clang-format off + boundingHex[0] = Point3D{0.0, -len, -len}; + boundingHex[1] = Point3D{len, -len, -len}; + boundingHex[2] = Point3D{len, len, -len}; + boundingHex[3] = Point3D{0.0, len, -len}; + boundingHex[4] = Point3D{0.0, -len, len}; + boundingHex[5] = Point3D{len, -len, len}; + boundingHex[6] = Point3D{len, len, len}; + boundingHex[7] = Point3D{0.0, len, len}; + // clang-format on + numerics::Matrix rotate = vorAxisRotMatrix(plane.getNormal()); + const auto translate = plane.getNormal() * plane.getOffset(); + for(int i = 0; i < 8; ++i) + { + Point3D newCoords; + numerics::matrix_vector_multiply(rotate, + boundingHex[i].data(), + newCoords.data()); + newCoords.array() += translate.array(); + boundingHex[i].array() = newCoords.array(); + } + + axom::StackArray tets; + boundingHex.triangulate(tets); + + const axom::IndexType tetCount = HexType::NUM_TRIANGULATE; + const axom::IndexType nodeCount = HexType::NUM_TRIANGULATE * 4; + axom::Array nodeCoords(nodeCount, 3); + auto nodeCoordsView = nodeCoords.view(); + axom::Array connectivity(tetCount, 4); + auto connectivityView = connectivity.view(); + // NOTE: This is not much computation, so just run on host. + axom::for_all( + tetCount, + AXOM_LAMBDA(axom::IndexType iTet) { + const auto& tet = tets[iTet]; + for(int i = 0; i < 4; ++i) + { + axom::IndexType iNode = iTet * 4 + i; + const auto& coords = tet[i]; + nodeCoordsView[iNode][0] = coords[0]; + nodeCoordsView[iNode][1] = coords[1]; + nodeCoordsView[iNode][2] = coords[2]; + connectivityView[iTet][i] = iNode; + } + }); + + TetMesh* tetMesh = nullptr; + if(m_sidreGroup != nullptr) + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + m_sidreGroup, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + else + { + tetMesh = new TetMesh(3, + axom::mint::CellType::TET, + nodeCoords.shape()[0], + connectivity.shape()[0]); + } + tetMesh->appendNodes((double*)nodeCoords.data(), nodeCoords.shape()[0]); + tetMesh->appendCells(connectivity.data(), connectivity.shape()[0]); + m_meshRep.reset(tetMesh); + + applyTransforms(); +} + +void DiscreteShape::createSphereRepresentation() +{ + const axom::klee::Geometry& geometry = m_shape.getGeometry(); + + /* + Discretize the sphere into m_meshRep and apply transforms: + 1. Discretize the sphere into a set of Octahedra. + 2. Split each Octahedron into 8 tets. + 3. Insert tets into tet mesh. + 4. Transform the tet mesh. + We can reduce memory use by eliminating repeated points if it + is a problem. + */ + const auto& sphere = geometry.getSphere(); + axom::Array octs; + int octCount = 0; + axom::quest::discretize(sphere, geometry.getLevelOfRefinement(), octs, octCount); + + constexpr int TETS_PER_OCT = 8; + constexpr int NODES_PER_TET = 4; + const axom::IndexType tetCount = octCount * TETS_PER_OCT; + const axom::IndexType nodeCount = tetCount * NODES_PER_TET; + + axom::Array nodeCoords(nodeCount, nodeCount); + axom::Array connectivity( + axom::StackArray {tetCount, NODES_PER_TET}); + + auto nodeCoordsView = nodeCoords.view(); + auto connectivityView = connectivity.view(); + axom::for_all( + octCount, + AXOM_LAMBDA(axom::IndexType octIdx) { + TetType tetsInOct[TETS_PER_OCT]; + axom::primal::split(octs[octIdx], tetsInOct); + for(int iTet = 0; iTet < TETS_PER_OCT; ++iTet) + { + axom::IndexType tetIdx = octIdx * TETS_PER_OCT + iTet; + for(int iNode = 0; iNode < NODES_PER_TET; ++iNode) + { + axom::IndexType nodeIdx = tetIdx * NODES_PER_TET + iNode; + nodeCoordsView[nodeIdx] = tetsInOct[iTet][iNode]; + connectivityView[tetIdx][iNode] = nodeIdx; + } + } + }); + + TetMesh* tetMesh = nullptr; + if(m_sidreGroup != nullptr) + { + tetMesh = + new TetMesh(3, axom::mint::CellType::TET, m_sidreGroup, nodeCount, tetCount); + } + else + { + tetMesh = new TetMesh(3, axom::mint::CellType::TET, nodeCount, tetCount); + } + tetMesh->appendNodes((double*)nodeCoords.data(), nodeCount); + tetMesh->appendCells(connectivity.data(), tetCount); + m_meshRep.reset(tetMesh); + + applyTransforms(); +} + +void DiscreteShape::createVORRepresentation() +{ + // 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 octs; + int octCount = 0; + axom::ArrayView polyline((Point2D*)discreteFcn.data(), + discreteFcn.shape()[0]); + const bool good = axom::quest::discretize( + polyline, + int(polyline.size()), + m_shape.getGeometry().getLevelOfRefinement(), + octs, + octCount); + SLIC_ASSERT(good); + + // Rotate to the VOR axis direction and translate to the base location. + numerics::Matrix rotate = vorAxisRotMatrix(vorGeom.getVorDirection()); + const auto& translate = vorGeom.getVorBaseCoords(); + auto octsView = octs.view(); + axom::for_all( + octCount, + AXOM_LAMBDA(axom::IndexType iOct) { + auto& oct = octsView[iOct]; + for(int iVert = 0; iVert < OctType::NUM_VERTS; ++iVert) + { + auto& newCoords = oct[iVert]; + auto oldCoords = newCoords; + numerics::matrix_vector_multiply(rotate, + oldCoords.data(), + newCoords.data()); + newCoords.array() += translate.array(); + } + }); + + // 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) + { + // If using sidre, copy the tetMesh into sidre. + auto* tetMesh = static_cast(mesh); + axom::mint::UnstructuredMesh* siMesh = + new axom::mint::UnstructuredMesh( + 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); + delete mesh; + mesh = nullptr; + } + else + { + m_meshRep.reset(mesh); + } + + // Transform the coordinates of the linearized mesh. + applyTransforms(); +} + +void DiscreteShape::applyTransforms() +{ + numerics::Matrix transformation = getTransforms(); + + // Apply transformation to coordinates of each vertex in mesh + if(!transformation.isIdentity()) + { + const int spaceDim = m_meshRep->getDimension(); + const int numSurfaceVertices = m_meshRep->getNumberOfNodes(); + double* x = m_meshRep->getCoordinateArray(mint::X_COORDINATE); + double* y = m_meshRep->getCoordinateArray(mint::Y_COORDINATE); + double* z = spaceDim > 2 ? m_meshRep->getCoordinateArray(mint::Z_COORDINATE) + : nullptr; + + double xformed[4]; + for(int i = 0; i < numSurfaceVertices; ++i) + { + double coords[4] = {x[i], y[i], (z == nullptr ? 0. : z[i]), 1.}; + numerics::matrix_vector_multiply(transformation, coords, xformed); + x[i] = xformed[0]; + y[i] = xformed[1]; + if(z != nullptr) + { + z[i] = xformed[2]; + } + } + } +} + +numerics::Matrix DiscreteShape::getTransforms() const +{ + const auto identity4x4 = numerics::Matrix::identity(4); + numerics::Matrix transformation(identity4x4); + auto& geometryOperator = m_shape.getGeometry().getGeometryOperator(); + if(geometryOperator) + { + auto composite = + std::dynamic_pointer_cast(geometryOperator); + if(composite) + { + // Concatenate the transformations + + // Why don't we multiply the matrices in CompositeOperator::addOperator()? + // Why keep the matrices factored and multiply them here repeatedly? + // Combining them would also avoid this if-else logic. BTNG + for(auto op : composite->getOperators()) + { + // Use visitor pattern to extract the affine matrix from supported operators + internal::AffineMatrixVisitor visitor; + op->accept(visitor); + if(!visitor.isValid()) + { + continue; + } + const auto& matrix = visitor.getMatrix(); + numerics::Matrix res(identity4x4); + numerics::matrix_multiply(matrix, transformation, res); + transformation = res; + } + } + else + { + internal::AffineMatrixVisitor visitor; + geometryOperator->accept(visitor); + if(visitor.isValid()) + { + transformation = visitor.getMatrix(); + } + } + } + return transformation; +} + +// Return a 3x3 matrix that rotate coordinates from the x-axis to the given direction. +numerics::Matrix DiscreteShape::vorAxisRotMatrix(const Vector3D& dir) +{ + // Note that the rotation matrix is not unique. + static const Vector3D x {1.0, 0.0, 0.0}; + Vector3D a = dir.unitVector(); + Vector3D u; // Rotation vector, the cross product of x and a. + numerics::cross_product(x.data(), a.data(), u.data()); + double sinT = u.norm(); + double cosT = numerics::dot_product(x.data(), a.data(), 3); + double ccosT = 1 - cosT; + + // Degenerate case with angle near 0 or pi. + if(utilities::isNearlyEqual(sinT, 0.0)) + { + if(cosT > 0) + { + return numerics::Matrix::identity(3); + } + else + { + // Give u a tiny component in any non-x direction + // so we can rotate around it. + u[1] = 1e-8; + } + } + + u = u.unitVector(); + numerics::Matrix rot(3, 3, 0.0); + rot(0, 0) = u[0] * u[0] * ccosT + cosT; + rot(0, 1) = u[0] * u[1] * ccosT - u[2] * sinT; + rot(0, 2) = u[0] * u[2] * ccosT + u[1] * sinT; + rot(1, 0) = u[1] * u[0] * ccosT + u[2] * sinT; + rot(1, 1) = u[1] * u[1] * ccosT + cosT; + rot(1, 2) = u[1] * u[2] * ccosT - u[0] * sinT; + rot(2, 0) = u[2] * u[0] * ccosT - u[1] * sinT; + rot(2, 1) = u[2] * u[1] * ccosT + u[0] * sinT; + rot(2, 2) = u[2] * u[2] * ccosT + cosT; + + return rot; +} + +void DiscreteShape::setSamplesPerKnotSpan(int nSamples) +{ + using axom::utilities::clampLower; + SLIC_WARNING_IF( + nSamples < 1, + axom::fmt::format( + "Samples per knot span must be at least 1. Provided value was {}", + nSamples)); + + m_samplesPerKnotSpan = clampLower(nSamples, 1); +} + +void DiscreteShape::setVertexWeldThreshold(double threshold) +{ + SLIC_WARNING_IF( + threshold <= 0., + axom::fmt::format( + "Vertex weld threshold should be positive Provided value was {}", + threshold)); + + m_vertexWeldThreshold = threshold; +} + +void DiscreteShape::setPercentError(double percent) +{ + using axom::utilities::clampVal; + SLIC_WARNING_IF( + percent <= MINIMUM_PERCENT_ERROR, + axom::fmt::format("Percent error must be greater than {}. Provided value " + "was {}. Dynamic refinement will not be used.", + MINIMUM_PERCENT_ERROR, + percent)); + SLIC_WARNING_IF(percent > MAXIMUM_PERCENT_ERROR, + axom::fmt::format( + "Percent error must be less than {}. Provided value was {}", + MAXIMUM_PERCENT_ERROR, + percent)); + if(percent <= MINIMUM_PERCENT_ERROR) + { + m_refinementType = DiscreteShape::RefinementUniformSegments; + } + m_percentError = + clampVal(percent, MINIMUM_PERCENT_ERROR, MAXIMUM_PERCENT_ERROR); +} + +std::string DiscreteShape::resolvePath() const +{ + const std::string& geomPath = m_shape.getGeometry().getPath(); + if(geomPath[0] == '/') + { + return geomPath; + } + if(m_prefixPath.empty()) + { + throw std::logic_error("Relative geometry path requires a parent path."); + } + std::string dir; + utilities::filesystem::getDirName(dir, m_prefixPath); + return utilities::filesystem::joinPath(dir, geomPath); +} + +void DiscreteShape::setParentGroup(axom::sidre::Group* parentGroup) +{ + if(parentGroup) + { + // Use object address to create a unique name for sidre group under parent. + std::string myGroupName; + int i = 0; + while(myGroupName.empty()) + { + std::ostringstream os; + os << "DiscreteShapeTemp-" << i; + if(!parentGroup->hasGroup(os.str())) + { + myGroupName = os.str(); + } + ++i; + } + m_sidreGroup = parentGroup->createGroup(myGroupName); + } +} + +} // namespace quest +} // namespace axom diff --git a/src/axom/quest/DiscreteShape.hpp b/src/axom/quest/DiscreteShape.hpp new file mode 100644 index 0000000000..d64ae81f39 --- /dev/null +++ b/src/axom/quest/DiscreteShape.hpp @@ -0,0 +1,198 @@ +// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#ifndef AXOM_QUEST_DISCRETE_SHAPE_HPP +#define AXOM_QUEST_DISCRETE_SHAPE_HPP + +#include +#include + +#include "axom/klee/Shape.hpp" +#include "axom/mint/mesh/UnstructuredMesh.hpp" + +#if defined(AXOM_USE_MPI) + #include "mpi.h" +#endif + +namespace axom +{ +namespace quest +{ +/*! + @brief Post-processed klee::Shape, with the geometry discretized and transformed + according to the Shape's operators. + + TODO: Move this class into internal namespace. +*/ +class DiscreteShape +{ +public: + /// Refinement type. + using RefinementType = enum { RefinementUniformSegments, RefinementDynamic }; + + using Point3D = axom::primal::Point; + using Vector3D = axom::primal::Vector; + using TetType = axom::primal::Tetrahedron; + using OctType = axom::primal::Octahedron; + using HexType = axom::primal::Hexahedron; + + // Common type for the mesh approximation of the shape. + using TetMesh = + axom::mint::UnstructuredMesh; + + static constexpr int DEFAULT_SAMPLES_PER_KNOT_SPAN {25}; + static constexpr double MINIMUM_PERCENT_ERROR {0.}; + static constexpr double MAXIMUM_PERCENT_ERROR {100.}; + static constexpr double DEFAULT_VERTEX_WELD_THRESHOLD {1e-9}; + + /*! + @brief Constructor. + + @param shape The Klee specifications for the shape. + @param parentGroup Group under which to put the discrete mesh + and support blueprint-tets shapes. + If null, don't use sidre and don't support blueprint-tets. + @param prefixPath Path prefix for shape files specified + with a relative path. + + Refinement type is set to DiscreteShape::RefinementUniformSegments + and percent erro is set to 0. See setPercentError() and + setRefinementType(). + */ + DiscreteShape(const axom::klee::Shape& shape, + axom::sidre::Group* parentGroup, + const std::string& prefixPath = {}); + + virtual ~DiscreteShape() { clearInternalData(); } + + //@{ + //! @name Functions to get and set shaping parameters + + void setPrefixPath(const std::string& prefixPath) + { + m_prefixPath = prefixPath; + } + + /*! + @brief Set the refinement type. + Refinement type is used for shaping with C2C contours. + */ + void setRefinementType(RefinementType refinementType) + { + m_refinementType = refinementType; + } + + void setSamplesPerKnotSpan(int nSamples); + void setVertexWeldThreshold(double threshold); + /*! + @brief Set the percentage error tolerance. + + If percent <= MINIMUM_PERCENT_ERROR, the refinement type + will change to DiscreteShape::RefinementUniformSegments. + */ + void setPercentError(double percent); + + //@} + +#if defined(AXOM_USE_MPI) + /** + @brief Set the MPI communicator used when reading C2C files. + */ + void setMPICommunicator(MPI_Comm comm) { m_comm = comm; } +#endif + + /*! + Get the name of this shape. + \return the shape's name + */ + const axom::klee::Shape& getShape() const { return m_shape; } + + /*! + \brief Get the discrete mesh representation. + + If the sidre parent group was used in the constructor, the + mesh data is stored under that group. + + If the discrete mesh isn't generated yet (for analytical shapes), + generate it. + */ + std::shared_ptr createMeshRepresentation(); + + //!@brief Get the discrete mesh representation. + std::shared_ptr getMeshRepresentation() const + { + return m_meshRep; + } + + /*! + \brief Get the revolved volume for volumes of revolution. + */ + double getRevolvedVolume() const { return m_revolvedVolume; } + +private: + const axom::klee::Shape& m_shape; + + /*! + \brief Discrete mesh representation. + + This is either an internally generated mesh representing a + discretized analytical shape or a modifiable copy of the + discrete input geometry for geometries specified as a + discrete mesh. + */ + std::shared_ptr m_meshRep; + + //!@brief Internal DataStore for working space + axom::sidre::DataStore m_dataStore; + + //!@brief Sidre store for m_meshRep. + axom::sidre::Group* m_sidreGroup {nullptr}; + + //!@brief Prefix for disc files with relative path. + std::string m_prefixPath; + + //@{ + //!@name Various parameters for discretization of analytical shapes. + RefinementType m_refinementType; + double m_percentError {MINIMUM_PERCENT_ERROR}; + int m_samplesPerKnotSpan {DEFAULT_SAMPLES_PER_KNOT_SPAN}; + double m_vertexWeldThreshold {DEFAULT_VERTEX_WELD_THRESHOLD}; + double m_revolvedVolume {0.0}; + //@} + +#if defined(AXOM_USE_MPI) + MPI_Comm m_comm {MPI_COMM_WORLD}; +#endif + + void applyTransforms(); + numerics::Matrix getTransforms() const; + + //! @brief Returns the full geometry path. + std::string resolvePath() const; + + /*! + @brief Set the parent group for this object to store data. + */ + void setParentGroup(axom::sidre::Group* parentGroup); + + //!@brief Return a 3x3 matrix that rotate coordinates from the x-axis to the given direction. + numerics::Matrix vorAxisRotMatrix(const Vector3D& dir); + + void clearInternalData(); + +public: + // These are public only for the CUDA device compiler. + void createBlueprintTetsRepresentation(); + void createTetRepresentation(); + void createHexRepresentation(); + void createPlaneRepresentation(); + void createSphereRepresentation(); + void createVORRepresentation(); +}; + +} // namespace quest +} // namespace axom + +#endif diff --git a/src/axom/quest/Discretize.cpp b/src/axom/quest/Discretize.cpp index 065a6aec2f..41abb9600a 100644 --- a/src/axom/quest/Discretize.cpp +++ b/src/axom/quest/Discretize.cpp @@ -208,7 +208,7 @@ double octPolyVolume(const OctType& o) } // namespace //------------------------------------------------------------------------------ -int mesh_from_discretized_polyline(axom::ArrayView& octs, +int mesh_from_discretized_polyline(const axom::ArrayView& octs, int octcount, int segcount, mint::Mesh*& mesh) diff --git a/src/axom/quest/Discretize.hpp b/src/axom/quest/Discretize.hpp index e97efafd4e..b95bcca47e 100644 --- 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 -bool discretize(axom::Array& polyline, +bool discretize(const axom::ArrayView& polyline, int len, int levels, axom::Array& out, @@ -107,7 +107,7 @@ bool discretize(axom::Array& polyline, * the caller is responsible for properly deallocating the mesh object that * the return mesh pointer points to. */ -int mesh_from_discretized_polyline(axom::ArrayView& octs, +int mesh_from_discretized_polyline(const axom::ArrayView& octs, int octcount, int segcount, mint::Mesh*& mesh); diff --git a/src/axom/quest/InOutOctree.hpp b/src/axom/quest/InOutOctree.hpp index ebd99fdf94..5b8992aee1 100644 --- a/src/axom/quest/InOutOctree.hpp +++ b/src/axom/quest/InOutOctree.hpp @@ -175,7 +175,8 @@ class InOutOctree : public spin::SpatialOctree * \note The InOutOctree modifies its mesh in an effort to repair common * problems. Please make sure to discard all old copies of the meshPtr. */ - InOutOctree(const GeometricBoundingBox& bb, SurfaceMesh*& meshPtr) + InOutOctree(const GeometricBoundingBox& bb, + std::shared_ptr& meshPtr) : SpatialOctreeType( GeometricBoundingBox(bb).scale(DEFAULT_BOUNDING_BOX_SCALE_FACTOR)) , m_meshWrapper(meshPtr) diff --git a/src/axom/quest/IntersectionShaper.hpp b/src/axom/quest/IntersectionShaper.hpp index 406c4a264b..573288eba6 100644 --- a/src/axom/quest/IntersectionShaper.hpp +++ b/src/axom/quest/IntersectionShaper.hpp @@ -41,7 +41,7 @@ #endif // clang-format off -#if defined (AXOM_USE_RAJA) && defined (AXOM_USE_UMPIRE) +#if defined (AXOM_USE_RAJA) using seq_exec = axom::SEQ_EXEC; #if defined(AXOM_USE_OPENMP) @@ -50,14 +50,14 @@ using omp_exec = seq_exec; #endif - #if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) && defined (AXOM_USE_UMPIRE) constexpr int CUDA_BLOCK_SIZE = 256; using cuda_exec = axom::CUDA_EXEC; #else using cuda_exec = seq_exec; #endif - #if defined(AXOM_USE_HIP) + #if defined(AXOM_USE_HIP) && defined (AXOM_USE_UMPIRE) constexpr int HIP_BLOCK_SIZE = 64; using hip_exec = axom::HIP_EXEC; #else @@ -299,6 +299,7 @@ class IntersectionShaper : public Shaper using Point3D = primal::Point; using TetrahedronType = primal::Tetrahedron; using SegmentMesh = mint::UnstructuredMesh; + using TetMesh = mint::UnstructuredMesh; using RuntimePolicy = axom::runtime_policy::Policy; @@ -354,7 +355,7 @@ class IntersectionShaper : public Shaper */ double getApproximateRevolvedVolume() const { - return volume(m_surfaceMesh, m_level); + return volume(m_surfaceMesh.get(), m_level); } virtual void loadShape(const klee::Shape& shape) override @@ -366,9 +367,8 @@ class IntersectionShaper : public Shaper if(shape.getGeometry().getFormat() == "c2c") { SegmentMesh* newm = - filterMesh(dynamic_cast(m_surfaceMesh)); - delete m_surfaceMesh; - m_surfaceMesh = newm; + filterMesh(dynamic_cast(m_surfaceMesh.get())); + m_surfaceMesh.reset(newm); } } @@ -376,11 +376,11 @@ class IntersectionShaper : public Shaper //@{ //! @name Functions related to the stages for a given shape -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) - // Prepares the ProE mesh cells for the spatial index + // Prepares the tet mesh mesh cells for the spatial index template - void prepareProECells() + void prepareTetCells() { const int host_allocator = axom::execution_space::allocatorID(); @@ -454,7 +454,7 @@ class IntersectionShaper : public Shaper num_degenerate.get())); // Dump tet mesh as a vtk mesh - axom::mint::write_vtk(m_surfaceMesh, "proe_tet.vtk"); + axom::mint::write_vtk(m_surfaceMesh.get(), "proe_tet.vtk"); } // end of verbose output for Pro/E } @@ -600,18 +600,16 @@ class IntersectionShaper : public Shaper AXOM_UNUSED_VAR(shapeDimension); AXOM_UNUSED_VAR(shapeName); - SLIC_INFO(axom::fmt::format("Current shape is {}", shapeName)); - std::string shapeFormat = shape.getGeometry().getFormat(); + // C2C mesh is not discretized into tets, but all others are. if(shapeFormat == "c2c") { prepareC2CCells(); } - - else if(shapeFormat == "proe") + else if(surfaceMeshIsTet()) { - prepareProECells(); + prepareTetCells(); } else { @@ -621,7 +619,7 @@ class IntersectionShaper : public Shaper } #endif -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) template void runShapeQueryImpl(const klee::Shape& shape, axom::Array& shapes, @@ -999,7 +997,7 @@ class IntersectionShaper : public Shaper } // end of runShapeQuery() function #endif -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) // These methods are private in support of replacement rules. private: /*! @@ -1410,7 +1408,7 @@ class IntersectionShaper : public Shaper switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: applyReplacementRulesImpl(shape); break; @@ -1419,12 +1417,12 @@ class IntersectionShaper : public Shaper applyReplacementRulesImpl(shape); break; #endif // AXOM_USE_OPENMP - #if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::cuda: applyReplacementRulesImpl(shape); break; #endif // AXOM_USE_CUDA - #if defined(AXOM_USE_HIP) + #if defined(AXOM_USE_HIP) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::hip: applyReplacementRulesImpl(shape); break; @@ -1439,9 +1437,7 @@ class IntersectionShaper : public Shaper AXOM_ANNOTATE_SCOPE("finalizeShapeQuery"); // Implementation here -- destroy BVH tree and other shape-based data structures - delete m_surfaceMesh; - - m_surfaceMesh = nullptr; + m_surfaceMesh.reset(); } //@} @@ -1468,7 +1464,7 @@ class IntersectionShaper : public Shaper // Now that the mesh is refined, dispatch to device implementations. switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: prepareShapeQueryImpl(shapeDimension, shape); break; @@ -1477,12 +1473,12 @@ class IntersectionShaper : public Shaper prepareShapeQueryImpl(shapeDimension, shape); break; #endif // AXOM_USE_OPENMP - #if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::cuda: prepareShapeQueryImpl(shapeDimension, shape); break; #endif // AXOM_USE_CUDA - #if defined(AXOM_USE_HIP) + #if defined(AXOM_USE_HIP) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::hip: prepareShapeQueryImpl(shapeDimension, shape); break; @@ -1504,12 +1500,12 @@ class IntersectionShaper : public Shaper AXOM_ANNOTATE_SCOPE("runShapeQuery"); const std::string shapeFormat = shape.getGeometry().getFormat(); - // Testing separate workflow for Pro/E - if(shapeFormat == "proe") + // C2C mesh is not discretized into tets, but all others are. + if(surfaceMeshIsTet()) { switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: runShapeQueryImpl(shape, m_tets, m_tetcount); break; @@ -1518,12 +1514,12 @@ class IntersectionShaper : public Shaper runShapeQueryImpl(shape, m_tets, m_tetcount); break; #endif // AXOM_USE_OPENMP - #if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::cuda: runShapeQueryImpl(shape, m_tets, m_tetcount); break; #endif // AXOM_USE_CUDA - #if defined(AXOM_USE_HIP) + #if defined(AXOM_USE_HIP) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::hip: runShapeQueryImpl(shape, m_tets, m_tetcount); break; @@ -1535,7 +1531,7 @@ class IntersectionShaper : public Shaper { switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: runShapeQueryImpl(shape, m_octs, m_octcount); break; @@ -1544,12 +1540,12 @@ class IntersectionShaper : public Shaper runShapeQueryImpl(shape, m_octs, m_octcount); break; #endif // AXOM_USE_OPENMP - #if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::cuda: runShapeQueryImpl(shape, m_octs, m_octcount); break; #endif // AXOM_USE_CUDA - #if defined(AXOM_USE_HIP) + #if defined(AXOM_USE_HIP) && defined(AXOM_USE_UMPIRE) case RuntimePolicy::hip: runShapeQueryImpl(shape, m_octs, m_octcount); break; @@ -1833,7 +1829,7 @@ class IntersectionShaper : public Shaper { // If we are not refining dynamically, return. if(m_percentError <= MINIMUM_PERCENT_ERROR || - m_refinementType != RefinementDynamic) + m_refinementType != DiscreteShape::RefinementDynamic) { return; } @@ -1921,7 +1917,7 @@ class IntersectionShaper : public Shaper for(int level = m_level; level < MAX_LEVELS; level++) { // Compute the revolved volume of the surface mesh at level. - currentVol = volume(m_surfaceMesh, level); + currentVol = volume(m_surfaceMesh.get(), level); pct = 100. * (1. - currentVol / revolvedVolume); SLIC_INFO( @@ -1994,8 +1990,7 @@ class IntersectionShaper : public Shaper curvePercentError = ce; // Free the previous surface mesh. - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + m_surfaceMesh.reset(); // Reload the shape using new curvePercentError. This will cause // a new m_surfaceMesh to be created. @@ -2008,9 +2003,8 @@ class IntersectionShaper : public Shaper // Filter the mesh, store in m_surfaceMesh. SegmentMesh* newm = - filterMesh(dynamic_cast(m_surfaceMesh)); - delete m_surfaceMesh; - m_surfaceMesh = newm; + filterMesh(dynamic_cast(m_surfaceMesh.get())); + m_surfaceMesh.reset(newm); } else { @@ -2043,6 +2037,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}; @@ -2052,7 +2054,7 @@ class IntersectionShaper : public Shaper axom::Array m_hex_volumes; axom::Array m_overlap_volumes; -#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) +#if defined(AXOM_USE_RAJA) double m_vertexWeldThreshold {1.e-10}; int m_octcount {0}; int m_tetcount {0}; diff --git a/src/axom/quest/MeshViewUtil.hpp b/src/axom/quest/MeshViewUtil.hpp index b8a633090f..03bfabef83 100644 --- a/src/axom/quest/MeshViewUtil.hpp +++ b/src/axom/quest/MeshViewUtil.hpp @@ -198,7 +198,7 @@ static void stridesAndOffsetsToShapes(const axom::StackArray& realSh This class recognizes potential ghost (a.k.a. phony, image) data layers around the domain. Some methods and paramenters names refer to the data with ghosts, while others refer to the data without - hosts (a.k.a. real data). + ghosts (a.k.a. real data). TODO: Figure out if there's a better place for this utility. It's only in axom/quest because the initial need was there. diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index 9c4d6dffcb..360f1f0d79 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -80,14 +80,15 @@ class InOutSampler * * \note Does not take ownership of the surface mesh */ - InOutSampler(const std::string& shapeName, mint::Mesh* surfaceMesh) + InOutSampler(const std::string& shapeName, + std::shared_ptr surfaceMesh) : m_shapeName(shapeName) , m_surfaceMesh(surfaceMesh) { } ~InOutSampler() { delete m_octree; } - mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh; } + std::shared_ptr getSurfaceMesh() const { return m_surfaceMesh; } /// Computes the bounding box of the surface mesh void computeBounds() @@ -303,9 +304,9 @@ class InOutSampler std::string m_shapeName; GeometricBoundingBox m_bbox; - mint::Mesh* m_surfaceMesh {nullptr}; + std::shared_ptr m_surfaceMesh {nullptr}; InOutOctreeType* m_octree {nullptr}; -}; +}; // class InOutSampler } // end namespace shaping @@ -415,14 +416,12 @@ class SamplingShaper : public Shaper m_inoutSampler2D = new shaping::InOutSampler<2>(shapeName, m_surfaceMesh); m_inoutSampler2D->computeBounds(); m_inoutSampler2D->initSpatialIndex(this->m_vertexWeldThreshold); - m_surfaceMesh = m_inoutSampler2D->getSurfaceMesh(); break; case klee::Dimensions::Three: m_inoutSampler3D = new shaping::InOutSampler<3>(shapeName, m_surfaceMesh); m_inoutSampler3D->computeBounds(); m_inoutSampler3D->initSpatialIndex(this->m_vertexWeldThreshold); - m_surfaceMesh = m_inoutSampler3D->getSurfaceMesh(); break; default: @@ -446,7 +445,7 @@ class SamplingShaper : public Shaper "After welding, surface mesh has {} vertices and {} triangles.", nVerts, nCells)); - mint::write_vtk(m_surfaceMesh, + mint::write_vtk(m_surfaceMesh.get(), axom::fmt::format("meldedTriMesh_{}.vtk", shapeName)); } } @@ -584,8 +583,7 @@ class SamplingShaper : public Shaper delete m_inoutSampler3D; m_inoutSampler3D = nullptr; - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + m_surfaceMesh.reset(); } //@} diff --git a/src/axom/quest/Shaper.cpp b/src/axom/quest/Shaper.cpp index cbf6e5b594..c3f0a9ca3e 100644 --- a/src/axom/quest/Shaper.cpp +++ b/src/axom/quest/Shaper.cpp @@ -7,7 +7,9 @@ #include "axom/config.hpp" #include "axom/core.hpp" +#include "axom/primal/operators/split.hpp" #include "axom/quest/interface/internal/QuestHelpers.hpp" +#include "axom/quest/DiscreteShape.hpp" #include "axom/fmt.hpp" @@ -21,60 +23,6 @@ namespace axom { namespace quest { -namespace internal -{ -/*! - * \brief Implementation of a GeometryOperatorVisitor for processing klee shape operators - * - * This class extracts the matrix form of supported operators and marks the operator as unvalid otherwise - * To use, check the \a isValid() function after visiting and then call the \a getMatrix() function. - */ -class AffineMatrixVisitor : public klee::GeometryOperatorVisitor -{ -public: - AffineMatrixVisitor() : m_matrix(4, 4) { } - - void visit(const klee::Translation& translation) override - { - m_matrix = translation.toMatrix(); - m_isValid = true; - } - void visit(const klee::Rotation& rotation) override - { - m_matrix = rotation.toMatrix(); - m_isValid = true; - } - void visit(const klee::Scale& scale) override - { - m_matrix = scale.toMatrix(); - m_isValid = true; - } - void visit(const klee::UnitConverter& converter) override - { - m_matrix = converter.toMatrix(); - m_isValid = true; - } - - void visit(const klee::CompositeOperator&) override - { - SLIC_WARNING("CompositeOperator not supported for Shaper query"); - m_isValid = false; - } - void visit(const klee::SliceOperator&) override - { - SLIC_WARNING("SliceOperator not yet supported for Shaper query"); - m_isValid = false; - } - - const numerics::Matrix& getMatrix() const { return m_matrix; } - bool isValid() const { return m_isValid; } - -private: - bool m_isValid {false}; - numerics::Matrix m_matrix; -}; - -} // end namespace internal // These were needed for linking - but why? They are constexpr. constexpr int Shaper::DEFAULT_SAMPLES_PER_KNOT_SPAN; @@ -130,7 +78,7 @@ void Shaper::setPercentError(double percent) percent)); if(percent <= MINIMUM_PERCENT_ERROR) { - m_refinementType = RefinementUniformSegments; + m_refinementType = DiscreteShape::RefinementUniformSegments; } m_percentError = clampVal(percent, MINIMUM_PERCENT_ERROR, MAXIMUM_PERCENT_ERROR); @@ -144,7 +92,9 @@ void Shaper::setRefinementType(Shaper::RefinementType t) bool Shaper::isValidFormat(const std::string& format) const { return (format == "stl" || format == "proe" || format == "c2c" || - format == "none"); + format == "blueprint-tets" || format == "tet3D" || + format == "hex3D" || format == "plane3D" || format == "sphere3D" || + format == "vor3D" || format == "none"); } void Shaper::loadShape(const klee::Shape& shape) @@ -160,8 +110,6 @@ void Shaper::loadShapeInternal(const klee::Shape& shape, double percentError, double& revolvedVolume) { - using axom::utilities::string::endsWith; - internal::ScopedLogLevelChanger logLevelChanger( this->isVerbose() ? slic::message::Debug : slic::message::Warning); @@ -169,157 +117,22 @@ void Shaper::loadShapeInternal(const klee::Shape& shape, "{:-^80}", axom::fmt::format(" Loading shape '{}' ", shape.getName()))); - const std::string& file_format = shape.getGeometry().getFormat(); + const axom::klee::Geometry& geometry = shape.getGeometry(); + const std::string& geometryFormat = geometry.getFormat(); SLIC_ASSERT_MSG( - this->isValidFormat(file_format), - axom::fmt::format("Shape has unsupported format: '{}", file_format)); - - if(!shape.getGeometry().hasGeometry()) - { - SLIC_DEBUG( - axom::fmt::format("Current shape '{}' of material '{}' has no geometry", - shape.getName(), - shape.getMaterial())); - return; - } - - std::string shapePath = m_shapeSet.resolvePath(shape.getGeometry().getPath()); - SLIC_INFO("Reading file: " << shapePath << "..."); - - // Initialize revolved volume. - revolvedVolume = 0.; - - if(endsWith(shapePath, ".stl")) - { - SLIC_ASSERT_MSG( - file_format == "stl", - axom::fmt::format(" '{}' format requires .stl file type", file_format)); - - quest::internal::read_stl_mesh(shapePath, m_surfaceMesh, m_comm); - // Transform the coordinates of the linearized mesh. - applyTransforms(shape); - } - else if(endsWith(shapePath, ".proe")) - { - SLIC_ASSERT_MSG( - file_format == "proe", - axom::fmt::format(" '{}' format requires .proe file type", file_format)); - - quest::internal::read_pro_e_mesh(shapePath, m_surfaceMesh, m_comm); - } -#ifdef AXOM_USE_C2C - else if(endsWith(shapePath, ".contour")) - { - SLIC_ASSERT_MSG( - file_format == "c2c", - axom::fmt::format(" '{}' format requires .contour file type", file_format)); - - // Get the transforms that are being applied to the mesh. Get them - // as a single concatenated matrix. - auto transform = getTransforms(shape); - - // Pass in the transform so any transformations can figure into - // computing the revolved volume. - if(m_refinementType == RefinementDynamic && - percentError > MINIMUM_PERCENT_ERROR) - { - quest::internal::read_c2c_mesh_non_uniform(shapePath, - transform, - percentError, - m_vertexWeldThreshold, - m_surfaceMesh, - revolvedVolume, // output arg - m_comm); - } - else - { - quest::internal::read_c2c_mesh_uniform(shapePath, - transform, - m_samplesPerKnotSpan, - m_vertexWeldThreshold, - m_surfaceMesh, - revolvedVolume, // output arg - m_comm); - } - - // Transform the coordinates of the linearized mesh. - applyTransforms(transform); - } -#endif - else - { - SLIC_ERROR( - axom::fmt::format("Unsupported filetype for this Axom configuration. " - "Provided file was '{}', with format '{}'", - shapePath, - file_format)); - } -} + this->isValidFormat(geometryFormat), + axom::fmt::format("Shape has unsupported format: '{}", geometryFormat)); -numerics::Matrix Shaper::getTransforms(const klee::Shape& shape) const -{ - const auto identity4x4 = numerics::Matrix::identity(4); - numerics::Matrix transformation(identity4x4); - auto& geometryOperator = shape.getGeometry().getGeometryOperator(); - auto composite = - std::dynamic_pointer_cast(geometryOperator); - if(composite) + // Code for discretizing shapes has been factored into DiscreteShape class. + DiscreteShape discreteShape(shape, m_dataStore.getRoot(), m_shapeSet.getPath()); + discreteShape.setVertexWeldThreshold(m_vertexWeldThreshold); + discreteShape.setRefinementType(m_refinementType); + if(percentError > 0) { - // Concatenate the transformations - for(auto op : composite->getOperators()) - { - // Use visitor pattern to extract the affine matrix from supported operators - internal::AffineMatrixVisitor visitor; - op->accept(visitor); - if(!visitor.isValid()) - { - continue; - } - const auto& matrix = visitor.getMatrix(); - numerics::Matrix res(identity4x4); - numerics::matrix_multiply(matrix, transformation, res); - transformation = res; - } - } - return transformation; -} - -void Shaper::applyTransforms(const klee::Shape& shape) -{ - // Concatenate the transformations - numerics::Matrix transformation = getTransforms(shape); - // Transform the surface mesh coordinates. - applyTransforms(transformation); -} - -void Shaper::applyTransforms(const numerics::Matrix& transformation) -{ - AXOM_ANNOTATE_SCOPE("applyTransforms"); - - // Apply transformation to coordinates of each vertex in mesh - if(!transformation.isIdentity()) - { - const int spaceDim = m_surfaceMesh->getDimension(); - const int numSurfaceVertices = m_surfaceMesh->getNumberOfNodes(); - double* x = m_surfaceMesh->getCoordinateArray(mint::X_COORDINATE); - double* y = m_surfaceMesh->getCoordinateArray(mint::Y_COORDINATE); - double* z = spaceDim > 2 - ? m_surfaceMesh->getCoordinateArray(mint::Z_COORDINATE) - : nullptr; - - double xformed[4]; - for(int i = 0; i < numSurfaceVertices; ++i) - { - double coords[4] = {x[i], y[i], (z == nullptr ? 0. : z[i]), 1.}; - numerics::matrix_vector_multiply(transformation, coords, xformed); - x[i] = xformed[0]; - y[i] = xformed[1]; - if(z != nullptr) - { - z[i] = xformed[2]; - } - } + discreteShape.setPercentError(percentError); } + m_surfaceMesh = discreteShape.createMeshRepresentation(); + revolvedVolume = discreteShape.getRevolvedVolume(); } // ---------------------------------------------------------------------------- diff --git a/src/axom/quest/Shaper.hpp b/src/axom/quest/Shaper.hpp index 01942de782..cf8ee7d27f 100644 --- a/src/axom/quest/Shaper.hpp +++ b/src/axom/quest/Shaper.hpp @@ -23,6 +23,7 @@ #include "axom/sidre.hpp" #include "axom/klee.hpp" #include "axom/mint.hpp" +#include "axom/quest/DiscreteShape.hpp" #include "axom/quest/interface/internal/mpicomm_wrapper.hpp" @@ -48,7 +49,7 @@ class Shaper static constexpr double DEFAULT_VERTEX_WELD_THRESHOLD {1e-9}; /// Refinement type. - using RefinementType = enum { RefinementUniformSegments, RefinementDynamic }; + using RefinementType = DiscreteShape::RefinementType; //@{ //! @name Functions to get and set shaping parameters @@ -64,7 +65,7 @@ class Shaper bool isVerbose() const { return m_verboseOutput; } sidre::MFEMSidreDataCollection* getDC() { return m_dc; } - mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh; } + mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh.get(); } /*! * \brief Predicate to determine if the specified format is valid @@ -108,8 +109,8 @@ class Shaper protected: /*! - * \brief Loads the shape from file into m_surfaceMesh and computes a revolvedVolume - * for the shape. + * \brief Loads the shape from file into m_surfaceMesh and, if its a C2C + * contour, computes a revolvedVolume for the shape. * \param shape The shape. * \param percentError A percent error to use when refining the shape. If it * positive then Axom will try to refine dynamically @@ -151,14 +152,16 @@ class Shaper int getRank() const; protected: + sidre::DataStore m_dataStore; + const klee::ShapeSet& m_shapeSet; sidre::MFEMSidreDataCollection* m_dc; - mint::Mesh* m_surfaceMesh {nullptr}; + std::shared_ptr m_surfaceMesh; int m_samplesPerKnotSpan {DEFAULT_SAMPLES_PER_KNOT_SPAN}; double m_percentError {MINIMUM_PERCENT_ERROR}; - RefinementType m_refinementType {RefinementUniformSegments}; + RefinementType m_refinementType {DiscreteShape::RefinementUniformSegments}; double m_vertexWeldThreshold {DEFAULT_VERTEX_WELD_THRESHOLD}; bool m_verboseOutput {false}; diff --git a/src/axom/quest/detail/Discretize_detail.hpp b/src/axom/quest/detail/Discretize_detail.hpp index 9fdb582ba2..e24521abf6 100644 --- 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 -bool discretize(axom::Array &polyline, +bool discretize(const axom::ArrayView &polyline, int pointcount, int levels, axom::Array &out, @@ -276,8 +276,8 @@ bool discretize(axom::Array &polyline, int segmentcount = pointcount - 1; for(int seg = 0; seg < segmentcount && stillValid; ++seg) { - Point2D &a = polyline[seg]; - Point2D &b = polyline[seg + 1]; + const Point2D &a = polyline[seg]; + const Point2D &b = polyline[seg + 1]; // invalid if a.x > b.x if(a[0] > b[0]) { diff --git a/src/axom/quest/detail/inout/MeshWrapper.hpp b/src/axom/quest/detail/inout/MeshWrapper.hpp index 1f72ab339d..72cf8244ed 100644 --- a/src/axom/quest/detail/inout/MeshWrapper.hpp +++ b/src/axom/quest/detail/inout/MeshWrapper.hpp @@ -81,7 +81,7 @@ class SimplexMeshWrapper static constexpr CellIndex NO_CELL = -1; protected: - SimplexMeshWrapper(SurfaceMesh*& meshPtr) + SimplexMeshWrapper(std::shared_ptr& meshPtr) : m_surfaceMesh(meshPtr) , m_vertexPositions(&m_vertexSet) , m_cellToVertexRelation() @@ -278,7 +278,7 @@ class SimplexMeshWrapper } protected: - SurfaceMesh*& m_surfaceMesh; // ref to pointer to allow changing the mesh + std::shared_ptr& m_surfaceMesh; // ref to pointer to allow changing the mesh MeshVertexSet m_vertexSet {0}; MeshElementSet m_elementSet {0}; @@ -312,7 +312,7 @@ class MeshWrapper<2> : public SimplexMeshWrapper<2, MeshWrapper<2>> public: /// \brief Constructor for a mesh wrapper */ - MeshWrapper(SurfaceMesh*& meshPtr) : Base(meshPtr) { } + MeshWrapper(std::shared_ptr& meshPtr) : Base(meshPtr) { } /** * \brief Helper function to compute the bounding box of an edge @@ -419,7 +419,7 @@ class MeshWrapper<2> : public SimplexMeshWrapper<2, MeshWrapper<2>> // Grab relation from mesh using UMesh = mint::UnstructuredMesh; axom::IndexType* vertIds = - static_cast(m_surfaceMesh)->getCellNodeIDs(i); + std::static_pointer_cast(m_surfaceMesh)->getCellNodeIDs(i); // Remap the vertex IDs for(int j = 0; j < NUM_EDGE_VERTS; ++j) @@ -442,9 +442,8 @@ class MeshWrapper<2> : public SimplexMeshWrapper<2, MeshWrapper<2>> m_cellToVertexRelation.bindIndices(static_cast(m_cv_data.size()), &m_cv_data); - // Delete old mesh, and NULL its pointer - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + // Delete old mesh + m_surfaceMesh.reset(); m_meshWasReindexed = true; } @@ -454,8 +453,7 @@ class MeshWrapper<2> : public SimplexMeshWrapper<2, MeshWrapper<2>> { if(m_surfaceMesh != nullptr) { - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + m_surfaceMesh.reset(); } using UMesh = mint::UnstructuredMesh; @@ -476,7 +474,7 @@ class MeshWrapper<2> : public SimplexMeshWrapper<2, MeshWrapper<2>> edgeMesh->appendCell(tv); } - m_surfaceMesh = edgeMesh; + m_surfaceMesh.reset(edgeMesh); } }; @@ -500,7 +498,7 @@ class MeshWrapper<3> : public SimplexMeshWrapper<3, MeshWrapper<3>> public: /// \brief Constructor for a mesh wrapper */ - MeshWrapper(SurfaceMesh*& meshPtr) : Base(meshPtr) { } + MeshWrapper(std::shared_ptr& meshPtr) : Base(meshPtr) { } /** * \brief Helper function to compute the bounding box of a triangle @@ -570,7 +568,7 @@ class MeshWrapper<3> : public SimplexMeshWrapper<3, MeshWrapper<3>> // Grab relation from mesh using UMesh = mint::UnstructuredMesh; axom::IndexType* vertIds = - static_cast(m_surfaceMesh)->getCellNodeIDs(i); + std::static_pointer_cast(m_surfaceMesh)->getCellNodeIDs(i); // Remap the vertex IDs for(int j = 0; j < NUM_TRI_VERTS; ++j) @@ -595,9 +593,8 @@ class MeshWrapper<3> : public SimplexMeshWrapper<3, MeshWrapper<3>> m_cellToVertexRelation.bindIndices(static_cast(m_cv_data.size()), &m_cv_data); - // Delete old mesh, and NULL its pointer - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + // Delete old mesh + m_surfaceMesh.reset(); m_meshWasReindexed = true; } @@ -607,8 +604,7 @@ class MeshWrapper<3> : public SimplexMeshWrapper<3, MeshWrapper<3>> { if(m_surfaceMesh != nullptr) { - delete m_surfaceMesh; - m_surfaceMesh = nullptr; + m_surfaceMesh.reset(); } using UMesh = mint::UnstructuredMesh; @@ -629,7 +625,7 @@ class MeshWrapper<3> : public SimplexMeshWrapper<3, MeshWrapper<3>> triMesh->appendCell(tv); } - m_surfaceMesh = triMesh; + m_surfaceMesh.reset(triMesh); } }; diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 4841fa0fe5..9835592167 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -264,6 +264,50 @@ if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI endif() endif() +# Shaping in-memory example ------------------------------------------------------ +if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI + AND AXOM_ENABLE_SIDRE AND AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION + AND AXOM_ENABLE_KLEE AND RAJA_FOUND) + axom_add_executable( + NAME quest_shape_in_memory_ex + SOURCES quest_shape_in_memory.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON ${quest_example_depends} mfem + FOLDER axom/quest/examples + ) + + if(AXOM_ENABLE_TESTS AND AXOM_DATA_DIR) + # 3D shaping tests. No 2D shaping in this feature, at least for now. + set(_nranks 1) + + # Run the in-memory shaping example on with each enabled policy + set(_policies "seq") + if(RAJA_FOUND) + blt_list_append(TO _policies ELEMENTS "omp" IF AXOM_ENABLE_OPENMP) + blt_list_append(TO _policies ELEMENTS "cuda" IF AXOM_ENABLE_CUDA) + blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) + endif() + + set(_testshapes "tetmesh" "tet" "hex" "sphere" "cyl" "cone" "vor" "all" "plane") + + foreach(_policy ${_policies}) + foreach(_testshape ${_testshapes}) + set(_testname "quest_shape_in_memory_${_policy}_${_testshape}") + axom_add_test( + NAME ${_testname} + COMMAND quest_shape_in_memory_ex + --policy ${_policy} + --testShape ${_testshape} + --refinements 2 + --scale 0.75 0.75 0.75 + --dir 0.2 0.4 0.8 + inline_mesh --min -2 -2 -2 --max 2 2 2 --res 30 30 30 + NUM_MPI_TASKS ${_nranks}) + endforeach() + endforeach() + endif() +endif() + # Distributed closest point example ------------------------------------------- if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) axom_add_executable( diff --git a/src/axom/quest/examples/containment_driver.cpp b/src/axom/quest/examples/containment_driver.cpp index 4606b33ab7..7ef006f68d 100644 --- a/src/axom/quest/examples/containment_driver.cpp +++ b/src/axom/quest/examples/containment_driver.cpp @@ -64,7 +64,7 @@ class ContainmentDriver ~ContainmentDriver() { - delete m_surfaceMesh; + m_surfaceMesh.reset(); delete m_octree; } @@ -77,8 +77,8 @@ class ContainmentDriver reader.read(); // Create surface mesh - m_surfaceMesh = new UMesh(2, mint::SEGMENT); - reader.getLinearMeshUniform(static_cast(m_surfaceMesh), + m_surfaceMesh.reset(new UMesh(2, mint::SEGMENT)); + reader.getLinearMeshUniform(static_cast(m_surfaceMesh.get()), segmentsPerKnotSpan); } #else @@ -100,11 +100,11 @@ class ContainmentDriver reader.read(); // Create surface mesh - m_surfaceMesh = new UMesh(3, mint::TRIANGLE); - reader.getMesh(static_cast(m_surfaceMesh)); + m_surfaceMesh.reset(new UMesh(3, mint::TRIANGLE)); + reader.getMesh(static_cast(m_surfaceMesh.get())); } - mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh; } + mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh.get(); } int dimension() const { return DIM; } @@ -482,7 +482,7 @@ class ContainmentDriver } private: - mint::Mesh* m_surfaceMesh {nullptr}; + std::shared_ptr m_surfaceMesh {nullptr}; InOutOctreeType* m_octree {nullptr}; GeometricBoundingBox m_meshBB; GeometricBoundingBox m_queryBB; diff --git a/src/axom/quest/examples/quest_shape_in_memory.cpp b/src/axom/quest/examples/quest_shape_in_memory.cpp new file mode 100644 index 0000000000..19ea3d7435 --- /dev/null +++ b/src/axom/quest/examples/quest_shape_in_memory.cpp @@ -0,0 +1,1410 @@ +// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file quest_shape_mint_mesh.cpp + * \brief Driver application for shaping material volume fractions onto a simulation mesh + * using a mint::UnstructuredMesh of tets. + * Modeled after shaping_driver.cpp. + */ + +// Axom includes +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/mint.hpp" +#include "axom/primal.hpp" +#include "axom/sidre.hpp" +#include "axom/klee.hpp" +#include "axom/quest.hpp" + +#include "axom/fmt.hpp" +#include "axom/CLI11.hpp" + +#include "conduit_relay_io_blueprint.hpp" + +// NOTE: The shaping driver requires Axom to be configured with mfem as well as +// the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION CMake option +#ifndef AXOM_USE_MFEM + #error Shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option +#endif + +#include "mfem.hpp" + +#ifdef AXOM_USE_MPI + #include "mpi.h" +#endif + +// RAJA +#ifdef AXOM_USE_RAJA + #include "RAJA/RAJA.hpp" +#endif + +// C/C++ includes +#include +#include +#include + +namespace klee = axom::klee; +namespace primal = axom::primal; +namespace quest = axom::quest; +namespace slic = axom::slic; +namespace sidre = axom::sidre; + +//------------------------------------------------------------------------------ + +using RuntimePolicy = axom::runtime_policy::Policy; + +/// Struct to parse and store the input parameters +// Some parameters are used to override defaults. +struct Input +{ +public: + std::string outputFile; + + std::vector center; + double radius {-1.0}; + double radius2 {-0.3}; + double length {-2.0}; + std::vector direction; + + // Shape transformation parameters + std::vector scaleFactors; + + // Mesh format: mfem or blueprint + std::string meshFormat = "mfem"; + + // Inline mesh parameters + std::vector boxMins {-2, -2, -2}; + std::vector boxMaxs {2, 2, 2}; + std::vector boxResolution {20, 20, 20}; + int getBoxDim() const + { + auto d = boxResolution.size(); + SLIC_ASSERT(boxMins.size() == d); + SLIC_ASSERT(boxMaxs.size() == d); + return int(d); + } + + // The shape to run. + std::string testShape {"tetmesh"}; + // The shapes this example is set up to run. + const std::set availableShapes {"tetmesh", + "sphere", + "cyl", + "cone", + "vor", + "tet", + "hex", + "plane", + "all"}; + + RuntimePolicy policy {RuntimePolicy::seq}; + int outputOrder {2}; + int refinementLevel {7}; + double weldThresh {1e-9}; + double percentError {-1.}; + std::string annotationMode {"none"}; + + std::string backgroundMaterial; + +private: + bool m_verboseOutput {false}; + +public: + bool isVerbose() const { return m_verboseOutput; } + + /// @brief Return volume of input box mesh + double boxMeshVolume() const + { + primal::Vector x {boxMaxs[0] - boxMins[0], 0, 0}; + primal::Vector y {0, boxMaxs[1] - boxMins[1], 0}; + primal::Vector z {0, 0, boxMaxs[2] - boxMins[2]}; + double volume = primal::Vector::scalar_triple_product(x, y, z); + return volume; + } + + /// Generate an mfem Cartesian mesh, scaled to the bounding box range + mfem::Mesh* createBoxMesh() + { + mfem::Mesh* mesh = nullptr; + + switch(getBoxDim()) + { + case 2: + { + using BBox2D = primal::BoundingBox; + using Pt2D = primal::Point; + auto res = primal::NumericArray(boxResolution.data()); + auto bbox = BBox2D(Pt2D(boxMins.data()), Pt2D(boxMaxs.data())); + + SLIC_INFO(axom::fmt::format( + "Creating inline box mesh of resolution {} and bounding box {}", + res, + bbox)); + + mesh = quest::util::make_cartesian_mfem_mesh_2D(bbox, res, outputOrder); + } + break; + case 3: + { + using BBox3D = primal::BoundingBox; + using Pt3D = primal::Point; + auto res = primal::NumericArray(boxResolution.data()); + auto bbox = BBox3D(Pt3D(boxMins.data()), Pt3D(boxMaxs.data())); + + SLIC_INFO(axom::fmt::format( + "Creating inline box mesh of resolution {} and bounding box {}", + res, + bbox)); + + mesh = quest::util::make_cartesian_mfem_mesh_3D(bbox, res, outputOrder); + } + break; + default: + SLIC_ERROR("Only 2D and 3D meshes are currently supported."); + break; + } + + // Handle conversion to parallel mfem mesh +#if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) + { + int* partitioning = nullptr; + int part_method = 0; + mfem::Mesh* parallelMesh = + new mfem::ParMesh(MPI_COMM_WORLD, *mesh, partitioning, part_method); + delete[] partitioning; + delete mesh; + mesh = parallelMesh; + } +#endif + + return mesh; + } + + std::unique_ptr loadComputationalMesh() + { + constexpr bool dc_owns_data = true; + mfem::Mesh* mesh = createBoxMesh(); + std::string name = "mesh"; + + auto dc = std::unique_ptr( + new sidre::MFEMSidreDataCollection(name, mesh, dc_owns_data)); + dc->SetComm(MPI_COMM_WORLD); + + return dc; + } + + void parse(int argc, char** argv, axom::CLI::App& app) + { + app.add_option("-o,--outputFile", outputFile) + ->description("Path to output file(s)"); + + app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) + ->description("Enable/disable verbose output") + ->capture_default_str(); + + app.add_option("-t,--weld-threshold", weldThresh) + ->description("Threshold for welding") + ->check(axom::CLI::NonNegativeNumber) + ->capture_default_str(); + + app.add_option("-e,--percent-error", percentError) + ->description( + "Percent error used for calculating curve refinement and revolved " + "volume.\n" + "If this value is provided then dynamic curve refinement will be used\n" + "instead of segment-based curve refinement.") + ->check(axom::CLI::PositiveNumber) + ->capture_default_str(); + + app.add_option("-s,--testShape", testShape) + ->description("The shape to run") + ->check(axom::CLI::IsMember(availableShapes)); + +#ifdef AXOM_USE_CALIPER + app.add_option("--caliper", annotationMode) + ->description( + "caliper annotation mode. Valid options include 'none' and 'report'. " + "Use 'help' to see full list.") + ->capture_default_str() + ->check(axom::utilities::ValidCaliperMode); +#endif + + app.add_option("--center", center) + ->description("Center of sphere or base of cone/cyl/VOR (x,y[,z]) shape") + ->expected(2, 3); + + app.add_option("--radius", radius) + ->description("Radius of sphere or cylinder shape") + ->check(axom::CLI::PositiveNumber); + + app.add_option("--length", length) + ->description("Length of cone/cyl/VOR shape") + ->check(axom::CLI::PositiveNumber); + + app.add_option("--dir", direction) + ->description( + "Direction of axis of rotation (cone/cyl/VOR (x,y[,z])), or rotated " + "x-axis (hex, tet, tetmesh, and sphere), or positive normal direction " + "(plane).") + ->expected(2, 3); + + app.add_option("--radius2", radius2) + ->description("Second radius of cone shape") + ->check(axom::CLI::PositiveNumber); + + app.add_option("--scale", scaleFactors) + ->description("Scale factor to apply to shape (x,y[,z])") + ->expected(2, 3) + ->check(axom::CLI::PositiveNumber); + + // use either an input mesh file or a simple inline Cartesian mesh + { + auto* inline_mesh_subcommand = + app.add_subcommand("inline_mesh") + ->description("Options for setting up a simple inline mesh") + ->fallthrough(); + + inline_mesh_subcommand->add_option("--min", boxMins) + ->description("Min bounds for box mesh (x,y[,z])") + ->expected(2, 3) + ->required(); + inline_mesh_subcommand->add_option("--max", boxMaxs) + ->description("Max bounds for box mesh (x,y[,z])") + ->expected(2, 3) + ->required(); + + inline_mesh_subcommand->add_option("--res", boxResolution) + ->description("Resolution of the box mesh (i,j[,k])") + ->expected(2, 3) + ->required(); + } + + app.add_option("--background-material", backgroundMaterial) + ->description("Sets the name of the background material"); + + // parameters that only apply to the intersection method + { + auto* intersection_options = + app.add_option_group("intersection", + "Options related to intersection-based queries"); + + intersection_options->add_option("-r, --refinements", refinementLevel) + ->description("Number of refinements to perform for revolved contour") + ->capture_default_str() + ->check(axom::CLI::NonNegativeNumber); + + std::stringstream pol_sstr; + pol_sstr << "Set runtime policy for intersection-based sampling method."; +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) + pol_sstr << "\nSet to 'seq' or 0 to use the RAJA sequential policy."; + #ifdef AXOM_USE_OPENMP + pol_sstr << "\nSet to 'omp' or 1 to use the RAJA OpenMP policy."; + #endif + #ifdef AXOM_USE_CUDA + pol_sstr << "\nSet to 'cuda' or 2 to use the RAJA CUDA policy."; + #endif + #ifdef AXOM_USE_HIP + pol_sstr << "\nSet to 'hip' or 3 to use the RAJA HIP policy."; + #endif +#endif + + intersection_options->add_option("-p, --policy", policy, pol_sstr.str()) + ->capture_default_str() + ->transform( + axom::CLI::CheckedTransformer(axom::runtime_policy::s_nameToPolicy)); + } + app.get_formatter()->column_width(50); + + // could throw an exception + app.parse(argc, argv); + + slic::setLoggingMsgLevel(m_verboseOutput ? slic::message::Debug + : slic::message::Info); + } +}; // struct Input +Input params; + +// Start property for all 3D shapes. +axom::klee::TransformableGeometryProperties startProp { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + +// Add scale operator if specified by input parameters. +void addScaleOperator(axom::klee::CompositeOperator& compositeOp) +{ + SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); + if(!params.scaleFactors.empty()) + { + std::shared_ptr scaleOp = + std::make_shared(params.scaleFactors[0], + params.scaleFactors[1], + params.scaleFactors[2], + startProp); + compositeOp.addOperator(scaleOp); + } +} + +// Add translate operator. +void addTranslateOperator(axom::klee::CompositeOperator& compositeOp, + double shiftx, + double shifty, + double shiftz) +{ + primal::Vector3D shift({shiftx, shifty, shiftz}); + auto translateOp = std::make_shared(shift, startProp); + compositeOp.addOperator(translateOp); +} + +// Add operator to rotate x-axis to params.direction, if it is given. +void addRotateOperator(axom::klee::CompositeOperator& compositeOp) +{ + if(!params.direction.empty()) + { + static const primal::Point3D center {0.0, 0.0, 0.0}; + static const primal::Vector3D x {1.0, 0.0, 0.0}; + primal::Vector3D rotateTo(params.direction.data()); + // Note that the rotation matrix is not unique. + primal::Vector3D a = rotateTo.unitVector(); + primal::Vector3D u; // Rotation vector, the cross product of x and a. + axom::numerics::cross_product(x.data(), a.data(), u.data()); + double angle = asin(u.norm()) * 180 / M_PI; + + auto rotateOp = + std::make_shared(angle, center, u, startProp); + compositeOp.addOperator(rotateOp); + } +} + +/** + * \brief Print some info about the mesh + * + * \note In MPI-based configurations, this is a collective call, but + * only prints on rank 0 + */ +void printMfemMeshInfo(mfem::Mesh* mesh, const std::string& prefixMessage = "") +{ + namespace primal = axom::primal; + + int myRank = 0; +#ifdef AXOM_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &myRank); +#endif + + int numElements = mesh->GetNE(); + + mfem::Vector mins, maxs; +#ifdef MFEM_USE_MPI + auto* parallelMesh = dynamic_cast(mesh); + if(parallelMesh != nullptr) + { + parallelMesh->GetBoundingBox(mins, maxs); + numElements = parallelMesh->ReduceInt(numElements); + myRank = parallelMesh->GetMyRank(); + } + else +#endif + { + mesh->GetBoundingBox(mins, maxs); + } + + if(myRank == 0) + { + switch(mesh->Dimension()) + { + case 2: + SLIC_INFO(axom::fmt::format( + axom::utilities::locale(), + "{} mesh has {:L} elements and (approximate) bounding box {}", + prefixMessage, + numElements, + primal::BoundingBox(primal::Point(mins.GetData()), + primal::Point(maxs.GetData())))); + break; + case 3: + SLIC_INFO(axom::fmt::format( + axom::utilities::locale(), + "{} mesh has {:L} elements and (approximate) bounding box {}", + prefixMessage, + numElements, + primal::BoundingBox(primal::Point(mins.GetData()), + primal::Point(maxs.GetData())))); + break; + } + } + + slic::flushStreams(); +} + +/// \brief Utility function to initialize the logger +void initializeLogger() +{ + // Initialize Logger + slic::initialize(); + slic::setLoggingMsgLevel(slic::message::Info); + + slic::LogStream* logStream {nullptr}; + +#ifdef AXOM_USE_MPI + int num_ranks = 1; + MPI_Comm_size(MPI_COMM_WORLD, &num_ranks); + if(num_ranks > 1) + { + std::string fmt = "[][]: \n"; + #ifdef AXOM_USE_LUMBERJACK + const int RLIMIT = 8; + logStream = + new slic::LumberjackStream(&std::cout, MPI_COMM_WORLD, RLIMIT, fmt); + #else + logStream = new slic::SynchronizedStream(&std::cout, MPI_COMM_WORLD, fmt); + #endif + } + else +#endif // AXOM_USE_MPI + { + std::string fmt = "[]: \n"; + logStream = new slic::GenericOutputStream(&std::cout, fmt); + } + + slic::addStreamToAllMsgLevels(logStream); +} + +/// \brief Utility function to finalize the logger +void finalizeLogger() +{ + if(slic::isInitialized()) + { + slic::flushStreams(); + slic::finalize(); + } +} + +// Single triangle ShapeSet. +axom::klee::ShapeSet create2DShapeSet(sidre::DataStore& ds) +{ + sidre::Group* meshGroup = ds.getRoot()->createGroup("triangleMesh"); + AXOM_UNUSED_VAR(meshGroup); // variable is only referenced in debug configs + const std::string topo = "mesh"; + const std::string coordset = "coords"; + axom::mint::UnstructuredMesh triangleMesh( + 2, + axom::mint::CellType::TRIANGLE, + meshGroup, + topo, + coordset); + + double lll = 2.0; + + // Insert tet at origin. + triangleMesh.appendNode(0.0, 0.0); + triangleMesh.appendNode(lll, 0.0); + triangleMesh.appendNode(0.0, lll); + axom::IndexType conn[3] = {0, 1, 2}; + triangleMesh.appendCell(conn); + + meshGroup->print(); + SLIC_ASSERT(axom::mint::blueprint::isValidRootGroup(meshGroup)); + + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Two, + axom::klee::LengthUnit::unspecified}; + axom::klee::Geometry triangleGeom(prop, + triangleMesh.getSidreGroup(), + topo, + nullptr); + + std::vector shapes; + axom::klee::Shape triangleShape( + "triangle", + "AL", + {}, + {}, + axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}); + shapes.push_back(axom::klee::Shape { + "triangle", + "AL", + {}, + {}, + axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}}); + + axom::klee::ShapeSet shapeSet; + shapeSet.setShapes(shapes); + shapeSet.setDimensions(axom::klee::Dimensions::Two); + + return shapeSet; +} + +axom::klee::Shape createShape_Sphere() +{ + Point3D center = + params.center.empty() ? Point3D {0, 0, 0} : Point3D {params.center.data()}; + double radius = params.radius < 0 ? 0.6 : params.radius; + axom::primal::Sphere sphere {center, radius}; + + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp, 1, 1, 1); + + const axom::IndexType levelOfRefinement = params.refinementLevel; + axom::klee::Geometry sphereGeometry(prop, sphere, levelOfRefinement, compositeOp); + axom::klee::Shape sphereShape("sphere", "AU", {}, {}, sphereGeometry); + + return sphereShape; +} + +axom::klee::Shape createShape_TetMesh(sidre::DataStore& ds) +{ + // Shape a tetrahedal mesh. + sidre::Group* meshGroup = ds.getRoot()->createGroup("tetMesh"); + AXOM_UNUSED_VAR(meshGroup); // variable is only referenced in debug configs + const std::string topo = "mesh"; + const std::string coordset = "coords"; + axom::mint::UnstructuredMesh tetMesh( + 3, + axom::mint::CellType::TET, + meshGroup, + topo, + coordset); + + double lll = params.length < 0 ? 0.7 : params.length; + + // Insert tets around origin. + tetMesh.appendNode(-lll, -lll, -lll); + tetMesh.appendNode(+lll, -lll, -lll); + tetMesh.appendNode(-lll, +lll, -lll); + tetMesh.appendNode(-lll, -lll, +lll); + tetMesh.appendNode(+lll, +lll, +lll); + tetMesh.appendNode(-lll, +lll, +lll); + tetMesh.appendNode(+lll, +lll, -lll); + tetMesh.appendNode(+lll, -lll, +lll); + axom::IndexType conn0[4] = {0, 1, 2, 3}; + tetMesh.appendCell(conn0); + axom::IndexType conn1[4] = {4, 5, 6, 7}; + tetMesh.appendCell(conn1); + + SLIC_ASSERT(axom::mint::blueprint::isValidRootGroup(meshGroup)); + + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp, -1, 1, 1); + + axom::klee::Geometry tetMeshGeometry(prop, + tetMesh.getSidreGroup(), + topo, + compositeOp); + axom::klee::Shape tetShape("tetmesh", "TETMESH", {}, {}, tetMeshGeometry); + + return tetShape; +} + +axom::klee::Geometry createGeometry_Vor( + axom::primal::Point& vorBase, + axom::primal::Vector& vorDirection, + axom::Array& discreteFunction, + std::shared_ptr& compositeOp) +{ + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); + + const axom::IndexType levelOfRefinement = params.refinementLevel; + axom::klee::Geometry vorGeometry(prop, + discreteFunction, + vorBase, + vorDirection, + levelOfRefinement, + compositeOp); + return vorGeometry; +} + +axom::klee::Shape createShape_Vor() +{ + Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + : Point3D {params.center.data()}; + axom::primal::Vector vorDirection = params.direction.empty() + ? primal::Vector3D {0.1, 0.2, 0.4} + : primal::Vector3D {params.direction.data()}; + const int numIntervals = 5; + axom::Array discreteFunction({numIntervals + 1, 2}, + axom::ArrayStrideOrder::ROW); + double zLen = params.length < 0 ? 1.6 : params.length; + double zShift = -zLen / 2; + double maxR = params.radius < 0 ? 0.75 : params.radius; + double dz = zLen / numIntervals; + discreteFunction[0][0] = 0 * dz + zShift; + discreteFunction[0][1] = 0.0 * maxR; + discreteFunction[1][0] = 1 * dz + zShift; + discreteFunction[1][1] = 0.8 * maxR; + discreteFunction[2][0] = 2 * dz + zShift; + discreteFunction[2][1] = 0.4 * maxR; + discreteFunction[3][0] = 3 * dz + zShift; + discreteFunction[3][1] = 0.5 * maxR; + discreteFunction[4][0] = 4 * dz + zShift; + discreteFunction[4][1] = 1.0 * maxR; + discreteFunction[5][0] = 5 * dz + zShift; + discreteFunction[5][1] = 0.0; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp, -1, -1, 1); + + axom::klee::Geometry vorGeometry = + createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + + axom::klee::Shape vorShape("vor", "VOR", {}, {}, vorGeometry); + + return vorShape; +} + +axom::klee::Shape createShape_Cylinder() +{ + Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + : Point3D {params.center.data()}; + axom::primal::Vector vorDirection = params.direction.empty() + ? primal::Vector3D {0.1, 0.2, 0.4} + : primal::Vector3D {params.direction.data()}; + axom::Array discreteFunction({2, 2}, axom::ArrayStrideOrder::ROW); + double radius = params.radius < 0 ? 0.5 : params.radius; + double height = params.length < 0 ? 1.2 : params.length; + discreteFunction[0][0] = -height / 2; + discreteFunction[0][1] = radius; + discreteFunction[1][0] = height / 2; + discreteFunction[1][1] = radius; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp, 1, -1, 1); + + axom::klee::Geometry vorGeometry = + createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + + axom::klee::Shape vorShape("cyl", "CYL", {}, {}, vorGeometry); + + return vorShape; +} + +axom::klee::Shape createShape_Cone() +{ + Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + : Point3D {params.center.data()}; + axom::primal::Vector vorDirection = params.direction.empty() + ? primal::Vector3D {0.1, 0.2, 0.4} + : primal::Vector3D {params.direction.data()}; + axom::Array discreteFunction({2, 2}, axom::ArrayStrideOrder::ROW); + double baseRadius = params.radius < 0 ? 0.7 : params.radius; + double topRadius = params.radius2 < 0 ? 0.1 : params.radius2; + double height = params.length < 0 ? 1.3 : params.length; + discreteFunction[0][0] = -height / 2; + discreteFunction[0][1] = baseRadius; + discreteFunction[1][0] = height / 2; + discreteFunction[1][1] = topRadius; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp, 1, 1, -1); + + axom::klee::Geometry vorGeometry = + createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + + axom::klee::Shape vorShape("cone", "CONE", {}, {}, vorGeometry); + + return vorShape; +} + +axom::klee::Shape createShape_Tet() +{ + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); + + // Tetrahedron at origin. + const double len = params.length < 0 ? 0.8 : params.length; + const Point3D a {-len, -len, -len}; + const Point3D b {+len, -len, -len}; + const Point3D c {+len, +len, -len}; + const Point3D d {-len, +len, +len}; + const primal::Tetrahedron tet {a, b, c, d}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp, -1, 1, -1); + + axom::klee::Geometry tetGeometry(prop, tet, compositeOp); + axom::klee::Shape tetShape("tet", "TET", {}, {}, tetGeometry); + + return tetShape; +} + +axom::klee::Shape createShape_Hex() +{ + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); + + const double md = params.length < 0 ? 0.6 : params.length; + const double lg = 1.2 * md; + const double sm = 0.8 * md; + const Point3D p {-lg, -md, -sm}; + const Point3D q {+lg, -md, -sm}; + const Point3D r {+lg, +md, -sm}; + const Point3D s {-lg, +md, -sm}; + const Point3D t {-lg, -md, +sm}; + const Point3D u {+lg, -md, +sm}; + const Point3D v {+lg, +md, +sm}; + const Point3D w {-lg, +md, +sm}; + const primal::Hexahedron hex {p, q, r, s, t, u, v, w}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp, -1, -1, -1); + + axom::klee::Geometry hexGeometry(prop, hex, compositeOp); + axom::klee::Shape hexShape("hex", "HEX", {}, {}, hexGeometry); + + return hexShape; +} + +axom::klee::Shape createShape_Plane() +{ + axom::klee::TransformableGeometryProperties prop { + axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); + std::shared_ptr scaleOp; + if(!params.scaleFactors.empty()) + { + scaleOp = std::make_shared(params.scaleFactors[0], + params.scaleFactors[1], + params.scaleFactors[2], + prop); + } + + // Create a plane crossing center of mesh. No matter the normal, + // it cuts the mesh in half. + Point3D center {0.5 * + (primal::NumericArray(params.boxMins.data()) + + primal::NumericArray(params.boxMaxs.data()))}; + primal::Vector normal = params.direction.empty() + ? primal::Vector3D {1.0, 0.0, 0.0} + : primal::Vector3D {params.direction.data()}.unitVector(); + const primal::Plane plane {normal, center, true}; + + axom::klee::Geometry planeGeometry(prop, plane, scaleOp); + axom::klee::Shape planeShape("plane", "PLANE", {}, {}, planeGeometry); + + return planeShape; +} + +//!@brief Create a ShapeSet with a single shape. +axom::klee::ShapeSet createShapeSet(const axom::klee::Shape& shape) +{ + axom::klee::ShapeSet shapeSet; + shapeSet.setShapes(std::vector {shape}); + shapeSet.setDimensions(axom::klee::Dimensions::Three); + + return shapeSet; +} + +double volumeOfTetMesh( + const axom::mint::UnstructuredMesh& tetMesh) +{ + using TetType = axom::primal::Tetrahedron; + if(0) + { + std::ofstream os("tets.js"); + tetMesh.getSidreGroup()->print(os); + } + axom::StackArray nodesShape {tetMesh.getNumberOfNodes()}; + axom::ArrayView x(tetMesh.getCoordinateArray(0), nodesShape); + axom::ArrayView y(tetMesh.getCoordinateArray(1), nodesShape); + axom::ArrayView z(tetMesh.getCoordinateArray(2), nodesShape); + const axom::IndexType cellCount = tetMesh.getNumberOfCells(); + axom::Array tetVolumes(cellCount, cellCount); + double meshVolume = 0.0; + for(axom::IndexType ic = 0; ic < cellCount; ++ic) + { + const axom::IndexType* nodeIds = tetMesh.getCellNodeIDs(ic); + TetType tet; + for(int j = 0; j < 4; ++j) + { + auto cornerNodeId = nodeIds[j]; + tet[j][0] = x[cornerNodeId]; + tet[j][1] = y[cornerNodeId]; + tet[j][2] = z[cornerNodeId]; + } + meshVolume += tet.volume(); + } + return meshVolume; +} + +/*! + @brief Return the element volumes as a sidre::View. + + If it doesn't exist, allocate and compute it. + \post The volume data is in \c dc->GetNamedBuffer(volFieldName). + + Most of this is lifted from IntersectionShaper::runShapeQueryImpl. +*/ +template +axom::sidre::View* getElementVolumes( + sidre::MFEMSidreDataCollection* dc, + const std::string& volFieldName = std::string("elementVolumes")) +{ + using HexahedronType = axom::primal::Hexahedron; + + axom::sidre::View* volSidreView = dc->GetNamedBuffer(volFieldName); + if(volSidreView == nullptr) + { + mfem::Mesh* mesh = dc->GetMesh(); + const axom::IndexType cellCount = mesh->GetNE(); + + constexpr int NUM_VERTS_PER_HEX = 8; + constexpr int NUM_COMPS_PER_VERT = 3; + constexpr double ZERO_THRESHOLD = 1.e-10; + + axom::Array vertCoords(cellCount * NUM_VERTS_PER_HEX, + cellCount * NUM_VERTS_PER_HEX); + auto vertCoordsView = vertCoords.view(); + + // This runs only only on host, because the mfem::Mesh only uses host memory, I think. + for(axom::IndexType cellIdx = 0; cellIdx < cellCount; ++cellIdx) + { + // Get the indices of this element's vertices + mfem::Array verts; + mesh->GetElementVertices(cellIdx, verts); + SLIC_ASSERT(verts.Size() == NUM_VERTS_PER_HEX); + + // Get the coordinates for the vertices + for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) + { + int vertIdx = cellIdx * NUM_VERTS_PER_HEX + j; + for(int k = 0; k < NUM_COMPS_PER_VERT; k++) + { + vertCoordsView[vertIdx][k] = (mesh->GetVertex(verts[j]))[k]; + } + } + } + + // Set vertex coords to zero if within threshold. + // (I don't know why we do this. I'm following examples.) + axom::ArrayView flatCoordsView( + (double*)vertCoords.data(), + vertCoords.size() * Point3D::dimension()); + assert(flatCoordsView.size() == cellCount * NUM_VERTS_PER_HEX * 3); + axom::for_all( + cellCount * 3, + AXOM_LAMBDA(axom::IndexType i) { + if(axom::utilities::isNearlyEqual(flatCoordsView[i], 0.0, ZERO_THRESHOLD)) + { + flatCoordsView[i] = 0.0; + } + }); + + // Initialize hexahedral elements. + axom::Array hexes(cellCount, cellCount); + auto hexesView = hexes.view(); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellIdx) { + // Set each hexahedral element vertices + hexesView[cellIdx] = HexahedronType(); + for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) + { + int vertIndex = (cellIdx * NUM_VERTS_PER_HEX) + j; + auto& hex = hexesView[cellIdx]; + hex[j] = vertCoordsView[vertIndex]; + } + }); // end of loop to initialize hexahedral elements and bounding boxes + + // Allocate and populate cell volumes. + volSidreView = dc->AllocNamedBuffer(volFieldName, cellCount); + axom::ArrayView volView(volSidreView->getData(), + volSidreView->getNumElements()); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellIdx) { + volView[cellIdx] = hexesView[cellIdx].volume(); + }); + } + + return volSidreView; +} + +/*! + @brief Return global sum of volume of the given material. +*/ +template +double sumMaterialVolumes(sidre::MFEMSidreDataCollection* dc, + const std::string& material) +{ + mfem::Mesh* mesh = dc->GetMesh(); + int const cellCount = mesh->GetNE(); + + // Get cell volumes from dc. + axom::sidre::View* elementVols = getElementVolumes(dc); + axom::ArrayView elementVolsView(elementVols->getData(), + elementVols->getNumElements()); + + // Get material volume fractions + const std::string materialFieldName = + axom::fmt::format("vol_frac_{}", material); + mfem::GridFunction* volFracGf = dc->GetField(materialFieldName); + axom::quest::GridFunctionView volFracView(volFracGf); + + using ReducePolicy = typename axom::execution_space::reduce_policy; + RAJA::ReduceSum localVol(0); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType i) { + localVol += volFracView[i] * elementVolsView[i]; + }); + + double globalVol = localVol.get(); +#ifdef AXOM_USE_MPI + MPI_Allreduce(MPI_IN_PLACE, &globalVol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + return globalVol; +} + +//------------------------------------------------------------------------------ +int main(int argc, char** argv) +{ + axom::utilities::raii::MPIWrapper mpi_raii_wrapper(argc, argv); + const int my_rank = mpi_raii_wrapper.my_rank(); + + initializeLogger(); + + //--------------------------------------------------------------------------- + // Set up and parse command line arguments + //--------------------------------------------------------------------------- + axom::CLI::App app {"Driver for Klee shaping query"}; + + try + { + params.parse(argc, argv, app); + } + catch(const axom::CLI::ParseError& e) + { + int retval = -1; + if(my_rank == 0) + { + retval = app.exit(e); + } + finalizeLogger(); + +#ifdef AXOM_USE_MPI + MPI_Bcast(&retval, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Finalize(); +#endif + exit(retval); + } + + axom::utilities::raii::AnnotationsWrapper annotations_raii_wrapper( + params.annotationMode); + + AXOM_ANNOTATE_BEGIN("quest example for shaping primals"); + AXOM_ANNOTATE_BEGIN("init"); + + // Storage for the shape geometry meshes. + sidre::DataStore ds; + + //--------------------------------------------------------------------------- + // Create simple ShapeSet for the example. + //--------------------------------------------------------------------------- + axom::klee::ShapeSet shapeSet; + SLIC_ERROR_IF(params.getBoxDim() != 3, "This example is only in 3D."); + switch(params.getBoxDim()) + { + case 2: + shapeSet = create2DShapeSet(ds); + break; + case 3: + if(params.testShape == "tetmesh") + { + shapeSet = createShapeSet(createShape_TetMesh(ds)); + } + else if(params.testShape == "tet") + { + shapeSet = createShapeSet(createShape_Tet()); + } + else if(params.testShape == "hex") + { + shapeSet = createShapeSet(createShape_Hex()); + } + else if(params.testShape == "sphere") + { + shapeSet = createShapeSet(createShape_Sphere()); + } + else if(params.testShape == "cyl") + { + shapeSet = createShapeSet(createShape_Cylinder()); + } + else if(params.testShape == "cone") + { + shapeSet = createShapeSet(createShape_Cone()); + } + else if(params.testShape == "vor") + { + shapeSet = createShapeSet(createShape_Vor()); + } + else if(params.testShape == "plane") + { + shapeSet = createShapeSet(createShape_Plane()); + } + else if(params.testShape == "all") + { + std::vector shapesVec; + shapesVec.push_back(createShape_TetMesh(ds)); + shapesVec.push_back(createShape_Tet()); + shapesVec.push_back(createShape_Hex()); + shapesVec.push_back(createShape_Sphere()); + shapesVec.push_back(createShape_Vor()); + shapesVec.push_back(createShape_Cylinder()); + shapesVec.push_back(createShape_Cone()); + shapeSet.setShapes(shapesVec); + shapeSet.setDimensions(axom::klee::Dimensions::Three); + } + break; + } + + // Save the discrete shapes for viz and testing. + auto* shapeMeshGroup = ds.getRoot()->createGroup("shapeMeshGroup"); + std::vector> discreteShapeMeshes; + for(const auto& shape : shapeSet.getShapes()) + { + axom::quest::DiscreteShape dShape(shape, shapeMeshGroup); + auto dMesh = + std::dynamic_pointer_cast>( + dShape.createMeshRepresentation()); + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Shape '{}' discrete geometry has {} cells", + shape.getName(), + dMesh->getNumberOfCells()))); + + discreteShapeMeshes.push_back(dMesh); + + if(!params.outputFile.empty()) + { + std::string shapeFileName = params.outputFile + ".shape"; + conduit::Node tmpNode, info; + dMesh->getSidreGroup()->createNativeLayout(tmpNode); + conduit::relay::io::blueprint::save_mesh(tmpNode, shapeFileName, "hdf5"); + } + } + + const klee::Dimensions shapeDim = shapeSet.getDimensions(); + + // Apply error checking +#ifndef AXOM_USE_C2C + SLIC_ERROR_IF(shapeDim == klee::Dimensions::Two, + "Shaping with contour files requires an Axom configured with " + "the C2C library"); +#endif + + AXOM_ANNOTATE_BEGIN("load mesh"); + //--------------------------------------------------------------------------- + // Load the computational mesh + // originalMeshDC is the input MFEM mesh + // It's converted to shapingDC below, and that is used for shaping calls. + //--------------------------------------------------------------------------- + auto originalMeshDC = params.loadComputationalMesh(); + + //--------------------------------------------------------------------------- + // Set up DataCollection for shaping + // shapingDC is a "copy" of originalMeshDC, the MFEM mesh. + // It's created empty then populated by the SetMesh call. + // shapingMesh and parallelMesh are some kind of temporary versions of originalMeshDC. + //--------------------------------------------------------------------------- + mfem::Mesh* shapingMesh = nullptr; + constexpr bool dc_owns_data = true; + sidre::MFEMSidreDataCollection shapingDC("shaping", shapingMesh, dc_owns_data); + { + shapingDC.SetMeshNodesName("positions"); + + // With MPI, loadComputationalMesh returns a parallel mesh. + mfem::ParMesh* parallelMesh = + dynamic_cast(originalMeshDC->GetMesh()); + shapingMesh = (parallelMesh != nullptr) + ? new mfem::ParMesh(*parallelMesh) + : new mfem::Mesh(*originalMeshDC->GetMesh()); + shapingDC.SetMesh(shapingMesh); + } + AXOM_ANNOTATE_END("load mesh"); + printMfemMeshInfo(shapingDC.GetMesh(), "After loading"); + + const axom::IndexType cellCount = shapingMesh->GetNE(); + + // TODO Port to GPUs. Shaper should be data-parallel, but data may not be on devices yet. + using ExecSpace = typename axom::SEQ_EXEC; + using ReducePolicy = typename axom::execution_space::reduce_policy; + + //--------------------------------------------------------------------------- + // Initialize the shaping query object + //--------------------------------------------------------------------------- + AXOM_ANNOTATE_BEGIN("setup shaping problem"); + quest::Shaper* shaper = nullptr; + shaper = new quest::IntersectionShaper(shapeSet, &shapingDC); + SLIC_ASSERT_MSG(shaper != nullptr, "Invalid shaping method selected!"); + + // Set generic parameters for the base Shaper instance + shaper->setVertexWeldThreshold(params.weldThresh); + shaper->setVerbosity(params.isVerbose()); + if(params.percentError > 0.) + { + shaper->setPercentError(params.percentError); + shaper->setRefinementType(quest::DiscreteShape::RefinementDynamic); + } + + // Associate any fields that begin with "vol_frac" with "material" so when + // the data collection is written, a matset will be created. + shaper->getDC()->AssociateMaterialSet("vol_frac", "material"); + + // Set specific parameters here for IntersectionShaper + if(auto* intersectionShaper = dynamic_cast(shaper)) + { + intersectionShaper->setLevel(params.refinementLevel); + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Setting IntersectionShaper policy to '{}'", + axom::runtime_policy::policyToName(params.policy)))); + intersectionShaper->setExecPolicy(params.policy); + + if(!params.backgroundMaterial.empty()) + { + intersectionShaper->setFreeMaterialName(params.backgroundMaterial); + } + } + + AXOM_ANNOTATE_END("setup shaping problem"); + AXOM_ANNOTATE_END("init"); + + //--------------------------------------------------------------------------- + // Process each of the shapes + //--------------------------------------------------------------------------- + SLIC_INFO(axom::fmt::format("{:=^80}", "Shaping loop")); + AXOM_ANNOTATE_BEGIN("shaping"); + for(const auto& shape : shapeSet.getShapes()) + { + const std::string shapeFormat = shape.getGeometry().getFormat(); + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Processing shape '{}' of material '{}' (format '{}')", + shape.getName(), + shape.getMaterial(), + shapeFormat))); + + // Load the shape from file. This also applies any transformations. + shaper->loadShape(shape); + slic::flushStreams(); + + // Generate a spatial index over the shape + shaper->prepareShapeQuery(shapeDim, shape); + slic::flushStreams(); + + // Query the mesh against this shape + shaper->runShapeQuery(shape); + slic::flushStreams(); + + // Apply the replacement rules for this shape against the existing materials + shaper->applyReplacementRules(shape); + slic::flushStreams(); + + // Finalize data structures associated with this shape and spatial index + shaper->finalizeShapeQuery(); + slic::flushStreams(); + } + AXOM_ANNOTATE_END("shaping"); + + //--------------------------------------------------------------------------- + // After shaping in all shapes, generate/adjust the material volume fractions + //--------------------------------------------------------------------------- + AXOM_ANNOTATE_BEGIN("adjust"); + SLIC_INFO( + axom::fmt::format("{:=^80}", + "Generating volume fraction fields for materials")); + + shaper->adjustVolumeFractions(); + + //--------------------------------------------------------------------------- + // Compute and print volumes of each material's volume fraction + //--------------------------------------------------------------------------- + using axom::utilities::string::startsWith; + auto dc = shaper->getDC(); + auto& fieldMap = dc->GetFieldMap(); + for(auto& kv : fieldMap) + { + if(startsWith(kv.first, "vol_frac_")) + { + const auto mat_name = kv.first.substr(9); + auto* gf = kv.second; + + mfem::ConstantCoefficient one(1.0); + mfem::LinearForm vol_form(gf->FESpace()); + vol_form.AddDomainIntegrator(new mfem::DomainLFIntegrator(one)); + vol_form.Assemble(); + + const double volume = shaper->allReduceSum(*gf * vol_form); + + SLIC_INFO(axom::fmt::format(axom::utilities::locale(), + "Volume of material '{}' is {:.6Lf}", + mat_name, + volume)); + } + } + AXOM_ANNOTATE_END("adjust"); + + int failCounts = 0; + + auto matsetsGrp = shapingDC.GetBPGroup()->getGroup("matsets"); + auto materialGrp = matsetsGrp->getGroup("material"); + auto volFracGroups = materialGrp->getGroup("volume_fractions"); + + //--------------------------------------------------------------------------- + // Correctness test: volume fractions should be in [0,1]. + //--------------------------------------------------------------------------- + RAJA::ReduceSum rangeViolationCount(0); + for(axom::sidre::Group& materialGroup : volFracGroups->groups()) + { + axom::sidre::View* values = materialGroup.getView("value"); + double* volFracData = values->getArray(); + axom::ArrayView volFracDataView(volFracData, cellCount); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType i) { + bool bad = volFracDataView[i] < 0.0 || volFracDataView[i] > 1.0; + rangeViolationCount += bad; + }); + } + + failCounts += (rangeViolationCount.get() != 0); + + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Count of volume fractions outside of [0,1]: {}.", + rangeViolationCount.get()))); + slic::flushStreams(); + + //--------------------------------------------------------------------------- + // Correctness test: volume fractions in each cell should sum to 1.0. + //--------------------------------------------------------------------------- + axom::Array volSums(cellCount); + volSums.fill(0.0); + axom::ArrayView volSumsView = volSums.view(); + for(axom::sidre::Group& materialGroup : volFracGroups->groups()) + { + axom::sidre::View* values = materialGroup.getView("value"); + double* volFracData = values->getArray(); + axom::ArrayView volFracDataView(volFracData, cellCount); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType i) { volSumsView[i] += volFracDataView[i]; }); + } + RAJA::ReduceSum nonUnitSums(0); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType i) { + bool bad = !axom::utilities::isNearlyEqual(volSums[i], 0.0); + nonUnitSums += bad; + }); + + failCounts += (nonUnitSums.get() != 0); + + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Count non-unit volume fraction sums: {}.", + nonUnitSums.get()))); + slic::flushStreams(); + + //--------------------------------------------------------------------------- + // Correctness test: shape volume in shapingMesh should match volume of the + // shape mesh for closes shape. As long as the shapes don't overlap, this + // should be a good correctness check. + //--------------------------------------------------------------------------- + auto* meshVerificationGroup = ds.getRoot()->createGroup("meshVerification"); + for(const auto& shape : shapeSet.getShapes()) + { + axom::quest::DiscreteShape dShape(shape, meshVerificationGroup); + auto shapeMesh = + std::dynamic_pointer_cast>( + dShape.createMeshRepresentation()); + double shapeMeshVol = volumeOfTetMesh(*shapeMesh); + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Shape '{}' discrete geometry has {} cells", + shape.getName(), + shapeMesh->getNumberOfCells()))); + + const std::string& materialName = shape.getMaterial(); + double shapeVol = + sumMaterialVolumes(&shapingDC, materialName); + double correctShapeVol = + params.testShape == "plane" ? params.boxMeshVolume() / 2 : shapeMeshVol; + double diff = shapeVol - correctShapeVol; + + bool err = !axom::utilities::isNearlyEqual(shapeVol, correctShapeVol); + failCounts += err; + + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format( + "Material '{}' in shape '{}' has volume {} vs {}, diff of {}, {}.", + materialName, + shape.getName(), + shapeVol, + correctShapeVol, + diff, + (err ? "ERROR" : "OK")))); + } + slic::flushStreams(); + ds.getRoot()->destroyGroupAndData("meshVerification"); + + //--------------------------------------------------------------------------- + // Save meshes and fields + //--------------------------------------------------------------------------- + +#ifdef MFEM_USE_MPI + if(!params.outputFile.empty()) + { + std::string fileName = params.outputFile + ".volfracs"; + shaper->getDC()->Save(fileName, sidre::Group::getDefaultIOProtocol()); + SLIC_INFO(axom::fmt::format("{:=^80}", "Wrote output mesh " + fileName)); + } +#endif + + delete shaper; + + //--------------------------------------------------------------------------- + // Cleanup and exit + //--------------------------------------------------------------------------- + SLIC_INFO(axom::fmt::format("{:-^80}", "")); + slic::flushStreams(); + + AXOM_ANNOTATE_END("quest shaping example"); + + finalizeLogger(); + + return failCounts; +} diff --git a/src/axom/quest/examples/shaping_driver.cpp b/src/axom/quest/examples/shaping_driver.cpp index 77badcd544..31abd26041 100644 --- a/src/axom/quest/examples/shaping_driver.cpp +++ b/src/axom/quest/examples/shaping_driver.cpp @@ -578,7 +578,7 @@ int main(int argc, char** argv) if(params.percentError > 0.) { shaper->setPercentError(params.percentError); - shaper->setRefinementType(quest::Shaper::RefinementDynamic); + shaper->setRefinementType(quest::DiscreteShape::RefinementDynamic); } // Associate any fields that begin with "vol_frac" with "material" so when diff --git a/src/axom/quest/interface/inout.cpp b/src/axom/quest/interface/inout.cpp index af2521b840..05316e103b 100644 --- a/src/axom/quest/interface/inout.cpp +++ b/src/axom/quest/interface/inout.cpp @@ -92,7 +92,7 @@ struct InOutHelper */ int initialize(const std::string& file, MPI_Comm comm) { - mint::Mesh* mesh = nullptr; + mint::Mesh* tmpMeshPtr = nullptr; m_params.m_dimension = getDimension(); // load the mesh @@ -106,7 +106,7 @@ struct InOutHelper numerics::Matrix::identity(4), m_params.m_segmentsPerKnotSpan, m_params.m_vertexWeldThreshold, - mesh, + tmpMeshPtr, revolvedVolume, comm); #else @@ -116,12 +116,14 @@ struct InOutHelper #endif break; case 3: - rc = internal::read_stl_mesh(file, mesh, comm); + rc = internal::read_stl_mesh(file, tmpMeshPtr, comm); break; default: // no-op break; } + std::shared_ptr mesh(tmpMeshPtr); + if(rc != QUEST_INOUT_SUCCESS) { SLIC_WARNING("reading mesh from [" << file << "] failed!"); @@ -138,7 +140,7 @@ struct InOutHelper * * \sa inout_init */ - int initialize(mint::Mesh*& mesh, MPI_Comm comm) + int initialize(std::shared_ptr& mesh, MPI_Comm comm) { // initialize logger, if necessary internal::logger_init(m_state.m_logger_is_initialized, @@ -221,7 +223,7 @@ struct InOutHelper // deal with mesh if(m_state.m_should_delete_mesh) { - delete m_surfaceMesh; + m_surfaceMesh.reset(); } m_surfaceMesh = nullptr; @@ -294,7 +296,7 @@ struct InOutHelper } private: - mint::Mesh* m_surfaceMesh {nullptr}; + std::shared_ptr m_surfaceMesh {nullptr}; InOutOctree* m_inoutTree {nullptr}; GeometricBoundingBox m_meshBoundingBox; SpacePt m_meshCenterOfMass; @@ -360,7 +362,7 @@ int inout_init(const std::string& file, MPI_Comm comm) return rc; } -int inout_init(mint::Mesh*& mesh, MPI_Comm comm) +int inout_init(std::shared_ptr& mesh, MPI_Comm comm) { const int dim = inout_get_dimension(); int rc = QUEST_INOUT_FAILED; diff --git a/src/axom/quest/interface/inout.hpp b/src/axom/quest/interface/inout.hpp index e2d8284a3b..20ecbace65 100644 --- a/src/axom/quest/interface/inout.hpp +++ b/src/axom/quest/interface/inout.hpp @@ -14,6 +14,7 @@ // C/C++ includes #include +#include /*! * \file inout.hpp @@ -78,7 +79,7 @@ int inout_init(const std::string& file, MPI_Comm comm = MPI_COMM_SELF); * by welding vertices) and updates the \a mesh pointer. It is the user's * responsibility to update any other pointers to this same mesh. */ -int inout_init(mint::Mesh*& mesh, MPI_Comm comm = MPI_COMM_SELF); +int inout_init(std::shared_ptr& mesh, MPI_Comm comm = MPI_COMM_SELF); /*! * \brief Finalizes the inout query diff --git a/src/axom/quest/tests/quest_initialize.cpp b/src/axom/quest/tests/quest_initialize.cpp index cc78cff609..ed18bc6310 100644 --- a/src/axom/quest/tests/quest_initialize.cpp +++ b/src/axom/quest/tests/quest_initialize.cpp @@ -29,7 +29,8 @@ TEST(quest_initialize, inout_pointer_initialize) { int rc = quest::QUEST_INOUT_SUCCESS; - mint::Mesh* input_mesh = axom::quest::utilities::make_tetrahedron_mesh(); + std::shared_ptr input_mesh { + axom::quest::utilities::make_tetrahedron_mesh()}; // Note: the following call updates the input_mesh pointer #ifdef AXOM_USE_MPI @@ -47,8 +48,6 @@ TEST(quest_initialize, inout_pointer_initialize) rc = quest::inout_finalize(); EXPECT_EQ(quest::QUEST_INOUT_SUCCESS, rc); - - delete input_mesh; } // Test initializing quest signed_distance from a preloaded mesh diff --git a/src/axom/quest/tests/quest_inout_interface.cpp b/src/axom/quest/tests/quest_inout_interface.cpp index 25d3e337d8..3643a2b7ed 100644 --- a/src/axom/quest/tests/quest_inout_interface.cpp +++ b/src/axom/quest/tests/quest_inout_interface.cpp @@ -137,7 +137,7 @@ TYPED_TEST(InOutInterfaceTest, initialize_from_mesh) EXPECT_TRUE(axom::utilities::filesystem::pathExists(this->meshfile)); - axom::mint::Mesh* mesh {nullptr}; + axom::mint::Mesh* tmpMeshPtr {nullptr}; int rc = failCode; @@ -152,15 +152,18 @@ TYPED_TEST(InOutInterfaceTest, initialize_from_mesh) identity, segmentsPerKnotSpan, weldThreshold, - mesh, + tmpMeshPtr, revolvedVolume); #endif // AXOM_USE_C2C } else // DIM == 3 { - rc = axom::quest::internal::read_stl_mesh(this->meshfile, mesh); + rc = axom::quest::internal::read_stl_mesh(this->meshfile, tmpMeshPtr); } + std::shared_ptr mesh {tmpMeshPtr}; + tmpMeshPtr = nullptr; + EXPECT_EQ(0, rc); ASSERT_NE(nullptr, mesh); @@ -179,8 +182,6 @@ TYPED_TEST(InOutInterfaceTest, initialize_from_mesh) // InOut should no longer be initialized EXPECT_FALSE(axom::quest::inout_initialized()); - - delete mesh; } TYPED_TEST(InOutInterfaceTest, query_properties) diff --git a/src/axom/quest/tests/quest_inout_octree.cpp b/src/axom/quest/tests/quest_inout_octree.cpp index 61a9f45a53..4a028360d7 100644 --- a/src/axom/quest/tests/quest_inout_octree.cpp +++ b/src/axom/quest/tests/quest_inout_octree.cpp @@ -46,18 +46,18 @@ using BlockIndex = Octree3D::BlockIndex; #endif /// Returns a SpacePt corresponding to the given vertex id \a vIdx in \a mesh -SpacePt getVertex(axom::mint::Mesh*& mesh, int vIdx) +SpacePt getVertex(axom::mint::Mesh& mesh, int vIdx) { SpacePt pt; - mesh->getNode(vIdx, pt.data()); + mesh.getNode(vIdx, pt.data()); return pt; } -GeometricBoundingBox computeBoundingBox(axom::mint::Mesh*& mesh) +GeometricBoundingBox computeBoundingBox(axom::mint::Mesh& mesh) { GeometricBoundingBox bbox; - for(int i = 0; i < mesh->getNumberOfNodes(); ++i) + for(int i = 0; i < mesh.getNumberOfNodes(); ++i) { bbox.addPoint(getVertex(mesh, i)); } @@ -66,7 +66,8 @@ GeometricBoundingBox computeBoundingBox(axom::mint::Mesh*& mesh) } /// Runs randomized inout queries on an octahedron mesh -void queryOctahedronMesh(axom::mint::Mesh*& mesh, const GeometricBoundingBox& bbox) +void queryOctahedronMesh(std::shared_ptr& mesh, + const GeometricBoundingBox& bbox) { const double bbMin = bbox.getMin()[0]; const double bbMax = bbox.getMax()[0]; @@ -96,7 +97,7 @@ void queryOctahedronMesh(axom::mint::Mesh*& mesh, const GeometricBoundingBox& bb case 3: case 4: case 5: - pt = getVertex(mesh, i); + pt = getVertex(*mesh, i); break; case 6: case 7: @@ -111,9 +112,9 @@ void queryOctahedronMesh(axom::mint::Mesh*& mesh, const GeometricBoundingBox& bb GridPt vertInds; mesh->getCellNodeIDs(tIdx, vertInds.data()); - pt = axom::quest::utilities::getCentroid(getVertex(mesh, vertInds[0]), - getVertex(mesh, vertInds[1]), - getVertex(mesh, vertInds[2])); + pt = axom::quest::utilities::getCentroid(getVertex(*mesh, vertInds[0]), + getVertex(*mesh, vertInds[1]), + getVertex(*mesh, vertInds[2])); } break; @@ -135,8 +136,8 @@ void queryOctahedronMesh(axom::mint::Mesh*& mesh, const GeometricBoundingBox& bb const int v2[] = {1, 2, 3, 4, 2, 3, 4, 1, 1, 2, 3, 4}; int eIdx = (i - 14); - pt = axom::quest::utilities::getCentroid(getVertex(mesh, v1[eIdx]), - getVertex(mesh, v2[eIdx])); + pt = axom::quest::utilities::getCentroid(getVertex(*mesh, v1[eIdx]), + getVertex(*mesh, v2[eIdx])); } break; @@ -182,7 +183,8 @@ TEST(quest_inout_octree, octahedron_mesh) << " and tests point containment.\n"); // Generate the InOutOctree - axom::mint::Mesh* mesh = axom::quest::utilities::make_octahedron_mesh(); + std::shared_ptr mesh( + axom::quest::utilities::make_octahedron_mesh()); // axom::mint::write_vtk(mesh, "octahedron.vtk"); /// @@ -203,9 +205,6 @@ TEST(quest_inout_octree, octahedron_mesh) SLIC_INFO("Testing InOutOctree on octahedron mesh with shifted bounding box " << bbox2); queryOctahedronMesh(mesh, bbox2); - - delete mesh; - mesh = nullptr; } TEST(quest_inout_octree, tetrahedron_mesh) @@ -219,24 +218,22 @@ TEST(quest_inout_octree, tetrahedron_mesh) for(auto thresh : thresholds) { - mint::Mesh* mesh = quest::utilities::make_tetrahedron_mesh(); - GeometricBoundingBox bbox = computeBoundingBox(mesh); + std::shared_ptr mesh {quest::utilities::make_tetrahedron_mesh()}; + GeometricBoundingBox bbox = computeBoundingBox(*mesh); Octree3D octree(bbox, mesh); octree.setVertexWeldThreshold(thresh); octree.generateIndex(); - SpacePt queryInside = quest::utilities::getCentroid(getVertex(mesh, 0), - getVertex(mesh, 1), - getVertex(mesh, 2), - getVertex(mesh, 3)); + SpacePt queryInside = quest::utilities::getCentroid(getVertex(*mesh, 0), + getVertex(*mesh, 1), + getVertex(*mesh, 2), + getVertex(*mesh, 3)); SpacePt queryOutside = SpacePt(2. * bbox.getMax().array()); EXPECT_TRUE(octree.within(queryInside)); EXPECT_FALSE(octree.within(queryOutside)); - - delete mesh; } } diff --git a/src/axom/quest/tests/quest_inout_quadtree.cpp b/src/axom/quest/tests/quest_inout_quadtree.cpp index d25b0338cc..bbf08c3ed9 100644 --- a/src/axom/quest/tests/quest_inout_quadtree.cpp +++ b/src/axom/quest/tests/quest_inout_quadtree.cpp @@ -40,18 +40,18 @@ using BlockIndex = Octree2D::BlockIndex; } // namespace /// Returns a SpacePt corresponding to the given vertex id \a vIdx in \a mesh -SpacePt getVertex(axom::mint::Mesh*& mesh, int vIdx) +SpacePt getVertex(axom::mint::Mesh& mesh, int vIdx) { SpacePt pt; - mesh->getNode(vIdx, pt.data()); + mesh.getNode(vIdx, pt.data()); return pt; } -GeometricBoundingBox computeBoundingBox(axom::mint::Mesh*& mesh) +GeometricBoundingBox computeBoundingBox(axom::mint::Mesh& mesh) { GeometricBoundingBox bbox; - for(int i = 0; i < mesh->getNumberOfNodes(); ++i) + for(int i = 0; i < mesh.getNumberOfNodes(); ++i) { bbox.addPoint(getVertex(mesh, i)); } @@ -71,9 +71,10 @@ TEST(quest_inout_quadtree, triangle_boundary_mesh) for(auto thresh : thresholds) { // create a simple mesh on the boundary of an equilateral triangle - mint::Mesh* mesh = [=]() { - auto* mesh = - new mint::UnstructuredMesh(DIM, mint::SEGMENT); + std::shared_ptr mesh = [=]() { + auto mesh = std::make_shared>( + DIM, + mint::SEGMENT); mesh->appendNode(1, 1); mesh->appendNode(4, 1); mesh->appendNode(2.5, 3 * sqrt(3) / 2.); @@ -86,22 +87,20 @@ TEST(quest_inout_quadtree, triangle_boundary_mesh) return mesh; }(); - GeometricBoundingBox bbox = computeBoundingBox(mesh); + GeometricBoundingBox bbox = computeBoundingBox(*mesh); Octree2D octree(bbox, mesh); octree.setVertexWeldThreshold(thresh); octree.generateIndex(); - SpacePt queryInside = quest::utilities::getCentroid(getVertex(mesh, 0), - getVertex(mesh, 1), - getVertex(mesh, 2)); + SpacePt queryInside = quest::utilities::getCentroid(getVertex(*mesh, 0), + getVertex(*mesh, 1), + getVertex(*mesh, 2)); SpacePt queryOutside = SpacePt(2. * bbox.getMax().array()); EXPECT_TRUE(octree.within(queryInside)); EXPECT_FALSE(octree.within(queryOutside)); - - delete mesh; } } @@ -117,11 +116,11 @@ TEST(quest_inout_quadtree, circle_mesh) for(double radius : {1. / 3., 1., sqrt(2.), 1234.5678}) { ASSERT_TRUE(num_segments >= 3); - mint::Mesh* mesh = - quest::utilities::make_circle_mesh_2d(radius, num_segments); + std::shared_ptr mesh { + quest::utilities::make_circle_mesh_2d(radius, num_segments)}; //mint::write_vtk(mesh,axom::fmt::format("circle_mesh_r{:.3f}_s{:06}.vtk", radius, num_segments)); - GeometricBoundingBox bbox = computeBoundingBox(mesh).scale(1.2); + GeometricBoundingBox bbox = computeBoundingBox(*mesh).scale(1.2); Octree2D octree(bbox, mesh); octree.generateIndex(); @@ -191,8 +190,6 @@ TEST(quest_inout_quadtree, circle_mesh) 100. * (static_cast(insideCount) / NUM_PT_TESTS), 100. * (static_cast(outsideCount) / NUM_PT_TESTS), 100. * (static_cast(uncertainCount) / NUM_PT_TESTS))); - - delete mesh; } } } diff --git a/src/axom/quest/tests/quest_intersection_shaper.cpp b/src/axom/quest/tests/quest_intersection_shaper.cpp index 95809d2d41..6d32f07c9f 100644 --- a/src/axom/quest/tests/quest_intersection_shaper.cpp +++ b/src/axom/quest/tests/quest_intersection_shaper.cpp @@ -462,7 +462,7 @@ void IntersectionWithErrorTolerances(const std::string &filebase, quest::IntersectionShaper shaper(shapeSet, &dc); shaper.setLevel(refinementLevel); shaper.setPercentError(targetPercentError); - shaper.setRefinementType(quest::Shaper::RefinementDynamic); + shaper.setRefinementType(quest::DiscreteShape::RefinementDynamic); shaper.setExecPolicy(policy); // Borrowed from shaping_driver (there should just be one shape) diff --git a/src/axom/sidre/core/Group.cpp b/src/axom/sidre/core/Group.cpp index f4314dfadc..3bb2ef67b3 100644 --- a/src/axom/sidre/core/Group.cpp +++ b/src/axom/sidre/core/Group.cpp @@ -963,7 +963,7 @@ View* Group::copyView(View* view) * ************************************************************************* */ -View* Group::deepCopyView(View* view, int allocID) +View* Group::deepCopyView(const View* view, int allocID) { allocID = getValidAllocatorID(allocID); @@ -1402,7 +1402,7 @@ Group* Group::copyGroup(Group* group) * ************************************************************************* */ -Group* Group::deepCopyGroup(Group* srcGroup, int allocID) +Group* Group::deepCopyGroup(const Group* srcGroup, int allocID) { allocID = getValidAllocatorID(allocID); diff --git a/src/axom/sidre/core/Group.hpp b/src/axom/sidre/core/Group.hpp index 690f48e7d4..7209f92928 100644 --- a/src/axom/sidre/core/Group.hpp +++ b/src/axom/sidre/core/Group.hpp @@ -862,13 +862,13 @@ class Group /*! * \brief Destroy View with given name or path owned by this Group, but leave - * its data intect. + * its data intact. */ void destroyView(const std::string& path); /*! * \brief Destroy View with given index owned by this Group, but leave - * its data intect. + * its data intact. */ void destroyView(IndexType idx); @@ -950,7 +950,7 @@ class Group * \return pointer to the new copied View object or nullptr if a View * is not copied into this Group. */ - View* deepCopyView(View* view, int allocID = INVALID_ALLOCATOR_ID); + View* deepCopyView(const View* view, int allocID = INVALID_ALLOCATOR_ID); //@} @@ -1277,7 +1277,7 @@ class Group /*! * \brief Create a (shallow) copy of Group hierarchy rooted at given - * Group and make it a child of this Group. + * Group and make the copy a child of this Group. * * Note that all Views in the Group hierarchy are copied as well. * @@ -1319,7 +1319,7 @@ class Group * \return pointer to the new copied Group object or nullptr if a Group * is not copied into this Group. */ - Group* deepCopyGroup(Group* srcGroup, int allocID = INVALID_ALLOCATOR_ID); + Group* deepCopyGroup(const Group* srcGroup, int allocID = INVALID_ALLOCATOR_ID); //@}