diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index be388661e1..1b19578f7a 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -20,6 +20,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ ## [Unreleased] - Release date yyyy-mm-dd ### Added +- Support in `quest::IntersectionShaper` for Blueprint mesh stored in a `conduit::Node` + or `sidre::Group`. - A number of new `klee::Geometry` constructors are added, for the different shapes now supported. This is a temporary change. The class will be subclassed in the future to support a diversity of geometries. - Support some analytical shapes in `IntersectionShaper`. @@ -32,6 +34,8 @@ to use Open Cascade's file I/O capabilities in support of Quest applications. - Adds a Quest example that reads in a STEP file using Open Cascade and processes its geometry ### Changed +- `quest::Shaper` and `quest::IntersectionShaper` constructors require a runtime policy. + Changing the policy after construction is no longer supported. - Importing Conduit array data into `sidre::View` now allocates destination data using the `View`'s parent's allocator ID, instead of always using host memory. This is consistent with the behavior of deep-copying data diff --git a/src/axom/core/CMakeLists.txt b/src/axom/core/CMakeLists.txt index f1dbf22b26..1a87c367dd 100644 --- a/src/axom/core/CMakeLists.txt +++ b/src/axom/core/CMakeLists.txt @@ -110,6 +110,7 @@ blt_list_append( TO core_depends ELEMENTS caliper IF CALIPER_FOUND ) blt_list_append( TO core_depends ELEMENTS camp IF CAMP_FOUND ) blt_list_append( TO core_depends ELEMENTS umpire IF UMPIRE_FOUND ) blt_list_append( TO core_depends ELEMENTS RAJA IF RAJA_FOUND ) +blt_list_append( TO core_depends ELEMENTS conduit::conduit IF CONDUIT_FOUND ) blt_list_append( TO core_depends ELEMENTS mpi IF AXOM_ENABLE_MPI ) # HACK: RAJA's dependencies are not getting added to core due to a bug in diff --git a/src/axom/core/Types.hpp b/src/axom/core/Types.hpp index 535c852570..d39789691f 100644 --- a/src/axom/core/Types.hpp +++ b/src/axom/core/Types.hpp @@ -15,6 +15,10 @@ // Axom includes #include "axom/config.hpp" +#if defined(AXOM_USE_CONDUIT) + #include "conduit/conduit_data_type.hpp" +#endif + // C/C++ includes #include @@ -62,8 +66,16 @@ using uint64 = std::uint64_t; /*!< 64-bit unsigned integer type */ #if defined(AXOM_USE_64BIT_INDEXTYPE) && !defined(AXOM_NO_INT64_T) using IndexType = std::int64_t; + #if defined(AXOM_USE_CONDUIT) +static constexpr conduit::DataType::TypeID conduitDataIdOfAxomIndexType = + conduit::DataType::INT64_ID; + #endif #else using IndexType = std::int32_t; + #if defined(AXOM_USE_CONDUIT) +static constexpr conduit::DataType::TypeID conduitDataIdOfAxomIndexType = + conduit::DataType::INT32_ID; + #endif #endif static constexpr IndexType InvalidIndex = -1; diff --git a/src/axom/core/execution/execution_space.hpp b/src/axom/core/execution/execution_space.hpp index 1c23a74d90..f619fc66df 100644 --- a/src/axom/core/execution/execution_space.hpp +++ b/src/axom/core/execution/execution_space.hpp @@ -8,6 +8,7 @@ #include "axom/config.hpp" #include "axom/core/memory_management.hpp" +#include "axom/core/execution/runtime_policy.hpp" /*! * \file execution_space.hpp @@ -88,6 +89,19 @@ struct execution_space static constexpr bool onDevice() noexcept { return false; } static constexpr char* name() noexcept { return (char*)"[UNDEFINED]"; } static int allocatorID() noexcept { return axom::INVALID_ALLOCATOR_ID; } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::seq; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; } // namespace axom @@ -109,4 +123,36 @@ struct execution_space #include "axom/core/execution/internal/hip_exec.hpp" #endif +namespace axom +{ + +/// \brief Return default allocator id for a runtime policy. +static const std::map s_policyToDefaultAllocatorID { + {axom::runtime_policy::Policy::seq, + axom::execution_space::allocatorID()} +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + , + {axom::runtime_policy::Policy::omp, + axom::execution_space::allocatorID()} +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + , + {axom::runtime_policy::Policy::cuda, + axom::execution_space>::allocatorID()} +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + , + {axom::runtime_policy::Policy::hip, + axom::execution_space>::allocatorID()} +#endif +}; + +/// \brief Return default allocator id for a runtime policy. +inline int policyToDefaultAllocatorID(axom::runtime_policy::Policy policy) +{ + return s_policyToDefaultAllocatorID.find(policy)->second; +} + +} // namespace axom + #endif // AXOM_EXECUTIONSPACE_HPP_ diff --git a/src/axom/core/execution/internal/cuda_exec.hpp b/src/axom/core/execution/internal/cuda_exec.hpp index e4ede04544..36ff1901fe 100644 --- a/src/axom/core/execution/internal/cuda_exec.hpp +++ b/src/axom/core/execution/internal/cuda_exec.hpp @@ -65,6 +65,19 @@ struct execution_space> { return axom::getUmpireResourceAllocatorID(umpire::resource::Device); } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::cuda; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space || m == MemorySpace::Unified; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; /*! @@ -95,6 +108,19 @@ struct execution_space> { return axom::getUmpireResourceAllocatorID(umpire::resource::Device); } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::cuda; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space || m == MemorySpace::Unified; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; } // namespace axom diff --git a/src/axom/core/execution/internal/hip_exec.hpp b/src/axom/core/execution/internal/hip_exec.hpp index 42b43a3eba..8d292f0fb0 100644 --- a/src/axom/core/execution/internal/hip_exec.hpp +++ b/src/axom/core/execution/internal/hip_exec.hpp @@ -63,6 +63,19 @@ struct execution_space> { return axom::getUmpireResourceAllocatorID(umpire::resource::Device); } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::hip; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space || m == MemorySpace::Unified; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; /*! @@ -90,6 +103,19 @@ struct execution_space> { return axom::getUmpireResourceAllocatorID(umpire::resource::Device); } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::hip; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space || m == MemorySpace::Unified; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; } // namespace axom diff --git a/src/axom/core/execution/internal/omp_exec.hpp b/src/axom/core/execution/internal/omp_exec.hpp index 448b0fc923..860bc7d677 100644 --- a/src/axom/core/execution/internal/omp_exec.hpp +++ b/src/axom/core/execution/internal/omp_exec.hpp @@ -60,6 +60,23 @@ struct execution_space return axom::getDefaultAllocatorID(); #endif } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::omp; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space +#ifdef AXOM_USE_UMPIRE + || m == MemorySpace::Unified +#endif + ; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; } // namespace axom diff --git a/src/axom/core/execution/internal/seq_exec.hpp b/src/axom/core/execution/internal/seq_exec.hpp index 9c50f80ee8..10833cb584 100644 --- a/src/axom/core/execution/internal/seq_exec.hpp +++ b/src/axom/core/execution/internal/seq_exec.hpp @@ -69,6 +69,23 @@ struct execution_space return axom::getDefaultAllocatorID(); #endif } + static constexpr runtime_policy::Policy runtimePolicy() noexcept + { + return runtime_policy::Policy::seq; + } + static bool usesMemorySpace(axom::MemorySpace m) noexcept + { + return m == memory_space +#ifdef AXOM_USE_UMPIRE + || m == MemorySpace::Unified +#endif + ; + } + static bool usesAllocId(int allocId) noexcept + { + return allocId == 0 || + usesMemorySpace(axom::detail::getAllocatorSpace(allocId)); + } }; } // namespace axom diff --git a/src/axom/core/execution/runtime_policy.hpp b/src/axom/core/execution/runtime_policy.hpp index 445bb538ff..246488ea67 100644 --- a/src/axom/core/execution/runtime_policy.hpp +++ b/src/axom/core/execution/runtime_policy.hpp @@ -6,7 +6,8 @@ #ifndef AXOM_CORE_EXECUTION_RUNTIME_POLICY_HPP_ #define AXOM_CORE_EXECUTION_RUNTIME_POLICY_HPP_ -#include "axom/config.hpp" /* for compile time defs. */ +#include "axom/config.hpp" /* for compile time defs. */ +#include #include "axom/fmt/format.h" /* for axom::fmt */ #include diff --git a/src/axom/core/memory_management.hpp b/src/axom/core/memory_management.hpp index 5a8cebb424..fe85f95590 100644 --- a/src/axom/core/memory_management.hpp +++ b/src/axom/core/memory_management.hpp @@ -110,6 +110,27 @@ inline int getDefaultAllocatorID() #endif } +/*! + * \brief Get the allocator id from which data has been allocated. + * \return Allocator id. If Umpire doesn't have an allocator for + * the pointer, or Axom wasn't configured with Umpire, return 0. + * + * \pre ptr has a valid pointer value. + */ +inline int getAllocatorIDFromPointer(const void* ptr) +{ +#ifdef AXOM_USE_UMPIRE + umpire::ResourceManager& rm = umpire::ResourceManager::getInstance(); + if(rm.hasAllocator(const_cast(ptr))) + { + umpire::Allocator allocator = rm.getAllocator(const_cast(ptr)); + return allocator.getId(); + } +#endif + AXOM_UNUSED_VAR(ptr); + return 0; +} + /*! * \brief Allocates a chunk of memory of type T. * diff --git a/src/axom/core/utilities/FileUtilities.cpp b/src/axom/core/utilities/FileUtilities.cpp index c359ed9dd6..81c1811468 100644 --- a/src/axom/core/utilities/FileUtilities.cpp +++ b/src/axom/core/utilities/FileUtilities.cpp @@ -106,6 +106,32 @@ int makeDirsForPath(const std::string& path) return err; } +//----------------------------------------------------------------------------- +std::string prefixRelativePath(const std::string& path, const std::string& prefix) +{ + if(path[0] == '/' || prefix.empty()) + { + return path; + } + return utilities::filesystem::joinPath(prefix, path); +} + +//----------------------------------------------------------------------------- +std::string getParentPath(const std::string& path) +{ + char separator = '/'; + std::size_t found = path.rfind(separator); + + std::string parent; + + if(found != std::string::npos) + { + parent = path.substr(0, found); + } + + return parent; +} + //----------------------------------------------------------------------------- void getDirName(std::string& dir, const std::string& path) { diff --git a/src/axom/core/utilities/FileUtilities.hpp b/src/axom/core/utilities/FileUtilities.hpp index 9f70620ab5..bfe1406f96 100644 --- a/src/axom/core/utilities/FileUtilities.hpp +++ b/src/axom/core/utilities/FileUtilities.hpp @@ -66,6 +66,28 @@ std::string joinPath(const std::string& fileDir, */ int makeDirsForPath(const std::string& path); +/*! + * \brief Add a prefix to a path if it is relative. + * + * \param [in] path string representing an absolute or relative path + * \param [in] prefix string representing a directory path + * + * \return \c prefix + \c path (with delimiter) if the path is + * relative or \c prefix is empty, or just \c path otherwise. + */ +std::string prefixRelativePath(const std::string& path, + const std::string& prefix); + +/*! + * \brief Get parent path name from a filesystem path. + * + * \param [in] path an absolute or relative filesystem path + * + * \return a directory path formed by removing the last part of the + * input path + */ +std::string getParentPath(const std::string& path); + /*! * \brief Get directory name from a path that contains a file name * diff --git a/src/axom/klee/Geometry.cpp b/src/axom/klee/Geometry.cpp index 939e542926..8ddcf51aca 100644 --- a/src/axom/klee/Geometry.cpp +++ b/src/axom/klee/Geometry.cpp @@ -86,19 +86,19 @@ Geometry::Geometry(const TransformableGeometryProperties& startProperties, Geometry::Geometry(const TransformableGeometryProperties& startProperties, const axom::Array& discreteFunction, - const Point3D& vorBase, - const Vector3D& vorDirection, + const Point3D& sorBase, + const Vector3D& sorDirection, axom::IndexType levelOfRefinement, std::shared_ptr operator_) : m_startProperties(startProperties) - , m_format("vor3D") + , m_format("sor3D") , m_path() , m_meshGroup(nullptr) , m_topology() , m_sphere() , m_discreteFunction(discreteFunction) - , m_vorBase(vorBase) - , m_vorDirection(vorDirection) + , m_sorBase(sorBase) + , m_sorDirection(sorDirection) , m_levelOfRefinement(levelOfRefinement) , m_operator(std::move(operator_)) { } diff --git a/src/axom/klee/Geometry.hpp b/src/axom/klee/Geometry.hpp index e23d095a99..c08812200b 100644 --- a/src/axom/klee/Geometry.hpp +++ b/src/axom/klee/Geometry.hpp @@ -132,27 +132,27 @@ class Geometry std::shared_ptr operator_); /** - * Create a volume-of-revolution (VOR) Geometry object. + * Create a volume-of-revolution (SOR) 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 sorBase Coordinates of the base of the SOR. + * \param sorDirection SOR axis, in the direction of increasing z. * \param levelOfRefinement Number of refinement levels to use for - * discretizing the VOR. + * discretizing the SOR. * \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. + * The \c sorBase coordinates corresponds to z=0. + * \c sorAxis should point in the direction of increasing z. */ Geometry(const TransformableGeometryProperties &startProperties, const axom::Array &discreteFunction, - const Point3D &vorBase, - const Vector3D &vorDirection, + const Point3D &sorBase, + const Vector3D &sorDirection, axom::IndexType levelOfRefinement, std::shared_ptr operator_); @@ -181,7 +181,7 @@ class Geometry * - "blueprint-tets" = Blueprint tetrahedral mesh in memory * - "tet3D" = 3D tetrahedron (4 points) * - "sphere3D" = 3D sphere, as \c primal::Sphere - * - "vor3D" = 3D volume of revolution. + * - "sor3D" = 3D volume of revolution. * - "cone3D" = 3D cone, as \c primal::Cone * - "cylinder3D" = 3D cylinder, as \c primal::Cylinder * - "hex3D" = 3D hexahedron (8 points) @@ -216,14 +216,14 @@ class Geometry const std::string &getBlueprintTopology() const; /** - @brief Return the VOR axis direction. + @brief Return the SOR axis direction. */ - const Vector3D getVorDirection() const { return m_vorDirection; } + const Vector3D getSorDirection() const { return m_sorDirection; } /** - @brief Return the 3D coordinates of the VOR base. + @brief Return the 3D coordinates of the SOR base. */ - const Point3D getVorBaseCoords() const { return m_vorBase; } + const Point3D getSorBaseCoords() const { return m_sorBase; } /*! @brief Predicate that returns true when the shape has an associated geometry @@ -336,11 +336,11 @@ class Geometry //! @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 The point corresponding to z=0 on the SOR axis. + Point3D m_sorBase; - //!@brief VOR axis in the direction of increasing z. - Vector3D m_vorDirection; + //!@brief SOR axis in the direction of increasing z. + Vector3D m_sorDirection; //!@brief Level of refinement for discretizing curved // analytical shapes and surfaces of revolutions. diff --git a/src/axom/klee/ShapeSet.cpp b/src/axom/klee/ShapeSet.cpp index 326000bce9..6a5d384f1c 100644 --- a/src/axom/klee/ShapeSet.cpp +++ b/src/axom/klee/ShapeSet.cpp @@ -21,21 +21,6 @@ void ShapeSet::setShapes(std::vector shapes) void ShapeSet::setPath(const std::string &path) { m_path = path; } -std::string ShapeSet::resolvePath(const std::string &filePath) const -{ - if(m_path.empty()) - { - throw std::logic_error("The ShapeSet's path has not been set"); - } - if(filePath[0] == '/') - { - return filePath; - } - std::string dir; - utilities::filesystem::getDirName(dir, m_path); - return utilities::filesystem::joinPath(dir, filePath); -} - void ShapeSet::setDimensions(Dimensions dimensions) { m_dimensions = dimensions; diff --git a/src/axom/klee/ShapeSet.hpp b/src/axom/klee/ShapeSet.hpp index f8c197b31a..602ada06b2 100644 --- a/src/axom/klee/ShapeSet.hpp +++ b/src/axom/klee/ShapeSet.hpp @@ -51,16 +51,6 @@ class ShapeSet */ const std::string &getPath() const { return m_path; } - /** - * Resolves a path relative to the path of this ShapeSet. - * - * \param filePath the path to resolve - * \return if the given path is absolute, then the given path. Otherwise, - * the path is interpreted as being relative to the directory containing - * this ShapeSet, and that is is returned. - */ - std::string resolvePath(const std::string &filePath) const; - /** * Sets the dimensions for all shapes in the ShapeSet. * diff --git a/src/axom/klee/tests/klee_shape_set.cpp b/src/axom/klee/tests/klee_shape_set.cpp index 40be40a136..9c7f4f3888 100644 --- a/src/axom/klee/tests/klee_shape_set.cpp +++ b/src/axom/klee/tests/klee_shape_set.cpp @@ -15,43 +15,6 @@ namespace klee { namespace { -TEST(ShapeSetTest, resolvePath_noPathSet) -{ - ShapeSet shapeSet; - EXPECT_THROW(shapeSet.resolvePath("anyPath"), std::logic_error); -} - -TEST(ShapeSetTest, resolvePath_startWithSimpleFileName) -{ - ShapeSet shapeSet; - shapeSet.setPath("file.yaml"); - EXPECT_EQ("newPath.txt", shapeSet.resolvePath("newPath.txt")); - EXPECT_EQ("d1/d2/newPath.txt", shapeSet.resolvePath("d1/d2/newPath.txt")); - EXPECT_EQ("/abs/path/newPath.txt", - shapeSet.resolvePath("/abs/path/newPath.txt")); -} - -TEST(ShapeSetTest, resolvePath_startWithRelativeFileName) -{ - ShapeSet shapeSet; - shapeSet.setPath("path/to/file.yaml"); - EXPECT_EQ("path/to/newPath.txt", shapeSet.resolvePath("newPath.txt")); - EXPECT_EQ("path/to/d1/d2/newPath.txt", - shapeSet.resolvePath("d1/d2/newPath.txt")); - EXPECT_EQ("/abs/path/newPath.txt", - shapeSet.resolvePath("/abs/path/newPath.txt")); -} - -TEST(ShapeSetTest, resolvePath_startWithAbsoluteFileName) -{ - ShapeSet shapeSet; - shapeSet.setPath("/path/to/file.yaml"); - EXPECT_EQ("/path/to/newPath.txt", shapeSet.resolvePath("newPath.txt")); - EXPECT_EQ("/path/to/d1/d2/newPath.txt", - shapeSet.resolvePath("d1/d2/newPath.txt")); - EXPECT_EQ("/other/abs/path/newPath.txt", - shapeSet.resolvePath("/other/abs/path/newPath.txt")); -} TEST(ShapeSetTest, dimensions_getAndSet) { diff --git a/src/axom/mint/mesh/Mesh.cpp b/src/axom/mint/mesh/Mesh.cpp index fc54c38395..e1df1fd075 100644 --- a/src/axom/mint/mesh/Mesh.cpp +++ b/src/axom/mint/mesh/Mesh.cpp @@ -123,22 +123,19 @@ Mesh::Mesh(sidre::Group* group, const std::string& topo) m_coordset = blueprint::getCoordsetGroup(m_group, getTopologyGroup())->getName(); - SLIC_ERROR_IF(!m_group->hasChildGroup("state"), - "root group does not have a state group."); - sidre::Group* state_group = m_group->getGroup("state"); - SLIC_ERROR_IF(!state_group->hasChildGroup(m_topology), - "state group has no " << m_topology << " child group."); - - state_group = state_group->getGroup(m_topology); - if(state_group->hasChildView("block_id")) + if(state_group != nullptr && state_group->hasChildGroup(m_topology)) { - m_block_idx = state_group->getView("block_id")->getScalar(); - } + state_group = state_group->getGroup(m_topology); + if(state_group->hasChildView("block_id")) + { + m_block_idx = state_group->getView("block_id")->getScalar(); + } - if(state_group->hasChildView("partition_id")) - { - m_part_idx = state_group->getView("partition_id")->getScalar(); + if(state_group->hasChildView("partition_id")) + { + m_part_idx = state_group->getView("partition_id")->getScalar(); + } } SLIC_ERROR_IF(!validMeshType(), "invalid mesh type=" << m_type); diff --git a/src/axom/primal/geometry/KnotVector.hpp b/src/axom/primal/geometry/KnotVector.hpp index bfc3ca2f3a..a82d230021 100644 --- a/src/axom/primal/geometry/KnotVector.hpp +++ b/src/axom/primal/geometry/KnotVector.hpp @@ -157,8 +157,8 @@ class KnotVector */ T& operator[](axom::IndexType i) { - return m_knots[i]; SLIC_ASSERT(isValid()); + return m_knots[i]; } /// \brief Return the degree of the knot vector diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index 3aabb7d98c..8e7308da15 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -1369,4 +1369,4 @@ template struct axom::fmt::formatter> : ostream_formatter { }; -#endif // AXOM_PRIMAL_NURBSCURVE_HPP_ \ No newline at end of file +#endif // AXOM_PRIMAL_NURBSCURVE_HPP_ diff --git a/src/axom/quest/CMakeLists.txt b/src/axom/quest/CMakeLists.txt index c9dd6f67da..bac906450f 100644 --- a/src/axom/quest/CMakeLists.txt +++ b/src/axom/quest/CMakeLists.txt @@ -130,16 +130,20 @@ blt_list_append( blt_list_append( TO quest_depends_on ELEMENTS conduit::conduit IF CONDUIT_FOUND ) -if(MFEM_FOUND AND AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE AND AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION) +if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE AND ((MFEM_FOUND AND AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION) OR CONDUIT_FOUND)) list(APPEND quest_headers Shaper.hpp - SamplingShaper.hpp - IntersectionShaper.hpp - DiscreteShape.hpp - detail/shaping/shaping_helpers.hpp) + DiscreteShape.hpp) list(APPEND quest_sources Shaper.cpp - DiscreteShape.cpp - detail/shaping/shaping_helpers.cpp) + DiscreteShape.cpp) list(APPEND quest_depends_on klee) + if(RAJA_FOUND) + list(APPEND quest_headers IntersectionShaper.hpp) + endif() +endif() + +if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE AND (MFEM_FOUND AND AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION)) + list(APPEND quest_headers SamplingShaper.hpp detail/shaping/shaping_helpers.hpp) + list(APPEND quest_sources detail/shaping/shaping_helpers.cpp) endif() if(C2C_FOUND) diff --git a/src/axom/quest/DiscreteShape.cpp b/src/axom/quest/DiscreteShape.cpp index 4fe0431d53..805d60ddb4 100644 --- a/src/axom/quest/DiscreteShape.cpp +++ b/src/axom/quest/DiscreteShape.cpp @@ -134,9 +134,9 @@ std::shared_ptr DiscreteShape::createMeshRepresentation() { createRepresentationOfSphere(); } - if(geometryFormat == "vor3D") + if(geometryFormat == "sor3D") { - createRepresentationOfVOR(); + createRepresentationOfSOR(); } if(m_meshRep) @@ -157,7 +157,9 @@ std::shared_ptr DiscreteShape::createMeshRepresentation() // We handled all the non-file formats. The rest are file formats. const std::string& file_format = geometryFormat; - std::string shapePath = resolvePath(); + std::string shapePath = axom::utilities::filesystem::prefixRelativePath( + m_shape.getGeometry().getPath(), + m_prefixPath); SLIC_INFO("Reading file: " << shapePath << "..."); // Initialize revolved volume. @@ -406,7 +408,7 @@ void DiscreteShape::createRepresentationOfPlane() // clang-format on // Rotate and translate boundingHex to align with the plane. - numerics::Matrix rotate = vorAxisRotMatrix(plane.getNormal()); + numerics::Matrix rotate = sorAxisRotMatrix(plane.getNormal()); const auto translate = plane.getNormal() * plane.getOffset(); for(int i = 0; i < 8; ++i) { @@ -529,11 +531,11 @@ void DiscreteShape::createRepresentationOfSphere() applyTransforms(); } -void DiscreteShape::createRepresentationOfVOR() +void DiscreteShape::createRepresentationOfSOR() { // Construct the tet m_meshRep from the volume-of-revolution. - auto& vorGeom = m_shape.getGeometry(); - const auto& discreteFcn = vorGeom.getDiscreteFunction(); + auto& sorGeom = m_shape.getGeometry(); + const auto& discreteFcn = sorGeom.getDiscreteFunction(); // Generate the Octahedra axom::Array octs; @@ -546,11 +548,12 @@ void DiscreteShape::createRepresentationOfVOR() m_shape.getGeometry().getLevelOfRefinement(), octs, octCount); + AXOM_UNUSED_VAR(good); 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(); + // Rotate to the SOR axis direction and translate to the base location. + numerics::Matrix rotate = sorAxisRotMatrix(sorGeom.getSorDirection()); + const auto& translate = sorGeom.getSorBaseCoords(); auto octsView = octs.view(); axom::for_all( octCount, @@ -568,6 +571,13 @@ void DiscreteShape::createRepresentationOfVOR() }); // Dump discretized octs as a tet mesh + // + // NOTE: mesh_from_discretized_polyline returns a tet mesh. + // IntersectionShaper can handle octahedra, which may be faster + // because it has 8x fewer elements. This method returns the + // discretization as a mesh. But mint doesn't support octs and + // Blueprint support octs only as polyhedra. We can always return + // an array of primal::Octahedron instead of a mesh. axom::mint::Mesh* mesh; axom::quest::mesh_from_discretized_polyline(octs.view(), octCount, @@ -680,7 +690,7 @@ numerics::Matrix DiscreteShape::getTransforms() const } // Return a 3x3 matrix that rotates coordinates from the x-axis to the given direction. -numerics::Matrix DiscreteShape::vorAxisRotMatrix(const Vector3D& dir) +numerics::Matrix DiscreteShape::sorAxisRotMatrix(const Vector3D& dir) { // Note that the rotation matrix is not unique. static const Vector3D x {1.0, 0.0, 0.0}; @@ -766,20 +776,12 @@ void DiscreteShape::setPercentError(double percent) clampVal(percent, MINIMUM_PERCENT_ERROR, MAXIMUM_PERCENT_ERROR); } -std::string DiscreteShape::resolvePath() const +void DiscreteShape::setPrefixPath(const std::string& prefixPath) { - 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); + SLIC_ERROR_IF( + !prefixPath.empty() && !axom::utilities::filesystem::pathExists(prefixPath), + "Path '" + prefixPath + "' does not exist."); + m_prefixPath = prefixPath; } void DiscreteShape::setParentGroup(axom::sidre::Group* parentGroup) diff --git a/src/axom/quest/DiscreteShape.hpp b/src/axom/quest/DiscreteShape.hpp index 0d43679726..5a3d23768d 100644 --- a/src/axom/quest/DiscreteShape.hpp +++ b/src/axom/quest/DiscreteShape.hpp @@ -54,11 +54,11 @@ class DiscreteShape @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 + @param prefixPath Path prefix for shape from a file specified with a relative path. Refinement type is set to DiscreteShape::RefinementUniformSegments - and percent erro is set to 0. See setPercentError() and + and percent error is set to 0. See setPercentError() and setRefinementType(). */ DiscreteShape(const axom::klee::Shape& shape, @@ -70,10 +70,8 @@ class DiscreteShape //@{ //! @name Functions to get and set shaping parameters - void setPrefixPath(const std::string& prefixPath) - { - m_prefixPath = prefixPath; - } + //! @brief Set prefix for shape files specified as relative path. + void setPrefixPath(const std::string& prefixPath); /*! @brief Set the refinement type. @@ -170,16 +168,13 @@ class DiscreteShape 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); + numerics::Matrix sorAxisRotMatrix(const Vector3D& dir); void clearInternalData(); @@ -195,8 +190,8 @@ class DiscreteShape void createRepresentationOfPlane(); //!@brief Create the internal mesh representation of the analytical sphere. void createRepresentationOfSphere(); - //!@brief Create the internal mesh representation of the analytical VOR. - void createRepresentationOfVOR(); + //!@brief Create the internal mesh representation of the analytical SOR. + void createRepresentationOfSOR(); }; } // namespace quest diff --git a/src/axom/quest/IntersectionShaper.hpp b/src/axom/quest/IntersectionShaper.hpp index ed2cb6ef49..3a942d5e4b 100644 --- a/src/axom/quest/IntersectionShaper.hpp +++ b/src/axom/quest/IntersectionShaper.hpp @@ -13,35 +13,37 @@ #define AXOM_QUEST_INTERSECTION_SHAPER__HPP_ #include "axom/config.hpp" -#include "axom/core.hpp" -#include "axom/slic.hpp" -#include "axom/slam.hpp" -#include "axom/primal.hpp" -#include "axom/mint.hpp" -#include "axom/spin.hpp" -#include "axom/klee.hpp" - -#ifndef AXOM_USE_MFEM - #error Shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option -#endif - -#include "axom/quest/Shaper.hpp" -#include "axom/spin/BVH.hpp" -#include "axom/quest/interface/internal/mpicomm_wrapper.hpp" -#include "axom/quest/interface/internal/QuestHelpers.hpp" -#include "axom/quest/detail/shaping/shaping_helpers.hpp" - -#include "mfem.hpp" +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) + + #include "axom/core.hpp" + #include "axom/slic.hpp" + #include "axom/slam.hpp" + #include "axom/primal.hpp" + #include "axom/sidre/core/Group.hpp" + #include "axom/sidre/core/View.hpp" + #include "axom/mint/mesh/UnstructuredMesh.hpp" + #include "axom/mint/utils/vtk_utils.hpp" + #include "axom/klee.hpp" + #include "axom/quest/Shaper.hpp" + #include "axom/quest/Discretize.hpp" + #include "axom/spin/BVH.hpp" + #include "axom/quest/interface/internal/mpicomm_wrapper.hpp" + #include "axom/quest/interface/internal/QuestHelpers.hpp" + #include "axom/fmt.hpp" + + #ifdef AXOM_USE_MFEM + #include "mfem.hpp" + #endif -#include "axom/fmt.hpp" + #include "axom/fmt.hpp" -// RAJA -#if defined(AXOM_USE_RAJA) - #include "RAJA/RAJA.hpp" -#endif + #if defined(AXOM_USE_CONDUIT) + #include "conduit_node.hpp" + #include "conduit_blueprint_mesh.hpp" + #include "conduit_blueprint_mcarray.hpp" + #endif // clang-format off -#if defined (AXOM_USE_RAJA) using seq_exec = axom::SEQ_EXEC; #if defined(AXOM_USE_OPENMP) @@ -63,7 +65,6 @@ #else using hip_exec = seq_exec; #endif -#endif // clang-format on namespace axom @@ -71,7 +72,7 @@ namespace axom namespace quest { /*! - * \class GridFunctionView + * \class TempArrayView * * \brief Provides a view over an MFEM grid function. MFEM grid functions are * assumed to live in host memory. This class performs data movement @@ -82,19 +83,30 @@ namespace quest * be accessed. */ template -class GridFunctionView +class TempArrayView { public: /*! - * \brief Host constructor that accepts the grid function. + * \brief Host constructor that accepts data in an Array. + * + * \param gf The data that will be accessed/modified by this view. + * \param _needResult Whether the data needs to be brought back to the host + * from the device. + */ + AXOM_HOST TempArrayView(axom::Array& gf, bool _needResult = true) + { + initialize(gf.data(), gf.size(), _needResult); + } + /*! + * \brief Host constructor that accepts data in an ArrayView. * - * \param gf The grid function that will be accessed/modified by the view. + * \param gf The data that will be accessed/modified by this view. * \param _needResult Whether the data needs to be brought back to the host * from the device. */ - AXOM_HOST GridFunctionView(mfem::GridFunction* gf, bool _needResult = true) + AXOM_HOST TempArrayView(axom::ArrayView& gf, bool _needResult = true) { - initialize(gf->GetData(), gf->Size(), _needResult); + initialize(gf.data(), gf.size(), _needResult); } /*! @@ -103,7 +115,7 @@ class GridFunctionView * happened in the host constructor. This version sets hostData to * nullptr so we know not to clean up in the destructor. */ - AXOM_HOST_DEVICE GridFunctionView(const GridFunctionView& obj) + AXOM_HOST_DEVICE TempArrayView(const TempArrayView& obj) : m_hostData(nullptr) , m_deviceData(obj.m_deviceData) , m_numElements(obj.m_numElements) @@ -114,7 +126,7 @@ class GridFunctionView * \brief Destructor. On the host, this method may move data from the device and deallocate device storage. */ - AXOM_HOST_DEVICE ~GridFunctionView() { finalize(); } + AXOM_HOST_DEVICE ~TempArrayView() { finalize(); } /*! * \brief Indexing operator for accessing the data. @@ -148,7 +160,7 @@ class GridFunctionView */ AXOM_HOST_DEVICE void finalize() { m_deviceData = nullptr; } -#if defined(AXOM_USE_CUDA) || defined(AXOM_USE_HIP) + #if defined(AXOM_USE_CUDA) || defined(AXOM_USE_HIP) /*! * \brief Initializes members using data from the grid function. This method * is called on the host and it copies data to the device. @@ -175,7 +187,7 @@ class GridFunctionView */ AXOM_HOST_DEVICE void finalizeDevice() { - #ifndef AXOM_DEVICE_CODE + #ifndef AXOM_DEVICE_CODE // Only the host will do this work. if(m_hostData != nullptr) { @@ -187,9 +199,9 @@ class GridFunctionView axom::deallocate(m_deviceData); m_deviceData = nullptr; } - #endif + #endif } -#endif + #endif private: double* m_hostData {nullptr}; @@ -198,7 +210,7 @@ class GridFunctionView bool m_needResult {false}; }; -#if defined(AXOM_USE_CUDA) + #if defined(AXOM_USE_CUDA) /*! * \brief CUDA specialization that calls initializeDevice to copy data * from the host to the device. @@ -208,9 +220,9 @@ class GridFunctionView * \param _needResult Whether any data are copied from device. */ template <> -AXOM_HOST inline void GridFunctionView::initialize(double* hostPtr, - int nElem, - bool _needResult) +AXOM_HOST inline void TempArrayView::initialize(double* hostPtr, + int nElem, + bool _needResult) { initializeDevice(hostPtr, nElem, _needResult); } @@ -220,12 +232,12 @@ AXOM_HOST inline void GridFunctionView::initialize(double* hostPtr, * and deallocate any associated device data. */ template <> -AXOM_HOST_DEVICE inline void GridFunctionView::finalize() +AXOM_HOST_DEVICE inline void TempArrayView::finalize() { finalizeDevice(); } -#endif -#if defined(AXOM_USE_HIP) + #endif + #if defined(AXOM_USE_HIP) /*! * \brief HIP specialization that calls initializeDevice to copy data * from the host to the device. @@ -235,9 +247,9 @@ AXOM_HOST_DEVICE inline void GridFunctionView::finalize() * \param _needResult Whether any data are copied from device. */ template <> -AXOM_HOST inline void GridFunctionView::initialize(double* hostPtr, - int nElem, - bool _needResult) +AXOM_HOST inline void TempArrayView::initialize(double* hostPtr, + int nElem, + bool _needResult) { initializeDevice(hostPtr, nElem, _needResult); } @@ -247,17 +259,17 @@ AXOM_HOST inline void GridFunctionView::initialize(double* hostPtr, * and deallocate any associated device data. */ template <> -AXOM_HOST_DEVICE inline void GridFunctionView::finalize() +AXOM_HOST_DEVICE inline void TempArrayView::finalize() { finalizeDevice(); } -#endif + #endif //--------------------------------------------------------------------------- /** * \class - * \brief Intersects a solid, created by revolving a c2c contour or by - * loading a Pro/E mesh, with an input mesh to produce volume fractions. + * \brief Intersects a solid with an input mesh to produce volume fractions. + * The solid may be any supported by the \c klee::Geometry class. * * The IntersectionShaper generates material volume fractions: * @@ -269,16 +281,30 @@ AXOM_HOST_DEVICE inline void GridFunctionView::finalize() * intersected with the mesh. The octahedra are intersected with the input * mesh to produce volume fractions. * - * - For Pro/E, an input mesh of 3D tetrahedra are loaded in. + * - For tetrahedral mesh (including Pro/E), an input mesh of 3D tetrahedra is loaded in. * Each tetrahedron has its own respective volume. The tetrahedra are * intersected with the input mesh to produce volume fractions. * - * Volume fractions are represented as a GridFunction with a special prefix, + * - For analytical geometries, the shapes are discretized into a tetrahedral + * mesh first. Sphere and surfaces-of-revolution discretization uses + * the refinement level specified in the \c Geometry. + * + * The input mesh can be an MFEM mesh stored as a \c + * sidre::MFEMSidreDataCollection or be a Blueprint mesh stored as a + * \c conduit::Node or a \c sidre::Group. + * + * IntersectionShaper requires Axom configured with RAJA and Umpire. + * + * Support for replacement rules exists for MFEM input meshes. + * Replacement rules for Blueprint meshes is not yet supported. + * The following comments apply to replacement rules. + * + * Volume fractions are represented in the input mesh as a GridFunction with a special prefix, * currently "vol_frac_", followed by a material name. Volume fractions * can be present in the input data collection prior to shaping and the * IntersectionShaper will augment them when changes are needed such as when - * a material overwrites them. If a new material is not yet represented by - * a grid function, one will be added. + * a material overwrites them. If a new material is not yet represented in + * the mesh, one will be added. * * In addition to user-specified materials, the IntersectionShaper creates * a "free" material that is used to account for volume fractions that are @@ -307,20 +333,58 @@ class IntersectionShaper : public Shaper static constexpr double DEFAULT_REVOLVED_VOLUME {0.}; public: - IntersectionShaper(const klee::ShapeSet& shapeSet, + #if defined(AXOM_USE_MFEM) + /*! + \brief Construct Shaper to operate on an MFEM mesh. + */ + IntersectionShaper(RuntimePolicy runtimePolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, sidre::MFEMSidreDataCollection* dc) - : Shaper(shapeSet, dc) + : Shaper(runtimePolicy, allocatorId, shapeSet, dc) { m_free_mat_name = "free"; } + #endif + + #if defined(AXOM_USE_CONDUIT) + /*! + \brief Construct Shaper to operate on a blueprint-formatted mesh + stored in a sidre Group. + \param [in] runtimePolicy A value from RuntimePolicy. + The simplest policy is RuntimePolicy::seq, which specifies + running sequentially on the CPU. + \param [in] allocatorID Data allocator ID. Choose something compatible + with \c runtimePolicy. See \c execution_space. + */ + IntersectionShaper(RuntimePolicy runtimePolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + sidre::Group* bpGrp, + const std::string& topo = "") + : Shaper(runtimePolicy, allocatorId, shapeSet, bpGrp, topo) + , m_free_mat_name("free") + { } + + /*! + \brief Construct Shaper to operate on a blueprint-formatted mesh + stored in a Conduit Node. + */ + IntersectionShaper(RuntimePolicy runtimePolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + conduit::Node* bpNode, + const std::string& topo = "") + : Shaper(runtimePolicy, allocatorId, shapeSet, bpNode, topo) + , m_free_mat_name("free") + { } + #endif //@{ //! @name Functions to get and set shaping parameters related to intersection; supplements parameters in base class void setLevel(int level) { m_level = level; } - void setExecPolicy(RuntimePolicy policy) { m_execPolicy = policy; } - /*! * \brief Set the name of the material used to account for free volume fractions. * \param name The new name of the material. This name cannot contain @@ -333,7 +397,7 @@ class IntersectionShaper : public Shaper { SLIC_ERROR("The free material name cannot contain underscores."); } - if(m_num_elements > 0) + if(m_cellCount > 0) { SLIC_ERROR( "The free material name cannot be set once shaping has occurred."); @@ -372,19 +436,17 @@ class IntersectionShaper : public Shaper } } - // The following private methods are made public due to CUDA compilers - // requirements for methods that call device functions. -#if defined(__CUDACC__) + // The following private methods are made public due to CUDA compilers + // requirements for methods that call device functions. + #if defined(__CUDACC__) public: -#else + #else private: -#endif + #endif //@{ //! @name Private functions related to the stages for a given shape -#if defined(AXOM_USE_RAJA) - // Prepares the tet mesh cells for the spatial index template void prepareTetCells() @@ -624,21 +686,31 @@ class IntersectionShaper : public Shaper axom::fmt::format("The shape format {} is unsupported", shapeFormat)); } } -#endif -#if defined(AXOM_USE_RAJA) + /*! + \tparam ShapeType either TetrahedeonType or OctahedronType. + depending on whether shape is tet or c2c. + \param shape the input shape to be superimposed on the mesh. + \param shapes either the array m_tets (if surface mesh is tet) + or m_octs (if surface mesh is c2c). + \param shape_count the count for the shapes array, maintained + separately from the array. + */ template void runShapeQueryImpl(const klee::Shape& shape, axom::Array& shapes, int shape_count) { + AXOM_ANNOTATE_SCOPE("IntersectionShaper::runShapeQueryImpl"); + // No need for shape, because it has been converted into m_tets or m_octs, + // which is what the parameter shapes is. + AXOM_UNUSED_VAR(shape); + const int host_allocator = axom::execution_space::allocatorID(); const int device_allocator = axom::execution_space::allocatorID(); - constexpr int NUM_VERTS_PER_HEX = 8; - constexpr int NUM_COMPS_PER_VERT = 3; constexpr int NUM_TETS_PER_HEX = 24; constexpr double ZERO_THRESHOLD = 1.e-10; @@ -646,7 +718,8 @@ class IntersectionShaper : public Shaper " Inserting shapes' bounding boxes into BVH ")); // Generate the BVH tree over the shapes - // Access-aligned bounding boxes + // Axis-aligned bounding boxes + // Does m_aabbs need to be a member? It's only used here. BTNG. m_aabbs = axom::Array(shape_count, shape_count, device_allocator); @@ -663,164 +736,55 @@ class IntersectionShaper : public Shaper }); // Insert shapes' Bounding Boxes into BVH. - //bvh.setAllocatorID(poolID); spin::BVH<3, ExecSpace, double> bvh; bvh.initialize(aabbs_device_view, shape_count); SLIC_INFO(axom::fmt::format("{:-^80}", " Querying the BVH tree ")); - mfem::Mesh* mesh = getDC()->GetMesh(); - - // Intersection algorithm only works on linear elements - SLIC_ASSERT(mesh != nullptr); - int const NE = mesh->GetNE(); - m_num_elements = NE; - - if(this->isVerbose()) + if(m_hexes.empty()) { - SLIC_INFO(axom::fmt::format( - "{:-^80}", - axom::fmt::format( - " Initializing {} hexahedral elements from given mesh ", - m_num_elements))); + // m_hexes depend only on mesh so should not change once computed. + populateHexesFromMesh(); } + axom::ArrayView hexes_device_view = m_hexes.view(); - if(NE > 0) + if(m_hex_bbs.empty()) { - SLIC_ASSERT(mesh->GetNodes() == nullptr || - mesh->GetNodes()->FESpace()->GetOrder(0)); - } - - // Create and register a scalar field for this shape's volume fractions - // The Degrees of Freedom will be in correspondence with the elements - auto* volFrac = this->newVolFracGridFunction(); - auto volFracName = axom::fmt::format("shape_vol_frac_{}", shape.getName()); - this->getDC()->RegisterField(volFracName, volFrac); - - // Initialize hexahedral elements - m_hexes = axom::Array(NE, NE, device_allocator); - axom::ArrayView hexes_device_view = m_hexes.view(); - - m_hex_bbs = axom::Array(NE, NE, device_allocator); - axom::ArrayView hex_bbs_device_view = m_hex_bbs.view(); - - // Initialize vertices from mfem mesh and - // set each shape volume fraction to 1 - // Allocation size is: - // # of elements * # of vertices per hex * # of components per vertex - axom::Array vertCoords_host( - NE * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, - NE * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, - host_allocator); - - for(int i = 0; i < NE; i++) - { - // Get the indices of this element's vertices - mfem::Array verts; - mesh->GetElementVertices(i, verts); - SLIC_ASSERT(verts.Size() == NUM_VERTS_PER_HEX); - - // Set each shape volume fraction to 1 - double vf = 1.0; - (*volFrac)(i) = vf; + // Compute m_hex_bbs only once. We assume that the mesh doesn't change + // once set. + m_hex_bbs = + axom::Array(m_cellCount, m_cellCount, device_allocator); - // Get the coordinates for the vertices - for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) - { - for(int k = 0; k < NUM_COMPS_PER_VERT; k++) - { - vertCoords_host[(i * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT) + - (j * NUM_COMPS_PER_VERT) + k] = - (mesh->GetVertex(verts[j]))[k]; - } - } + // Get bounding boxes for hexahedral elements + axom::ArrayView hex_bbs_device_view = m_hex_bbs.view(); + axom::for_all( + m_cellCount, + AXOM_LAMBDA(axom::IndexType i) { + hex_bbs_device_view[i] = + primal::compute_bounding_box(hexes_device_view[i]); + }); // end of loop to initialize hexahedral elements and bounding boxes } - axom::Array vertCoords_device = - axom::Array(vertCoords_host, device_allocator); - auto vertCoords_device_view = vertCoords_device.view(); - - // Initialize each hexahedral element and its bounding box - axom::for_all( - NE, - AXOM_LAMBDA(axom::IndexType i) { - // Set each hexahedral element vertices - hexes_device_view[i] = HexahedronType(); - for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) - { - int vertIndex = (i * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT) + - j * NUM_COMPS_PER_VERT; - hexes_device_view[i][j] = - Point3D({vertCoords_device_view[vertIndex], - vertCoords_device_view[vertIndex + 1], - vertCoords_device_view[vertIndex + 2]}); - - // Set hexahedra components to zero if within threshold - if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][0], - 0.0, - ZERO_THRESHOLD)) - { - hexes_device_view[i][j][0] = 0.0; - } - - if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][1], - 0.0, - ZERO_THRESHOLD)) - { - hexes_device_view[i][j][1] = 0.0; - } - - if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][2], - 0.0, - ZERO_THRESHOLD)) - { - hexes_device_view[i][j][2] = 0.0; - } - } - - // Get bounding box for hexahedral element - hex_bbs_device_view[i] = - primal::compute_bounding_box(hexes_device_view[i]); - }); // end of loop to initialize hexahedral elements and bounding boxes - // Set shape components to zero if within threshold - axom::for_all( - shape_count, - AXOM_LAMBDA(axom::IndexType i) { - for(int j = 0; j < ShapeType::NUM_VERTS; j++) - { - if(axom::utilities::isNearlyEqual(shapes_device_view[i][j][0], - 0.0, - ZERO_THRESHOLD)) - { - shapes_device_view[i][j][0] = 0.0; - } - - if(axom::utilities::isNearlyEqual(shapes_device_view[i][j][1], - 0.0, - ZERO_THRESHOLD)) - { - shapes_device_view[i][j][1] = 0.0; - } - - if(axom::utilities::isNearlyEqual(shapes_device_view[i][j][2], - 0.0, - ZERO_THRESHOLD)) - { - shapes_device_view[i][j][2] = 0.0; - } - } - }); + snapShapeVerticesToZero(shapes, + shape_count, + ZERO_THRESHOLD); // Find which shape bounding boxes intersect hexahedron bounding boxes SLIC_INFO(axom::fmt::format( "{:-^80}", " Finding shape candidates for each hexahedral element ")); - axom::Array offsets(NE, NE, device_allocator); - axom::Array counts(NE, NE, device_allocator); + axom::ArrayView hex_bbs_device_view = m_hex_bbs.view(); + + axom::Array offsets(m_cellCount, m_cellCount, device_allocator); + axom::Array counts(m_cellCount, m_cellCount, device_allocator); axom::Array candidates; - bvh.findBoundingBoxes(offsets, counts, candidates, NE, hex_bbs_device_view); + bvh.findBoundingBoxes(offsets, + counts, + candidates, + m_cellCount, + hex_bbs_device_view); // Get the total number of candidates using REDUCE_POL = typename axom::execution_space::reduce_policy; @@ -829,7 +793,7 @@ class IntersectionShaper : public Shaper const auto counts_device_view = counts.view(); RAJA::ReduceSum totalCandidates(0); axom::for_all( - NE, + m_cellCount, AXOM_LAMBDA(axom::IndexType i) { totalCandidates += counts_device_view[i]; }); @@ -848,11 +812,24 @@ class IntersectionShaper : public Shaper auto shape_candidates_device_view = shape_candidates_device.view(); // Tetrahedrons from hexes (24 for each hex) - axom::Array tets_from_hexes_device(NE * NUM_TETS_PER_HEX, - NE * NUM_TETS_PER_HEX, - device_allocator); + // NOTE: Why break hexes into tets if we can hex-tet clip (clip_impl.hpp:502). + // Possibly because ShapeType can also be oct, and we don't currently do hex-oct clip. + bool doInitializeTetsFromHexes = false; + if(m_tets_from_hexes_device.empty()) + { + m_tets_from_hexes_device = + axom::Array(m_cellCount * NUM_TETS_PER_HEX, + m_cellCount * NUM_TETS_PER_HEX, + device_allocator); + doInitializeTetsFromHexes = true; + } + else + { + SLIC_ASSERT(m_tets_from_hexes_device.size() == + m_cellCount * NUM_TETS_PER_HEX); + } axom::ArrayView tets_from_hexes_device_view = - tets_from_hexes_device.view(); + m_tets_from_hexes_device.view(); // Index into 'tets' axom::Array tet_indices_device( @@ -868,25 +845,28 @@ class IntersectionShaper : public Shaper axom::Array(newTotalCandidates_host, device_allocator); auto newTotalCandidates_device_view = newTotalCandidates_device.view(); - SLIC_INFO(axom::fmt::format( - "{:-^80}", - " Decomposing each hexahedron element into 24 tetrahedrons ")); + if(doInitializeTetsFromHexes) + { + SLIC_INFO(axom::fmt::format( + "{:-^80}", + " Decomposing each hexahedron element into 24 tetrahedrons ")); - using TetHexArray = axom::StackArray; + using TetHexArray = axom::StackArray; - { - AXOM_ANNOTATE_SCOPE("init_tets"); - axom::for_all( - NE, - AXOM_LAMBDA(axom::IndexType i) { - TetHexArray cur_tets; - hexes_device_view[i].triangulate(cur_tets); + { + AXOM_ANNOTATE_SCOPE("init_tets"); + axom::for_all( + m_cellCount, + AXOM_LAMBDA(axom::IndexType i) { + TetHexArray cur_tets; + hexes_device_view[i].triangulate(cur_tets); - for(int j = 0; j < NUM_TETS_PER_HEX; j++) - { - tets_from_hexes_device_view[i * NUM_TETS_PER_HEX + j] = cur_tets[j]; - } - }); + for(int j = 0; j < NUM_TETS_PER_HEX; j++) + { + tets_from_hexes_device_view[i * NUM_TETS_PER_HEX + j] = cur_tets[j]; + } + }); + } } SLIC_INFO( @@ -898,7 +878,7 @@ class IntersectionShaper : public Shaper { AXOM_ANNOTATE_SCOPE("init_candidates"); axom::for_all( - NE, + m_cellCount, AXOM_LAMBDA(axom::IndexType i) { for(int j = 0; j < counts_device_view[i]; j++) { @@ -919,30 +899,32 @@ class IntersectionShaper : public Shaper // Overlap volume is the volume of clip(oct,tet) for c2c // or clip(tet,tet) for Pro/E meshes - m_overlap_volumes = axom::Array(NE, NE, device_allocator); + if(m_overlap_volumes.empty()) + { + m_overlap_volumes = + axom::Array(m_cellCount, m_cellCount, device_allocator); + } + m_overlap_volumes.fill(0.0); // Hex volume is the volume of the hexahedron element - m_hex_volumes = axom::Array(NE, NE, device_allocator); + if(m_hex_volumes.empty()) + { + m_hex_volumes = + axom::Array(m_cellCount, m_cellCount, device_allocator); + m_hex_volumes.fill(0.0); + } axom::ArrayView overlap_volumes_device_view = m_overlap_volumes.view(); axom::ArrayView hex_volumes_device_view = m_hex_volumes.view(); - // Set initial values to 0 - axom::for_all( - NE, - AXOM_LAMBDA(axom::IndexType i) { - overlap_volumes_device_view[i] = 0; - hex_volumes_device_view[i] = 0; - }); - SLIC_INFO( axom::fmt::format("{:-^80}", " Calculating hexahedron element volume ")); { AXOM_ANNOTATE_SCOPE("hex_volume"); axom::for_all( - NE, + m_cellCount, AXOM_LAMBDA(axom::IndexType i) { hex_volumes_device_view[i] = hexes_device_view[i].volume(); }); @@ -955,13 +937,13 @@ class IntersectionShaper : public Shaper constexpr bool tryFixOrientation = true; { + AXOM_ANNOTATE_SCOPE("IntersectionShaper::clipLoop"); // Copy calculated total back to host axom::Array newTotalCandidates_calc_host = axom::Array(newTotalCandidates_device, host_allocator); - AXOM_ANNOTATE_SCOPE("tet_shape_volume"); axom::for_all( - newTotalCandidates_calc_host[0], + newTotalCandidates_calc_host[0], // Number of candidates found. AXOM_LAMBDA(axom::IndexType i) { const int index = hex_indices_device_view[i]; const int shapeIndex = shape_candidates_device_view[i]; @@ -989,7 +971,7 @@ class IntersectionShaper : public Shaper RAJA::ReduceSum totalHex(0); axom::for_all( - NE, + m_cellCount, AXOM_LAMBDA(axom::IndexType i) { totalOverlap += overlap_volumes_device_view[i]; totalHex += hex_volumes_device_view[i]; @@ -1001,10 +983,8 @@ class IntersectionShaper : public Shaper SLIC_INFO(axom::fmt::format(axom::utilities::locale(), "Total mesh volume is {:.3Lf}", this->allReduceSum(totalHex))); - } // end of runShapeQuery() function -#endif + } // end of runShapeQueryImpl() function -#if defined(AXOM_USE_RAJA) // These methods are private in support of replacement rules. private: /*! @@ -1039,48 +1019,6 @@ class IntersectionShaper : public Shaper return name; } - /*! - * \brief Gets the grid function and material number for a material name. - * - * \param materialName The name of the material. - * - * \return A pair containing the associated grid function and material - * number (its order in the list). - */ - std::pair getMaterial(const std::string& materialName) - { - // If we already know about the material, return it. - for(size_t i = 0; i < m_vf_material_names.size(); i++) - { - if(m_vf_material_names[i] == materialName) - { - return std::make_pair(m_vf_grid_functions[i], static_cast(i)); - } - } - - // Get or create the volume fraction field for this shape's material - auto materialVolFracName = materialNameToFieldName(materialName); - mfem::GridFunction* matVolFrac = nullptr; - if(this->getDC()->HasField(materialVolFracName)) - { - matVolFrac = this->getDC()->GetField(materialVolFracName); - } - else - { - matVolFrac = newVolFracGridFunction(); - this->getDC()->RegisterField(materialVolFracName, matVolFrac); - // Zero out the volume fractions (on host). - memset(matVolFrac->begin(), 0, matVolFrac->Size() * sizeof(double)); - } - - // Add the material to our vectors. - int idx = static_cast(m_vf_grid_functions.size()); - m_vf_grid_functions.push_back(matVolFrac); - m_vf_material_names.push_back(materialName); - - return std::make_pair(matVolFrac, idx); - } - /*! * \brief Scans the grid functions in the data collection and creates * a material entry for any that do not already exist. We maintain @@ -1091,15 +1029,7 @@ class IntersectionShaper : public Shaper */ void populateMaterials() { - std::vector materialNames; - for(auto it : this->getDC()->GetFieldMap()) - { - std::string materialName = fieldNameToMaterialName(it.first); - if(!materialName.empty()) - { - materialNames.emplace_back(materialName); - } - } + std::vector materialNames = getMaterialNames(); // Add any of these existing fields to this class' bookkeeping. for(const auto& materialName : materialNames) { @@ -1124,34 +1054,31 @@ class IntersectionShaper : public Shaper * free space in each zone. */ template - mfem::GridFunction* getCompletelyFree() + axom::ArrayView getCompletelyFree() { - // Add the material prefix so the MFEMSidreDataCollection will automatically + // Add the material prefix so the mesh will automatically // consider the free material something it needs to write as a matset. const std::string fieldName(materialNameToFieldName(m_free_mat_name)); - mfem::GridFunction* cfgf = nullptr; - if(this->getDC()->HasField(fieldName)) - { - cfgf = this->getDC()->GetField(fieldName); - } - else - { - // Make the new grid function. - cfgf = newVolFracGridFunction(); - this->getDC()->RegisterField(fieldName, cfgf); + bool makeNewData = !hasData(fieldName); + axom::ArrayView cfgf = getScalarCellData(fieldName); + SLIC_ASSERT(!cfgf.empty()); + + if(makeNewData) + { AXOM_ANNOTATE_SCOPE("compute_free"); - int dataSize = cfgf->Size(); - GridFunctionView cfView(cfgf); + int dataSize = cfgf.size(); + TempArrayView cfView(cfgf, true); + axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { cfView[i] = 1.; }); // Iterate over all materials and subtract off their VFs from cfgf. - for(auto& gf : m_vf_grid_functions) + for(axom::ArrayView& gf : m_vf_grid_functions) { - GridFunctionView matVFView(gf, false); + TempArrayView matVFView(gf, false); axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { @@ -1160,6 +1087,7 @@ class IntersectionShaper : public Shaper }); } } + return cfgf; } @@ -1220,17 +1148,24 @@ class IntersectionShaper : public Shaper populateMaterials(); // Get the free material so it is created first. - mfem::GridFunction* freeMat = getCompletelyFree(); + // mfem::GridFunction* freeMat = getCompletelyFree(); + axom::ArrayView freeMat = getCompletelyFree(); // Get this shape's material, creating the GridFunction if needed. auto matVF = getMaterial(shape.getMaterial()); - int dataSize = matVF.first->Size(); + // int dataSize = matVF.first->Size(); + int dataSize = matVF.first.size(); // Get this shape's array. auto shapeVolFracName = axom::fmt::format("shape_vol_frac_{}", shape.getName()); - auto* shapeVolFrac = this->getDC()->GetField(shapeVolFracName); - SLIC_ASSERT(shapeVolFrac != nullptr); + // auto* shapeVolFrac = this->getDC()->GetField(shapeVolFracName); + auto shapeVolFrac = getScalarCellData(shapeVolFracName); + SLIC_ERROR_IF(shapeVolFrac.empty(), + "Field '" + shapeVolFracName + + "' must be pre-allocated" + " in the Conduit Node computational mesh before using" + " IntersectionShaper::applyReplacementRules."); // Allocate some memory for the replacement rule data arrays. int execSpaceAllocatorID = axom::execution_space::allocatorID(); @@ -1241,8 +1176,10 @@ class IntersectionShaper : public Shaper ArrayView vf_writable(vf_writable_array); // Determine which grid functions need to be considered for VF updates. - std::vector> gf_order_by_matnumber; - std::vector updateVFs, excludeVFs; + // std::vector> gf_order_by_matnumber; + // std::vector updateVFs, excludeVFs; + std::vector, int>> gf_order_by_matnumber; + std::vector> updateVFs, excludeVFs; if(!shape.getMaterialsReplaced().empty()) { // Include materials replaced in updateVFs. @@ -1256,12 +1193,12 @@ class IntersectionShaper : public Shaper // Include all materials except those in "does_not_replace". // We'll also sort them by material number since the field map // sorts them by name rather than order added. - for(auto it : this->getDC()->GetFieldMap()) + std::vector materialNames = getMaterialNames(); + for(auto name : materialNames) { - // Check whether the field name looks like a VF field (and is not the + // Check whether the field name is not the // "free" field, which we handle specially) - std::string name = fieldNameToMaterialName(it.first); - if(!name.empty() && name != m_free_mat_name) + if(name != m_free_mat_name) { // See if the field is in the exclusion list. For the normal // case, the list is empty so we'd add the material. @@ -1291,8 +1228,8 @@ class IntersectionShaper : public Shaper // Sort eligible update materials by material number. std::sort(gf_order_by_matnumber.begin(), gf_order_by_matnumber.end(), - [&](const std::pair& lhs, - const std::pair& rhs) { + [&](const std::pair, int>& lhs, + const std::pair, int>& rhs) { return lhs.second < rhs.second; }); @@ -1320,10 +1257,11 @@ class IntersectionShaper : public Shaper axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { vf_writable[i] = 0.; }); + for(const auto& name : shape.getMaterialsReplaced()) { auto mat = getMaterial(name); - GridFunctionView matVFView(mat.first, false); + TempArrayView matVFView(mat.first, false); axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { @@ -1342,7 +1280,7 @@ class IntersectionShaper : public Shaper for(auto& gf : excludeVFs) { - GridFunctionView matVFView(gf, false); + TempArrayView matVFView(gf, false); axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { @@ -1356,8 +1294,8 @@ class IntersectionShaper : public Shaper { AXOM_ANNOTATE_SCOPE("compute_vf"); - GridFunctionView matVFView(matVF.first); - GridFunctionView shapeVFView(shapeVolFrac); + TempArrayView matVFView(matVF.first, true); + TempArrayView shapeVFView(shapeVolFrac, true); axom::ArrayView overlap_volumes_view = m_overlap_volumes.view(); axom::ArrayView hex_volumes_view = m_hex_volumes.view(); @@ -1387,7 +1325,7 @@ class IntersectionShaper : public Shaper AXOM_ANNOTATE_SCOPE("update_vf"); for(auto& gf : updateVFs) { - GridFunctionView matVFView(gf); + TempArrayView matVFView(gf, true); axom::for_all( dataSize, AXOM_LAMBDA(axom::IndexType i) { @@ -1403,7 +1341,6 @@ class IntersectionShaper : public Shaper } } } -#endif /*! * \brief Apply material replacement rules for the current shape, using @@ -1415,7 +1352,6 @@ class IntersectionShaper : public Shaper switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: applyReplacementRulesImpl(shape); break; @@ -1434,7 +1370,6 @@ class IntersectionShaper : public Shaper applyReplacementRulesImpl(shape); break; #endif // AXOM_USE_HIP -#endif // AXOM_USE_RAJA && AXOM_USE_UMPIRE } AXOM_UNUSED_VAR(shape); } @@ -1450,8 +1385,13 @@ class IntersectionShaper : public Shaper //@} public: - // Prepares the shaping query, based on the policy member set - // (default is sequential) + /*! + \brief Prepares the shaping query, based on the policy member set + (default is sequential) + + \internal This method populates m_tets or m_octs from the given \c + shape. These arrays are used in runShapeQuery. + */ void prepareShapeQuery(klee::Dimensions shapeDimension, const klee::Shape& shape) override { @@ -1471,7 +1411,6 @@ class IntersectionShaper : public Shaper // Now that the mesh is refined, dispatch to device implementations. switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: prepareShapeQueryImpl(shapeDimension, shape); break; @@ -1490,7 +1429,6 @@ class IntersectionShaper : public Shaper prepareShapeQueryImpl(shapeDimension, shape); break; #endif // AXOM_USE_HIP -#endif // AXOM_USE_RAJA && AXOM_USE_UMPIRE } AXOM_UNUSED_VAR(shapeDimension); AXOM_UNUSED_VAR(shape); @@ -1502,6 +1440,10 @@ class IntersectionShaper : public Shaper // Runs the shaping query, based on the policy member and shape format set // (default is sequential) + // Fills m_overlap_volumes and m_hex_volumes, whose data + // will be in the default memory space for m_execPolicy. + // The data will be used in applyReplacementRules and can + // also be accessed by getOverlapVolumes() and getHexVolumes(). void runShapeQuery(const klee::Shape& shape) override { AXOM_ANNOTATE_SCOPE("runShapeQuery"); @@ -1512,7 +1454,6 @@ class IntersectionShaper : public Shaper { switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: runShapeQueryImpl(shape, m_tets, m_tetcount); break; @@ -1531,14 +1472,12 @@ class IntersectionShaper : public Shaper runShapeQueryImpl(shape, m_tets, m_tetcount); break; #endif // AXOM_USE_HIP -#endif // AXOM_USE_RAJA && AXOM_USE_UMPIRE } } else if(shapeFormat == "c2c") { switch(m_execPolicy) { -#if defined(AXOM_USE_RAJA) case RuntimePolicy::seq: runShapeQueryImpl(shape, m_octs, m_octcount); break; @@ -1557,7 +1496,6 @@ class IntersectionShaper : public Shaper runShapeQueryImpl(shape, m_octs, m_octcount); break; #endif // AXOM_USE_HIP -#endif // AXOM_USE_RAJA && AXOM_USE_UMPIRE } } else @@ -1567,11 +1505,177 @@ class IntersectionShaper : public Shaper } } + axom::ArrayView getOverlapVolumes() const + { + return m_overlap_volumes.view(); + } + + axom::ArrayView getHexVolumes() const + { + return m_hex_volumes.view(); + } + + double sumOverlapVolumes(bool global = true) const + { + double overlapVol = 0.0; + switch(m_execPolicy) + { + #if defined(AXOM_USE_OPENMP) + case RuntimePolicy::omp: + overlapVol = + sumArray(m_overlap_volumes.data(), m_overlap_volumes.size()); + break; + #endif // AXOM_USE_OPENMP + #if defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE) + case RuntimePolicy::cuda: + overlapVol = + sumArray(m_overlap_volumes.data(), m_overlap_volumes.size()); + break; + #endif // AXOM_USE_CUDA + #if defined(AXOM_USE_HIP) && defined(AXOM_USE_UMPIRE) + case RuntimePolicy::hip: + overlapVol = + sumArray(m_overlap_volumes.data(), m_overlap_volumes.size()); + break; + #endif // AXOM_USE_HIP + case RuntimePolicy::seq: + default: + overlapVol = + sumArray(m_overlap_volumes.data(), m_overlap_volumes.size()); + break; + } + + if(global) + { + overlapVol = this->allReduceSum(overlapVol); + } + return overlapVol; + } + + template + Summable sumArray(const Summable* a, axom::IndexType count) const + { + using LoopPolicy = typename axom::execution_space::loop_policy; + using ReducePolicy = typename axom::execution_space::reduce_policy; + RAJA::ReduceSum vsum {0}; + RAJA::forall( + RAJA::RangeSegment(0, count), + AXOM_LAMBDA(RAJA::Index_type i) { vsum += a[i]; }); + Summable sum = static_cast(vsum.get()); + return sum; + } + void adjustVolumeFractions() override { // Implementation here -- not sure if this will require anything for intersection-based shaping } + std::vector getMaterialNames() const + { + std::vector materialNames; + #if defined(AXOM_USE_MFEM) + if(m_dc) + { + for(auto it : this->getDC()->GetFieldMap()) + { + std::string materialName = fieldNameToMaterialName(it.first); + if(!materialName.empty()) + { + materialNames.emplace_back(materialName); + } + } + } + #endif + #if defined(AXOM_USE_CONDUIT) + if(m_bpGrp) + { + auto fieldsGrp = m_bpGrp->getGroup("fields"); + SLIC_ERROR_IF(fieldsGrp == nullptr, + "Input blueprint mesh lacks the 'fields' Group/Node."); + for(auto& group : fieldsGrp->groups()) + { + std::string materialName = fieldNameToMaterialName(group.getName()); + if(!materialName.empty()) + { + materialNames.emplace_back(materialName); + } + } + } + #endif + return materialNames; + } + + /*! + * \brief Gets the grid function and material number for a material name. + * + * \param materialName The name of the material. + * + * \return A pair containing the associated grid function (an element + * in the m_vf_grid_functions array) and material + * number (its index in the array). + */ + std::pair, int> getMaterial(const std::string& materialName) + { + // If we already know about the material, return it. + for(size_t i = 0; i < m_vf_material_names.size(); i++) + { + if(m_vf_material_names[i] == materialName) + { + return std::make_pair(m_vf_grid_functions[i], static_cast(i)); + } + } + + // Get or create the volume fraction field for this shape's material + auto materialVolFracName = materialNameToFieldName(materialName); + + bool makeNewData = !hasData(materialVolFracName); + auto matVolFrac = getScalarCellData(materialVolFracName); + SLIC_ASSERT(!matVolFrac.empty()); + if(makeNewData) + { + // Zero out the volume fractions (on host). + #ifdef AXOM_USE_UMPIRE + auto allocId = matVolFrac.getAllocatorID(); + const axom::MemorySpace memorySpace = + axom::detail::getAllocatorSpace(allocId); + const bool onDevice = memorySpace == axom::MemorySpace::Device || + memorySpace == axom::MemorySpace::Unified; + #else + const bool onDevice = false; + #endif + if(onDevice) + { + #if defined(AXOM_USE_CUDA) + if(m_execPolicy == RuntimePolicy::cuda) + { + axom::for_all>( + matVolFrac.size(), + AXOM_LAMBDA(axom::IndexType i) { matVolFrac[i] = 0.0; }); + } + #endif + #if defined(AXOM_USE_HIP) + if(m_execPolicy == RuntimePolicy::hip) + { + axom::for_all>( + matVolFrac.size(), + AXOM_LAMBDA(axom::IndexType i) { matVolFrac[i] = 0.0; }); + } + #endif + } + else + { + memset(matVolFrac.data(), 0, matVolFrac.size() * sizeof(double)); + } + } + + // Add the material to our vectors. + int idx = static_cast(m_vf_grid_functions.size()); + m_vf_grid_functions.push_back(matVolFrac); + m_vf_material_names.push_back(materialName); + + return std::make_pair(matVolFrac, idx); + } + private: /*! * \brief Filter the input mesh so it does not have any degenerate segments @@ -2027,6 +2131,381 @@ class IntersectionShaper : public Shaper m_level = circleLevel; } + /// Whether the given field exists in the mesh. + bool hasData(const std::string& fieldName) + { + bool has = false; + #if defined(AXOM_USE_MFEM) + if(m_dc != nullptr) + { + has = m_dc->HasField(fieldName); + } + #endif + #if defined(AXOM_USE_CONDUIT) + if(m_bpGrp != nullptr) + { + std::string fieldPath = axom::fmt::format("fields/{}", fieldName); + has = m_bpGrp->hasGroup(fieldPath); + } + #endif + return has; + } + + /*! \brief Get a scalar double-type field data from the mesh, + "fields/fieldName/values", creating it if it doesn't exist. + + Also add the corresponding entry in the blueprint field + "matsets/fieldName/volume_fractions". + + If mesh is in an external Conduit Node, and the field + doesn't exist, emit a warning and return an empty ArrayView. + We don't add fields to the external Conduit Node. + \see Shaper::Shaper(). + */ + axom::ArrayView getScalarCellData(const std::string& fieldName, + bool volumeDependent = false) + { + axom::ArrayView rval; + + #if defined(AXOM_USE_MFEM) + if(m_dc != nullptr) + { + mfem::GridFunction* gridFunc = nullptr; + if(m_dc->HasField(fieldName)) + { + gridFunc = m_dc->GetField(fieldName); + } + else + { + gridFunc = newVolFracGridFunction(); + m_dc->RegisterField(fieldName, gridFunc); + } + rval = axom::ArrayView(gridFunc->GetData(), gridFunc->Size()); + } + #endif + + #if defined(AXOM_USE_CONDUIT) + if(m_bpGrp != nullptr) + { + std::string fieldPath = "fields/" + fieldName; + auto dtype = conduit::DataType::float64(m_cellCount); + axom::sidre::View* valuesView = nullptr; + if(m_bpGrp->hasGroup(fieldPath)) + { + auto* fieldGrp = m_bpGrp->getGroup(fieldPath); + valuesView = fieldGrp->getView("values"); + SLIC_ASSERT(fieldGrp->getView("association")->getString() == + std::string("element")); + SLIC_ASSERT(fieldGrp->getView("topology")->getString() == m_bpTopo); + SLIC_ASSERT(valuesView->getNumElements() == m_cellCount); + SLIC_ASSERT(valuesView->getNode().dtype().id() == dtype.id()); + } + else + { + if(m_bpNodeExt != nullptr) + { + /* + If the computational mesh is an external conduit::Node, it + must have all necessary fields. We will only generate + fields for meshes in sidre::Group, where the user can set + the allocator id for only array data. conduit::Node doesn't + have this capability. + */ + SLIC_WARNING_IF(m_bpNodeExt != nullptr, + "For a computational mesh in a conduit::Node, all" + " output fields must be preallocated before shaping." + " IntersectionShaper will NOT contravene the user's" + " memory management. The cell-centered field '" + + fieldPath + + "' is missing. Please pre-allocate" + " this output memory, or to have IntersectionShaper" + " allocate it, construct the IntersectionShaper" + " with the mesh as a sidre::Group with your" + " specific allocator id."); + } + else + { + constexpr axom::IndexType componentCount = 1; + axom::IndexType shape[2] = {m_cellCount, componentCount}; + auto* fieldGrp = m_bpGrp->createGroup(fieldPath); + // valuesView = fieldGrp->createView("values"); + valuesView = + fieldGrp->createViewWithShape("values", + axom::sidre::DataTypeId::FLOAT64_ID, + 2, + shape); + fieldGrp->createView("association")->setString("element"); + fieldGrp->createView("topology")->setString(m_bpTopo); + fieldGrp->createView("volume_dependent") + ->setString(std::string(volumeDependent ? "true" : "false")); + valuesView->allocate(); + } + } + + rval = + axom::ArrayView(static_cast(valuesView->getVoidPtr()), + m_cellCount); + } + #endif + return rval; + } + + #if defined(__CUDACC__) +public: + // These methods should be private, but NVCC complains unless they're public. + #endif + + template + void populateHexesFromMesh() + { + using XS = axom::execution_space; + + constexpr int NUM_VERTS_PER_HEX = 8; + constexpr int NUM_COMPS_PER_VERT = 3; + const int allocId = XS::allocatorID(); + + axom::Array vertCoords( + m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, + m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, + allocId); + + #if defined(AXOM_USE_MFEM) + if(m_dc != nullptr) + { + populateVertCoordsFromMFEMMesh(vertCoords); + } + #endif + #if defined(AXOM_USE_CONDUIT) + if(m_bpGrp != nullptr) + { + populateVertCoordsFromBlueprintMesh(vertCoords); + } + #endif + + auto vertCoords_device_view = vertCoords.view(); + + SLIC_ASSERT(m_hexes.empty()); // No need to repeat this work. + m_hexes = axom::Array(m_cellCount, m_cellCount, allocId); + axom::ArrayView hexes_device_view = m_hexes.view(); + constexpr double ZERO_THRESHOLD = 1.e-10; + axom::for_all( + m_cellCount, + AXOM_LAMBDA(axom::IndexType i) { + // Set each hexahedral element vertices + hexes_device_view[i] = HexahedronType(); + for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) + { + int vertIndex = (i * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT) + + j * NUM_COMPS_PER_VERT; + hexes_device_view[i][j] = + Point3D({vertCoords_device_view[vertIndex], + vertCoords_device_view[vertIndex + 1], + vertCoords_device_view[vertIndex + 2]}); + + // Set hexahedra components to zero if within threshold + if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][0], + 0.0, + ZERO_THRESHOLD)) + { + hexes_device_view[i][j][0] = 0.0; + } + + if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][1], + 0.0, + ZERO_THRESHOLD)) + { + hexes_device_view[i][j][1] = 0.0; + } + + if(axom::utilities::isNearlyEqual(hexes_device_view[i][j][2], + 0.0, + ZERO_THRESHOLD)) + { + hexes_device_view[i][j][2] = 0.0; + } + } + }); // end of loop to initialize hexahedral elements and bounding boxes + } + + //!\brief Set shape vertices to zero of within threshold. + template + void snapShapeVerticesToZero(axom::Array& shapes, + axom::IndexType shapeCount, + double zeroThreshold) + { + axom::ArrayView shapesView = shapes.view(); + axom::for_all( + shapeCount, + AXOM_LAMBDA(axom::IndexType i) { + for(int j = 0; j < ShapeType::NUM_VERTS; j++) + { + if(axom::utilities::isNearlyEqual(shapesView[i][j][0], 0.0, zeroThreshold)) + { + shapesView[i][j][0] = 0.0; + } + + if(axom::utilities::isNearlyEqual(shapesView[i][j][1], 0.0, zeroThreshold)) + { + shapesView[i][j][1] = 0.0; + } + + if(axom::utilities::isNearlyEqual(shapesView[i][j][2], 0.0, zeroThreshold)) + { + shapesView[i][j][2] = 0.0; + } + } + }); + } + + #if defined(AXOM_USE_CONDUIT) + template + void populateVertCoordsFromBlueprintMesh(axom::Array& vertCoords) + { + using XS = axom::execution_space; + + const int allocId = XS::allocatorID(); + + // Initialize vertices from blueprint mesh and + // set each shape volume fraction to 1 + // Allocation size is: + // # of elements * # of vertices per hex * # of components per vertex + constexpr int NUM_VERTS_PER_HEX = 8; + constexpr int NUM_COMPS_PER_VERT = 3; + + // Put mesh in Node so we can use conduit::blueprint utilities. + // conduit::Node meshNode; + // m_bpGrp->createNativeLayout(m_bpNodeInt); + + const conduit::Node& topoNode = + m_bpNodeInt.fetch_existing("topologies").fetch_existing(m_bpTopo); + const std::string coordsetName = + topoNode.fetch_existing("coordset").as_string(); + + // Assume unstructured and hexahedral + SLIC_ERROR_IF(topoNode["type"].as_string() != "unstructured", + "topology type must be 'unstructured'"); + SLIC_ERROR_IF(topoNode["elements/shape"].as_string() != "hex", + "element shape must be 'hex'"); + + const auto& connNode = topoNode["elements/connectivity"]; + SLIC_ERROR_IF( + !XS::usesAllocId(axom::getAllocatorIDFromPointer(connNode.data_ptr())), + std::string(XS::name()) + + axom::fmt::format( + " execution space cannot use the connectivity allocator id {}", + axom::getAllocatorIDFromPointer(connNode.data_ptr()))); + SLIC_ERROR_IF(connNode.dtype().id() != conduitDataIdOfAxomIndexType, + "IntersectionShaper error: connectivity data type must be " + "axom::IndexType."); + const auto* connPtr = + static_cast(connNode.data_ptr()); + axom::ArrayView conn(connPtr, + m_cellCount, + NUM_VERTS_PER_HEX); + + const conduit::Node& coordNode = m_bpNodeInt["coordsets"][coordsetName]; + const conduit::Node& coordValues = coordNode.fetch_existing("values"); + axom::IndexType vertexCount = coordValues["x"].dtype().number_of_elements(); + bool isInterleaved = conduit::blueprint::mcarray::is_interleaved(coordValues); + int stride = isInterleaved ? NUM_COMPS_PER_VERT : 1; + + axom::StackArray, 3> coordArrays { + axom::ArrayView(coordValues["x"].as_double_ptr(), + {vertexCount}, + stride), + axom::ArrayView(coordValues["y"].as_double_ptr(), + {vertexCount}, + stride), + axom::ArrayView(coordValues["z"].as_double_ptr(), + {vertexCount}, + stride)}; + + vertCoords = + axom::Array(m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, + m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, + allocId); + auto vertCoordsView = vertCoords.view(); + + axom::for_all( + m_cellCount, + AXOM_LAMBDA(axom::IndexType i) { + // Get the indices of this element's vertices + auto hexVerts = conn[i]; + + // Get the coordinates for the vertices + for(int j = 0; j < NUM_VERTS_PER_HEX; ++j) + { + auto vertId = hexVerts[j]; + for(int k = 0; k < NUM_COMPS_PER_VERT; k++) + { + vertCoordsView[(i * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT) + + (j * NUM_COMPS_PER_VERT) + k] = coordArrays[k][vertId]; + } + } + }); + } + #endif // AXOM_USE_CONDUIT + + #if defined(AXOM_USE_MFEM) + template + void populateVertCoordsFromMFEMMesh(axom::Array& vertCoords) + { + mfem::Mesh* mesh = getDC()->GetMesh(); + // Intersection algorithm only works on linear elements + SLIC_ASSERT(mesh != nullptr); + + if(m_cellCount > 0) + { + SLIC_ASSERT(mesh->GetNodes() == nullptr || + mesh->GetNodes()->FESpace()->GetOrder(0)); + } + + // Allocation size is: + // # of elements * # of vertices per hex * # of components per vertex + constexpr int NUM_VERTS_PER_HEX = 8; + constexpr int NUM_COMPS_PER_VERT = 3; + + // The MFEM mesh interface works only on host. + // If on device, fill temporary host array then copy to device. + axom::Array tmpVertCoords; + + axom::Array& fillVertCoords = + axom::execution_space::onDevice() ? tmpVertCoords : vertCoords; + fillVertCoords = + axom::Array(m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT, + m_cellCount * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT); + + // Initialize vertices from mfem mesh and + // set each shape volume fraction to 1 + + auto fillVertCoordsView = fillVertCoords.view(); + for(int i = 0; i < m_cellCount; i++) + { + // Get the indices of this element's vertices + mfem::Array verts; + mesh->GetElementVertices(i, 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) + { + for(int k = 0; k < NUM_COMPS_PER_VERT; k++) + { + fillVertCoordsView[(i * NUM_VERTS_PER_HEX * NUM_COMPS_PER_VERT) + + (j * NUM_COMPS_PER_VERT) + k] = + (mesh->GetVertex(verts[j]))[k]; + } + } + } + if(vertCoords.data() != fillVertCoords.data()) + { + axom::copy(vertCoords.data(), + fillVertCoords.data(), + sizeof(double) * vertCoords.size()); + } + } + +private: /// Create and return a new volume fraction grid function for the current mesh mfem::GridFunction* newVolFracGridFunction() { @@ -2043,6 +2522,7 @@ class IntersectionShaper : public Shaper return volFrac; } + #endif // AXOM_USE_MFEM bool surfaceMeshIsTet() const { @@ -2053,16 +2533,18 @@ class IntersectionShaper : public Shaper } private: - RuntimePolicy m_execPolicy {RuntimePolicy::seq}; int m_level {DEFAULT_CIRCLE_REFINEMENT_LEVEL}; double m_revolvedVolume {DEFAULT_REVOLVED_VOLUME}; - int m_num_elements {0}; std::string m_free_mat_name; + //! \brief Volumes of cells in the computational mesh. axom::Array m_hex_volumes; + + //! \brief Overlap volumes of cells in the computational mesh for the last shape. axom::Array m_overlap_volumes; -#if defined(AXOM_USE_RAJA) + double m_vertexWeldThreshold {1.e-10}; + // Guard these to prevent warnings. int m_octcount {0}; int m_tetcount {0}; @@ -2072,13 +2554,16 @@ class IntersectionShaper : public Shaper axom::Array m_aabbs; axom::Array m_hexes; axom::Array m_hex_bbs; + axom::Array m_tets_from_hexes_device; - std::vector m_vf_grid_functions; + // Views of volume-fraction data owned by grid. + std::vector> m_vf_grid_functions; std::vector m_vf_material_names; -#endif }; } // end namespace quest } // end namespace axom +#endif // AXOM_USE_RAJA && AXOM_USE_UMPIRE + #endif // AXOM_QUEST_INTERSECTION_SHAPER__HPP_ diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index 360f1f0d79..63e2c0bdad 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -318,7 +318,10 @@ class SamplingShaper : public Shaper public: SamplingShaper(const klee::ShapeSet& shapeSet, sidre::MFEMSidreDataCollection* dc) - : Shaper(shapeSet, dc) + : Shaper(Shaper::RuntimePolicy::seq, + axom::policyToDefaultAllocatorID(axom::runtime_policy::Policy::seq), + shapeSet, + dc) { } ~SamplingShaper() diff --git a/src/axom/quest/Shaper.cpp b/src/axom/quest/Shaper.cpp index c3f0a9ca3e..b3de404cd4 100644 --- a/src/axom/quest/Shaper.cpp +++ b/src/axom/quest/Shaper.cpp @@ -9,16 +9,13 @@ #include "axom/core.hpp" #include "axom/primal/operators/split.hpp" #include "axom/quest/interface/internal/QuestHelpers.hpp" +#include "axom/quest/Shaper.hpp" #include "axom/quest/DiscreteShape.hpp" +#include "axom/quest/util/mesh_helpers.hpp" +#include "conduit_blueprint_mesh.hpp" #include "axom/fmt.hpp" -#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" - namespace axom { namespace quest @@ -30,13 +27,127 @@ constexpr double Shaper::MINIMUM_PERCENT_ERROR; constexpr double Shaper::MAXIMUM_PERCENT_ERROR; constexpr double Shaper::DEFAULT_VERTEX_WELD_THRESHOLD; -Shaper::Shaper(const klee::ShapeSet& shapeSet, sidre::MFEMSidreDataCollection* dc) - : m_shapeSet(shapeSet) +#if defined(AXOM_USE_MFEM) +Shaper::Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + sidre::MFEMSidreDataCollection* dc) + : m_execPolicy(execPolicy) + , m_allocatorId(allocatorId) + , m_shapeSet(shapeSet) , m_dc(dc) + #if defined(AXOM_USE_CONDUIT) + , m_bpGrp(nullptr) + , m_bpTopo() + , m_bpNodeExt(nullptr) + , m_bpNodeInt() + #endif { -#if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) + #if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) m_comm = m_dc->GetComm(); + #endif + m_cellCount = m_dc->GetMesh()->GetNE(); + + setFilePath(shapeSet.getPath()); +} +#endif + +Shaper::Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + sidre::Group* bpGrp, + const std::string& topo) + : m_execPolicy(execPolicy) + , m_allocatorId(allocatorId) + , m_shapeSet(shapeSet) +#if defined(AXOM_USE_CONDUIT) + , m_bpGrp(bpGrp) + , m_bpTopo(topo.empty() ? bpGrp->getGroup("topologies")->getGroupName(0) : topo) + , m_bpNodeExt(nullptr) + , m_bpNodeInt() +#endif + , m_comm(MPI_COMM_WORLD) +{ + SLIC_ASSERT(m_bpTopo != sidre::InvalidName); + + // This may take too long if there are repeated construction. + m_bpGrp->createNativeLayout(m_bpNodeInt); + +#if defined(AXOM_DEBUG) + std::string whyBad; + bool goodMesh = verifyInputMesh(whyBad); + SLIC_ASSERT_MSG(goodMesh, whyBad); +#endif + + m_cellCount = conduit::blueprint::mesh::topology::length( + m_bpNodeInt.fetch_existing("topologies").fetch_existing(m_bpTopo)); + + setFilePath(shapeSet.getPath()); +} + +Shaper::Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + conduit::Node* bpNode, + const std::string& topo) + : m_execPolicy(execPolicy) + , m_allocatorId(allocatorId) + , m_shapeSet(shapeSet) +#if defined(AXOM_USE_CONDUIT) + , m_bpGrp(nullptr) + , m_bpTopo(topo.empty() ? bpNode->fetch_existing("topologies").child(0).name() + : topo) + , m_bpNodeExt(bpNode) + , m_bpNodeInt() #endif + , m_comm(MPI_COMM_WORLD) +{ + AXOM_ANNOTATE_SCOPE("Shaper::Shaper_Node"); + m_bpGrp = m_ds.getRoot()->createGroup("internalGrp"); + m_bpGrp->setDefaultAllocator(policyToDefaultAllocatorID(m_execPolicy)); + + m_bpGrp->importConduitTreeExternal(*bpNode); + + // We want unstructured topo but can accomodate structured. + const std::string topoType = bpNode->fetch_existing("topologies") + .fetch_existing(m_bpTopo) + .fetch_existing("type") + .as_string(); + if(topoType == "structured") + { + AXOM_ANNOTATE_SCOPE("Shaper::convertStructured"); + axom::quest::util::convert_blueprint_structured_explicit_to_unstructured( + m_bpGrp, + m_bpTopo, + m_execPolicy); + } + + m_bpGrp->createNativeLayout(m_bpNodeInt); + +#if defined(AXOM_DEBUG) + std::string whyBad; + bool goodMesh = verifyInputMesh(whyBad); + SLIC_ASSERT_MSG(goodMesh, whyBad); +#endif + + m_cellCount = conduit::blueprint::mesh::topology::length( + bpNode->fetch_existing("topologies").fetch_existing(m_bpTopo)); + + setFilePath(shapeSet.getPath()); +} + +Shaper::~Shaper() { } + +void Shaper::setFilePath(const std::string& filePath) +{ + if(filePath.empty()) + { + m_prefixPath.clear(); + } + else + { + m_prefixPath = utilities::filesystem::getParentPath(filePath); + } } void Shaper::setSamplesPerKnotSpan(int nSamples) @@ -94,7 +205,7 @@ bool Shaper::isValidFormat(const std::string& format) const return (format == "stl" || format == "proe" || format == "c2c" || format == "blueprint-tets" || format == "tet3D" || format == "hex3D" || format == "plane3D" || format == "sphere3D" || - format == "vor3D" || format == "none"); + format == "sor3D" || format == "none"); } void Shaper::loadShape(const klee::Shape& shape) @@ -110,6 +221,7 @@ void Shaper::loadShapeInternal(const klee::Shape& shape, double percentError, double& revolvedVolume) { + AXOM_ANNOTATE_SCOPE("Shaper::loadShapeInternal"); internal::ScopedLogLevelChanger logLevelChanger( this->isVerbose() ? slic::message::Debug : slic::message::Warning); @@ -117,14 +229,12 @@ void Shaper::loadShapeInternal(const klee::Shape& shape, "{:-^80}", axom::fmt::format(" Loading shape '{}' ", shape.getName()))); - const axom::klee::Geometry& geometry = shape.getGeometry(); - const std::string& geometryFormat = geometry.getFormat(); - SLIC_ASSERT_MSG( - this->isValidFormat(geometryFormat), - axom::fmt::format("Shape has unsupported format: '{}", geometryFormat)); + SLIC_ASSERT_MSG(this->isValidFormat(shape.getGeometry().getFormat()), + axom::fmt::format("Shape has unsupported format: '{}", + shape.getGeometry().getFormat())); // Code for discretizing shapes has been factored into DiscreteShape class. - DiscreteShape discreteShape(shape, m_dataStore.getRoot(), m_shapeSet.getPath()); + DiscreteShape discreteShape(shape, m_dataStore.getRoot(), m_prefixPath); discreteShape.setVertexWeldThreshold(m_vertexWeldThreshold); discreteShape.setRefinementType(m_refinementType); if(percentError > 0) @@ -135,29 +245,62 @@ void Shaper::loadShapeInternal(const klee::Shape& shape, revolvedVolume = discreteShape.getRevolvedVolume(); } +bool Shaper::verifyInputMesh(std::string& whyBad) const +{ + bool rval = true; + +#if defined(AXOM_USE_CONDUIT) + if(m_bpGrp != nullptr) + { + conduit::Node info; + rval = conduit::blueprint::mesh::verify(m_bpNodeInt, info); + if(rval) + { + std::string topoType = + m_bpNodeInt.fetch("topologies")[m_bpTopo]["type"].as_string(); + rval = topoType == "unstructured"; + info[0].set_string("Topology is not unstructured."); + } + if(rval) + { + std::string elemShape = + m_bpNodeInt.fetch("topologies")[m_bpTopo]["elements"]["shape"].as_string(); + rval = elemShape == "hex"; + info[0].set_string("Topology elements are not hex."); + } + whyBad = info.to_summary_string(); + } +#endif + +#if defined(AXOM_USE_MFEM) + if(m_dc != nullptr) + { + // No specific requirements for MFEM mesh. + } +#endif + + return rval; +} + // ---------------------------------------------------------------------------- int Shaper::getRank() const { -#if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) - if(auto* pmesh = static_cast(m_dc->GetMesh())) - { - return pmesh->GetMyRank(); - } +#if defined(AXOM_USE_MPI) + int rank = -1; + MPI_Comm_rank(m_comm, &rank); + return rank; #endif return 0; } double Shaper::allReduceSum(double val) const { - if(m_dc != nullptr && m_dc->GetNumProcs() > 1) - { #if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) - double global; - MPI_Allreduce(&val, &global, 1, MPI_DOUBLE, MPI_SUM, m_comm); - return global; + double global; + MPI_Allreduce(&val, &global, 1, MPI_DOUBLE, MPI_SUM, m_comm); + return global; #endif - } return val; } diff --git a/src/axom/quest/Shaper.hpp b/src/axom/quest/Shaper.hpp index 93b2b54d3e..88ed87cb79 100644 --- a/src/axom/quest/Shaper.hpp +++ b/src/axom/quest/Shaper.hpp @@ -16,14 +16,23 @@ #ifndef AXOM_USE_KLEE #error Shaping functionality requires Axom to be configured with the Klee component #endif -#ifndef AXOM_USE_MFEM - #error Shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option + +#if !defined(AXOM_USE_MFEM) && !defined(AXOM_USE_CONDUIT) + #error Shaping functionality requires Axom to be configured with Conduit or MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option #endif #include "axom/sidre.hpp" #include "axom/klee.hpp" #include "axom/mint.hpp" #include "axom/quest/DiscreteShape.hpp" +#include "axom/core/execution/runtime_policy.hpp" + +#if defined(AXOM_USE_MFEM) + #include "mfem.hpp" +#endif +#if defined(AXOM_USE_CONDUIT) + #include "conduit_node.hpp" +#endif #include "axom/quest/interface/internal/mpicomm_wrapper.hpp" @@ -33,13 +42,51 @@ namespace quest { /** * Abstract base class for shaping material volume fractions + * + * Shaper requires Axom to be configured with Conduit or MFEM + * or both. */ class Shaper { public: - Shaper(const klee::ShapeSet& shapeSet, sidre::MFEMSidreDataCollection* dc); + using RuntimePolicy = axom::runtime_policy::Policy; - virtual ~Shaper() = default; +#if defined(AXOM_USE_MFEM) + /*! + @brief Construct Shaper to operate on an MFEM mesh. + */ + Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + sidre::MFEMSidreDataCollection* dc); +#endif + + /*! + @brief Construct Shaper to operate on a blueprint-formatted mesh + stored in a sidre Group. + */ + Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + sidre::Group* bpMesh, + const std::string& topo = ""); + + /*! + @brief Construct Shaper to operate on a blueprint-formatted mesh + stored in a conduit Node. + + Because \c conduit::Node doesn't support application-specified + allocator id for (only) arrays, the incoming \c bpNode must have + all arrays pre-allocated in a space accessible by the runtime + policy. Any needed-but-missing space would lead to an exception. + */ + Shaper(RuntimePolicy execPolicy, + int allocatorId, + const klee::ShapeSet& shapeSet, + conduit::Node* bpNode, + const std::string& topo = ""); + + virtual ~Shaper(); public: // Some default values. @@ -51,6 +98,9 @@ class Shaper /// Refinement type. using RefinementType = DiscreteShape::RefinementType; + //! @brief Verify the input mesh is okay for this class to work with. + bool verifyInputMesh(std::string& whyBad) const; + //@{ //! @name Functions to get and set shaping parameters @@ -62,10 +112,22 @@ class Shaper //@} + /*! + @brief Set path of shape input file. + + The path is used to resolve relative paths that may have been + specified in the file. + */ + void setFilePath(const std::string& filePath); + + mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh.get(); } + bool isVerbose() const { return m_verboseOutput; } +#ifdef AXOM_USE_MFEM sidre::MFEMSidreDataCollection* getDC() { return m_dc; } - mint::Mesh* getSurfaceMesh() const { return m_surfaceMesh.get(); } + const sidre::MFEMSidreDataCollection* getDC() const { return m_dc; } +#endif /*! * \brief Predicate to determine if the specified format is valid @@ -152,10 +214,35 @@ class Shaper int getRank() const; protected: + RuntimePolicy m_execPolicy; + int m_allocatorId; + sidre::DataStore m_dataStore; const klee::ShapeSet& m_shapeSet; - sidre::MFEMSidreDataCollection* m_dc; + + //! \brief Prefix path for shape file names with relative path. + std::string m_prefixPath; + +#if defined(AXOM_USE_MFEM) + // For mesh represented as MFEMSidreDataCollection + sidre::MFEMSidreDataCollection* m_dc {nullptr}; +#endif + +#if defined(AXOM_USE_CONDUIT) + // For mesh represented in Conduit or sidre + sidre::DataStore m_ds; + //! @brief Version of the mesh for computations. + axom::sidre::Group* m_bpGrp; + const std::string m_bpTopo; + //! @brief Mesh in an external Node, when provided as a Node. + conduit::Node* m_bpNodeExt; + //! @brief Initial copy of mesh in an internal Node storage. + conduit::Node m_bpNodeInt; +#endif + + //! @brief Number of cells in computational mesh (m_dc or m_bpGrp). + axom::IndexType m_cellCount; std::shared_ptr m_surfaceMesh; diff --git a/src/axom/quest/detail/shaping/shaping_helpers.cpp b/src/axom/quest/detail/shaping/shaping_helpers.cpp index 0cac0b8958..735e91787c 100644 --- a/src/axom/quest/detail/shaping/shaping_helpers.cpp +++ b/src/axom/quest/detail/shaping/shaping_helpers.cpp @@ -11,10 +11,6 @@ #include "axom/fmt.hpp" -#ifndef AXOM_USE_MFEM - #error Shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option -#endif - namespace axom { namespace quest diff --git a/src/axom/quest/detail/shaping/shaping_helpers.hpp b/src/axom/quest/detail/shaping/shaping_helpers.hpp index 4e223e3ab4..8e9550315a 100644 --- a/src/axom/quest/detail/shaping/shaping_helpers.hpp +++ b/src/axom/quest/detail/shaping/shaping_helpers.hpp @@ -14,8 +14,8 @@ #include "axom/config.hpp" -#ifndef AXOM_USE_MFEM - #error Shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option +#if !defined(AXOM_USE_MFEM) + #error Sampling-shaping functionality requires Axom to be configured with MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option #endif #include "mfem.hpp" diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 747f74936e..22b0a52dd6 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -265,14 +265,19 @@ if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI 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) +if((CONDUIT_FOUND OR + (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 UMPIRE_FOUND AND RAJA_FOUND) + if(MFEM_FOUND) + set(optional_dependency, "mfem") + endif() 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 + DEPENDS_ON ${quest_example_depends} ${optional_dependency} FOLDER axom/quest/examples ) @@ -288,8 +293,7 @@ if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) endif() - set(_testshapes "tetmesh" "tet" "hex" "sphere" "cyl" "cone" "vor" "all" "plane") - + set(_testshapes "tetmesh" "tet" "hex" "sphere" "cyl" "cone" "sor" "all" "plane") foreach(_policy ${_policies}) foreach(_testshape ${_testshapes}) set(_testname "quest_shape_in_memory_${_policy}_${_testshape}") @@ -301,10 +305,34 @@ if(AXOM_ENABLE_MPI AND MFEM_FOUND AND MFEM_USE_MPI --refinements 2 --scale 0.75 0.75 0.75 --dir 0.2 0.4 0.8 + --meshType bpSidre inline_mesh --min -2 -2 -2 --max 2 2 2 --res 30 30 30 NUM_MPI_TASKS ${_nranks}) endforeach() endforeach() + + # Test Conduit and MFEM mesh types with a subset of shapes (so it doesn't take too long) + set(_testshapes "tetmesh" "sphere" "cyl") + set(_testMeshTypes "bpConduit") + blt_list_append(TO _testMeshTypes ELEMENTS "mfem" IF MFEM_FOUND) + foreach(_policy ${_policies}) + foreach(_testMeshType ${_testMeshTypes}) + foreach(_testshape ${_testshapes}) + set(_testname "quest_shape_in_memory_${_policy}_${_testMeshType}_${_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 + --meshType ${_testMeshType} + inline_mesh --min -2 -2 -2 --max 2 2 2 --res 8 8 8 + NUM_MPI_TASKS ${_nranks}) + endforeach() + endforeach() + endforeach() endif() endif() @@ -372,6 +400,7 @@ if(AXOM_ENABLE_MPI AND AXOM_ENABLE_SIDRE AND HDF5_FOUND) endforeach() endforeach() + unset(optional_dependency) unset(_nranks) unset(_policies) unset(_test) diff --git a/src/axom/quest/examples/quest_marching_cubes_example.cpp b/src/axom/quest/examples/quest_marching_cubes_example.cpp index 3c1223533e..0143d4300f 100644 --- a/src/axom/quest/examples/quest_marching_cubes_example.cpp +++ b/src/axom/quest/examples/quest_marching_cubes_example.cpp @@ -30,7 +30,9 @@ #include "axom/core/MDMapping.hpp" #include "axom/quest/MarchingCubes.hpp" #include "axom/quest/MeshViewUtil.hpp" -#include "axom/sidre.hpp" +#if defined(AXOM_USE_SIDRE) + #include "axom/sidre.hpp" +#endif #include "axom/core/Types.hpp" #include "axom/core/numerics/floating_point_limits.hpp" @@ -56,7 +58,9 @@ namespace quest = axom::quest; namespace slic = axom::slic; +#if defined(AXOM_USE_SIDRE) namespace sidre = axom::sidre; +#endif namespace primal = axom::primal; namespace mint = axom::mint; namespace numerics = axom::numerics; @@ -667,6 +671,7 @@ void saveMesh(const conduit::Node& mesh, const std::string& filename) #endif } +#if defined(AXOM_USE_SIDRE) /// Write blueprint mesh to disk void saveMesh(const sidre::Group& mesh, const std::string& filename) { @@ -676,11 +681,11 @@ void saveMesh(const sidre::Group& mesh, const std::string& filename) mesh.createNativeLayout(tmpMesh); { conduit::Node info; -#ifdef AXOM_USE_MPI + #ifdef AXOM_USE_MPI if(!conduit::blueprint::mpi::verify("mesh", tmpMesh, info, MPI_COMM_WORLD)) -#else + #else if(!conduit::blueprint::verify("mesh", tmpMesh, info)) -#endif + #endif { SLIC_INFO("Invalid blueprint for mesh: \n" << info.to_yaml()); slic::flushStreams(); @@ -690,6 +695,7 @@ void saveMesh(const sidre::Group& mesh, const std::string& filename) } saveMesh(tmpMesh, filename); } +#endif template T product(const axom::StackArray& a) @@ -914,10 +920,10 @@ struct ContourTestBase AXOM_ANNOTATE_BEGIN("error checking"); AXOM_ANNOTATE_BEGIN("convert to mint mesh"); - std::string sidreGroupName = "contour_mesh"; - sidre::DataStore objectDS; #ifdef AXOM_MINT_USE_SIDRE + std::string sidreGroupName = "contour_mesh"; + sidre::DataStore objectDS; auto* meshGroup = objectDS.getRoot()->createGroup(sidreGroupName); axom::mint::UnstructuredMesh contourMesh( DIM, @@ -959,7 +965,7 @@ struct ContourTestBase checkCellsContainingContour(computationalMesh, contourMesh); } -#ifdef AXOM_MINT_USE_SIDRE +#if defined(AXOM_MINT_USE_SIDRE) if(contourMesh.hasSidreGroup()) { assert(contourMesh.getSidreGroup() == meshGroup); @@ -968,10 +974,10 @@ struct ContourTestBase saveMesh(*contourMesh.getSidreGroup(), outputName); SLIC_INFO(axom::fmt::format("Wrote contour mesh to {}", outputName)); } + objectDS.getRoot()->destroyGroupAndData(sidreGroupName); #endif - AXOM_ANNOTATE_END("error checking"); - objectDS.getRoot()->destroyGroupAndData(sidreGroupName); + AXOM_ANNOTATE_END("error checking"); return localErrCount; } @@ -1670,6 +1676,7 @@ int allocatorIdToTest(axom::runtime_policy::Policy policy) #endif axom::INVALID_ALLOCATOR_ID; #else + AXOM_UNUSED_VAR(policy); int allocatorID = axom::getDefaultAllocatorID(); #endif return allocatorID; diff --git a/src/axom/quest/examples/quest_shape_in_memory.cpp b/src/axom/quest/examples/quest_shape_in_memory.cpp index efa026ec14..4fda13f9f7 100644 --- a/src/axom/quest/examples/quest_shape_in_memory.cpp +++ b/src/axom/quest/examples/quest_shape_in_memory.cpp @@ -23,24 +23,25 @@ #include "axom/fmt.hpp" #include "axom/CLI11.hpp" +#if !defined(AXOM_USE_CONDUIT) + #error Shaping functionality requires Axom to be configured with Conduit +#endif + #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 +#if defined(AXOM_USE_MFEM) + #include "mfem.hpp" #endif -#include "mfem.hpp" - #ifdef AXOM_USE_MPI #include "mpi.h" #endif // RAJA -#ifdef AXOM_USE_RAJA - #include "RAJA/RAJA.hpp" +#if !defined(AXOM_USE_RAJA) + #error quest_shape_in_memory example and IntersectionShaper require RAJA #endif +#include "RAJA/RAJA.hpp" // C/C++ includes #include @@ -73,9 +74,6 @@ struct Input // 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}; @@ -87,6 +85,10 @@ struct Input SLIC_ASSERT(boxMaxs.size() == d); return int(d); } + int getBoxCellCount() const + { + return boxResolution[0] * boxResolution[1] * boxResolution[2]; + } // The shape to run. std::string testShape {"tetmesh"}; @@ -95,7 +97,7 @@ struct Input "sphere", "cyl", "cone", - "vor", + "sor", "tet", "hex", "plane", @@ -110,6 +112,16 @@ struct Input std::string backgroundMaterial; + // clang-format off + enum class MeshType { bpSidre = 0, bpConduit = 1, mfem = 2 }; + const std::map meshTypeChoices + { {"bpSidre", MeshType::bpSidre} , {"bpConduit", MeshType::bpConduit}, {"mfem", MeshType::mfem} }; + // clang-format on + MeshType meshType {MeshType::bpSidre}; + bool useMfem() { return meshType == MeshType::mfem; } + bool useBlueprintSidre() { return meshType == MeshType::bpSidre; } + bool useBlueprintConduit() { return meshType == MeshType::bpConduit; } + private: bool m_verboseOutput {false}; @@ -126,6 +138,7 @@ struct Input return volume; } +#if defined(AXOM_USE_MFEM) /// Generate an mfem Cartesian mesh, scaled to the bounding box range mfem::Mesh* createBoxMesh() { @@ -169,7 +182,7 @@ struct Input } // Handle conversion to parallel mfem mesh -#if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) + #if defined(AXOM_USE_MPI) && defined(MFEM_USE_MPI) { int* partitioning = nullptr; int part_method = 0; @@ -179,7 +192,7 @@ struct Input delete mesh; mesh = parallelMesh; } -#endif + #endif return mesh; } @@ -196,6 +209,7 @@ struct Input return dc; } +#endif void parse(int argc, char** argv, axom::CLI::App& app) { @@ -206,6 +220,11 @@ struct Input ->description("Enable/disable verbose output") ->capture_default_str(); + app.add_option("--meshType", meshType) + ->description("Type of computational mesh to shape on") + ->capture_default_str() + ->transform(axom::CLI::CheckedTransformer(meshTypeChoices)); + app.add_option("-t,--weld-threshold", weldThresh) ->description("Threshold for welding") ->check(axom::CLI::NonNegativeNumber) @@ -234,7 +253,7 @@ struct Input #endif app.add_option("--center", center) - ->description("Center of sphere or base of cone/cyl/VOR (x,y[,z]) shape") + ->description("Center of sphere or base of cone/cyl/SOR (x,y[,z]) shape") ->expected(2, 3); app.add_option("--radius", radius) @@ -242,12 +261,12 @@ struct Input ->check(axom::CLI::PositiveNumber); app.add_option("--length", length) - ->description("Length of cone/cyl/VOR shape") + ->description("Length of cone/cyl/SOR shape, avg length of hex.") ->check(axom::CLI::PositiveNumber); app.add_option("--dir", direction) ->description( - "Direction of axis of rotation (cone/cyl/VOR (x,y[,z])), or rotated " + "Direction of axis of rotation (cone/cyl/SOR (x,y[,z])), or rotated " "x-axis (hex, tet, tetmesh, and sphere), or positive normal direction " "(plane).") ->expected(2, 3); @@ -379,6 +398,7 @@ void addRotateOperator(axom::klee::CompositeOperator& compositeOp) } } +#if defined(AXOM_USE_MFEM) /** * \brief Print some info about the mesh * @@ -390,14 +410,14 @@ void printMfemMeshInfo(mfem::Mesh* mesh, const std::string& prefixMessage = "") namespace primal = axom::primal; int myRank = 0; -#ifdef AXOM_USE_MPI + #ifdef AXOM_USE_MPI MPI_Comm_rank(MPI_COMM_WORLD, &myRank); -#endif + #endif int numElements = mesh->GetNE(); mfem::Vector mins, maxs; -#ifdef MFEM_USE_MPI + #ifdef MFEM_USE_MPI auto* parallelMesh = dynamic_cast(mesh); if(parallelMesh != nullptr) { @@ -406,7 +426,7 @@ void printMfemMeshInfo(mfem::Mesh* mesh, const std::string& prefixMessage = "") myRank = parallelMesh->GetMyRank(); } else -#endif + #endif { mesh->GetBoundingBox(mins, maxs); } @@ -438,6 +458,42 @@ void printMfemMeshInfo(mfem::Mesh* mesh, const std::string& prefixMessage = "") slic::flushStreams(); } +#endif + +const std::string topoName = "mesh"; +const std::string coordsetName = "coords"; +int cellCount = -1; + +// Computational mesh in different forms, initialized in main +#if defined(AXOM_USE_MFEM) +std::shared_ptr shapingDC; +#endif +axom::sidre::Group* compMeshGrp = nullptr; +std::shared_ptr compMeshNode; + +axom::sidre::Group* createBoxMesh(axom::sidre::Group* meshGrp) +{ + using BBox3D = primal::BoundingBox; + using Pt3D = primal::Point; + auto res = primal::NumericArray(params.boxResolution.data()); + auto bbox = BBox3D(Pt3D(params.boxMins.data()), Pt3D(params.boxMaxs.data())); + axom::quest::util::make_unstructured_blueprint_box_mesh(meshGrp, + bbox, + res, + topoName, + coordsetName, + params.policy); +#if defined(AXOM_DEBUG) + conduit::Node meshNode, info; + meshGrp->createNativeLayout(meshNode); + SLIC_ASSERT(conduit::blueprint::mesh::verify(meshNode, info)); +#endif + + // State group is optional to blueprint, and we don't use it, but mint checks for it. + meshGrp->createGroup("state"); + + return meshGrp; +} /// \brief Utility function to initialize the logger void initializeLogger() @@ -505,7 +561,6 @@ axom::klee::ShapeSet create2DShapeSet(sidre::DataStore& ds) axom::IndexType conn[3] = {0, 1, 2}; triangleMesh.appendCell(conn); - meshGroup->print(); SLIC_ASSERT(axom::mint::blueprint::isValidRootGroup(meshGroup)); axom::klee::TransformableGeometryProperties prop { @@ -610,9 +665,9 @@ axom::klee::Shape createShape_TetMesh(sidre::DataStore& ds) return tetShape; } -axom::klee::Geometry createGeometry_Vor( - axom::primal::Point& vorBase, - axom::primal::Vector& vorDirection, +axom::klee::Geometry createGeometry_Sor( + axom::primal::Point& sorBase, + axom::primal::Vector& sorDirection, axom::Array& discreteFunction, std::shared_ptr& compositeOp) { @@ -623,20 +678,20 @@ axom::klee::Geometry createGeometry_Vor( SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); const axom::IndexType levelOfRefinement = params.refinementLevel; - axom::klee::Geometry vorGeometry(prop, + axom::klee::Geometry sorGeometry(prop, discreteFunction, - vorBase, - vorDirection, + sorBase, + sorDirection, levelOfRefinement, compositeOp); - return vorGeometry; + return sorGeometry; } -axom::klee::Shape createShape_Vor() +axom::klee::Shape createShape_Sor() { - Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; - axom::primal::Vector vorDirection = params.direction.empty() + axom::primal::Vector sorDirection = params.direction.empty() ? primal::Vector3D {0.1, 0.2, 0.4} : primal::Vector3D {params.direction.data()}; const int numIntervals = 5; @@ -665,19 +720,19 @@ axom::klee::Shape createShape_Vor() addScaleOperator(*compositeOp); addTranslateOperator(*compositeOp, -1, -1, 1); - axom::klee::Geometry vorGeometry = - createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + axom::klee::Geometry sorGeometry = + createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape vorShape("vor", "VOR", {}, {}, vorGeometry); + axom::klee::Shape sorShape("sor", "SOR", {}, {}, sorGeometry); - return vorShape; + return sorShape; } axom::klee::Shape createShape_Cylinder() { - Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; - axom::primal::Vector vorDirection = params.direction.empty() + axom::primal::Vector sorDirection = params.direction.empty() ? primal::Vector3D {0.1, 0.2, 0.4} : primal::Vector3D {params.direction.data()}; // discreteFunction are discrete z-r pairs describing the function @@ -694,19 +749,19 @@ axom::klee::Shape createShape_Cylinder() addScaleOperator(*compositeOp); addTranslateOperator(*compositeOp, 1, -1, 1); - axom::klee::Geometry vorGeometry = - createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + axom::klee::Geometry sorGeometry = + createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape vorShape("cyl", "CYL", {}, {}, vorGeometry); + axom::klee::Shape sorShape("cyl", "CYL", {}, {}, sorGeometry); - return vorShape; + return sorShape; } axom::klee::Shape createShape_Cone() { - Point3D vorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; - axom::primal::Vector vorDirection = params.direction.empty() + axom::primal::Vector sorDirection = params.direction.empty() ? primal::Vector3D {0.1, 0.2, 0.4} : primal::Vector3D {params.direction.data()}; // discreteFunction are discrete z-r pairs describing the function @@ -724,12 +779,12 @@ axom::klee::Shape createShape_Cone() addScaleOperator(*compositeOp); addTranslateOperator(*compositeOp, 1, 1, -1); - axom::klee::Geometry vorGeometry = - createGeometry_Vor(vorBase, vorDirection, discreteFunction, compositeOp); + axom::klee::Geometry sorGeometry = + createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape vorShape("cone", "CONE", {}, {}, vorGeometry); + axom::klee::Shape sorShape("cone", "CONE", {}, {}, sorGeometry); - return vorShape; + return sorShape; } axom::klee::Shape createShape_Tet() @@ -767,7 +822,7 @@ axom::klee::Shape createShape_Hex() SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); - const double md = params.length < 0 ? 0.6 : params.length; + const double md = params.length < 0 ? 0.6 : params.length / 2; const double lg = 1.2 * md; const double sm = 0.8 * md; const Point3D p {-lg, -md, -sm}; @@ -846,10 +901,10 @@ double volumeOfTetMesh( 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); + const axom::IndexType tetCount = tetMesh.getNumberOfCells(); + axom::Array tetVolumes(tetCount, tetCount); double meshVolume = 0.0; - for(axom::IndexType ic = 0; ic < cellCount; ++ic) + for(axom::IndexType ic = 0; ic < tetCount; ++ic) { const axom::IndexType* nodeIds = tetMesh.getCellNodeIDs(ic); TetType tet; @@ -865,6 +920,7 @@ double volumeOfTetMesh( return meshVolume; } +#if defined(AXOM_USE_MFEM) /*! @brief Return the element volumes as a sidre::View. @@ -884,7 +940,6 @@ axom::sidre::View* getElementVolumes( 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; @@ -957,7 +1012,157 @@ axom::sidre::View* getElementVolumes( return volSidreView; } +#endif +/*! + @brief Return the element volumes as a sidre::View containing + the volumes in an array. + + If it doesn't exist, allocate and compute it. + \post The volume data is in the blueprint field \c volFieldName. + + Most of this is lifted from IntersectionShaper::runShapeQueryImpl. +*/ +template +axom::sidre::View* getElementVolumes( + sidre::Group* meshGrp, + const std::string& volFieldName = std::string("elementVolumes")) +{ + using XS = axom::execution_space; + using HexahedronType = axom::primal::Hexahedron; + + axom::sidre::View* volSidreView = nullptr; + + const auto fieldPath = axom::fmt::format("fields/{}", volFieldName); + if(meshGrp->hasGroup(fieldPath)) + { + sidre::Group* fieldGrp = meshGrp->getGroup(fieldPath); + volSidreView = fieldGrp->getView("values"); + } + else + { + axom::mint::UnstructuredMesh mesh(meshGrp, + topoName); + + constexpr int NUM_VERTS_PER_HEX = 8; + constexpr int NUM_COMPS_PER_VERT = 3; + constexpr double ZERO_THRESHOLD = 1.e-10; + + /* + Get vertex coordinates. We use UnstructuredMesh for this, + so get it on host first then transfer to device if needed. + */ + auto* connData = meshGrp->getGroup("topologies") + ->getGroup(topoName) + ->getGroup("elements") + ->getView("connectivity"); + SLIC_ASSERT(connData->getNode().dtype().id() == + axom::conduitDataIdOfAxomIndexType); + + conduit::Node coordNode; + meshGrp->getGroup("coordsets") + ->getGroup(coordsetName) + ->createNativeLayout(coordNode); + const conduit::Node& coordValues = coordNode.fetch_existing("values"); + axom::IndexType vertexCount = coordValues["x"].dtype().number_of_elements(); + bool isInterleaved = conduit::blueprint::mcarray::is_interleaved(coordValues); + int stride = isInterleaved ? NUM_COMPS_PER_VERT : 1; + axom::StackArray, 3> coordArrays { + axom::ArrayView(coordValues["x"].as_double_ptr(), + {vertexCount}, + stride), + axom::ArrayView(coordValues["y"].as_double_ptr(), + {vertexCount}, + stride), + axom::ArrayView(coordValues["z"].as_double_ptr(), + {vertexCount}, + stride)}; + + const axom::IndexType* connPtr = connData->getArray(); + SLIC_ASSERT(connPtr != nullptr); + axom::ArrayView conn(connPtr, + cellCount, + NUM_VERTS_PER_HEX); + axom::Array vertCoords(cellCount * NUM_VERTS_PER_HEX, + cellCount * NUM_VERTS_PER_HEX, + XS::allocatorID()); + auto vertCoordsView = vertCoords.view(); + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellIdx) { + // Get the indices of this element's vertices + auto verts = conn[cellIdx]; + + // 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] = coordArrays[k][verts[j]]; + // vertCoordsView[vertIdx][k] = mesh.getNodeCoordinate(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, + meshGrp->getDefaultAllocatorID()); + 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. + axom::sidre::Group* fieldGrp = meshGrp->createGroup(fieldPath); + fieldGrp->createViewString("topology", topoName); + fieldGrp->createViewString("association", "element"); + fieldGrp->createViewString("volume_dependent", "true"); + volSidreView = + fieldGrp->createViewAndAllocate("values", + axom::sidre::detail::SidreTT::id, + cellCount); + axom::IndexType shape2d[] = {cellCount, 1}; + volSidreView->reshapeArray(2, shape2d); + axom::ArrayView volView(volSidreView->getData(), + volSidreView->getNumElements()); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellIdx) { + volView[cellIdx] = hexesView[cellIdx].volume(); + }); + } + + return volSidreView; +} + +#if defined(AXOM_USE_MFEM) /*! @brief Return global sum of volume of the given material. */ @@ -965,9 +1170,6 @@ 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(), @@ -977,7 +1179,60 @@ double sumMaterialVolumes(sidre::MFEMSidreDataCollection* dc, const std::string materialFieldName = axom::fmt::format("vol_frac_{}", material); mfem::GridFunction* volFracGf = dc->GetField(materialFieldName); - axom::quest::GridFunctionView volFracView(volFracGf); + axom::ArrayView volFracGfArrayView(volFracGf->GetData(), + volFracGf->Size()); + axom::Array volFracGfArray( + volFracGfArrayView, + axom::execution_space::allocatorID()); + auto volFracView = volFracGfArray.view(); + + 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; +} +#endif + +template +double sumMaterialVolumesImpl(sidre::Group* meshGrp, const std::string& material) +{ + conduit::Node meshNode; + meshGrp->createNativeLayout(meshNode); +#if defined(AXOM_DEBUG) + // Conduit can verify Blueprint mesh, but only if data is on host. + if(axom::execution_space::usesAllocId( + meshGrp->getDefaultAllocatorID())) + { + conduit::Node info; + conduit::blueprint::mesh::verify(meshNode, info); + SLIC_ASSERT(conduit::blueprint::mesh::verify(meshNode, info)); + } +#endif + std::string topoPath = "topologies/" + topoName; + conduit::Node& topoNode = meshNode.fetch_existing(topoPath); + const int cellCount = conduit::blueprint::mesh::topology::length(topoNode); + + // Get cell volumes from meshGrp. + const std::string volsName = "vol_" + material; + axom::sidre::View* elementVols = + getElementVolumes(meshGrp, volsName); + axom::ArrayView elementVolsView(elementVols->getData(), + elementVols->getNumElements()); + + // Get material volume fractions + const auto vfFieldName = "vol_frac_" + material; + const auto vfFieldValuesPath = "fields/" + vfFieldName + "/values"; + axom::sidre::View* volFrac = meshGrp->getView(vfFieldValuesPath); + axom::ArrayView volFracView(volFrac->getArray(), cellCount); using ReducePolicy = typename axom::execution_space::reduce_policy; RAJA::ReduceSum localVol(0); @@ -994,6 +1249,152 @@ double sumMaterialVolumes(sidre::MFEMSidreDataCollection* dc, return globalVol; } +double sumMaterialVolumes(sidre::Group* meshGrp, const std::string& material) +{ + double rval = 0.0; + if(params.policy == RuntimePolicy::seq) + { + rval = sumMaterialVolumesImpl(meshGrp, material); + } +#if defined(AXOM_USE_OPENMP) + if(params.policy == RuntimePolicy::omp) + { + rval = sumMaterialVolumesImpl(meshGrp, material); + } +#endif +#if defined(AXOM_USE_CUDA) + if(params.policy == RuntimePolicy::cuda) + { + rval = sumMaterialVolumesImpl>(meshGrp, material); + } +#endif +#if defined(AXOM_USE_HIP) + if(params.policy == RuntimePolicy::hip) + { + rval = sumMaterialVolumesImpl>(meshGrp, material); + } +#endif + return rval; +} + +/// Write blueprint mesh to disk +void saveMesh(const conduit::Node& mesh, const std::string& filename) +{ + AXOM_ANNOTATE_SCOPE("save mesh (conduit)"); + +#ifdef AXOM_USE_MPI + conduit::relay::mpi::io::blueprint::save_mesh(mesh, + filename, + "hdf5", + MPI_COMM_WORLD); +#else + conduit::relay::io::blueprint::save_mesh(mesh, filename, "hdf5"); +#endif +} + +/// Write blueprint mesh to disk +void saveMesh(const sidre::Group& mesh, const std::string& filename) +{ + AXOM_ANNOTATE_SCOPE("save mesh (sidre)"); + + axom::sidre::DataStore ds; + const sidre::Group* meshOnHost = &mesh; + if(mesh.getDefaultAllocatorID() != + axom::execution_space::allocatorID()) + { + meshOnHost = ds.getRoot()->deepCopyGroup( + &mesh, + axom::execution_space::allocatorID()); + } + conduit::Node tmpMesh; + meshOnHost->createNativeLayout(tmpMesh); + { + conduit::Node info; +#ifdef AXOM_USE_MPI + if(!conduit::blueprint::mpi::verify("mesh", tmpMesh, info, MPI_COMM_WORLD)) +#else + if(!conduit::blueprint::verify("mesh", tmpMesh, info)) +#endif + { + SLIC_INFO("Invalid blueprint for mesh: \n" << info.to_yaml()); + slic::flushStreams(); + assert(false); + } + // info.print(); + } + saveMesh(tmpMesh, filename); +} + +axom::ArrayView getFieldAsArrayView(const std::string& fieldName) +{ + axom::ArrayView fieldDataView; + if(params.useBlueprintSidre()) + { + std::string valuesPath = "fields/" + fieldName + "/values"; + axom::sidre::View* fieldValues = compMeshGrp->getView(valuesPath); + double* fieldData = fieldValues->getArray(); + fieldDataView = + axom::ArrayView(fieldData, fieldValues->getNumElements()); + } + if(params.useBlueprintConduit()) + { + std::string valuesPath = "fields/" + fieldName + "/values"; + conduit::Node& fieldValues = compMeshNode->fetch_existing(valuesPath); + double* fieldData = fieldValues.as_double_ptr(); + fieldDataView = + axom::ArrayView(fieldData, + fieldValues.dtype().number_of_elements()); + } +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + // auto* mfemMesh = shapingDC->GetMesh(); + mfem::GridFunction* gridFunc = shapingDC->GetField(fieldName); + fieldDataView = + axom::ArrayView(gridFunc->GetData(), gridFunc->Size()); + } +#endif + return fieldDataView; +} + +//!@brief Fill a sidre array View with a value. +// No error checking. +template +void fillSidreViewData(axom::sidre::View* view, const T& value) +{ + double* valuesPtr = view->getData(); + switch(params.policy) + { +#if defined(AXOM_USE_CUDA) + case RuntimePolicy::cuda: + axom::for_all>( + view->getNumElements(), + AXOM_LAMBDA(axom::IndexType i) { valuesPtr[i] = value; }); + break; +#endif +#if defined(AXOM_USE_HIP) + case RuntimePolicy::hip: + axom::for_all>( + view->getNumElements(), + AXOM_LAMBDA(axom::IndexType i) { valuesPtr[i] = value; }); + break; +#endif +#if defined(AXOM_USE_OMP) + case RuntimePolicy::omp: + axom::for_all( + view->getNumElements(), + AXOM_LAMBDA(axom::IndexType i) { valuesPtr[i] = value; }); + break; +#endif + case RuntimePolicy::seq: + default: + axom::for_all( + view->getNumElements(), + AXOM_LAMBDA(axom::IndexType i) { valuesPtr[i] = value; }); + break; + } +} + //------------------------------------------------------------------------------ int main(int argc, char** argv) { @@ -1030,6 +1431,8 @@ int main(int argc, char** argv) axom::utilities::raii::AnnotationsWrapper annotations_raii_wrapper( params.annotationMode); + const int allocatorId = axom::policyToDefaultAllocatorID(params.policy); + AXOM_ANNOTATE_BEGIN("quest example for shaping primals"); AXOM_ANNOTATE_BEGIN("init"); @@ -1041,58 +1444,50 @@ int main(int argc, char** argv) //--------------------------------------------------------------------------- axom::klee::ShapeSet shapeSet; SLIC_ERROR_IF(params.getBoxDim() != 3, "This example is only in 3D."); - switch(params.getBoxDim()) + if(params.testShape == "tetmesh") { - 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; + 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 == "sor") + { + shapeSet = createShapeSet(createShape_Sor()); + } + 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_Sor()); + shapesVec.push_back(createShape_Cylinder()); + shapesVec.push_back(createShape_Cone()); + shapeSet.setShapes(shapesVec); + shapeSet.setDimensions(axom::klee::Dimensions::Three); } // Save the discrete shapes for viz and testing. @@ -1121,59 +1516,125 @@ int main(int argc, char** argv) } } - const klee::Dimensions shapeDim = shapeSet.getDimensions(); +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + 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; + shapingDC = std::make_shared("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(), "Loaded MFEM mesh"); - // 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"); + cellCount = shapingMesh->GetNE(); + } +#else + SLIC_ERROR_IF(params.useMfem(), + "Cannot use MFEM mesh due to Axom configuration. Please use " + "Blueprint mesh or configure with MFEM and " + "-DAXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION."); #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); + if(params.useBlueprintSidre() || params.useBlueprintConduit()) { - 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"); + compMeshGrp = ds.getRoot()->createGroup("compMesh"); + compMeshGrp->setDefaultAllocator(allocatorId); - const axom::IndexType cellCount = shapingMesh->GetNE(); + createBoxMesh(compMeshGrp); - // 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; + if(params.useBlueprintConduit()) + { + // Intersection requires conduit mesh to have array data pre-allocated. + auto makeField = [&](const std::string& fieldName, double initValue) { + auto fieldGrp = compMeshGrp->createGroup("fields/" + fieldName); + axom::IndexType shape[] = {params.getBoxCellCount(), 1}; + fieldGrp->createViewString("association", "element"); + fieldGrp->createViewString("topology", topoName); + fieldGrp->createViewString("volume_dependent", "true"); + axom::sidre::View* valuesView = fieldGrp->createViewWithShapeAndAllocate( + "values", + axom::sidre::detail::SidreTT::id, + 2, + shape); + fillSidreViewData(valuesView, initValue); + }; + makeField("vol_frac_free", 1.0); + const auto& shapes = shapeSet.getShapes(); + for(const auto& shape : shapes) + { + makeField("vol_frac_" + shape.getMaterial(), + 0.0); // Used in volume fraction computation + makeField("shape_vol_frac_" + shape.getName(), + 0.0); // Used in applyReplacementRules + } + } + + /* + Shallow-copy compMeshGrp into compMeshNode, + so that any change in one is reflected in the other. + */ + compMeshNode = std::make_shared(); + compMeshGrp->createNativeLayout(*compMeshNode); + SLIC_INFO(axom::fmt::format("{:-^80}", "Generated Blueprint mesh")); + cellCount = params.getBoxCellCount(); + } //--------------------------------------------------------------------------- // 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!"); + std::shared_ptr shaper = nullptr; + if(params.useBlueprintSidre()) + { + shaper = std::make_shared(params.policy, + allocatorId, + shapeSet, + compMeshGrp); + } + if(params.useBlueprintConduit()) + { + shaper = std::make_shared(params.policy, + allocatorId, + shapeSet, + compMeshNode.get()); + } +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + shaper = std::make_shared(params.policy, + allocatorId, + shapeSet, + shapingDC.get()); + } +#endif + SLIC_ASSERT(shaper != nullptr); + // shaper->setExecPolicy(params.policy); // Set generic parameters for the base Shaper instance shaper->setVertexWeldThreshold(params.weldThresh); @@ -1186,22 +1647,23 @@ int main(int argc, char** argv) // 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"); +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + shaper->getDC()->AssociateMaterialSet("vol_frac", "material"); + } +#endif // 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); + shaper->setLevel(params.refinementLevel); + SLIC_INFO(axom::fmt::format( + "{:-^80}", + axom::fmt::format("Setting IntersectionShaper policy to '{}'", + axom::runtime_policy::policyToName(params.policy)))); - if(!params.backgroundMaterial.empty()) - { - intersectionShaper->setFreeMaterialName(params.backgroundMaterial); - } + if(!params.backgroundMaterial.empty()) + { + shaper->setFreeMaterialName(params.backgroundMaterial); } AXOM_ANNOTATE_END("setup shaping problem"); @@ -1227,7 +1689,7 @@ int main(int argc, char** argv) slic::flushStreams(); // Generate a spatial index over the shape - shaper->prepareShapeQuery(shapeDim, shape); + shaper->prepareShapeQuery(shapeSet.getDimensions(), shape); slic::flushStreams(); // Query the mesh against this shape @@ -1258,49 +1720,80 @@ int main(int argc, char** argv) // 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(params.useBlueprintSidre() || params.useBlueprintConduit()) { - if(startsWith(kv.first, "vol_frac_")) +#if defined(AXOM_DEBUG) + conduit::Node info; + SLIC_ASSERT(conduit::blueprint::mesh::verify(*compMeshNode, info)); +#endif + std::vector materialNames = shaper->getMaterialNames(); + for(const auto& materialName : materialNames) { - 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); - + // Compute and print volume of material. + const double volume = sumMaterialVolumes(compMeshGrp, materialName); SLIC_INFO(axom::fmt::format(axom::utilities::locale(), "Volume of material '{}' is {:.6Lf}", - mat_name, + materialName, volume)); } } +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + for(auto& kv : shaper->getDC()->GetFieldMap()) + { + 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)); + } + } + } +#endif AXOM_ANNOTATE_END("adjust"); int failCounts = 0; - auto matsetsGrp = shapingDC.GetBPGroup()->getGroup("matsets"); - auto materialGrp = matsetsGrp->getGroup("material"); - auto volFracGroups = materialGrp->getGroup("volume_fractions"); + std::vector allVfNames( + 1, + "free"); // All volume fraction names plus "free". + for(const auto& shape : shapeSet.getShapes()) + { + allVfNames.push_back(shape.getMaterial()); + } + + // For error checking, work on host. + using ExecSpace = typename axom::SEQ_EXEC; + using ReducePolicy = typename axom::execution_space::reduce_policy; //--------------------------------------------------------------------------- // Correctness test: volume fractions should be in [0,1]. //--------------------------------------------------------------------------- RAJA::ReduceSum rangeViolationCount(0); - for(axom::sidre::Group& materialGroup : volFracGroups->groups()) + for(const auto& vfName : allVfNames) { - axom::sidre::View* values = materialGroup.getView("value"); - double* volFracData = values->getArray(); - axom::ArrayView volFracDataView(volFracData, cellCount); + std::string fieldName = "vol_frac_" + vfName; + axom::ArrayView vfView = getFieldAsArrayView(fieldName); + axom::Array vfHostArray( + vfView, + axom::execution_space::allocatorID()); + vfView = vfHostArray.view(); axom::for_all( cellCount, AXOM_LAMBDA(axom::IndexType i) { - bool bad = volFracDataView[i] < 0.0 || volFracDataView[i] > 1.0; + bool bad = vfView[i] < 0.0 || vfView[i] > 1.0; rangeViolationCount += bad; }); } @@ -1319,20 +1812,23 @@ int main(int argc, char** argv) axom::Array volSums(cellCount); volSums.fill(0.0); axom::ArrayView volSumsView = volSums.view(); - for(axom::sidre::Group& materialGroup : volFracGroups->groups()) + RAJA::ReduceSum nonUnitSums(0); + for(const auto& vfName : allVfNames) { - axom::sidre::View* values = materialGroup.getView("value"); - double* volFracData = values->getArray(); - axom::ArrayView volFracDataView(volFracData, cellCount); + std::string fieldName = "vol_frac_" + vfName; + axom::ArrayView vfView = getFieldAsArrayView(fieldName); + axom::Array vfHostArray( + vfView, + axom::execution_space::allocatorID()); + vfView = vfHostArray.view(); axom::for_all( cellCount, - AXOM_LAMBDA(axom::IndexType i) { volSumsView[i] += volFracDataView[i]; }); + AXOM_LAMBDA(axom::IndexType i) { volSumsView[i] += vfView[i]; }); } - RAJA::ReduceSum nonUnitSums(0); axom::for_all( cellCount, AXOM_LAMBDA(axom::IndexType i) { - bool bad = !axom::utilities::isNearlyEqual(volSums[i], 0.0); + bool bad = !axom::utilities::isNearlyEqual(volSumsView[i], 1.0); nonUnitSums += bad; }); @@ -1364,13 +1860,27 @@ int main(int argc, char** argv) shapeMesh->getNumberOfCells()))); const std::string& materialName = shape.getMaterial(); - double shapeVol = - sumMaterialVolumes(&shapingDC, materialName); + double shapeVol = -1; + if(params.useBlueprintSidre() || params.useBlueprintConduit()) + { + shapeVol = sumMaterialVolumes(compMeshGrp, materialName); + } +#if defined(AXOM_USE_MFEM) + if(params.useMfem()) + { + shapeVol = + sumMaterialVolumes(shapingDC.get(), materialName); + } +#endif double correctShapeVol = params.testShape == "plane" ? params.boxMeshVolume() / 2 : shapeMeshVol; + SLIC_ASSERT(correctShapeVol > 0.0); // Indicates error in the test setup. double diff = shapeVol - correctShapeVol; - bool err = !axom::utilities::isNearlyEqual(shapeVol, correctShapeVol); + bool err = !axom::utilities::isNearlyEqualRelative(shapeVol, + correctShapeVol, + 1e-6, + 1e-8); failCounts += err; SLIC_INFO(axom::fmt::format( @@ -1391,16 +1901,26 @@ int main(int argc, char** argv) // 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)); - } + if(params.useBlueprintSidre()) + { + saveMesh(*compMeshGrp, fileName); + SLIC_INFO(axom::fmt::format("{:-^80}", "Wrote output mesh " + fileName)); + } +#if defined(AXOM_USE_MFEM) + else + { + shaper->getDC()->Save(fileName, sidre::Group::getDefaultIOProtocol()); + } #endif + } - delete shaper; +#if defined(AXOM_USE_MFEM) + shapingDC.reset(); +#endif + shaper.reset(); //--------------------------------------------------------------------------- // Cleanup and exit diff --git a/src/axom/quest/examples/shaping_driver.cpp b/src/axom/quest/examples/shaping_driver.cpp index 31abd26041..85ff4a7d23 100644 --- a/src/axom/quest/examples/shaping_driver.cpp +++ b/src/axom/quest/examples/shaping_driver.cpp @@ -20,10 +20,10 @@ #include "axom/fmt.hpp" #include "axom/CLI11.hpp" -// NOTE: The shaping driver requires Axom to be configured with mfem as well as +// NOTE: The shaping driver requires Axom to be configured with conduit or mfem and // 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 +#if !defined(AXOM_USE_MFEM) && !defined(AXOM_USE_CONDUIT) + #error Shaping functionality requires Axom to be configured with Conduit or MFEM and the AXOM_ENABLE_MFEM_SIDRE_DATACOLLECTION option #endif #include "mfem.hpp" @@ -566,7 +566,16 @@ int main(int argc, char** argv) shaper = new quest::SamplingShaper(params.shapeSet, &shapingDC); break; case ShapingMethod::Intersection: - shaper = new quest::IntersectionShaper(params.shapeSet, &shapingDC); +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) + shaper = new quest::IntersectionShaper( + params.policy, + axom::policyToDefaultAllocatorID(axom::runtime_policy::Policy::seq), + params.shapeSet, + &shapingDC); +#else + SLIC_ERROR( + "IntersectionShaper requires Axom to be configured with Umpire."); +#endif break; } SLIC_ASSERT_MSG(shaper != nullptr, "Invalid shaping method selected!"); @@ -604,17 +613,18 @@ int main(int argc, char** argv) } } +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) // Set specific parameters here for IntersectionShaper if(auto* intersectionShaper = dynamic_cast(shaper)) { intersectionShaper->setLevel(params.refinementLevel); - intersectionShaper->setExecPolicy(params.policy); if(!params.backgroundMaterial.empty()) { intersectionShaper->setFreeMaterialName(params.backgroundMaterial); } } +#endif //--------------------------------------------------------------------------- // Project initial volume fractions, if applicable diff --git a/src/axom/quest/interface/inout.cpp b/src/axom/quest/interface/inout.cpp index 07aa782e35..c3e0a992b1 100644 --- a/src/axom/quest/interface/inout.cpp +++ b/src/axom/quest/interface/inout.cpp @@ -97,7 +97,9 @@ struct InOutHelper // load the mesh int rc = QUEST_INOUT_FAILED; +#ifdef AXOM_USE_C2C double revolvedVolume = 0.; +#endif switch(DIM) { case 2: @@ -110,7 +112,6 @@ struct InOutHelper revolvedVolume, comm); #else - AXOM_UNUSED_VAR(revolvedVolume); SLIC_WARNING(fmt::format( "Cannot read contour file: C2C not enabled in this configuration.", file)); diff --git a/src/axom/quest/tests/quest_intersection_shaper.cpp b/src/axom/quest/tests/quest_intersection_shaper.cpp index 6d32f07c9f..f7d474d962 100644 --- a/src/axom/quest/tests/quest_intersection_shaper.cpp +++ b/src/axom/quest/tests/quest_intersection_shaper.cpp @@ -284,6 +284,7 @@ bool loadBaseline(const std::string &filename, conduit::Node &n) return loaded; } +#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE) void replacementRuleTest(const std::string &shapeFile, const std::string &policyName, RuntimePolicy policy, @@ -316,16 +317,19 @@ void replacementRuleTest(const std::string &shapeFile, // Need to do the pipeline of the shaping driver. SLIC_INFO(axom::fmt::format("Shaping materials...")); const int refinementLevel = 7; -#ifdef AXOM_USE_MPI + #ifdef AXOM_USE_MPI // This has to happen here because the shaper gets its communicator from it. // If we do it before the mfem mesh is added to the data collection then the // data collection communicator gets set to MPI_COMM_NULL, which is bad for // the C2C reader. dc.SetComm(MPI_COMM_WORLD); -#endif - quest::IntersectionShaper shaper(shapeSet, &dc); + #endif + quest::IntersectionShaper shaper( + policy, + axom::policyToDefaultAllocatorID(axom::runtime_policy::Policy::seq), + shapeSet, + &dc); shaper.setLevel(refinementLevel); - shaper.setExecPolicy(policy); // Borrowed from shaping_driver. const klee::Dimensions shapeDim = shapeSet.getDimensions(); @@ -360,10 +364,10 @@ void replacementRuleTest(const std::string &shapeFile, conduit::Node current; dcToConduit(dc, current); -#ifdef VISUALIZE_DATASETS + #ifdef VISUALIZE_DATASETS saveVisIt("", baselineName, dc); -#endif -#ifdef GENERATE_BASELINES + #endif + #ifdef GENERATE_BASELINES for(const auto &path : baselinePaths) { SLIC_INFO(axom::fmt::format("Saving baseline to {}", path)); @@ -371,7 +375,7 @@ void replacementRuleTest(const std::string &shapeFile, std::string filename(pjoin(path, baselineName)); saveBaseline(filename, current); } -#endif + #endif // TODO: I might want an auto compare for generating baselines so I know if I need a policy-specific baseline. @@ -452,18 +456,21 @@ void IntersectionWithErrorTolerances(const std::string &filebase, // Need to do the pipeline of the shaping driver. SLIC_INFO(axom::fmt::format("Shaping materials...")); -#ifdef AXOM_USE_MPI + #ifdef AXOM_USE_MPI // This has to happen here because the shaper gets its communicator from it. // If we do it before the mfem mesh is added to the data collection then the // data collection communicator gets set to MPI_COMM_NULL, which is bad for // the C2C reader. dc.SetComm(MPI_COMM_WORLD); -#endif - quest::IntersectionShaper shaper(shapeSet, &dc); + #endif + quest::IntersectionShaper shaper( + policy, + axom::policyToDefaultAllocatorID(axom::runtime_policy::Policy::seq), + shapeSet, + &dc); shaper.setLevel(refinementLevel); shaper.setPercentError(targetPercentError); shaper.setRefinementType(quest::DiscreteShape::RefinementDynamic); - shaper.setExecPolicy(policy); // Borrowed from shaping_driver (there should just be one shape) const klee::Dimensions shapeDim = shapeSet.getDimensions(); @@ -771,6 +778,7 @@ dimensions: 3 policy); } } +#endif //--------------------------------------------------------------------------- // Define testing functions for different modes. diff --git a/src/axom/quest/util/mesh_helpers.cpp b/src/axom/quest/util/mesh_helpers.cpp index 0bdac364b8..aea68cd79c 100644 --- a/src/axom/quest/util/mesh_helpers.cpp +++ b/src/axom/quest/util/mesh_helpers.cpp @@ -4,6 +4,15 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include "mesh_helpers.hpp" +#if defined(AXOM_USE_SIDRE) + #include +#endif +#include +#include "axom/mint/mesh/UnstructuredMesh.hpp" +#if defined(AXOM_USE_CONDUIT) + #include +#endif +#include namespace axom { @@ -83,6 +92,421 @@ mfem::Mesh* make_cartesian_mfem_mesh_3D(const primal::BoundingBox& bb #endif // AXOM_USE_MFEM +#if defined(AXOM_USE_SIDRE) + +axom::sidre::Group* make_structured_blueprint_box_mesh( + axom::sidre::Group* meshGrp, + const primal::BoundingBox& bbox, + const primal::NumericArray& res, + const std::string& topologyName, + const std::string& coordsetName, + axom::runtime_policy::Policy runtimePolicy) +{ + auto* topoGrp = meshGrp->createGroup("topologies")->createGroup(topologyName); + SLIC_ERROR_IF(topoGrp == nullptr, + "Cannot allocate topology '" + topologyName + + "' in blueprint mesh '" + meshGrp->getName() + + "'. It already exists."); + + auto* coordsetGrp = + meshGrp->createGroup("coordsets")->createGroup(coordsetName); + SLIC_ERROR_IF(coordsetGrp == nullptr, + "Cannot allocate coordset '" + coordsetName + + "' in blueprint mesh '" + meshGrp->getName() + + "'. It already exists."); + + topoGrp->createView("type")->setString("structured"); + topoGrp->createView("coordset")->setString(coordsetName); + auto* dimsGrp = topoGrp->createGroup("elements/dims"); + + constexpr int DIM = 3; + + auto ni = res[0]; + auto nj = res[1]; + auto nk = res[2]; + + const axom::StackArray vertsShape {res[0] + 1, + res[1] + 1, + res[2] + 1}; + const auto numVerts = vertsShape[0] * vertsShape[1] * vertsShape[2]; + + dimsGrp->createViewScalar("i", ni); + dimsGrp->createViewScalar("j", nj); + dimsGrp->createViewScalar("k", nk); + + coordsetGrp->createView("type")->setString("explicit"); + auto* valuesGrp = coordsetGrp->createGroup("values"); + auto* xVu = + valuesGrp->createViewAndAllocate("x", + axom::sidre::DataTypeId::FLOAT64_ID, + numVerts); + auto* yVu = + valuesGrp->createViewAndAllocate("y", + axom::sidre::DataTypeId::FLOAT64_ID, + numVerts); + auto* zVu = + valuesGrp->createViewAndAllocate("z", + axom::sidre::DataTypeId::FLOAT64_ID, + numVerts); + #ifdef AXOM_USE_UMPIRE + SLIC_ASSERT(axom::getAllocatorIDFromPointer(xVu->getVoidPtr()) == + meshGrp->getDefaultAllocatorID()); + SLIC_ASSERT(axom::getAllocatorIDFromPointer(yVu->getVoidPtr()) == + meshGrp->getDefaultAllocatorID()); + SLIC_ASSERT(axom::getAllocatorIDFromPointer(zVu->getVoidPtr()) == + meshGrp->getDefaultAllocatorID()); + #endif + + const axom::MDMapping vertMapping(vertsShape, + axom::ArrayStrideOrder::COLUMN); + axom::ArrayView xView(xVu->getData(), + vertsShape, + vertMapping.strides()); + axom::ArrayView yView(yVu->getData(), + vertsShape, + vertMapping.strides()); + axom::ArrayView zView(zVu->getData(), + vertsShape, + vertMapping.strides()); + + fill_cartesian_coords_3d(runtimePolicy, bbox, xView, yView, zView); + + #if defined(AXOM_DEBUG) && defined(AXOM_USE_CONDUIT) + { + conduit::Node info; + bool isValid = verifyBlueprintMesh(meshGrp, info); + SLIC_ASSERT_MSG(isValid, "Internal error: Generated mesh is invalid."); + } + #endif + + return meshGrp; +} + +axom::sidre::Group* make_unstructured_blueprint_box_mesh( + axom::sidre::Group* meshGrp, + const primal::BoundingBox& bbox, + const primal::NumericArray& res, + const std::string& topologyName, + const std::string& coordsetName, + axom::runtime_policy::Policy runtimePolicy) +{ + make_structured_blueprint_box_mesh(meshGrp, + bbox, + res, + topologyName, + coordsetName, + runtimePolicy); + convert_blueprint_structured_explicit_to_unstructured(meshGrp, + topologyName, + runtimePolicy); + return meshGrp; +} + +void convert_blueprint_structured_explicit_to_unstructured( + axom::sidre::Group* meshGrp, + const std::string& topoName, + axom::runtime_policy::Policy runtimePolicy) +{ + if(runtimePolicy == axom::runtime_policy::Policy::seq) + { + convert_blueprint_structured_explicit_to_unstructured_impl( + meshGrp, + topoName); + } + #if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + if(runtimePolicy == axom::runtime_policy::Policy::omp) + { + convert_blueprint_structured_explicit_to_unstructured_impl( + meshGrp, + topoName); + } + #endif + #if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + if(runtimePolicy == axom::runtime_policy::Policy::cuda) + { + convert_blueprint_structured_explicit_to_unstructured_impl>( + meshGrp, + topoName); + } + #endif + #if defined(AXOM_RUNTIME_POLICY_USE_HIP) + if(runtimePolicy == axom::runtime_policy::Policy::hip) + { + convert_blueprint_structured_explicit_to_unstructured_impl>( + meshGrp, + topoName); + } + #endif +} + +template +void convert_blueprint_structured_explicit_to_unstructured_impl( + axom::sidre::Group* meshGrp, + const std::string& topoName) +{ + AXOM_ANNOTATE_SCOPE("convert_to_unstructured"); + constexpr int DIM = 3; + + const std::string& coordsetName = + meshGrp->getView("topologies/" + topoName + "/coordset")->getString(); + + sidre::Group* coordsetGrp = + meshGrp->getGroup("coordsets")->getGroup(coordsetName); + SLIC_ASSERT(std::string(coordsetGrp->getView("type")->getString()) == + "explicit"); + + axom::sidre::Group* topoGrp = + meshGrp->getGroup("topologies")->getGroup(topoName); + axom::sidre::View* topoTypeView = topoGrp->getView("type"); + SLIC_ASSERT(std::string(topoTypeView->getString()) == "structured"); + topoTypeView->setString("unstructured"); + topoGrp->createView("elements/shape")->setString("hex"); + + axom::sidre::Group* topoElemGrp = topoGrp->getGroup("elements"); + axom::sidre::Group* topoDimsGrp = topoElemGrp->getGroup("dims"); + + // Assuming no ghost, but we should eventually support ghosts. + SLIC_ASSERT(!topoGrp->hasGroup("elements/offsets")); + SLIC_ASSERT(!topoGrp->hasGroup("elements/strides")); + + const axom::StackArray cShape { + axom::IndexType(topoDimsGrp->getView("i")->getNode().to_value()), + axom::IndexType(topoDimsGrp->getView("j")->getNode().to_value()), + axom::IndexType(topoDimsGrp->getView("k")->getNode().to_value())}; + const axom::StackArray vShape {cShape[0] + 1, + cShape[1] + 1, + cShape[2] + 1}; + + const axom::IndexType cCount = cShape[0] * cShape[1] * cShape[2]; + + const axom::StackArray connShape {cCount, 8}; + axom::sidre::View* connView = topoGrp->createViewWithShapeAndAllocate( + "elements/connectivity", + axom::sidre::detail::SidreTT::id, + 2, + connShape.begin()); + axom::ArrayView connArrayView( + static_cast(connView->getVoidPtr()), + connShape); + + axom::MDMapping cIdMapping(cShape, axom::ArrayStrideOrder::COLUMN); + axom::MDMapping vIdMapping(vShape, axom::ArrayStrideOrder::COLUMN); + + const axom::StackArray, 8> cornerOffsets { + {{0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1}, + {0, 1, 1}}}; + axom::for_all( + cCount, + AXOM_LAMBDA(axom::IndexType iCell) { + axom::StackArray cIdx = + cIdMapping.toMultiIndex(iCell); + for(int n = 0; n < 8; ++n) + { + const auto& cornerOffset = cornerOffsets[n]; + axom::StackArray vIdx; + for(int d = 0; d < DIM; ++d) + { + vIdx[d] = cIdx[d] + cornerOffset[d]; + } + axom::IndexType iVert = vIdMapping.toFlatIndex(vIdx); + connArrayView(iCell, n) = iVert; + } + }); + + const bool addExtraDataForMint = true; + if(addExtraDataForMint) + { + AXOM_ANNOTATE_BEGIN("add_extra"); + /* + Constructing a mint mesh from meshGrp fails unless we add some + extra data. Blueprint doesn't require this extra data. (The mesh + passes conduit's Blueprint verification.) This should be fixed, + or we should write better blueprint support utilities. + */ + /* + Make the coordinate arrays 2D to use mint::Mesh. + For some reason, mint::Mesh requires the arrays to be + 2D, even though the second dimension is always 1. + */ + axom::IndexType curShape[2]; + int curDim; + auto* valuesGrp = coordsetGrp->getGroup("values"); + curDim = valuesGrp->getView("x")->getShape(2, curShape); + assert(curDim == 1); + const axom::IndexType vertsShape[2] = {curShape[0], 1}; + valuesGrp->getView("x")->reshapeArray(2, vertsShape); + valuesGrp->getView("y")->reshapeArray(2, vertsShape); + valuesGrp->getView("z")->reshapeArray(2, vertsShape); + + // Make connectivity array 2D for the same reason. + auto* elementsGrp = topoGrp->getGroup("elements"); + auto* connView = elementsGrp->getView("connectivity"); + curDim = connView->getShape(2, curShape); + constexpr axom::IndexType NUM_VERTS_PER_HEX = 8; + SLIC_ASSERT(curDim == 1 || curDim == 2); + if(curDim == 1) + { + SLIC_ASSERT(curShape[0] % NUM_VERTS_PER_HEX == 0); + axom::IndexType connShape[2] = {curShape[0] / NUM_VERTS_PER_HEX, + NUM_VERTS_PER_HEX}; + connView->reshapeArray(2, connShape); + } + + // mint::Mesh requires connectivity strides, even though Blueprint doesn't. + elementsGrp->createViewScalar("stride", NUM_VERTS_PER_HEX); + + // mint::Mesh requires field group, even though Blueprint doesn't. + meshGrp->createGroup("fields"); + AXOM_ANNOTATE_END("add_extra"); + } + + #if defined(AXOM_DEBUG) && defined(AXOM_USE_CONDUIT) + AXOM_ANNOTATE_BEGIN("validate_post"); + conduit::Node info; + bool isValid = verifyBlueprintMesh(meshGrp, info); + SLIC_ASSERT_MSG(isValid, "Internal error: Generated mesh is invalid."); + AXOM_ANNOTATE_END("validate_post"); + #endif + + return; +} + + #if defined(AXOM_USE_CONDUIT) +bool verifyBlueprintMesh(const axom::sidre::Group* meshGrp, conduit::Node info) +{ + conduit::Node meshNode; + meshGrp->createNativeLayout(meshNode); + bool isValid = conduit::blueprint::mesh::verify(meshNode, info); + if(!isValid) info.print(); + return isValid; +} + #endif + +#endif // AXOM_USE_SIDRE + +void fill_cartesian_coords_3d(axom::runtime_policy::Policy runtimePolicy, + const primal::BoundingBox& domainBox, + axom::ArrayView& xView, + axom::ArrayView& yView, + axom::ArrayView& zView) +{ + if(runtimePolicy == axom::runtime_policy::Policy::seq) + { + fill_cartesian_coords_3d_impl(domainBox, xView, yView, zView); + } +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + if(runtimePolicy == axom::runtime_policy::Policy::omp) + { + fill_cartesian_coords_3d_impl(domainBox, xView, yView, zView); + } +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + if(runtimePolicy == axom::runtime_policy::Policy::cuda) + { + fill_cartesian_coords_3d_impl>(domainBox, + xView, + yView, + zView); + } +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + if(runtimePolicy == axom::runtime_policy::Policy::hip) + { + fill_cartesian_coords_3d_impl>(domainBox, + xView, + yView, + zView); + } +#endif +} + +template +void fill_cartesian_coords_3d_impl(const primal::BoundingBox& domainBox, + axom::ArrayView& xView, + axom::ArrayView& yView, + axom::ArrayView& zView) +{ +#if defined(AXOM_DEBUG) + using XS = axom::execution_space; + SLIC_ASSERT_MSG(XS::usesAllocId(xView.getAllocatorID()) && + XS::usesAllocId(yView.getAllocatorID()) && + XS::usesAllocId(zView.getAllocatorID()), + std::string("fill_cartesian_coords_3d_impl: alloc ids ") + + std::to_string(xView.getAllocatorID()) + ", " + + std::to_string(yView.getAllocatorID()) + " and " + + std::to_string(zView.getAllocatorID()) + + " are not all compatible with execution space " + XS::name()); +#endif + + const auto& shape = xView.shape(); + const auto& mapping = xView.mapping(); + + SLIC_ASSERT(shape == yView.shape()); + SLIC_ASSERT(shape == zView.shape()); + SLIC_ASSERT(mapping == yView.mapping()); + SLIC_ASSERT(mapping == zView.mapping()); + + // Mesh resolution + const axom::primal::NumericArray res {shape[0] - 1, + shape[1] - 1, + shape[2] - 1}; + + // Mesh spacings. + double dx = (domainBox.getMax()[0] - domainBox.getMin()[0]) / res[0]; + double dy = (domainBox.getMax()[1] - domainBox.getMin()[1]) / res[1]; + double dz = (domainBox.getMax()[2] - domainBox.getMin()[2]) / res[2]; + +#if defined(AXOM_USE_RAJA) + RAJA::RangeSegment kRange(0, shape[2]); + RAJA::RangeSegment jRange(0, shape[1]); + RAJA::RangeSegment iRange(0, shape[0]); + using EXEC_POL = + typename axom::internal::nested_for_exec::loop3d_policy; + auto order = mapping.getStrideOrder(); + if(int(order) & int(axom::ArrayStrideOrder::COLUMN)) + { + RAJA::kernel( + RAJA::make_tuple(iRange, jRange, kRange), + AXOM_LAMBDA(axom::IndexType i, axom::IndexType j, axom::IndexType k) { + xView(i, j, k) = domainBox.getMin()[0] + i * dx; + yView(i, j, k) = domainBox.getMin()[1] + j * dy; + zView(i, j, k) = domainBox.getMin()[2] + k * dz; + }); + } + else + { + RAJA::kernel( + RAJA::make_tuple(kRange, jRange, iRange), + AXOM_LAMBDA(axom::IndexType k, axom::IndexType j, axom::IndexType i) { + xView(i, j, k) = domainBox.getMin()[0] + i * dx; + yView(i, j, k) = domainBox.getMin()[1] + j * dy; + zView(i, j, k) = domainBox.getMin()[2] + k * dz; + }); + } +#else + for(int k = 0; k < res[2]; ++k) + { + for(int j = 0; j < res[1]; ++j) + { + for(int i = 0; i < res[0]; ++i) + { + xView(i, j, k) = domainBox.getMin()[0] + i * dx; + yView(i, j, k) = domainBox.getMin()[1] + j * dy; + zView(i, j, k) = domainBox.getMin()[2] + k * dz; + } + } + } +#endif + + return; +} + } // namespace util } // namespace quest } // namespace axom diff --git a/src/axom/quest/util/mesh_helpers.hpp b/src/axom/quest/util/mesh_helpers.hpp index 9bafb4375b..aeccb68bd5 100644 --- a/src/axom/quest/util/mesh_helpers.hpp +++ b/src/axom/quest/util/mesh_helpers.hpp @@ -8,6 +8,12 @@ #include "axom/config.hpp" #include "axom/primal.hpp" +#if defined(AXOM_USE_SIDRE) + #include "axom/sidre.hpp" +#endif +#if defined(AXOM_USE_CONDUIT) + #include +#endif #ifdef AXOM_USE_MFEM #include "mfem.hpp" @@ -55,9 +61,99 @@ mfem::Mesh* make_cartesian_mfem_mesh_3D(const primal::BoundingBox& bb const primal::NumericArray& res, int polynomial_order, bool reorder_space_filling = true); +#endif + +#if defined(AXOM_USE_SIDRE) +/*! + /brief Creates a 3D Cartesian mfem mesh lying within a given bounding box + + \param meshGrp Put the mesh in this Group + \param bbox The bounding box for the mesh + \param res The resolution of the mesh + \param topologyName Name of the blueprint topoloyy to use. + \param coordsetName Name of the blueprint coordset to use. + \param runtimePolicy Runtime policy, see axom::runtime_policy. + Memory in \c meshGrp must be compatible with the + specified policy. + + \return The meshGrp pointer +*/ +axom::sidre::Group* make_structured_blueprint_box_mesh( + axom::sidre::Group* meshGrp, + const primal::BoundingBox& bbox, + const primal::NumericArray& res, + const std::string& topologyName = "mesh", + const std::string& coordsetName = "coords", + axom::runtime_policy::Policy runtimePolicy = axom::runtime_policy::Policy::seq); + +axom::sidre::Group* make_unstructured_blueprint_box_mesh( + axom::sidre::Group* meshGrp, + const primal::BoundingBox& bbox, + const primal::NumericArray& res, + const std::string& topologyName = "mesh", + const std::string& coordsetName = "coords", + axom::runtime_policy::Policy runtimePolicy = axom::runtime_policy::Policy::seq); + +/*! + \brief Convert a structured explicit blueprint mesh to unstructured. + \param meshGrp Put the mesh in this Group + \param topologyName Name of the blueprint topoloyy to use. + \param runtimePolicy Runtime policy, see axom::runtime_policy. + Memory in \c meshGrp must be compatible with the + specified policy. + + All input mesh data are expected to have the allocator id of + meshGrp->getDefaultAllocatorID(). On output, they will also have + the same allocator id, even if the intermediate steps have to + transfers memory to and from another space. +*/ +void convert_blueprint_structured_explicit_to_unstructured( + axom::sidre::Group* meshGrp, + const std::string& topoName, + axom::runtime_policy::Policy runtimePolicy); + +template +void convert_blueprint_structured_explicit_to_unstructured_impl( + axom::sidre::Group* meshGrp, + const std::string& topoName); + + #if defined(AXOM_USE_CONDUIT) +/*! + \brief Check if blueprint mesh is valid. +*/ +bool verifyBlueprintMesh(const axom::sidre::Group* meshGrp, conduit::Node info); + #endif #endif +/*! + @brief Fill in structured mesh cartesian coordinates. + \param xView Vertex x-values array + \param yView Vertex y-values array + \param zView Vertex z-values array + \param domainBox Physical domain +*/ +void fill_cartesian_coords_3d(axom::runtime_policy::Policy runtimePolicy, + const primal::BoundingBox& domainBox, + axom::ArrayView& xView, + axom::ArrayView& yView, + axom::ArrayView& zView); + +/*! + @brief Fill in structured mesh cartesian coordinates. + \tparam ExecSpace Execution space, e.g. \c axom::SEQ_EXEC. + @see axom::execution_space + \param xView Vertex x-values array + \param yView Vertex y-values array + \param zView Vertex z-values array + \param domainBox Physical domain +*/ +template +void fill_cartesian_coords_3d_impl(const primal::BoundingBox& domainBox, + axom::ArrayView& xView, + axom::ArrayView& yView, + axom::ArrayView& zView); + } // namespace util } // namespace quest } // namespace axom diff --git a/src/axom/sidre/core/Group.hpp b/src/axom/sidre/core/Group.hpp index 7209f92928..da64965dea 100644 --- a/src/axom/sidre/core/Group.hpp +++ b/src/axom/sidre/core/Group.hpp @@ -1346,7 +1346,7 @@ class Group * \brief Print given number of levels of Group sub-tree * starting at this Group object to an output stream. */ - void printTree(const int nlevels, std::ostream& os) const; + void printTree(const int nlevels, std::ostream& os = std::cout) const; //@} @@ -1802,11 +1802,13 @@ class Group /*! * \brief Detach Child Group with given name from this Group. + * \return the detached Group. */ Group* detachGroup(const std::string& name); /*! * \brief Detach Child Group with given index from this Group. + * \return the detached Group. */ Group* detachGroup(IndexType idx); diff --git a/src/axom/sidre/core/View.hpp b/src/axom/sidre/core/View.hpp index d451bc18e8..02a4f1b8c1 100644 --- a/src/axom/sidre/core/View.hpp +++ b/src/axom/sidre/core/View.hpp @@ -827,6 +827,9 @@ class View return getData(); } + /// \overload + Node::ConstValue getArray() const { return getData(); } + /*! * \brief Returns a pointer to the string contained in the view. *