diff --git a/src/axom/core/utilities/Utilities.hpp b/src/axom/core/utilities/Utilities.hpp index df17f8d362..510d505f65 100644 --- a/src/axom/core/utilities/Utilities.hpp +++ b/src/axom/core/utilities/Utilities.hpp @@ -348,6 +348,16 @@ inline AXOM_HOST_DEVICE bool isNearlyEqualRelative(RealType a, // return abs(a-b) <= max(absThresh, relThresh * maxFabs ); } +/*! + * \brief Sign of a value of any type that supports comparison and + * negation operators. + */ +template +inline int sign_of(const T& v, const T& eps = {0}) +{ + return v > eps ? 1 : v < -eps ? -1 : 0; +} + /*! * \brief Insertion sort of an array. * \accelerated diff --git a/src/axom/klee/Geometry.cpp b/src/axom/klee/Geometry.cpp index 09c8b03e6e..ff8a0f8c06 100644 --- a/src/axom/klee/Geometry.cpp +++ b/src/axom/klee/Geometry.cpp @@ -154,7 +154,7 @@ void Geometry::populateGeomInfo() m_discreteFunction = axom::Array(2, 2); m_discreteFunction(0, 0) = 0.0; m_discreteFunction(0, 1) = cone.getBaseRadius(); - m_discreteFunction(1, 1) = cone.getLength(); + m_discreteFunction(1, 0) = cone.getLength(); m_discreteFunction(1, 1) = cone.getTopRadius(); m_geomInfo["discreteFunction"].set(m_discreteFunction.data(), m_discreteFunction.size()); m_geomInfo["sorOrigin"].set(cone.getBaseCenter().data(), 3); diff --git a/src/axom/primal/geometry/Cone.hpp b/src/axom/primal/geometry/Cone.hpp index 8eb2e827d4..4584512364 100644 --- a/src/axom/primal/geometry/Cone.hpp +++ b/src/axom/primal/geometry/Cone.hpp @@ -143,7 +143,8 @@ class Cone } /*! - * \brief Returns the algebraic volume of the cone + * \brief Returns the algebraic volume of the cone, + * which is negative if the length is negative. * * Volume is only defined when NDIMS == 3. */ diff --git a/src/axom/primal/geometry/CoordinateTransformer.hpp b/src/axom/primal/geometry/CoordinateTransformer.hpp index ef983c28f7..798350fb37 100644 --- a/src/axom/primal/geometry/CoordinateTransformer.hpp +++ b/src/axom/primal/geometry/CoordinateTransformer.hpp @@ -66,6 +66,16 @@ class CoordinateTransformer */ CoordinateTransformer(const numerics::Matrix& matrix) { setMatrix(matrix); } + /*! + * @brief Contruct transformer that moves 4 starting points to 4 + * destination points. + */ + AXOM_HOST_DEVICE CoordinateTransformer(const primal::Point* startPts, + const primal::Point* destPts) + { + setByTerminusPts(startPts, destPts); + } + /*! * @brief Set the matrix, discarding the current transformation. * @param matrix [in] The transformation matrix for homogeneous diff --git a/src/axom/primal/geometry/Sphere.hpp b/src/axom/primal/geometry/Sphere.hpp index 9f89e37a32..0b52f3cbc3 100644 --- a/src/axom/primal/geometry/Sphere.hpp +++ b/src/axom/primal/geometry/Sphere.hpp @@ -178,6 +178,18 @@ class Sphere AXOM_HOST_DEVICE inline bool intersectsWith(const Sphere& sphere, double TOL = 1.e-9) const; + /*! + * \brief Tests if this sphere completely contains another sphere. + * + * \param [in] other The sphere object to check for containment + * \param [in] margin Amount that this sphere must contain the other sphere by. + * Positive means that the other sphere is "more inside". + * + * \return true if this sphere contains the other, false otherwise. + */ + AXOM_HOST_DEVICE + inline bool contains(const Sphere& other, double margin = 0.0) const; + /*! * \brief Prints the Sphere information in the given output stream. * \param [in,out] os the output stream to write to. @@ -233,6 +245,14 @@ AXOM_HOST_DEVICE inline bool Sphere::intersectsWith(const Sphere +AXOM_HOST_DEVICE inline bool Sphere::contains(const Sphere& sphere, double TOL) const +{ + const T center_sep = VectorType(sphere.getCenter(), m_center).norm(); + return (m_radius > center_sep + sphere.getRadius() + TOL); +} + //------------------------------------------------------------------------------ // implementation of free functions //------------------------------------------------------------------------------ diff --git a/src/axom/primal/tests/primal_sphere.cpp b/src/axom/primal/tests/primal_sphere.cpp index 3777b75f80..4c56570477 100644 --- a/src/axom/primal/tests/primal_sphere.cpp +++ b/src/axom/primal/tests/primal_sphere.cpp @@ -150,6 +150,31 @@ void check_sphere_intersection() EXPECT_FALSE(S0.intersectsWith(S3)); } +//------------------------------------------------------------------------------ +template +void check_sphere_containment() +{ + using PointType = primal::Point; + using SphereType = primal::Sphere; + + PointType center {0.0, 0.0, 0.0}; + double tol = 1e-12; + + // STEP 0: test fully containing + SphereType S0; + EXPECT_TRUE(S0.contains(S0, -tol)); + + // STEP 1: test barely not containing. + center[0] = tol; + SphereType S1(center); + EXPECT_FALSE(S0.contains(S1)); + + // STEP 2: test partial containment. + center[0] = 0.5; + SphereType S3(center); + EXPECT_FALSE(S0.contains(S3)); +} + //------------------------------------------------------------------------------ template void check_copy_constructor() @@ -244,6 +269,13 @@ TEST(primal_sphere, sphere_sphere_intersection) check_sphere_intersection<3>(); } +//------------------------------------------------------------------------------ +TEST(primal_sphere, sphere_sphere_containment) +{ + check_sphere_containment<2>(); + check_sphere_containment<3>(); +} + //------------------------------------------------------------------------------ int main(int argc, char* argv[]) { diff --git a/src/axom/quest/CMakeLists.txt b/src/axom/quest/CMakeLists.txt index a78b1f052c..64648e21db 100644 --- a/src/axom/quest/CMakeLists.txt +++ b/src/axom/quest/CMakeLists.txt @@ -150,11 +150,14 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE) if(RAJA_FOUND) list(APPEND quest_headers MeshClipper.hpp MeshClipperStrategy.hpp + detail/clipping/Plane3DClipper.hpp detail/clipping/TetClipper.hpp detail/clipping/MeshClipperImpl.hpp) list(APPEND quest_sources MeshClipper.cpp + MeshClipperStrategy.cpp + detail/clipping/Plane3DClipper.cpp detail/clipping/TetClipper.cpp - MeshClipperStrategy.cpp) + detail/clipping/MeshClipperImpl.hpp) endif() endif() diff --git a/src/axom/quest/Discretize.cpp b/src/axom/quest/Discretize.cpp index fcd317ba3e..448fdc1246 100644 --- a/src/axom/quest/Discretize.cpp +++ b/src/axom/quest/Discretize.cpp @@ -108,6 +108,8 @@ OctType new_inscribed_oct(const SphereType& sphere, OctType& o, int s, int t, in * * This routine allocates an array pointed to by \a out. The caller is responsible * to free the array. + * + * TODO: If possible, port to GPU (rewrite to be data-parallel). */ bool discretize(const SphereType& sphere, int levels, axom::Array& out, int& octcount) { @@ -125,7 +127,7 @@ bool discretize(const SphereType& sphere, int levels, axom::Array& out, octcount = count_sphere_octahedra(levels); - out = axom::Array(octcount, octcount); + out.resize(octcount); // index points to an octahedron of the last generation. We'll generate // new octahedra based on out[index]. diff --git a/src/axom/quest/Discretize.hpp b/src/axom/quest/Discretize.hpp index a815dff9e9..06396161d7 100644 --- a/src/axom/quest/Discretize.hpp +++ b/src/axom/quest/Discretize.hpp @@ -50,7 +50,8 @@ bool discretize(const SphereType& s, int levels, axom::Array& out, int& /*! * \brief Given a 2D polyline revolved around the positive X-axis, allocate * and return a list of Octahedra approximating the shape. - * \param [in] polyline The polyline to revolve around the X-axis + * \param [in] polyline The polyline to revolve around the X-axis. + * Data should be in a host-accessible memory space. * \param [in] len The number of points in \a polyline * \param [in] levels The number of refinements to perform, in addition to * a central level-zero octahedron in each segment diff --git a/src/axom/quest/InOutOctree.hpp b/src/axom/quest/InOutOctree.hpp index 4bfe21fbcb..a37fa786db 100644 --- a/src/axom/quest/InOutOctree.hpp +++ b/src/axom/quest/InOutOctree.hpp @@ -17,7 +17,6 @@ #include "axom/slic.hpp" #include "axom/slam.hpp" #include "axom/primal.hpp" -#include "axom/mint.hpp" #include "axom/spin.hpp" #include "detail/inout/BlockData.hpp" diff --git a/src/axom/quest/IntersectionShaper.hpp b/src/axom/quest/IntersectionShaper.hpp index e3a626d7d0..b7c5dd42fd 100644 --- a/src/axom/quest/IntersectionShaper.hpp +++ b/src/axom/quest/IntersectionShaper.hpp @@ -825,6 +825,7 @@ class IntersectionShaper : public Shaper // Generate the Octahedra // (octahedra m_octs will be on device) + m_octs = axom::Array(0, 0, axom::execution_space::allocatorID()); const bool disc_status = axom::quest::discretize(polyline, polyline_size, m_level, m_octs, m_octcount); @@ -1967,13 +1968,15 @@ class IntersectionShaper : public Shaper 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()) + if(fieldsGrp != nullptr) { - std::string materialName = fieldNameToMaterialName(group.getName()); - if(!materialName.empty()) + for(auto& group : fieldsGrp->groups()) { - materialNames.emplace_back(materialName); + std::string materialName = fieldNameToMaterialName(group.getName()); + if(!materialName.empty()) + { + materialNames.emplace_back(materialName); + } } } } diff --git a/src/axom/quest/MeshClipper.cpp b/src/axom/quest/MeshClipper.cpp index b2e5b7910d..4e1193d746 100644 --- a/src/axom/quest/MeshClipper.cpp +++ b/src/axom/quest/MeshClipper.cpp @@ -25,7 +25,18 @@ MeshClipper::MeshClipper(quest::experimental::ShapeMesh& shapeMesh, , m_strategy(strategy) , m_impl(newImpl()) , m_verbose(false) -{ } + , m_screenLevel(3) +{ + // Initialize statistics. + m_counterStats["cellsIn"].set_int64(0); + m_counterStats["cellsOn"].set_int64(0); + m_counterStats["cellsOut"].set_int64(0); + m_counterStats["tetsIn"].set_int64(0); + m_counterStats["tetsOn"].set_int64(0); + m_counterStats["tetsOut"].set_int64(0); + m_counterStats["clips"].set_int64(0); + m_counterStats["contribs"].set_int64(0); +} void MeshClipper::clip(axom::Array& ovlap) { @@ -35,7 +46,7 @@ void MeshClipper::clip(axom::Array& ovlap) // Resize output array and use appropriate allocator. if(ovlap.size() < cellCount || ovlap.getAllocatorID() != allocId) { - AXOM_ANNOTATE_SCOPE("MeshClipper::clip_alloc"); + AXOM_ANNOTATE_SCOPE("MeshClipper:clip_alloc"); ovlap = axom::Array(ArrayOptions::Uninitialized(), cellCount, cellCount, allocId); } clip(ovlap.view()); @@ -45,8 +56,8 @@ void MeshClipper::clip(axom::Array& ovlap) * @brief Orchestrates the geometry clipping by using the capabilities of the * MeshClipperStrategy implementation. * - * If the strategy can label cells as inside/on/outside geometry - * boundary, use that to reduce reliance on expensive clipping methods. + * If the strategy can label cells/tets as inside/on/outside geometry + * boundary, use that to reduce use of expensive primitive clipping methods. * * Regardless of labeling, try to use specialized clipping first. * If specialized methods aren't implemented, resort to discretizing @@ -58,50 +69,86 @@ void MeshClipper::clip(axom::ArrayView ovlap) const axom::IndexType cellCount = m_shapeMesh.getCellCount(); SLIC_ASSERT(ovlap.size() == cellCount); + auto& cellsInCount = *m_counterStats["cellsIn"].as_int64_ptr(); + auto& cellsOnCount = *m_counterStats["cellsOn"].as_int64_ptr(); + auto& cellsOutCount = *m_counterStats["cellsOut"].as_int64_ptr(); + auto& tetsInCount = *m_counterStats["tetsIn"].as_int64_ptr(); + auto& tetsOnCount = *m_counterStats["tetsOn"].as_int64_ptr(); + auto& tetsOutCount = *m_counterStats["tetsOut"].as_int64_ptr(); + // Try to label cells as inside, outside or on shape boundary axom::Array cellLabels; - bool withCellInOut = m_strategy->labelCellsInOut(m_shapeMesh, cellLabels); + bool withCellInOut = false; + if(m_screenLevel >= 1) + { + AXOM_ANNOTATE_BEGIN("MeshClipper:label_cells"); + withCellInOut = m_strategy->labelCellsInOut(m_shapeMesh, cellLabels); + AXOM_ANNOTATE_END("MeshClipper:label_cells"); + } bool done = false; if(withCellInOut) { - SLIC_ERROR_IF( - cellLabels.size() != m_shapeMesh.getCellCount(), - axom::fmt::format("MeshClipperStrategy '{}' did not return the correct array size of {}", - m_strategy->name(), - m_shapeMesh.getCellCount())); + SLIC_ERROR_IF(cellLabels.size() != m_shapeMesh.getCellCount(), + axom::fmt::format("MeshClipperStrategy '{}' did not return the correct" + " cell label array size of {}", + m_strategy->name(), + m_shapeMesh.getCellCount())); SLIC_ERROR_IF(cellLabels.getAllocatorID() != allocId, - axom::fmt::format("MeshClipperStrategy '{}' failed to provide cellLabels data " - "with the required allocator id {}", + axom::fmt::format("MeshClipperStrategy '{}' failed to provide" + " cellLabels data with the required allocator id {}", m_strategy->name(), allocId)); + // Counting labels is non-essential and presumed to be relatively very fast. + getLabelCounts(cellLabels, cellsInCount, cellsOnCount, cellsOutCount); if(m_verbose) { - logLabelStats(cellLabels, "cells"); + logClippingStats(); } - AXOM_ANNOTATE_BEGIN("MeshClipper::processInOut"); + AXOM_ANNOTATE_BEGIN("MeshClipper:process_in_out"); m_impl->initVolumeOverlaps(cellLabels.view(), ovlap); axom::Array cellsOnBdry; m_impl->collectOnIndices(cellLabels.view(), cellsOnBdry); + SLIC_ASSERT(cellsOnBdry.size() == cellsOnCount); axom::Array tetLabels; - bool withTetInOut = m_strategy->labelTetsInOut(m_shapeMesh, cellsOnBdry.view(), tetLabels); + bool withTetInOut = false; + if(m_screenLevel >= 2) + { + AXOM_ANNOTATE_BEGIN("MeshClipper:label_tets"); + withTetInOut = m_strategy->labelTetsInOut(m_shapeMesh, cellsOnBdry.view(), tetLabels); + AXOM_ANNOTATE_END("MeshClipper:label_tets"); + } axom::Array tetsOnBdry; if(withTetInOut) { + SLIC_ERROR_IF(tetLabels.size() != NUM_TETS_PER_HEX * cellsOnBdry.size(), + axom::fmt::format("MeshClipperStrategy '{}' did not return the correct" + " tet label array size of {}", + m_strategy->name(), + NUM_TETS_PER_HEX * cellsOnBdry.size())); + SLIC_ERROR_IF(tetLabels.getAllocatorID() != allocId, + axom::fmt::format("MeshClipperStrategy '{}' failed to provide" + "tetLabels data with the required allocator id {}", + m_strategy->name(), + allocId)); + + // Counting labels is non-essential and presumed to be very fast. + getLabelCounts(tetLabels, tetsInCount, tetsOnCount, tetsOutCount); if(m_verbose) { - logLabelStats(tetLabels, "tets"); + logClippingStats(); } + m_impl->collectOnIndices(tetLabels.view(), tetsOnBdry); - m_impl->remapTetIndices(tetsOnBdry, cellsOnBdry); + m_impl->remapTetIndices(cellsOnBdry, tetsOnBdry); SLIC_ASSERT(tetsOnBdry.getAllocatorID() == m_shapeMesh.getAllocatorID()); SLIC_ASSERT(tetsOnBdry.size() <= cellsOnBdry.size() * NUM_TETS_PER_HEX); @@ -109,51 +156,47 @@ void MeshClipper::clip(axom::ArrayView ovlap) m_impl->addVolumesOfInteriorTets(cellsOnBdry.view(), tetLabels.view(), ovlap); } - AXOM_ANNOTATE_END("MeshClipper::processInOut"); + AXOM_ANNOTATE_END("MeshClipper:process_in_out"); // // If implementation has a specialized clip, use it. // + AXOM_ANNOTATE_BEGIN("MeshClipper:specialized_clip"); if(withTetInOut) { - done = m_strategy->specializedClipTets(m_shapeMesh, ovlap, tetsOnBdry); + done = m_strategy->specializedClipTets(m_shapeMesh, ovlap, tetsOnBdry, m_counterStats); } else { - done = m_strategy->specializedClipCells(m_shapeMesh, ovlap, cellsOnBdry); + done = m_strategy->specializedClipCells(m_shapeMesh, ovlap, cellsOnBdry, m_counterStats); } + AXOM_ANNOTATE_END("MeshClipper:specialized_clip"); if(!done) { - AXOM_ANNOTATE_SCOPE("MeshClipper::clip3D_limited"); + AXOM_ANNOTATE_SCOPE("MeshClipper:clip_fcn"); if(withTetInOut) { - m_impl->computeClipVolumes3DTets(tetsOnBdry.view(), ovlap); + m_impl->computeClipVolumes3DTets(tetsOnBdry.view(), ovlap, m_counterStats); } else { - m_impl->computeClipVolumes3D(cellsOnBdry.view(), ovlap); + m_impl->computeClipVolumes3D(cellsOnBdry.view(), ovlap, m_counterStats); } } - - m_localCellInCount = cellsOnBdry.size(); } else // !withCellInOut { - std::string msg = - axom::fmt::format("MeshClipper strategy '{}' did not provide in/out cell labels.\n", - m_strategy->name()); - SLIC_INFO(msg); m_impl->zeroVolumeOverlaps(ovlap); - done = m_strategy->specializedClipCells(m_shapeMesh, ovlap); + AXOM_ANNOTATE_BEGIN("MeshClipper:specialized_clip"); + done = m_strategy->specializedClipCells(m_shapeMesh, ovlap, m_counterStats); + AXOM_ANNOTATE_END("MeshClipper:specialized_clip"); if(!done) { - AXOM_ANNOTATE_SCOPE("MeshClipper::clip3D"); - m_impl->computeClipVolumes3D(ovlap); + AXOM_ANNOTATE_SCOPE("MeshClipper:clip_fcn"); + m_impl->computeClipVolumes3D(ovlap, m_counterStats); } - - m_localCellInCount = cellCount; } } @@ -197,52 +240,90 @@ std::unique_ptr MeshClipper::newImpl() return impl; } -void MeshClipper::logLabelStats(axom::ArrayView labels, const std::string& labelType) +#if defined(AXOM_USE_MPI) +template +void globalReduce(axom::Array& values, int reduceOp) { - axom::IndexType countsa[4]; - axom::IndexType countsb[4]; - getLabelCounts(labels, countsa[0], countsa[1], countsa[2]); - countsa[3] = m_shapeMesh.getCellCount(); -#ifdef AXOM_USE_MPI - MPI_Reduce(countsa, countsb, 4, axom::mpi_traits::type, MPI_SUM, 0, MPI_COMM_WORLD); + axom::Array localValues(values); + MPI_Allreduce(localValues.data(), + values.data(), + values.size(), + axom::mpi_traits::type, + reduceOp, + MPI_COMM_WORLD); #endif - std::string msg = axom::fmt::format( - "MeshClipper strategy '{}' globally labeled {} {} inside, {} on and {} outside, for mesh with " - "{} cells ({} tets)\n", - m_strategy->name(), - labelType, - countsb[0], - countsb[1], - countsb[2], - countsb[3], - countsb[3] * NUM_TETS_PER_HEX); - SLIC_INFO(msg); } -void MeshClipper::getClippingStats(axom::IndexType& localCellInCount, - axom::IndexType& globalCellInCount, - axom::IndexType& maxLocalCellInCount) const +void MeshClipper::accumulateClippingStats(conduit::Node& curStats, const conduit::Node& newStats) { - localCellInCount = m_localCellInCount; -#ifdef AXOM_USE_MPI - MPI_Reduce(&localCellInCount, - &globalCellInCount, - 1, - axom::mpi_traits::type, - MPI_SUM, - 0, - MPI_COMM_WORLD); - MPI_Reduce(&localCellInCount, - &maxLocalCellInCount, - 1, - axom::mpi_traits::type, - MPI_MAX, - 0, - MPI_COMM_WORLD); -#else - maxLocalCellInCount = localCellInCount; - globalCellInCount = localCellInCount; + for(int i = 0; i < newStats.number_of_children(); ++i) + { + const auto& newStat = newStats.child(i); + SLIC_ERROR_IF(!newStat.dtype().is_integer(), + "MeshClipper statistic must be integer" + " (at least until a need for floats arises)."); + auto& currentStat = curStats[newStat.name()]; + if(currentStat.dtype().is_empty()) + { + currentStat.set_int64(newStat.as_int64()); + } + else + { + *currentStat.as_int64_ptr() += newStat.as_int64(); + } + } +} + +conduit::Node MeshClipper::getGlobalClippingStats() const +{ + conduit::Node stats; + auto& locNode = stats["loc"]; + auto& maxNode = stats["max"]; + auto& sumNode = stats["sum"]; + + locNode.set(m_counterStats); + sumNode.set(m_counterStats); + maxNode.set(m_counterStats); + +#if defined(AXOM_USE_MPI) + // Do sum and max reductions. + axom::Array sums(0, sumNode.number_of_children()); + for(int i = 0; i < sumNode.number_of_children(); ++i) + { + sums.push_back(locNode.child(i).as_int64()); + } + axom::Array maxs(sums); + globalReduce(maxs, MPI_MAX); + globalReduce(sums, MPI_SUM); + + for(int i = 0; i < sumNode.number_of_children(); ++i) + { + *maxNode.child(i).as_int64_ptr() = maxs[i]; + *sumNode.child(i).as_int64_ptr() = sums[i]; + } #endif + + return stats; +} + +void MeshClipper::logClippingStats(bool local, bool sum, bool max) const +{ + conduit::Node stats = getGlobalClippingStats(); + if(local) + { + SLIC_INFO(std::string("MeshClipper loc-stats: ") + + stats["loc"].to_string("yaml", 2, 0, "", " ")); + } + if(sum) + { + SLIC_INFO(std::string("MeshClipper sum-stats: ") + + stats["sum"].to_string("yaml", 2, 0, "", " ")); + } + if(max) + { + SLIC_INFO(std::string("MeshClipper max-stats: ") + + stats["max"].to_string("yaml", 2, 0, "", " ")); + } } } // namespace experimental diff --git a/src/axom/quest/MeshClipper.hpp b/src/axom/quest/MeshClipper.hpp index d91df1e275..4d0e9fcf5d 100644 --- a/src/axom/quest/MeshClipper.hpp +++ b/src/axom/quest/MeshClipper.hpp @@ -11,6 +11,7 @@ #include "axom/klee/Geometry.hpp" #include "axom/quest/MeshClipperStrategy.hpp" #include "axom/quest/ShapeMesh.hpp" +#include "conduit/conduit_node.hpp" namespace axom { @@ -38,7 +39,7 @@ class MeshClipper //!@brief Whether an element is in, out or on shape boundary. using LabelType = MeshClipperStrategy::LabelType; - static constexpr axom::IndexType NUM_TETS_PER_HEX = MeshClipperStrategy::NUM_TETS_PER_HEX; + static constexpr axom::IndexType NUM_TETS_PER_HEX = ShapeMesh::NUM_TETS_PER_HEX; /*! * @brief Construct a shape clipper @@ -81,6 +82,62 @@ class MeshClipper //!@brief Dimension of the shape (2 or 3) int dimension() const { return m_shapeMesh.dimension(); } + //@{ + + /*! + * @brief Return the number of times primitive clipping was used, + * intended for developer use. + */ + std::int64_t getClipCount() const { return m_counterStats["clips"].as_int64(); } + + /*! + * @brief Return the number of times primitive clipping contributed to the mesh clip volume, + * intended for developer use. + */ + std::int64_t getContribCount() const { return m_counterStats["contribs"].as_int64(); } + + /*! + * @brief Log clipping statistics. + * Intended for developer use. + * + * This is a collective method if MPI-parallel. + */ + void logClippingStats(bool local = false, bool sum = true, bool max = false) const; + + /*! + * @brief Get local assorted clipping statistics, + * intended for developer use. + */ + const conduit::Node& getClippingStats() const { return m_counterStats; } + + /*! + * @brief Get global assorted clipping statistics, + * intended for developer use. + * + * This is a collective method if MPI-parallel. + */ + conduit::Node getGlobalClippingStats() const; + + /*! + * @brief Set the level of screening, + * intended for developer use. + */ + void setScreenLevel(int screenLevel) { m_screenLevel = screenLevel; } + + /*! + * @brief Get the level of screening, + * intended for developer use. + */ + int getScreenLevel() const { return m_screenLevel; } + + /*! + * @brief Add new stats to current stats, + * intended for developer use. + */ + static void accumulateClippingStats(conduit::Node& curStats, const conduit::Node& newStats); + + //@} + /*! * @brief Single interface for methods implemented with * execution space templates. @@ -136,11 +193,12 @@ class MeshClipper axom::ArrayView ovlap) = 0; //!@brief Compute clip volumes for every cell. - virtual void computeClipVolumes3D(axom::ArrayView ovlap) = 0; + virtual void computeClipVolumes3D(axom::ArrayView ovlap, conduit::Node& statistics) = 0; //!@brief Compute clip volumes for cell in an index list. virtual void computeClipVolumes3D(const axom::ArrayView& cellIndices, - axom::ArrayView ovlap) = 0; + axom::ArrayView ovlap, + conduit::Node& statistics) = 0; /*! * @brief Compute clip volumes for cell tets in an index list. @@ -149,28 +207,24 @@ class MeshClipper * NUM_TETS_PER_HEX tets and stored consecutively. */ virtual void computeClipVolumes3DTets(const axom::ArrayView& tetIndices, - axom::ArrayView ovlap) = 0; + axom::ArrayView ovlap, + conduit::Node& statistics) = 0; //!@brief Count the number of labels of each type. virtual void getLabelCounts(axom::ArrayView labels, - axom::IndexType& inCount, - axom::IndexType& onCount, - axom::IndexType& outCount) = 0; + std::int64_t& inCount, + std::int64_t& onCount, + std::int64_t& outCount) = 0; ShapeMesh& getShapeMesh() { return m_myClipper.m_shapeMesh; } MeshClipperStrategy& getStrategy() { return *m_myClipper.m_strategy; } - private: + protected: //!@brief The MeshClipper that owns this Impl. MeshClipper& m_myClipper; }; - //! @brief For assessments, not general use. - void getClippingStats(axom::IndexType& localCellInCount, - axom::IndexType& globalCellInCount, - axom::IndexType& maxLocalCellInCount) const; - private: friend Impl; @@ -187,13 +241,13 @@ class MeshClipper * for multiple execution spaces. */ - ///@{ - //! @name Statistics - axom::IndexType m_localCellInCount {0}; - ///@} + //! @brief Statistics + conduit::Node m_counterStats; bool m_verbose; + int m_screenLevel; + #if defined(__CUDACC__) public: #endif @@ -204,14 +258,12 @@ class MeshClipper //!@name Convenience methods //!@brief Count the number of labels of each type. void getLabelCounts(const axom::Array& labels, - axom::IndexType& inCount, - axom::IndexType& onCount, - axom::IndexType& outCount) + std::int64_t& inCount, + std::int64_t& onCount, + std::int64_t& outCount) { m_impl->getLabelCounts(labels, inCount, onCount, outCount); } - - void logLabelStats(axom::ArrayView labels, const std::string& labelType); //@} }; diff --git a/src/axom/quest/MeshClipperStrategy.hpp b/src/axom/quest/MeshClipperStrategy.hpp index 48de9661ae..420e6f088d 100644 --- a/src/axom/quest/MeshClipperStrategy.hpp +++ b/src/axom/quest/MeshClipperStrategy.hpp @@ -53,12 +53,13 @@ namespace experimental * false if it was a no-op. * Subclasses of MeshClipperStrategy must implement either - * - a @c specializedClipCells method or + * - a @c specializedClipCells or @c specializedClipTets method or * - one of the @c getShapesAs...() methods. * The former is prefered if the use of geometry-specific information - * can make it faster. @c labelCellsInOut is optional but if provided, - * it can improve performance by limiting the slower clipping steps - * to a subset of cells. @c getBoundingBox2D or @c getBoundingBox3D + * can make it faster. @c labelCellsInOut and @c labelTetsInOut + * are optional but if provided, + * they can improve performance by limiting the slower clipping steps + * to a smaller subset. @c getBoundingBox2D or @c getBoundingBox3D * can also improve performance by reducing computation. */ class MeshClipperStrategy @@ -95,9 +96,8 @@ class MeshClipperStrategy using Ray2DType = axom::primal::Ray; using Segment2DType = axom::primal::Segment; - //!@brief Number of tetrahedra per hexahedron decomposes into - // @internal We could use a more efficient 18-tet decomposition in the future. static constexpr axom::IndexType NUM_TETS_PER_HEX = ShapeMesh::NUM_TETS_PER_HEX; + static constexpr axom::IndexType NUM_VERTS_PER_CELL_3D = ShapeMesh::NUM_VERTS_PER_CELL_3D; /*! * @brief Construct a strategy for the given klee::Geometry object. @@ -164,6 +164,9 @@ class MeshClipperStrategy * skip if it's not. It's safe to label cells as on the boundary if * it can't be positively determined as inside or outside. * + * Degenerate cells have zero volume and should be labeled outside + * for best clipping performance. + * * @return Whether the operation was done. (A false means * not done.) * @@ -191,10 +194,14 @@ class MeshClipperStrategy * * Indices [i*NUM_TETS_PER_HEX, (i+1)*NUM_TETS_PER_HEX) in \c tetLabels * correspond to parent cell index \c c = \c cellIds[i]. - * The \c NUM_TETS_PER_HEX tets in cell \c cid have indices + * + * The \c NUM_TETS_PER_HEX tets in cell \c c have indices * [c*NUM_TETS_PER_HEX, (c+1)*NUM_TETS_PER_HEX). * in \c shapeMesh.getCellsAsTets(). * + * Degenerate tets have zero volume and should be labeled outside + * for best clipping performance. + * * If implementation returns true, it should ensure these * post-conditions hold: * @post tetLabels.size() == NUM_TETS_PER_HEX * cellIds.size() @@ -219,6 +226,8 @@ class MeshClipperStrategy * @param [in] shapeMesh Blueprint mesh to shape into. * @param [out] ovlap Shape overlap volume of each cell * in the \c shapeMesh. It's initialized to zeros. + * @param [out] statistics Optional statistics to record + * consisting of child nodes with integer values. * * The default implementation has no specialized method, * so it's a no-op and returns false. @@ -231,16 +240,21 @@ class MeshClipperStrategy * This method need not be implemented if labelCellsInOut() * returns true. * + * Setting the statistics is not required except for getting + * accurate statistics. + * * If implementation returns true, it should ensure these * post-conditions hold: * @post ovlap.size() == shapeMesh.getCellCount() * @post ovlap.getAllocatorID() == shapeMesh.getAllocatorId() */ virtual bool specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, - axom::ArrayView ovlap) + axom::ArrayView ovlap, + conduit::Node& statistics) { AXOM_UNUSED_VAR(shapeMesh); AXOM_UNUSED_VAR(ovlap); + AXOM_UNUSED_VAR(statistics); return false; } @@ -253,18 +267,24 @@ class MeshClipperStrategy * in \c shapeMesh, initialized to the cell volumes * for cell inside the shape and zero for other cells. * @param [in] cellIds Limit computation to these cell ids. + * @param [out] statistics Optional statistics to record + * consisting of child nodes with integer values. * * The default implementation has no specialized method, * so it's a no-op and returns false. * - * If this method returns false, then exactly one of the - * shape discretization methods must be provided. + * If this method returns false, then exactly one of + * getGeometryAsTets() or getGeometryAsOcts() methods must be + * provided so MeshClipper can use the general clipping methods. * * @return True if clipping was done and false if a no-op. * * This method need not be implemented if labelCellsInOut() * returns false. * + * Setting the statistics is not required except for getting + * accurate statistics. + * * @pre @c ovlap is pre-initialized for the implementation * to add or subtract partial volumes to individual cells. * @@ -275,11 +295,13 @@ class MeshClipperStrategy */ virtual bool specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, axom::ArrayView ovlap, - const axom::ArrayView& cellIds) + const axom::ArrayView& cellIds, + conduit::Node& statistics) { AXOM_UNUSED_VAR(shapeMesh); AXOM_UNUSED_VAR(ovlap); AXOM_UNUSED_VAR(cellIds); + AXOM_UNUSED_VAR(statistics); return false; } @@ -293,19 +315,27 @@ class MeshClipperStrategy * done so far. Clip volumes computed by this method should * be added to the current values in this array. * + * @param [out] statistics Optional statistics to record + * consisting of child nodes with integer values. + * * @param [in] tetIds Indices of tets to clip, referring to the * shapeMesh.getCellsAsTets() array. tetIds[i] is the * \c (tetIds[i]%NUM_TETS_PER_HEX)-th tetrahedron of cell * \c = \c tetIds[i]/NUM_TETS_PER_HEX. Its overlap volume should * be added to \c ovlap[c]. + * + * Setting the statistics is not required except for getting + * accurate statistics. */ virtual bool specializedClipTets(quest::experimental::ShapeMesh& shapeMesh, axom::ArrayView ovlap, - const axom::ArrayView& tetIds) + const axom::ArrayView& tetIds, + conduit::Node& statistics) { AXOM_UNUSED_VAR(shapeMesh); AXOM_UNUSED_VAR(ovlap); AXOM_UNUSED_VAR(tetIds); + AXOM_UNUSED_VAR(statistics); return false; } @@ -313,8 +343,8 @@ class MeshClipperStrategy * @brief Get the geometry as discrete tetrahedra, or return false. * * @param [in] shapeMesh Blueprint mesh to shape into. - * @param [out] tets Array of tetrahedra filling the space of the shape, - * fully transformed. + * @param [out] tets Array of non-degenerate tetrahedra filling the + * space of the shape, fully transformed. * * Subclasses implementing this routine should snap to zero any * output vertex coordinate that is close to zero. @@ -338,8 +368,8 @@ class MeshClipperStrategy * @brief Get the geometry as discrete octahedra, or return false. * * @param [in] shapeMesh Blueprint mesh to shape into. - * @param [out] octs Array of octahedra filling the space of the shape, - * fully transformed. + * @param [out] tets Array of non-degenerate octahedra filling the + * space of the shape, fully transformed. * * Subclasses implementing this routine should snap to zero any * output vertex coordinate that is close to zero. diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index fe5f0c1cf3..f9a9e290af 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -15,10 +15,8 @@ #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 diff --git a/src/axom/quest/ShapeMesh.cpp b/src/axom/quest/ShapeMesh.cpp index 7dffd70eae..4a9cb10e97 100644 --- a/src/axom/quest/ShapeMesh.cpp +++ b/src/axom/quest/ShapeMesh.cpp @@ -184,6 +184,7 @@ void ShapeMesh::precomputeMeshData() getCellsAsHexes(); getCellsAsTets(); getCellVolumes(); + getTetVolumes(); getCellBoundingBoxes(); getCellLengths(); getCellNodeConnectivity(); @@ -217,6 +218,15 @@ axom::ArrayView ShapeMesh::getCellVolumes() return m_hexVolumes.view(); } +axom::ArrayView ShapeMesh::getTetVolumes() +{ + if(m_tetVolumes.size() != m_cellCount * NUM_TETS_PER_HEX) + { + computeTetVolumes(); + } + return m_tetVolumes.view(); +} + axom::ArrayView ShapeMesh::getCellBoundingBoxes() { if(m_hexBbs.size() != m_cellCount) @@ -624,6 +634,34 @@ void ShapeMesh::computeHexVolumes() } } +void ShapeMesh::computeTetVolumes() +{ + AXOM_ANNOTATE_SCOPE("ShapeMesh::computeTetVolumes"); + switch(m_runtimePolicy) + { + case RuntimePolicy::seq: + computeTetVolumesImpl(); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case RuntimePolicy::omp: + computeTetVolumesImpl(); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case RuntimePolicy::cuda: + computeTetVolumesImpl>(); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case RuntimePolicy::hip: + computeTetVolumesImpl>(); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } +} + void ShapeMesh::computeHexBbs() { AXOM_ANNOTATE_SCOPE("ShapeMesh::computeHexBoundingBoxes"); @@ -734,10 +772,10 @@ void ShapeMesh::computeCellsAsHexesImpl() SLIC_ASSERT(m_dim == NDIM); // or we shouldn't be here. - auto vertexCoords = getVertexCoords3D(); - const axom::ArrayView& vX = vertexCoords[0]; - const axom::ArrayView& vY = vertexCoords[1]; - const axom::ArrayView& vZ = vertexCoords[2]; + const auto& vertexCoords = getVertexCoords3D(); + const auto& vX = vertexCoords[0]; + const auto& vY = vertexCoords[1]; + const auto& vZ = vertexCoords[2]; axom::ArrayView connView = getCellNodeConnectivity(); @@ -750,7 +788,6 @@ void ShapeMesh::computeCellsAsHexesImpl() axom::for_all( m_cellCount, AXOM_LAMBDA(axom::IndexType cellId) { - // Set each hexahedral element vertices auto& hex = cellsAsHexesView[cellId]; for(int vi = 0; vi < NUM_VERTS_PER_HEX; ++vi) @@ -791,7 +828,7 @@ void ShapeMesh::computeCellsAsTetsImpl() AXOM_LAMBDA(axom::IndexType cellId) { const auto& hex = cellsAsHexesView[cellId]; auto* firstTetPtr = &cellsAsTetsView[cellId * NUM_TETS_PER_HEX]; - hex.triangulate(firstTetPtr); + hexToTets(hex, firstTetPtr); }); } @@ -809,6 +846,20 @@ void ShapeMesh::computeHexVolumesImpl() AXOM_LAMBDA(axom::IndexType i) { hexVolumesView[i] = cellsAsHexes[i].volume(); }); } +template +void ShapeMesh::computeTetVolumesImpl() +{ + axom::IndexType tetCount = m_cellCount * NUM_TETS_PER_HEX; + m_tetVolumes = axom::Array(ArrayOptions::Uninitialized(), tetCount, tetCount, m_allocId); + + auto cellsAsTets = getCellsAsTets(); + + auto tetVolumesView = m_tetVolumes.view(); + axom::for_all( + tetCount, + AXOM_LAMBDA(axom::IndexType i) { tetVolumesView[i] = cellsAsTets[i].volume(); }); +} + template void ShapeMesh::computeHexBbsImpl() { diff --git a/src/axom/quest/ShapeMesh.hpp b/src/axom/quest/ShapeMesh.hpp index b57ced3d15..63faa25dd9 100644 --- a/src/axom/quest/ShapeMesh.hpp +++ b/src/axom/quest/ShapeMesh.hpp @@ -58,11 +58,23 @@ class ShapeMesh using Point3DType = primal::Point; using TetrahedronType = primal::Tetrahedron; using HexahedronType = primal::Hexahedron; + using Plane3DType = axom::primal::Plane; using BoundingBox3DType = primal::BoundingBox; - //!@brief Number of tetrahedra per hexahedron decomposes into - // @internal We could use a more efficient 18-tet decomposition in the future. - static constexpr axom::IndexType NUM_TETS_PER_HEX = HexahedronType::NUM_TRIANGULATE; + /*! + * @brief Number of tetrahedra per hexahedron decomposes into + * @see hexToTets() + * + * @internal Values of 24 and 18 are valid. 18 is likely more + * performant because it generates fewer tets. + * + * @internal The code branches on the value at a low level, + * but this should be optimized out by the compiler. + */ + static constexpr axom::IndexType NUM_TETS_PER_HEX = 18; + + //!@brief Number of vertices per cell. + static constexpr axom::IndexType NUM_VERTS_PER_CELL_3D = 8; /*! * @brief Constructor with computational mesh in a conduit::Node. @@ -153,6 +165,21 @@ class ShapeMesh // zero to zero. Default threshold is 1e-10. void setZeroThreshold(double threshold) { m_zeroThreshold = threshold; } + /*! + * @brief Decompose a hexahedron into NUM_TETS_PER_HEX tetrahedra. + * @param hex [in] The hexahedron + * @param tets [out] Pointer to space for NUM_TETS_PER_HEX tetrahedra. + * + * To avoid ambiguity due to the choice of 2 diagonals for dividing + * each hex face into 2 triangles, we introduce a face-centered + * point at the average of the face vertices and decompose the face + * into 4 triangles. + * + * It is expected that this method will be used in long inner + * loops, so it is bare-bones for best performance. + */ + AXOM_HOST_DEVICE inline static void hexToTets(const HexahedronType& hex, TetrahedronType* tets); + //@{ //!@name Accessors to mesh data. //@} @@ -172,6 +199,8 @@ class ShapeMesh axom::ArrayView getCellsAsHexes(); //!@brief Get volume of mesh cells. axom::ArrayView getCellVolumes(); + //!@brief Get volumes of tets in getCellsAsTets(). + axom::ArrayView getTetVolumes(); //!@brief Get characteristic lengths of mesh cells. axom::ArrayView getCellLengths(); axom::ArrayView getCellBoundingBoxes(); @@ -276,6 +305,9 @@ class ShapeMesh //!@brief Volumes of hex cells. axom::Array m_hexVolumes; + //!@brief Volumes of cell tets. + axom::Array m_tetVolumes; + //!@brief Characteristic lengths of cells. axom::Array m_cellLengths; @@ -288,6 +320,7 @@ class ShapeMesh void computeCellsAsHexes(); void computeCellsAsTets(); void computeHexVolumes(); + void computeTetVolumes(); void computeHexBbs(); void computeCellLengths(); void computeVertPoints(); @@ -306,6 +339,9 @@ class ShapeMesh template void computeHexVolumesImpl(); + template + void computeTetVolumesImpl(); + template void computeHexBbsImpl(); @@ -345,6 +381,118 @@ class ShapeMesh const conduit::DataType& dtype); }; +AXOM_HOST_DEVICE inline void ShapeMesh::hexToTets(const HexahedronType& hex, TetrahedronType* tets) +{ + AXOM_STATIC_ASSERT(NUM_TETS_PER_HEX == 24 || NUM_TETS_PER_HEX == 18); + + if(NUM_TETS_PER_HEX == 24) + { + hex.triangulate(tets); + } + else + { + // Tets sharing the axis between hex vertices 4 and 2. + tets[0][0] = hex[4]; + tets[0][1] = hex[2]; + tets[0][2] = hex[1]; + tets[0][3] = hex[0]; + + tets[1][0] = hex[4]; + tets[1][1] = hex[2]; + tets[1][2] = hex[0]; + tets[1][3] = hex[3]; + + tets[2][0] = hex[4]; + tets[2][1] = hex[2]; + tets[2][2] = hex[3]; + tets[2][3] = hex[7]; + + tets[3][0] = hex[4]; + tets[3][1] = hex[2]; + tets[3][2] = hex[7]; + tets[3][3] = hex[6]; + + tets[4][0] = hex[4]; + tets[4][1] = hex[2]; + tets[4][2] = hex[6]; + tets[4][3] = hex[5]; + + tets[5][0] = hex[4]; + tets[5][1] = hex[2]; + tets[5][2] = hex[5]; + tets[5][3] = hex[1]; + + // Centroids of the 6 faces. + Point3DType mp0473 = Point3DType::midpoint(Point3DType::midpoint(hex[0], hex[4]), + Point3DType::midpoint(hex[7], hex[3])); + Point3DType mp1562 = Point3DType::midpoint(Point3DType::midpoint(hex[1], hex[5]), + Point3DType::midpoint(hex[6], hex[2])); + Point3DType mp0451 = Point3DType::midpoint(Point3DType::midpoint(hex[0], hex[4]), + Point3DType::midpoint(hex[5], hex[1])); + Point3DType mp3762 = Point3DType::midpoint(Point3DType::midpoint(hex[3], hex[7]), + Point3DType::midpoint(hex[6], hex[2])); + Point3DType mp0123 = Point3DType::midpoint(Point3DType::midpoint(hex[0], hex[1]), + Point3DType::midpoint(hex[2], hex[3])); + Point3DType mp4567 = Point3DType::midpoint(Point3DType::midpoint(hex[4], hex[5]), + Point3DType::midpoint(hex[6], hex[7])); + + // Tets from the 6 faces (two per face). + tets[6][0] = hex[4]; + tets[6][1] = hex[6]; + tets[6][2] = hex[7]; + tets[6][3] = mp4567; + tets[7][0] = hex[4]; + tets[7][1] = hex[5]; + tets[7][2] = hex[6]; + tets[7][3] = mp4567; + + tets[8][0] = hex[0]; + tets[8][1] = hex[2]; + tets[8][2] = hex[3]; + tets[8][3] = mp0123; + tets[9][0] = hex[0]; + tets[9][1] = hex[1]; + tets[9][2] = hex[2]; + tets[9][3] = mp0123; + + tets[10][0] = hex[4]; + tets[10][1] = hex[0]; + tets[10][2] = hex[3]; + tets[10][3] = mp0473; + tets[11][0] = hex[4]; + tets[11][1] = hex[3]; + tets[11][2] = hex[7]; + tets[11][3] = mp0473; + + tets[12][0] = hex[5]; + tets[12][1] = hex[1]; + tets[12][2] = hex[2]; + tets[12][3] = mp1562; + tets[13][0] = hex[5]; + tets[13][1] = hex[1]; + tets[13][2] = hex[2]; + tets[13][3] = mp1562; + + tets[14][0] = hex[4]; + tets[14][1] = hex[5]; + tets[14][2] = hex[1]; + tets[14][3] = mp0451; + tets[15][0] = hex[4]; + tets[15][1] = hex[1]; + tets[15][2] = hex[0]; + tets[15][3] = mp0451; + + tets[16][0] = hex[7]; + tets[16][1] = hex[6]; + tets[16][2] = hex[2]; + tets[16][3] = mp3762; + tets[17][0] = hex[7]; + tets[17][1] = hex[6]; + tets[17][2] = hex[3]; + tets[17][3] = mp3762; + } +} + } // namespace experimental } // namespace quest } // namespace axom diff --git a/src/axom/quest/detail/Discretize_detail.hpp b/src/axom/quest/detail/Discretize_detail.hpp index a1ea0c6815..587e1e18bb 100644 --- a/src/axom/quest/detail/Discretize_detail.hpp +++ b/src/axom/quest/detail/Discretize_detail.hpp @@ -265,7 +265,11 @@ bool discretize(const axom::ArrayView &polyline, axom::Array &out, int &octcount) { - int allocId = axom::execution_space::allocatorID(); + SLIC_ERROR_IF(!axom::execution_space::usesAllocId(out.getAllocatorID()), + axom::fmt::format("Execution space {} cannot access allocator id {}", + axom::execution_space::name(), + out.getAllocatorID())); + // Check for invalid input. If any segment is invalid, exit returning false. bool stillValid = true; int segmentcount = pointcount - 1; @@ -293,7 +297,8 @@ bool discretize(const axom::ArrayView &polyline, // That was the octahedron count for one segment. Multiply by the number // of segments we will compute. int totaloctcount = segoctcount * segmentcount; - out = axom::Array(totaloctcount, totaloctcount, allocId); + out.empty(); + out.resize(axom::ArrayOptions::Uninitialized(), totaloctcount); axom::ArrayView out_view = out.view(); octcount = 0; diff --git a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp index fd352aa602..6c32b117bf 100644 --- a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp +++ b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp @@ -6,13 +6,17 @@ #ifndef AXOM_MESHCLIPPERIMPL_HPP_ #define AXOM_MESHCLIPPERIMPL_HPP_ +#include "axom/config.hpp" + #ifndef AXOM_USE_RAJA #error "quest::MeshClipper requires RAJA." #endif +#include "axom/core/numerics/matvecops.hpp" #include "axom/quest/MeshClipperStrategy.hpp" #include "axom/quest/MeshClipper.hpp" #include "axom/spin/BVH.hpp" +#include "axom/primal/geometry/CoordinateTransformer.hpp" #include "RAJA/RAJA.hpp" namespace axom @@ -39,6 +43,10 @@ class MeshClipperImpl : public MeshClipper::Impl { public: using LabelType = MeshClipper::LabelType; + using Point3DType = primal::Point; + using BoundingBoxType = primal::BoundingBox; + using TetrahedronType = primal::Tetrahedron; + using OctahedronType = primal::Octahedron; MeshClipperImpl(MeshClipper& clipper) : MeshClipper::Impl(clipper) { } @@ -109,7 +117,7 @@ class MeshClipperImpl : public MeshClipper::Impl return; }; - AXOM_ANNOTATE_SCOPE("MeshClipper::collect_unlabeleds"); + AXOM_ANNOTATE_SCOPE("MeshClipper:collect_unlabeleds"); /*! * 1. Generate tmpLabels, having a value of 1 where labels is LABEL_ON and zero elsewhere. * 2. Inclusive scan on tmpLabels to generate values that step up at LABEL_ON cells. @@ -190,93 +198,35 @@ class MeshClipperImpl : public MeshClipper::Impl * and modified to work both tet and oct representations of the * geometry. */ - void computeClipVolumes3D(axom::ArrayView ovlap) override + void computeClipVolumes3D(axom::ArrayView ovlap, conduit::Node& statistics) override { - AXOM_ANNOTATE_SCOPE("MeshClipper::computeClipVolumes3D"); - - using BoundingBoxType = primal::BoundingBox; + using ATOMIC_POL = typename axom::execution_space::atomic_policy; ShapeMesh& shapeMesh = getShapeMesh(); - const int allocId = shapeMesh.getAllocatorID(); - const IndexType cellCount = shapeMesh.getCellCount(); - SLIC_INFO(axom::fmt::format( - "MeshClipper::computeClipVolumes3D: Getting discrete geometry for shape '{}'", - getStrategy().name())); - - // - // Get the geometry in discrete pieces, which can be tets or octs. - // - auto& strategy = getStrategy(); + /* + * Geometry as discrete tets or octs, and their bounding boxes. + */ axom::Array> geomAsTets; axom::Array> geomAsOcts; - const bool useOcts = strategy.getGeometryAsOcts(shapeMesh, geomAsOcts); - const bool useTets = strategy.getGeometryAsTets(shapeMesh, geomAsTets); - SLIC_ASSERT(useOcts || geomAsOcts.empty()); - SLIC_ASSERT(useTets || geomAsTets.empty()); - if(useTets == useOcts) - { - SLIC_ERROR( - axom::fmt::format("Problem with MeshClipperStrategy implementation '{}'." - " Implementations that don't provide a specializedClip function" - " must provide exactly one getGeometryAsOcts() or getGeometryAsTets()." - " This implementation provides {}.", - strategy.name(), - int(useOcts) + int(useTets))); - } - + axom::Array pieceBbs; + spin::BVH<3, ExecSpace, double> bvh; + bool useTets = getDiscreteGeometry(geomAsTets, geomAsOcts, pieceBbs, bvh); auto geomTetsView = geomAsTets.view(); auto geomOctsView = geomAsOcts.view(); - SLIC_INFO(axom::fmt::format("{:-^80}", " Inserting shapes' bounding boxes into BVH ")); - - // - // Generate the BVH over the (bounding boxes of the) discretized geometry - // - const axom::IndexType bbCount = useTets ? geomAsTets.size() : geomAsOcts.size(); - axom::Array pieceBbs(bbCount, bbCount, allocId); - axom::ArrayView pieceBbsView = pieceBbs.view(); - - // Get the bounding boxes for the shapes - if(useTets) - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomTetsView[i]); - }); - } - else - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomOctsView[i]); - }); - } - - spin::BVH<3, ExecSpace, double> bvh; - bvh.initialize(pieceBbsView, pieceBbsView.size()); - - SLIC_INFO(axom::fmt::format("{:-^80}", " Querying the BVH tree ")); - - axom::ArrayView cellBbsView = shapeMesh.getCellBoundingBoxes(); - - // Find which shape bounding boxes intersect hexahedron bounding boxes - SLIC_INFO( - axom::fmt::format("{:-^80}", " Finding shape candidates for each hexahedral element ")); + /* + * Find which shape bounding boxes intersect hexahedron bounding boxes + */ - axom::Array offsets(cellCount, cellCount, allocId); + AXOM_ANNOTATE_BEGIN("MeshClipper:find_candidates"); axom::Array counts(cellCount, cellCount, allocId); + axom::Array offsets(cellCount, cellCount, allocId); axom::Array candidates; - AXOM_ANNOTATE_BEGIN("bvh.findBoundingBoxes"); - bvh.findBoundingBoxes(offsets, counts, candidates, cellCount, cellBbsView); - AXOM_ANNOTATE_END("bvh.findBoundingBoxes"); - - // Get the total number of candidates - using ATOMIC_POL = typename axom::execution_space::atomic_policy; + bvh.findBoundingBoxes(offsets, counts, candidates, cellCount, shapeMesh.getCellBoundingBoxes()); + AXOM_ANNOTATE_END("MeshClipper:find_candidates"); const auto countsView = counts.view(); const int candidateCount = candidates.size(); @@ -293,8 +243,9 @@ class MeshClipperImpl : public MeshClipper::Impl allocId); auto shapeCandidatesView = shapeCandidates.view(); - // Tetrahedra from hexes (24 for each hex) + // Tetrahedra from hexes auto cellsAsTets = shapeMesh.getCellsAsTets(); + auto meshTetVolumes = getShapeMesh().getTetVolumes(); // Index into 'tets' axom::Array tetIndices(candidateCount * NUM_TETS_PER_HEX, @@ -337,77 +288,109 @@ class MeshClipperImpl : public MeshClipper::Impl }); } - constexpr double EPS = 1e-10; - constexpr bool tryFixOrientation = false; - - { - tetCandidatesCount = NUM_TETS_PER_HEX * candidates.size(); - AXOM_ANNOTATE_SCOPE("MeshClipper::clipLoop"); + tetCandidatesCount = NUM_TETS_PER_HEX * candidates.size(); #if defined(AXOM_DEBUG) - // Verifying: this should always pass. - if(tetCandidatesCountPtr != &tetCandidatesCount) - { - axom::copy(&tetCandidatesCount, tetCandidatesCountPtr, sizeof(IndexType)); - } - SLIC_ASSERT(tetCandidatesCount == candidateCount * NUM_TETS_PER_HEX); + // Verifying: this should always pass. + if(tetCandidatesCountPtr != &tetCandidatesCount) + { + axom::copy(&tetCandidatesCount, tetCandidatesCountPtr, sizeof(IndexType)); + } + SLIC_ASSERT(tetCandidatesCount == candidateCount * NUM_TETS_PER_HEX); #endif - SLIC_INFO( - axom::fmt::format("Running clip loop on {} candidate tets for of all {} hexes in the mesh", - tetCandidatesCount, - cellCount)); + SLIC_DEBUG( + axom::fmt::format("Running clip loop on {} candidate pieces for of all {} hexes in the mesh", + tetCandidatesCount, + cellCount)); - if(useTets) - { - axom::for_all( - tetCandidatesCount, - AXOM_LAMBDA(axom::IndexType i) { - const int index = hexIndicesView[i]; - const int shapeIndex = shapeCandidatesView[i]; - const int tetIndex = tetIndicesView[i]; - - const auto poly = primal::clip(geomTetsView[shapeIndex], - cellsAsTets[tetIndex], - EPS, - tryFixOrientation); - - // Poly is valid - if(poly.numVertices() >= 4) - { - // Workaround - intermediate volume variable needed for - // CUDA Pro/E test case correctness - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - RAJA::atomicAdd(ovlap.data() + index, volume); - } - }); - } - else // useOcts - { - axom::for_all( - tetCandidatesCount, - AXOM_LAMBDA(axom::IndexType i) { - const int index = hexIndicesView[i]; - const int shapeIndex = shapeCandidatesView[i]; - const int tetIndex = tetIndicesView[i]; - - const auto poly = primal::clip(geomOctsView[shapeIndex], - cellsAsTets[tetIndex], - EPS, - tryFixOrientation); - - // Poly is valid - if(poly.numVertices() >= 4) - { - // Workaround - intermediate volume variable needed for - // CUDA Pro/E test case correctness - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - RAJA::atomicAdd(ovlap.data() + index, volume); - } - }); - } + /* + * Initialize statistics and copy some to allocId space for computing. + */ + + const std::int64_t zero = 0; + std::int64_t& clipCount = *(statistics["clip"] = zero).as_int64_ptr(); + std::int64_t& contribCount = *(statistics["contribs"] = zero).as_int64_ptr(); + statistics["candidate"].set_int64(candidateCount); + + std::int64_t* clipCountPtr = axom::allocate(1, allocId); + std::int64_t* contribCountPtr = axom::allocate(1, allocId); + axom::copy(clipCountPtr, &clipCount, sizeof(zero)); + axom::copy(contribCountPtr, &contribCount, sizeof(zero)); + + const auto screenLevel = m_myClipper.getScreenLevel(); + constexpr double EPS1 = EPS; + + AXOM_ANNOTATE_BEGIN("MeshClipper:clipLoop_notScreened"); + if(useTets) + { + axom::for_all( + tetCandidatesCount, + AXOM_LAMBDA(axom::IndexType i) { + const int tetIndex = tetIndicesView[i]; + + // Skip degenerate mesh tets. + // Tet screening can filter out degenerate tets, but this method + // assumes no tet screening. + if(axom::utilities::isNearlyEqual(meshTetVolumes[tetIndex], 0.0, 1e-10)) + { + return; + } + + const int cellId = hexIndicesView[i]; + const int pieceId = shapeCandidatesView[i]; + const auto& cellTet = cellsAsTets[tetIndex]; + const TetrahedronType& geomPiece = geomTetsView[pieceId]; + + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(cellTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS1)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); + } + }); } + else // useOcts + { + axom::for_all( + tetCandidatesCount, + AXOM_LAMBDA(axom::IndexType i) { + const int tetIndex = tetIndicesView[i]; + + // Skip degenerate mesh tets. + if(axom::utilities::isNearlyEqual(meshTetVolumes[tetIndex], 0.0, 1e-10)) + { + return; + } + + const int cellId = hexIndicesView[i]; + const auto& cellTet = cellsAsTets[tetIndex]; + const int pieceId = shapeCandidatesView[i]; + const OctahedronType& geomPiece = geomOctsView[pieceId]; + + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(cellTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS1)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); + } + }); + } + AXOM_ANNOTATE_END("MeshClipper:clipLoop_notScreened"); + clipCount = tetCandidatesCount; + axom::copy(&contribCount, contribCountPtr, sizeof(contribCount)); + axom::deallocate(contribCountPtr); if(tetCandidatesCountPtr != &tetCandidatesCount) { @@ -422,99 +405,48 @@ class MeshClipperImpl : public MeshClipper::Impl * on the boundary. */ void computeClipVolumes3D(const axom::ArrayView& cellIndices, - axom::ArrayView ovlap) override + axom::ArrayView ovlap, + conduit::Node& statistics) override { - AXOM_ANNOTATE_SCOPE("MeshClipper::computeClipVolumes3D:limited"); - - using BoundingBoxType = primal::BoundingBox; + using ATOMIC_POL = typename axom::execution_space::atomic_policy; ShapeMesh& shapeMesh = getShapeMesh(); - const int allocId = shapeMesh.getAllocatorID(); - const IndexType cellCount = shapeMesh.getCellCount(); - SLIC_INFO(axom::fmt::format( - "MeshClipper::computeClipVolumes3D: Getting discrete geometry for shape '{}'", - getStrategy().name())); - - auto& strategy = getStrategy(); + /* + * Geometry as discrete tets or octs, and their bounding boxes. + */ axom::Array> geomAsTets; axom::Array> geomAsOcts; - const bool useOcts = strategy.getGeometryAsOcts(shapeMesh, geomAsOcts); - const bool useTets = strategy.getGeometryAsTets(shapeMesh, geomAsTets); - SLIC_ASSERT(useOcts || geomAsOcts.empty()); - SLIC_ASSERT(useTets || geomAsTets.empty()); - if(useTets == useOcts) - { - SLIC_ERROR( - axom::fmt::format("Problem with MeshClipperStrategy implementation '{}'." - " Implementations that don't provide a specializedClip function" - " must provide exactly one getGeometryAsOcts() or getGeometryAsTets()." - " This implementation provides {}.", - strategy.name(), - int(useOcts) + int(useTets))); - } - + axom::Array pieceBbs; + spin::BVH<3, ExecSpace, double> bvh; + bool useTets = getDiscreteGeometry(geomAsTets, geomAsOcts, pieceBbs, bvh); auto geomTetsView = geomAsTets.view(); auto geomOctsView = geomAsOcts.view(); - SLIC_INFO(axom::fmt::format("{:-^80}", " Inserting shapes' bounding boxes into BVH ")); - - // Generate the BVH tree over the shape's discretized geometry - // axis-aligned bounding boxes. "pieces" refers to tets or octs. - const axom::IndexType bbCount = useTets ? geomAsTets.size() : geomAsOcts.size(); - axom::Array pieceBbs(bbCount, bbCount, allocId); - axom::ArrayView pieceBbsView = pieceBbs.view(); - - // Get the bounding boxes for the shapes - if(useTets) - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomTetsView[i]); - }); - } - else - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomOctsView[i]); - }); - } - - // Insert shapes' Bounding Boxes into BVH. - spin::BVH<3, ExecSpace, double> bvh; - bvh.initialize(pieceBbsView, pieceBbsView.size()); - - SLIC_INFO(axom::fmt::format("{:-^80}", " Querying the BVH tree ")); + /* + * Find which shape bounding boxes intersect hexahedron bounding boxes + */ + AXOM_ANNOTATE_BEGIN("MeshClipper:find_candidates"); + // Find which shape bounding boxes intersect hexahedron bounding boxes // Create a temporary subset of cell bounding boxes, // containing only those listed in cellIndices. - const axom::IndexType limitedCellCount = cellIndices.size(); axom::ArrayView cellBbsView = shapeMesh.getCellBoundingBoxes(); + const axom::IndexType limitedCellCount = cellIndices.size(); axom::Array limitedCellBbs(limitedCellCount, limitedCellCount, allocId); axom::ArrayView limitedCellBbsView = limitedCellBbs.view(); axom::for_all( limitedCellBbsView.size(), AXOM_LAMBDA(axom::IndexType i) { limitedCellBbsView[i] = cellBbsView[cellIndices[i]]; }); - // 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(limitedCellCount, limitedCellCount, allocId); axom::Array counts(limitedCellCount, limitedCellCount, allocId); + axom::Array offsets(limitedCellCount, limitedCellCount, allocId); axom::Array candidates; - AXOM_ANNOTATE_BEGIN("bvh.findBoundingBoxes"); bvh.findBoundingBoxes(offsets, counts, candidates, limitedCellCount, limitedCellBbsView); - AXOM_ANNOTATE_END("bvh.findBoundingBoxes"); - - // Get the total number of candidates - using ATOMIC_POL = typename axom::execution_space::atomic_policy; + AXOM_ANNOTATE_END("MeshClipper:find_candidates"); const auto countsView = counts.view(); const int candidateCount = candidates.size(); @@ -531,8 +463,9 @@ class MeshClipperImpl : public MeshClipper::Impl allocId); auto shapeCandidatesView = shapeCandidates.view(); - // Tetrahedrons from hexes (24 for each hex) + // Tetrahedrons from hexes auto cellsAsTets = shapeMesh.getCellsAsTets(); + auto meshTetVolumes = getShapeMesh().getTetVolumes(); // Index into 'tets' axom::Array tetIndices(candidateCount * NUM_TETS_PER_HEX, @@ -575,92 +508,124 @@ class MeshClipperImpl : public MeshClipper::Impl }); } - SLIC_INFO(axom::fmt::format( - "Running clip loop on {} candidate tets for the select {} hexes of the full {} cells", + SLIC_DEBUG(axom::fmt::format( + "Running clip loop on {} candidate pieces for the select {} hexes of the full {} mesh cells", tetCandidatesCount, cellIndices.size(), cellCount)); - constexpr double EPS = 1e-10; - constexpr bool tryFixOrientation = false; - - { - tetCandidatesCount = NUM_TETS_PER_HEX * candidates.size(); - AXOM_ANNOTATE_SCOPE("MeshClipper::clipLoop_limited"); + tetCandidatesCount = NUM_TETS_PER_HEX * candidates.size(); #if defined(AXOM_DEBUG) - // Verifying: this should always pass. - if(tetCandidatesCountPtr != &tetCandidatesCount) - { - axom::copy(&tetCandidatesCount, tetCandidatesCountPtr, sizeof(IndexType)); - } - SLIC_ASSERT(tetCandidatesCount == candidateCount * NUM_TETS_PER_HEX); + // Verifying: this should always pass. + if(tetCandidatesCountPtr != &tetCandidatesCount) + { + axom::copy(&tetCandidatesCount, tetCandidatesCountPtr, sizeof(IndexType)); + } + SLIC_ASSERT(tetCandidatesCount == candidateCount * NUM_TETS_PER_HEX); #endif - if(useTets) - { - axom::for_all( - tetCandidatesCount, - AXOM_LAMBDA(axom::IndexType i) { - int index = hexIndicesView[i]; // index into limited mesh hex array - index = cellIndices[index]; // Now, it indexes into the full hex array. - - const int shapeIndex = shapeCandidatesView[i]; // index into pieces array - int tetIndex = - tetIndicesView[i]; // index into BVH results, implicit because BVH results specify hexes, not tets. - int tetIndex1 = tetIndex / NUM_TETS_PER_HEX; - int tetIndex2 = tetIndex % NUM_TETS_PER_HEX; - tetIndex = cellIndices[tetIndex1] * NUM_TETS_PER_HEX + - tetIndex2; // Now it indexes into the full tets-from-hexes array. - - const auto poly = primal::clip(geomTetsView[shapeIndex], - cellsAsTets[tetIndex], - EPS, - tryFixOrientation); - - // Poly is valid - if(poly.numVertices() >= 4) - { - // Workaround - intermediate volume variable needed for - // CUDA Pro/E test case correctness - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - RAJA::atomicAdd(ovlap.data() + index, volume); - } - }); - } - else // useOcts - { - axom::for_all( - tetCandidatesCount, - AXOM_LAMBDA(axom::IndexType i) { - int index = hexIndicesView[i]; // index into limited mesh hex array - index = cellIndices[index]; // Now, it indexes into the full hex array. - - const int shapeIndex = shapeCandidatesView[i]; // index into pieces array - int tetIndex = - tetIndicesView[i]; // index into BVH results, implicit because BVH results specify hexes, not tets. - int tetIndex1 = tetIndex / NUM_TETS_PER_HEX; - int tetIndex2 = tetIndex % NUM_TETS_PER_HEX; - tetIndex = cellIndices[tetIndex1] * NUM_TETS_PER_HEX + - tetIndex2; // Now it indexes into the full tets-from-hexes array. - - const auto poly = primal::clip(geomOctsView[shapeIndex], - cellsAsTets[tetIndex], - EPS, - tryFixOrientation); - - // Poly is valid - if(poly.numVertices() >= 4) - { - // Workaround - intermediate volume variable needed for - // CUDA Pro/E test case correctness - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - RAJA::atomicAdd(ovlap.data() + index, volume); - } - }); - } + /* + * Initialize statistics and copy some to allocId space for computing. + */ + + const std::int64_t zero = 0; + std::int64_t& clipCount = *(statistics["clip"] = zero).as_int64_ptr(); + std::int64_t& contribCount = *(statistics["contribs"] = zero).as_int64_ptr(); + statistics["candidate"].set_int64(candidateCount); + + std::int64_t* clipCountPtr = axom::allocate(1, allocId); + std::int64_t* contribCountPtr = axom::allocate(1, allocId); + axom::copy(clipCountPtr, &clipCount, sizeof(zero)); + axom::copy(contribCountPtr, &contribCount, sizeof(zero)); + + const auto screenLevel = m_myClipper.getScreenLevel(); + + AXOM_ANNOTATE_BEGIN("MeshClipper:clipLoop_hexScreened"); + if(useTets) + { + axom::for_all( + tetCandidatesCount, + AXOM_LAMBDA(axom::IndexType i) { + int tetIndex = + tetIndicesView[i]; // index into BVH results, implicit because BVH results specify hexes, not tets. + int tetIndex1 = tetIndex / NUM_TETS_PER_HEX; + int tetIndex2 = tetIndex % NUM_TETS_PER_HEX; + tetIndex = cellIndices[tetIndex1] * NUM_TETS_PER_HEX + + tetIndex2; // Now it indexes into the full tets-from-hexes array. + + // Skip degenerate mesh tets. + // Tet screening can filter out degenerate tets, but this method + // assumes no tet screening. + if(axom::utilities::isNearlyEqual(meshTetVolumes[tetIndex], 0.0, 1e-10)) + { + return; + } + + int cellId = hexIndicesView[i]; // index into limited mesh hex array + cellId = cellIndices[cellId]; // Now, it indexes into the full hex array. + + const int pieceId = shapeCandidatesView[i]; // index into pieces array + const auto& cellTet = cellsAsTets[tetIndex]; + const TetrahedronType& geomPiece = geomTetsView[pieceId]; + + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(cellTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); + } + }); + } + else // useOcts + { + axom::for_all( + tetCandidatesCount, + AXOM_LAMBDA(axom::IndexType i) { + int tetIndex = + tetIndicesView[i]; // index into BVH results, implicit because BVH results specify hexes, not tets. + int tetIndex1 = tetIndex / NUM_TETS_PER_HEX; + int tetIndex2 = tetIndex % NUM_TETS_PER_HEX; + tetIndex = cellIndices[tetIndex1] * NUM_TETS_PER_HEX + + tetIndex2; // Now it indexes into the full tets-from-hexes array. + + // Skip degenerate mesh tets. + if(axom::utilities::isNearlyEqual(meshTetVolumes[tetIndex], 0.0, 1e-10)) + { + return; + } + + int cellId = hexIndicesView[i]; // index into limited mesh hex array + cellId = cellIndices[cellId]; // Now, it indexes into the full hex array. + + const int pieceId = shapeCandidatesView[i]; // index into pieces array + const OctahedronType& geomPiece = geomOctsView[pieceId]; + const auto& cellTet = cellsAsTets[tetIndex]; + + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(cellTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); + } + }); } + AXOM_ANNOTATE_END("MeshClipper:clipLoop_hexScreened"); + + clipCount = tetCandidatesCount; + axom::copy(&contribCount, contribCountPtr, sizeof(contribCount)); + axom::deallocate(contribCountPtr); if(tetCandidatesCountPtr != &tetCandidatesCount) { @@ -675,79 +640,36 @@ class MeshClipperImpl : public MeshClipper::Impl * potentially on the boundary. */ void computeClipVolumes3DTets(const axom::ArrayView& tetIndices, - axom::ArrayView ovlap) override + axom::ArrayView ovlap, + conduit::Node& statistics) override { - AXOM_ANNOTATE_SCOPE("MeshClipper::computeClipVolumes3D:limited"); - - using BoundingBoxType = primal::BoundingBox; - using TetrahedronType = primal::Tetrahedron; - using OctahedronType = primal::Octahedron; + using ATOMIC_POL = typename axom::execution_space::atomic_policy; ShapeMesh& shapeMesh = getShapeMesh(); auto meshTets = getShapeMesh().getCellsAsTets(); const int allocId = shapeMesh.getAllocatorID(); - SLIC_INFO(axom::fmt::format( - "MeshClipper::computeClipVolumes3D: Getting discrete geometry for shape '{}'", - getStrategy().name())); - - auto& strategy = getStrategy(); - axom::Array geomAsTets; - axom::Array geomAsOcts; - const bool useOcts = strategy.getGeometryAsOcts(shapeMesh, geomAsOcts); - const bool useTets = strategy.getGeometryAsTets(shapeMesh, geomAsTets); - SLIC_ASSERT(useOcts || geomAsOcts.empty()); - SLIC_ASSERT(useTets || geomAsTets.empty()); - if(useTets == useOcts) - { - SLIC_ERROR( - axom::fmt::format("Problem with MeshClipperStrategy implementation '{}'." - " Implementations that don't provide a specializedClip function" - " must provide exactly one getGeometryAsOcts() or getGeometryAsTets()." - " This implementation provides {}.", - strategy.name(), - int(useOcts) + int(useTets))); - } - + /* + * Geometry as discrete tets or octs, and their bounding boxes. + */ + axom::Array> geomAsTets; + axom::Array> geomAsOcts; + axom::Array pieceBbs; + spin::BVH<3, ExecSpace, double> bvh; + bool useTets = getDiscreteGeometry(geomAsTets, geomAsOcts, pieceBbs, bvh); auto geomTetsView = geomAsTets.view(); auto geomOctsView = geomAsOcts.view(); - SLIC_INFO(axom::fmt::format("{:-^80}", " Inserting shapes' bounding boxes into BVH ")); - - // Generate the BVH tree over the shape's discretized geometry - // axis-aligned bounding boxes. "pieces" refers to tets or octs. - const axom::IndexType bbCount = useTets ? geomAsTets.size() : geomAsOcts.size(); - axom::Array pieceBbs(bbCount, bbCount, allocId); - axom::ArrayView pieceBbsView = pieceBbs.view(); - - // Get the bounding boxes for the shapes - if(useTets) - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomTetsView[i]); - }); - } - else - { - axom::for_all( - pieceBbsView.size(), - AXOM_LAMBDA(axom::IndexType i) { - pieceBbsView[i] = primal::compute_bounding_box(geomOctsView[i]); - }); - } - - // Insert shapes' Bounding Boxes into BVH. - spin::BVH<3, ExecSpace, double> bvh; - bvh.initialize(pieceBbsView, pieceBbsView.size()); - - SLIC_INFO(axom::fmt::format("{:-^80}", " Querying the BVH tree ")); + /* + * Find which shape bounding boxes intersect hexahedron bounding boxes + */ + AXOM_ANNOTATE_BEGIN("MeshClipper:find_candidates"); // Create a temporary subset of tet bounding boxes, // containing only those listed in tetIndices. + // The BVH searches on this array. const axom::IndexType tetCount = tetIndices.size(); axom::Array tetBbs(tetCount, tetCount, allocId); axom::ArrayView tetBbsView = tetBbs.view(); @@ -762,15 +684,11 @@ class MeshClipperImpl : public MeshClipper::Impl axom::Array counts(tetCount, tetCount, allocId); axom::Array offsets(tetCount, tetCount, allocId); - + axom::Array candidates; auto countsView = counts.view(); auto offsetsView = offsets.view(); - - AXOM_ANNOTATE_BEGIN("bvh.findBoundingBoxes"); - axom::Array candidates; bvh.findBoundingBoxes(offsets, counts, candidates, tetBbsView.size(), tetBbsView); auto candidatesView = candidates.view(); - AXOM_ANNOTATE_END("bvh.findBoundingBoxes"); axom::Array candToTetIdId(candidates.size(), candidates.size(), allocId); auto candToTetIdIdView = candToTetIdId.view(); @@ -781,35 +699,56 @@ class MeshClipperImpl : public MeshClipper::Impl auto offset = offsetsView[i]; for(int j = 0; j < count; ++j) candToTetIdIdView[offset + j] = i; }); + AXOM_ANNOTATE_END("MeshClipper:find_candidates"); - // Find which shape bounding boxes intersect hexahedron bounding boxes - SLIC_INFO(axom::fmt::format("Finding shape candidates for {} tet elements ", tetIndices.size())); + SLIC_DEBUG(axom::fmt::format( + "Running clip loop on {} candidate pieces for the select {} tets of the full {} mesh cells", + candidates.size(), + tetCount, + shapeMesh.getCellCount())); - using ATOMIC_POL = typename axom::execution_space::atomic_policy; - constexpr double EPS = 1e-10; - constexpr bool tryFixOrientation = false; + /* + * Initialize statistics and copy some to allocId space for computing. + */ + + const std::int64_t zero = 0; + std::int64_t& clipCount = *(statistics["clip"] = zero).as_int64_ptr(); + std::int64_t& contribCount = *(statistics["contribs"] = zero).as_int64_ptr(); + std::int64_t& candidateCount = *(statistics["candidate"] = zero).as_int64_ptr(); + candidateCount = candidates.size(); + std::int64_t* clipCountPtr = axom::allocate(1, allocId); + std::int64_t* contribCountPtr = axom::allocate(1, allocId); + axom::copy(clipCountPtr, &clipCount, sizeof(zero)); + axom::copy(contribCountPtr, &contribCount, sizeof(zero)); + + const auto screenLevel = m_myClipper.getScreenLevel(); + + AXOM_ANNOTATE_BEGIN("MeshClipper:clipLoop_tetScreened"); if(useTets) { axom::for_all( candidates.size(), AXOM_LAMBDA(axom::IndexType iCand) { - auto pieceId = candidatesView[iCand]; - const axom::primal::Tetrahedron& geomPiece = geomTetsView[pieceId]; - auto tetIdId = candToTetIdIdView[iCand]; auto tetId = tetIndices[tetIdId]; - const auto& tet = meshTets[tetId]; - const auto poly = primal::clip(tet, geomPiece, EPS, tryFixOrientation); + auto cellId = tetId / NUM_TETS_PER_HEX; + auto pieceId = candidatesView[iCand]; + const auto& meshTet = meshTets[tetId]; + const TetrahedronType& geomPiece = geomTetsView[pieceId]; - if(poly.numVertices() >= 4) + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(meshTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) { - // Poly is valid - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - auto cellId = tetId / NUM_TETS_PER_HEX; RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); } }); } @@ -818,39 +757,266 @@ class MeshClipperImpl : public MeshClipper::Impl axom::for_all( candidates.size(), AXOM_LAMBDA(axom::IndexType iCand) { - auto pieceId = candidatesView[iCand]; - const axom::primal::Octahedron& geomPiece = geomOctsView[pieceId]; - auto tetIdId = candToTetIdIdView[iCand]; auto tetId = tetIndices[tetIdId]; - const auto& tet = meshTets[tetId]; - const auto poly = primal::clip(tet, geomPiece, EPS, tryFixOrientation); + auto cellId = tetId / NUM_TETS_PER_HEX; + auto pieceId = candidatesView[iCand]; + const auto& meshTet = meshTets[tetId]; + const OctahedronType& geomPiece = geomOctsView[pieceId]; - if(poly.numVertices() >= 4) + double volume = 0.0; + LabelType tmpLabel = + computeMeshTetGeomPieceOverlap(meshTet, geomPiece, volume, screenLevel); + if(tmpLabel == LabelType::LABEL_IN || tmpLabel == LabelType::LABEL_ON) { - // Poly is valid - double volume = poly.volume(); - SLIC_ASSERT(volume >= 0); - auto cellId = tetId / NUM_TETS_PER_HEX; RAJA::atomicAdd(ovlap.data() + cellId, volume); + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS)); + } + if(tmpLabel == LabelType::LABEL_ON) + { + RAJA::atomicAdd(clipCountPtr, std::int64_t(1)); } }); } + AXOM_ANNOTATE_END("MeshClipper:clipLoop_tetScreened"); + + axom::copy(&contribCount, contribCountPtr, sizeof(contribCount)); + axom::copy(&clipCount, clipCountPtr, sizeof(clipCount)); + axom::deallocate(contribCountPtr); + axom::deallocate(clipCountPtr); + SLIC_DEBUG(axom::fmt::format("")); } // end of computeClipVolumes3DTets() function + /*! + * @brief Get the geometry in discrete pieces, + * which can be tets or octs, and place them in a search tree. + * @return whether geometry are tetrahedra instead of octahedra. + */ + bool getDiscreteGeometry(axom::Array>& geomAsTets, + axom::Array>& geomAsOcts, + axom::Array& pieceBbs, + spin::BVH<3, ExecSpace, double>& bvh) + { + auto& strategy = getStrategy(); + ShapeMesh& shapeMesh = getShapeMesh(); + int allocId = shapeMesh.getAllocatorID(); + + AXOM_ANNOTATE_BEGIN("MeshClipper:get_geometry"); + const bool useOcts = strategy.getGeometryAsOcts(shapeMesh, geomAsOcts); + const bool useTets = strategy.getGeometryAsTets(shapeMesh, geomAsTets); + AXOM_ANNOTATE_END("MeshClipper:get_geometry"); + + if(useTets) + { + SLIC_ASSERT(geomAsTets.getAllocatorID() == allocId); + } + else + { + SLIC_ASSERT(geomAsOcts.getAllocatorID() == allocId); + } + if(useTets == useOcts) + { + SLIC_ERROR( + axom::fmt::format("Problem with MeshClipperStrategy implementation '{}'." + " Implementations that don't provide a specializedClip function" + " must provide exactly one of either getGeometryAsOcts() or" + " getGeometryAsTets(). This implementation provides {}.", + strategy.name(), + int(useOcts) + int(useTets))); + } + + SLIC_DEBUG(axom::fmt::format("Geometry {} has {} discrete {}s", + strategy.name(), + useTets ? geomAsTets.size() : geomAsOcts.size(), + useTets ? "tet" : "oct")); + + /* + * Get the bounding boxes for the discrete geometry pieces. + * If debug build, check for degenerate pieces. + */ + const axom::IndexType bbCount = useTets ? geomAsTets.size() : geomAsOcts.size(); + pieceBbs = axom::Array(bbCount, bbCount, allocId); + axom::ArrayView pieceBbsView = pieceBbs.view(); + + if(useTets) + { + auto geomTetsView = geomAsTets.view(); + axom::for_all( + pieceBbsView.size(), + AXOM_LAMBDA(axom::IndexType i) { + pieceBbsView[i] = primal::compute_bounding_box(geomTetsView[i]); +#if defined(AXOM_DEBUG) + SLIC_ASSERT(!geomTetsView[i].degenerate()); +#endif + }); + } + else + { + auto geomOctsView = geomAsOcts.view(); + axom::for_all( + pieceBbsView.size(), + AXOM_LAMBDA(axom::IndexType i) { + pieceBbsView[i] = primal::compute_bounding_box(geomOctsView[i]); + }); + } + + bvh.setAllocatorID(allocId); + bvh.setTolerance(EPS); + bvh.setScaleFactor(BVH_SCALE_FACTOR); + bvh.initialize(pieceBbsView, pieceBbsView.size()); + + return useTets; + } + + /*! + * @brief Volume of a tetrahedron from discretized geometry. + */ + AXOM_HOST_DEVICE inline double geomPieceVolume(const TetrahedronType& tet) + { + return tet.volume(); + } + + /*! + * @brief Volume of a octahedron from discretized geometry. + * + * Assumes octahedron is convex. + */ + AXOM_HOST_DEVICE inline double geomPieceVolume(const OctahedronType& oct) + { + TetrahedronType tets[] = {TetrahedronType(oct[0], oct[3], oct[1], oct[2]), + TetrahedronType(oct[0], oct[3], oct[2], oct[4]), + TetrahedronType(oct[0], oct[3], oct[4], oct[5]), + TetrahedronType(oct[0], oct[3], oct[5], oct[1])}; + double octVol = 0.0; + for(int i = 0; i < 4; ++i) + { + double tetVol = tets[i].volume(); + SLIC_ASSERT(tetVol >= -EPS); // Tet may be degenerate but not inverted. + octVol += tetVol; + } + return octVol; + } + + /*! + * @brief Compute overlap volume between a tet (from the shape mesh) + * and a piece (tet or oct) of the discretized geometry. + * + * Becaue primal::clip is so expensive, we do a conservative + * overlap check on @c meshTet and @c geomPiece to avoid clipping. + * + * @return results of check whether the piece is IN/ON/OUT of the tet. + * + * @tparam TetOrOctType either a TetrahedronType or OctahedronType, + * the two types a geometry can be discretized into. + */ + template + AXOM_HOST_DEVICE inline LabelType computeMeshTetGeomPieceOverlap(const TetrahedronType& meshTet, + const TetOrOctType& geomPiece, + double& overlapVolume, + int screenLevel) + { + constexpr bool tryFixOrientation = false; + if(screenLevel >= 3) + { + LabelType geomLabel = labelPieceInOutOfTet(meshTet, geomPiece); + if(geomLabel == LabelType::LABEL_OUT) + { + overlapVolume = 0.0; + return geomLabel; + } + if(geomLabel == LabelType::LABEL_IN) + { + overlapVolume = geomPieceVolume(geomPiece); + return geomLabel; + } + } + + const auto poly = primal::clip(meshTet, geomPiece, EPS, tryFixOrientation); + if(poly.numVertices() >= 4) + { + // Poly is valid + overlapVolume = poly.volume(); + SLIC_ASSERT(overlapVolume >= 0); + } + else + { + overlapVolume = 0.0; + } + + return LabelType::LABEL_ON; + } + + /*! + * @brief Compute whether a tetrahedron or octhedron is inside, + * outside or on the boundary of a reference tetrahedron, + * and conservatively label it as on, if not known. + * + * @internal To reduce repeatedly computing toUnitTet for + * the same tet, precompute that in the calling function + * and use it instead of tet. + */ + template + AXOM_HOST_DEVICE inline LabelType labelPieceInOutOfTet(const TetrahedronType& tet, + const TetOrOctType& piece) + { + Point3DType unitTet[] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; + axom::primal::experimental::CoordinateTransformer toUnitTet(&tet[0], unitTet); + + /* + * Count (transformed) piece vertices above/below unitTet as unitTet + * rests on its 4 sides. Sides 0-2 are perpendicular to the axes. + * Side 3 is the diagonal side. + */ + int vsAbove[4] = {0, 0, 0, 0}; + int vsBelow[4] = {0, 0, 0, 0}; + int vsTetSide[4] = {0, 0, 0, 0}; + for(int i = 0; i < TetOrOctType::NUM_VERTS; ++i) + { + auto pVert = toUnitTet.getTransformed(piece[i]); + // h4 is height of pVert above the diagonal face, scaled by sqrt(3). + // h4 of 1 coresponds to the unitTet's height of sqrt(3). + double h4 = 1 - (pVert[0] + pVert[1] + pVert[2]); + vsAbove[0] += pVert[0] >= 1; + vsAbove[1] += pVert[1] >= 1; + vsAbove[2] += pVert[2] >= 1; + vsAbove[3] += h4 >= 1; + vsBelow[0] += pVert[0] <= 0; + vsBelow[1] += pVert[1] <= 0; + vsBelow[2] += pVert[2] <= 0; + vsBelow[3] += h4 <= 0; + vsTetSide[0] += pVert[0] >= 0; + vsTetSide[1] += pVert[1] >= 0; + vsTetSide[2] += pVert[2] >= 0; + vsTetSide[3] += h4 >= 0; + } + if(vsAbove[0] == TetOrOctType::NUM_VERTS || vsAbove[1] == TetOrOctType::NUM_VERTS || + vsAbove[2] == TetOrOctType::NUM_VERTS || vsAbove[3] == TetOrOctType::NUM_VERTS || + vsBelow[0] == TetOrOctType::NUM_VERTS || vsBelow[1] == TetOrOctType::NUM_VERTS || + vsBelow[2] == TetOrOctType::NUM_VERTS || vsBelow[3] == TetOrOctType::NUM_VERTS) + { + return LabelType::LABEL_OUT; + } + if(vsTetSide[0] == TetOrOctType::NUM_VERTS && vsTetSide[1] == TetOrOctType::NUM_VERTS && + vsTetSide[2] == TetOrOctType::NUM_VERTS && vsTetSide[3] == TetOrOctType::NUM_VERTS) + { + return LabelType::LABEL_IN; + } + return LabelType::LABEL_ON; + } + void getLabelCounts(axom::ArrayView labels, - axom::IndexType& inCount, - axom::IndexType& onCount, - axom::IndexType& outCount) override + std::int64_t& inCount, + std::int64_t& onCount, + std::int64_t& outCount) override { - AXOM_ANNOTATE_SCOPE("MeshClipper::getLabelCounts"); + AXOM_ANNOTATE_SCOPE("MeshClipper:count_labels"); using ReducePolicy = typename axom::execution_space::reduce_policy; using LoopPolicy = typename execution_space::loop_policy; - RAJA::ReduceSum inSum(0); - RAJA::ReduceSum onSum(0); - RAJA::ReduceSum outSum(0); + RAJA::ReduceSum inSum(0); + RAJA::ReduceSum onSum(0); + RAJA::ReduceSum outSum(0); RAJA::forall( RAJA::RangeSegment(0, labels.size()), AXOM_LAMBDA(axom::IndexType cellId) { @@ -868,12 +1034,14 @@ class MeshClipperImpl : public MeshClipper::Impl onSum += 1; } }); - inCount = static_cast(inSum.get()); - onCount = static_cast(onSum.get()); - outCount = static_cast(outSum.get()); + inCount = static_cast(inSum.get()); + onCount = static_cast(onSum.get()); + outCount = static_cast(outSum.get()); } private: + static constexpr double EPS = 1e-10; + static constexpr double BVH_SCALE_FACTOR = 1.0; static constexpr int MAX_VERTS_FOR_TET_CLIPPING = 32; static constexpr int MAX_NBRS_PER_VERT_FOR_TET_CLIPPING = 8; static constexpr int MAX_VERTS_FOR_OCT_CLIPPING = 32; diff --git a/src/axom/quest/detail/clipping/Plane3DClipper.cpp b/src/axom/quest/detail/clipping/Plane3DClipper.cpp new file mode 100644 index 0000000000..28bd071cb8 --- /dev/null +++ b/src/axom/quest/detail/clipping/Plane3DClipper.cpp @@ -0,0 +1,415 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "axom/config.hpp" + +#include "axom/quest/detail/clipping/Plane3DClipper.hpp" + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +Plane3DClipper::Plane3DClipper(const klee::Geometry& kGeom, const std::string& name) + : MeshClipperStrategy(kGeom) + , m_name(name.empty() ? std::string("Plane3D") : name) +{ + extractClipperInfo(); +} + +bool Plane3DClipper::labelCellsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& labels) +{ + int allocId = shapeMesh.getAllocatorID(); + auto cellCount = shapeMesh.getCellCount(); + if(labels.size() < cellCount || labels.getAllocatorID() != shapeMesh.getAllocatorID()) + { + labels = axom::Array(ArrayOptions::Uninitialized(), cellCount, cellCount, allocId); + } + + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + labelCellsInOutImpl(shapeMesh, labels.view()); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + labelCellsInOutImpl(shapeMesh, labels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + labelCellsInOutImpl>(shapeMesh, labels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + labelCellsInOutImpl>(shapeMesh, labels.view()); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +bool Plane3DClipper::labelTetsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::Array& tetLabels) +{ + int allocId = shapeMesh.getAllocatorID(); + const auto cellCount = cellIds.size(); + const auto tetCount = cellCount * NUM_TETS_PER_HEX; + if(tetLabels.size() < tetCount || tetLabels.getAllocatorID() != allocId) + { + tetLabels = axom::Array(ArrayOptions::Uninitialized(), tetCount, tetCount, allocId); + } + + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + labelTetsInOutImpl(shapeMesh, cellIds, tetLabels.view()); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + labelTetsInOutImpl(shapeMesh, cellIds, tetLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + labelTetsInOutImpl>(shapeMesh, cellIds, tetLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + labelTetsInOutImpl>(shapeMesh, cellIds, tetLabels.view()); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +bool Plane3DClipper::specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics) +{ + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + specializedClipCellsImpl(shapeMesh, ovlap, statistics); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + specializedClipCellsImpl(shapeMesh, ovlap, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + specializedClipCellsImpl>(shapeMesh, ovlap, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + specializedClipCellsImpl>(shapeMesh, ovlap, statistics); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +bool Plane3DClipper::specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& cellIds, + conduit::Node& statistics) +{ + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + specializedClipCellsImpl(shapeMesh, ovlap, cellIds, statistics); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + specializedClipCellsImpl(shapeMesh, ovlap, cellIds, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + specializedClipCellsImpl>(shapeMesh, ovlap, cellIds, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + specializedClipCellsImpl>(shapeMesh, ovlap, cellIds, statistics); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +bool Plane3DClipper::specializedClipTets(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& tetIds, + conduit::Node& statistics) +{ + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + specializedClipTetsImpl(shapeMesh, ovlap, tetIds, statistics); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + specializedClipTetsImpl(shapeMesh, ovlap, tetIds, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + specializedClipTetsImpl>(shapeMesh, ovlap, tetIds, statistics); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + specializedClipTetsImpl>(shapeMesh, ovlap, tetIds, statistics); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +template +void Plane3DClipper::labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView labels) +{ + int allocId = shapeMesh.getAllocatorID(); + auto cellCount = shapeMesh.getCellCount(); + auto vertCount = shapeMesh.getVertexCount(); + auto cellVolumes = shapeMesh.getCellVolumes(); + constexpr double EPS = 1e-10; + + const auto& vertCoords = shapeMesh.getVertexCoords3D(); + const auto& vX = vertCoords[0]; + const auto& vY = vertCoords[1]; + const auto& vZ = vertCoords[2]; + + /* + Compute whether vertices are inside shape. + */ + axom::Array vertIsInside {ArrayOptions::Uninitialized(), vertCount, vertCount, allocId}; + auto vertIsInsideView = vertIsInside.view(); + SLIC_ASSERT(axom::execution_space::usesAllocId(vX.getAllocatorID())); + SLIC_ASSERT(axom::execution_space::usesAllocId(vY.getAllocatorID())); + SLIC_ASSERT(axom::execution_space::usesAllocId(vZ.getAllocatorID())); + SLIC_ASSERT(axom::execution_space::usesAllocId(vertIsInsideView.getAllocatorID())); + + auto plane = m_plane; + axom::for_all( + vertCount, + AXOM_LAMBDA(axom::IndexType vertId) { + primal::Point3D vert {vX[vertId], vY[vertId], vZ[vertId]}; + double signedDist = plane.signedDistance(vert); + vertIsInsideView[vertId] = signedDist > 0; + }); + + /* + * Label cell by whether it has vertices inside, outside or both. + */ + axom::ArrayView connView = shapeMesh.getCellNodeConnectivity(); + SLIC_ASSERT(connView.shape()[1] == NUM_VERTS_PER_CELL_3D); + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellId) { + if(axom::utilities::isNearlyEqual(cellVolumes[cellId], 0.0, EPS)) + { + labels[cellId] = LabelType::LABEL_OUT; + return; + } + auto cellVertIds = connView[cellId]; + bool hasIn = vertIsInsideView[cellVertIds[0]]; + bool hasOut = !hasIn; + for(int vi = 0; vi < NUM_VERTS_PER_CELL_3D; ++vi) + { + int vertId = cellVertIds[vi]; + bool isIn = vertIsInsideView[vertId]; + hasIn |= isIn; + hasOut |= !isIn; + } + labels[cellId] = !hasOut ? LabelType::LABEL_IN + : !hasIn ? LabelType::LABEL_OUT + : LabelType::LABEL_ON; + }); + + return; +} + +template +void Plane3DClipper::labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::ArrayView tetLabels) +{ + auto cellCount = cellIds.size(); + auto meshTets = shapeMesh.getCellsAsTets(); + auto meshTetVolumes = shapeMesh.getTetVolumes(); + + auto plane = m_plane; + + constexpr double EPS = 1e-10; + + /* + * Label tet by whether it has vertices inside, outside or both. + * Degenerate tets as outside, because they contribute no volume. + */ + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType ci) { + axom::IndexType cellId = cellIds[ci]; + + const TetrahedronType* tetsForCell = &meshTets[cellId * NUM_TETS_PER_HEX]; + const double* tetVolumesForCell = &meshTetVolumes[cellId * NUM_TETS_PER_HEX]; + + for(IndexType ti = 0; ti < NUM_TETS_PER_HEX; ++ti) + { + const auto& tet = tetsForCell[ti]; + LabelType& tetLabel = tetLabels[ci * NUM_TETS_PER_HEX + ti]; + + if(axom::utilities::isNearlyEqual(tetVolumesForCell[ti], 0.0, EPS)) + { + tetLabel = LabelType::LABEL_OUT; + continue; + } + + bool hasIn = false; + bool hasOut = false; + for(int vi = 0; vi < TetrahedronType::NUM_VERTS; ++vi) + { + const auto& vert = tet[vi]; + double signedDist = plane.signedDistance(vert); + hasIn |= signedDist > 0; + hasOut |= signedDist < 0; + } + tetLabel = !hasOut ? LabelType::LABEL_IN + : !hasIn ? LabelType::LABEL_OUT + : LabelType::LABEL_ON; + } + }); + + return; +} + +template +void Plane3DClipper::specializedClipCellsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics) +{ + axom::IndexType cellCount = shapeMesh.getCellCount(); + axom::Array cellIds(cellCount, cellCount, shapeMesh.getAllocatorID()); + auto cellIdsView = cellIds.view(); + axom::for_all(cellCount, AXOM_LAMBDA(axom::IndexType i) { cellIdsView[i] = i; }); + specializedClipCellsImpl(shapeMesh, ovlap, cellIds, statistics); +} + +template +void Plane3DClipper::specializedClipCellsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& cellIds, + conduit::Node& statistics) +{ + using ATOMIC_POL = typename axom::execution_space::atomic_policy; + constexpr double EPS = 1e-10; + + int allocId = shapeMesh.getAllocatorID(); + + auto cellsAsTets = shapeMesh.getCellsAsTets(); + + auto plane = m_plane; + + const std::int64_t zero = 0; + std::int64_t& contribCount = *(statistics["contribs"] = zero).as_int64_ptr(); + std::int64_t* contribCountPtr = axom::allocate(1, allocId); + axom::copy(contribCountPtr, &contribCount, sizeof(zero)); + + axom::for_all( + cellIds.size(), + AXOM_LAMBDA(axom::IndexType i) { + axom::IndexType cellId = cellIds[i]; + const TetrahedronType* tetsInHex = cellsAsTets.data() + cellId * NUM_TETS_PER_HEX; + double vol = 0.0; + for(int ti = 0; ti < NUM_TETS_PER_HEX; ++ti) + { + const auto& tet = tetsInHex[ti]; + primal::Polyhedron overlap = primal::clip(tet, plane, EPS); + auto volume = overlap.volume(); + vol += volume; + RAJA::atomicAdd(contribCountPtr, std::int64_t(volume >= EPS)); + } + ovlap[cellId] = vol; + }); + + statistics["clip"].set_int64(cellIds.size() * NUM_TETS_PER_HEX); + axom::copy(&contribCount, contribCountPtr, sizeof(contribCount)); + axom::deallocate(contribCountPtr); +} + +template +void Plane3DClipper::specializedClipTetsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& tetIds, + conduit::Node& statistics) +{ + constexpr double EPS = 1e-10; + using ATOMIC_POL = typename axom::execution_space::atomic_policy; + + int allocId = shapeMesh.getAllocatorID(); + + auto meshTets = shapeMesh.getCellsAsTets(); + IndexType tetCount = tetIds.size(); + auto plane = m_plane; + + const std::int64_t zero = 0; + std::int64_t& contribCount = *(statistics["contribs"] = zero).as_int64_ptr(); + std::int64_t* contribCountPtr = axom::allocate(1, allocId); + axom::copy(contribCountPtr, &contribCount, sizeof(zero)); + + axom::for_all( + tetCount, + AXOM_LAMBDA(axom::IndexType ti) { + axom::IndexType tetId = tetIds[ti]; + axom::IndexType cellId = tetId / NUM_TETS_PER_HEX; + const auto& tet = meshTets[tetId]; + primal::Polyhedron overlap = primal::clip(tet, plane, EPS); + double vol = overlap.volume(); + RAJA::atomicAdd(ovlap.data() + cellId, vol); + RAJA::atomicAdd(contribCountPtr, std::int64_t(vol >= EPS)); + }); + + statistics["clip"].set_int64(tetIds.size()); + axom::copy(&contribCount, contribCountPtr, sizeof(contribCount)); + axom::deallocate(contribCountPtr); +} + +void Plane3DClipper::extractClipperInfo() +{ + const auto normal = m_info.fetch_existing("normal").as_double_array(); + const double offset = m_info.fetch_existing("offset").as_double(); + Vector3DType nVec; + for(int d = 0; d < 3; ++d) + { + nVec[d] = normal[d]; + } + m_plane = Plane3DType(nVec, offset); +} + +} // namespace experimental +} // end namespace quest +} // end namespace axom diff --git a/src/axom/quest/detail/clipping/Plane3DClipper.hpp b/src/axom/quest/detail/clipping/Plane3DClipper.hpp new file mode 100644 index 0000000000..37acda4502 --- /dev/null +++ b/src/axom/quest/detail/clipping/Plane3DClipper.hpp @@ -0,0 +1,102 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#ifndef AXOM_QUEST_PLANE3DCLIPPER_HPP +#define AXOM_QUEST_PLANE3DCLIPPER_HPP + +#include "axom/klee/Geometry.hpp" +#include "axom/quest/MeshClipperStrategy.hpp" + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +/*! + * @brief Geometry clipping operations for plane geometries. +*/ +class Plane3DClipper : public MeshClipperStrategy +{ +public: + /*! + * @brief Constructor. + + * @param [in] kGeom Describes the shape to place + * into the mesh. + * @param [in] name To override the default strategy name + + * Clipping operations for a semi-infinite half-space + * on the positive normal direction of a plane. + */ + Plane3DClipper(const klee::Geometry& kGeom, const std::string& name = ""); + + virtual ~Plane3DClipper() = default; + + const std::string& name() const override { return m_name; } + + bool labelCellsInOut(quest::experimental::ShapeMesh& shappeMesh, + axom::Array& label) override; + + bool labelTetsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::Array& tetLabels) override; + + bool specializedClipCells(quest::experimental::ShapeMesh& shappeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics) override; + + bool specializedClipCells(quest::experimental::ShapeMesh& shappeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& cellIds, + conduit::Node& statistics) override; + + bool specializedClipTets(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& tetIds, + conduit::Node& statistics) override; + +#if !defined(__CUDACC__) +private: +#endif + std::string m_name; + + axom::primal::Plane m_plane; + + template + void labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView label); + + template + void labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::ArrayView label); + + template + void specializedClipCellsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics); + + template + void specializedClipCellsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& cellIds, + conduit::Node& statistics); + + template + void specializedClipTetsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + const axom::ArrayView& tetIds, + conduit::Node& statistics); + + void extractClipperInfo(); +}; + +} // namespace experimental +} // namespace quest +} // namespace axom + +#endif // AXOM_QUEST_PLANE3DCLIPPER_HPP diff --git a/src/axom/quest/detail/clipping/TetClipper.cpp b/src/axom/quest/detail/clipping/TetClipper.cpp index 564f60cbb3..3eee2cb5a6 100644 --- a/src/axom/quest/detail/clipping/TetClipper.cpp +++ b/src/axom/quest/detail/clipping/TetClipper.cpp @@ -17,19 +17,306 @@ namespace experimental TetClipper::TetClipper(const klee::Geometry& kGeom, const std::string& name) : MeshClipperStrategy(kGeom) , m_name(name.empty() ? std::string("Tet") : name) - , m_transformer(m_extTrans) + , m_extTransformer(m_extTrans) { extractClipperInfo(); for(int i = 0; i < TetrahedronType::NUM_VERTS; ++i) { - m_tet[i] = m_transformer.getTransformed(m_tetBeforeTrans[i]); + m_tet[i] = m_extTransformer.getTransformed(m_tetBeforeTrans[i]); } - for(int i = 0; i < TetrahedronType::NUM_VERTS; ++i) + /* + * Compute the transformation from m_tet to a unit tet. Location of + * points w.r.t. the tet are done in the unit tet space. The unit + * tet has heights 1, 1, 1 and 1/sqrt(3) as it rests on each of its 4 + * faces. Height w.r.t. the 4th face are scaled by sqrt(3) to make + * it come out to 1. So height < 0 means means below the tet and + * height > 1 means above the tet. + * + * Points with any height outside [0,1] cannot possibly be in the tet. + * Points with all 4 heights in [0,1] must be in the tet. + */ + Point3DType unitTetPts[] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; + m_toUnitTet.setByTerminusPts(&m_tet[0], unitTetPts); +} + +bool TetClipper::labelCellsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& cellLabels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "FSorClipper requires a 3D mesh."); + + const int allocId = shapeMesh.getAllocatorID(); + const auto cellCount = shapeMesh.getCellCount(); + if(cellLabels.size() < cellCount || cellLabels.getAllocatorID() != allocId) + { + cellLabels = axom::Array(ArrayOptions::Uninitialized(), cellCount, cellCount, allocId); + } + + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + labelCellsInOutImpl(shapeMesh, cellLabels.view()); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + labelCellsInOutImpl(shapeMesh, cellLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + labelCellsInOutImpl>(shapeMesh, cellLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + labelCellsInOutImpl>(shapeMesh, cellLabels.view()); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +template +void TetClipper::labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellLabels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "TetClipper requires a 3D mesh."); + /* + * Compute whether the mesh vertices are above the tet or below as + * the tet rests on its 4 facets. + * + * - For any facet the tet rests on, if all cell vertices are + * above the tet or all are below, cell is OUT. + * + * - If all cell vertices are on the tet side of all facets, + * the cell is IN. + * + * - Otherwise, cell is ON. + */ + + int allocId = shapeMesh.getAllocatorID(); + auto vertCount = shapeMesh.getVertexCount(); + auto cellCount = shapeMesh.getCellCount(); + auto meshCellVolumes = shapeMesh.getCellVolumes(); + + const auto& vertCoords = shapeMesh.getVertexCoords3D(); + const auto& vX = vertCoords[0]; + const auto& vY = vertCoords[1]; + const auto& vZ = vertCoords[2]; + SLIC_ASSERT(axom::execution_space::usesAllocId(vX.getAllocatorID())); + SLIC_ASSERT(axom::execution_space::usesAllocId(vY.getAllocatorID())); + SLIC_ASSERT(axom::execution_space::usesAllocId(vZ.getAllocatorID())); + + axom::ArrayView connView = shapeMesh.getCellNodeConnectivity(); + SLIC_ASSERT(connView.shape() == + (axom::StackArray {cellCount, HexahedronType::NUM_HEX_VERTS})); + + axom::Array below[4]; + axom::Array above[4]; + axom::ArrayView belowView[4]; + axom::ArrayView aboveView[4]; + for(IndexType p = 0; p < 4; ++p) + { + below[p] = axom::Array(ArrayOptions::Uninitialized(), vertCount, vertCount, allocId); + above[p] = axom::Array(ArrayOptions::Uninitialized(), vertCount, vertCount, allocId); + belowView[p] = below[p].view(); + aboveView[p] = above[p].view(); + } + + auto toUnitTet = m_toUnitTet; + + axom::for_all( + vertCount, + AXOM_LAMBDA(axom::IndexType vertId) { + // vh is the heights of the vertex in the space of the unit tet. + // See comment on m_toUnitTet. + axom::NumericArray vh({vX[vertId], vY[vertId], vZ[vertId], 0}); + toUnitTet.transform(vh[0], vh[1], vh[2]); + vh[3] = 1 - vh[0] - vh[1] - vh[2]; + + for(int p = 0; p < 4; ++p) + { + belowView[p][vertId] = vh[p] < 0; + aboveView[p][vertId] = vh[p] > 1; + } + }); + + constexpr double EPS = 1e-10; + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellId) { + if(axom::utilities::isNearlyEqual(meshCellVolumes[cellId], 0.0, EPS)) + { + cellLabels[cellId] = LabelType::LABEL_OUT; + return; + } + + LabelType& cellLabel = cellLabels[cellId]; + auto cellVertIds = connView[cellId]; + + cellLabel = LabelType::LABEL_ON; + bool vertsAreOnTetSideOfAllPlanes = true; + for(IndexType p = 0; p < 4; ++p) + { + bool allVertsBelow = true; + bool allVertsAbove = true; + for(int vi = 0; vi < HexahedronType::NUM_HEX_VERTS; ++vi) + { + int vertId = cellVertIds[vi]; + auto vertIsBelow = belowView[p][vertId]; + auto vertIsAbove = aboveView[p][vertId]; + allVertsBelow &= vertIsBelow; + allVertsAbove &= vertIsAbove; + vertsAreOnTetSideOfAllPlanes &= !vertIsBelow; + } + if(allVertsBelow || allVertsAbove) + { + cellLabel = LabelType::LABEL_OUT; + break; + } + } + if(cellLabel != LabelType::LABEL_OUT && vertsAreOnTetSideOfAllPlanes) + { + cellLabel = LabelType::LABEL_IN; + } + }); + + return; +} + +bool TetClipper::labelTetsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::Array& tetLabels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "TetClipper requires a 3D mesh."); + + const int allocId = shapeMesh.getAllocatorID(); + const auto cellCount = cellIds.size(); + const auto tetCount = cellCount * NUM_TETS_PER_HEX; + if(tetLabels.size() < tetCount || tetLabels.getAllocatorID() != allocId) + { + tetLabels = axom::Array(ArrayOptions::Uninitialized(), tetCount, tetCount, allocId); + } + + switch(shapeMesh.getRuntimePolicy()) { - m_bb.addPoint(m_tet[i]); + case axom::runtime_policy::Policy::seq: + labelTetsInOutImpl(shapeMesh, cellIds, tetLabels.view()); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + labelTetsInOutImpl(shapeMesh, cellIds, tetLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + labelTetsInOutImpl>(shapeMesh, cellIds, tetLabels.view()); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + labelTetsInOutImpl>(shapeMesh, cellIds, tetLabels.view()); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); } + return true; +} + +template +void TetClipper::labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::ArrayView tetLabels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "TetClipper requires a 3D mesh."); + /* + * Compute whether the tets in hexes listed in cellIds + * are in, out or on the boundary. + * + * Picture the tet resting on each of of its 4 faces. + * - For any facet the tet rests on, if all cell vertices are + * above the tet or all are below, cell is OUT. + * + * - If all cell vertices are on the tet side of all facets, + * the cell is IN. + * + * - Otherwise, cell is ON. + */ + + auto meshTets = shapeMesh.getCellsAsTets(); + auto tetVolumes = shapeMesh.getTetVolumes(); + + const axom::IndexType cellCount = cellIds.size(); + + auto toUnitTet = m_toUnitTet; + constexpr double EPS = 1e-10; + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType ci) { + axom::IndexType cellId = cellIds[ci]; + + const TetrahedronType* tetsForCell = &meshTets[cellId * NUM_TETS_PER_HEX]; + + for(IndexType ti = 0; ti < NUM_TETS_PER_HEX; ++ti) + { + const TetrahedronType& cellTet = tetsForCell[ti]; + LabelType& tetLabel = tetLabels[ci * NUM_TETS_PER_HEX + ti]; + const axom::IndexType tetId = cellId * NUM_TETS_PER_HEX + ti; + + if(axom::utilities::isNearlyEqual(tetVolumes[tetId], 0.0, EPS)) + { + tetLabel = LabelType::LABEL_OUT; + continue; + } + + tetLabel = LabelType::LABEL_ON; + + bool allVertsBelow = true; + bool allVertsAbove = true; + bool vertsAreOnTetSideOfAllPlanes = true; + + for(IndexType vi = 0; vi < 4; ++vi) + { + const auto& vert = cellTet[vi]; + + // vh is the heights of vert in the space of the unit tet. + // See comment on m_toUnitTet. + axom::NumericArray vh({vert[0], vert[1], vert[2], 0}); + toUnitTet.transform(vh[0], vh[1], vh[2]); + vh[3] = 1 - vh[0] - vh[1] - vh[2]; + + // Where vertex vi is w.r.t. the tet resting on side pj. + for(int pj = 0; pj < 4; ++pj) + { + bool vertIsBelow = vh[pj] < 0; + bool vertIsAbove = vh[pj] > 1; + + allVertsBelow &= vertIsBelow; + allVertsAbove &= vertIsAbove; + vertsAreOnTetSideOfAllPlanes &= !vertIsBelow; + } + + if(allVertsBelow || allVertsAbove) + { + tetLabel = LabelType::LABEL_OUT; + break; + } + } + + if(tetLabel != LabelType::LABEL_OUT && vertsAreOnTetSideOfAllPlanes) + { + tetLabel = LabelType::LABEL_IN; + } + } + }); + + return; } bool TetClipper::getGeometryAsTets(quest::experimental::ShapeMesh& shapeMesh, @@ -63,6 +350,30 @@ void TetClipper::extractClipperInfo() m_tetBeforeTrans[2][d] = v2[d]; m_tetBeforeTrans[3][d] = v3[d]; } + + bool fixOrientation = false; + if(m_info.has_child("fixOrientation")) + { + fixOrientation = bool(m_info.fetch_existing("fixOrientation").as_int()); + } + + if(fixOrientation) + { + m_tetBeforeTrans.checkAndFixOrientation(); + } + else + { + constexpr double EPS = 1e-10; + double signedVol = m_tetBeforeTrans.signedVolume(); + if(signedVol < -EPS) + { + SLIC_ERROR( + axom::fmt::format("TetClipper tet {} has negative volume {}.:" + " (See TetClipper's 'fixOrientation' flag.)", + m_tetBeforeTrans, + signedVol)); + } + } } } // namespace experimental diff --git a/src/axom/quest/detail/clipping/TetClipper.hpp b/src/axom/quest/detail/clipping/TetClipper.hpp index c0062ea498..530ce312ff 100644 --- a/src/axom/quest/detail/clipping/TetClipper.hpp +++ b/src/axom/quest/detail/clipping/TetClipper.hpp @@ -33,6 +33,9 @@ class TetClipper : public MeshClipperStrategy * \c kGeom.asHierarchy() must contain the following data: * - v0, v1, v2, v3: each contains a 3D coordinates of the * tetrahedron vertices, in the order used by primal::Tetrahedron. + * The tet may be degenerate, but not inverted (negative volume). + * - "fixOrientation": Whether to fix inverted tetrahedra + * instead of aborting. */ TetClipper(const klee::Geometry& kGeom, const std::string& name = ""); @@ -40,13 +43,20 @@ class TetClipper : public MeshClipperStrategy const std::string& name() const override { return m_name; } + bool labelCellsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& cellLabels) override; + + bool labelTetsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::Array& tetLabels) override; + /*! * @copydoc MeshClipperStrategy::getGeometryAsTets() * * \c tets will have length one, because the geometry for this * class is a single tetrahedron. */ - bool getGeometryAsTets(quest::experimental::ShapeMesh& shappeMesh, + bool getGeometryAsTets(quest::experimental::ShapeMesh& shapeMesh, axom::Array& tets) override; #if !defined(__CUDACC__) @@ -54,15 +64,26 @@ class TetClipper : public MeshClipperStrategy #endif std::string m_name; - //!@brief Tetrahedron before transformation. + //!@brief Tetrahedron before external transformation. TetrahedronType m_tetBeforeTrans; - //!@brief Tetrahedron after transformation. + //!@brief Tetrahedron after external transformation. TetrahedronType m_tet; - axom::primal::BoundingBox m_bb; + //!@brief External transformation. + axom::primal::experimental::CoordinateTransformer m_extTransformer; + + //!@brief Transformation of m_tet to unit tet. + axom::primal::experimental::CoordinateTransformer m_toUnitTet; + + template + void labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellLabel); - axom::primal::experimental::CoordinateTransformer m_transformer; + template + void labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellsOnBdry, + axom::ArrayView tetLabels); // Extract clipper info from MeshClipperStrategy::m_info. void extractClipperInfo(); diff --git a/src/axom/quest/docs/sphinx/point_mesh_query_cpp.rst b/src/axom/quest/docs/sphinx/point_mesh_query_cpp.rst index 8663f160ab..e1479e01bf 100644 --- a/src/axom/quest/docs/sphinx/point_mesh_query_cpp.rst +++ b/src/axom/quest/docs/sphinx/point_mesh_query_cpp.rst @@ -56,7 +56,7 @@ Signed Distance --------------- The C++ signed distance query is provided by the ``quest::SignedDistance`` class, -which wraps an instance of ``primal::BVHTree``. +which wraps an instance of ``spin::BVH``. Examples from ``/src/axom/quest/tests/quest_signed_distance.cpp``. Class header: diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index ffa3e47f48..8c059b279b 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -337,6 +337,7 @@ endif() 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 RAJA_FOUND AND AXOM_ENABLE_KLEE) if(MFEM_FOUND) set(optional_dependency, "mfem") @@ -368,7 +369,7 @@ if((CONDUIT_FOUND OR blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) endif() - set(_testshapes "tetmesh" "tet" "hex" "sphere" "cyl" "cone" "sor" "all" "plane") + set(_testshapes "tetmesh" "tet" "hex" "sphere" "cyl" "cone" "sor" "plane" "tetmesh,tet,hex") foreach(_policy ${_policies}) foreach(_testshape ${_testshapes}) set(_testname "quest_shape_in_memory_${_policy}_${_testshape}") @@ -377,11 +378,11 @@ if((CONDUIT_FOUND OR 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 + --refinements 5 + --scale .99 .99 .99 + --dir 8 4 2 --meshType bpSidre - inline_mesh --min -2 -2 -2 --max 2 2 2 --res 30 30 30 + inline_mesh --min -2 -2 -2 --max 2 2 2 --res 16 16 16 NUM_MPI_TASKS ${_nranks}) endforeach() endforeach() @@ -399,9 +400,9 @@ if((CONDUIT_FOUND OR 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 + --refinements 5 + --scale .99 .99 .99 + --dir 8 4 2 --meshType ${_testMeshType} inline_mesh --min -2 -2 -2 --max 2 2 2 --res 8 8 8 NUM_MPI_TASKS ${_nranks}) diff --git a/src/axom/quest/examples/quest_shape_in_memory.cpp b/src/axom/quest/examples/quest_shape_in_memory.cpp index f3feb070de..bd67b0ac55 100644 --- a/src/axom/quest/examples/quest_shape_in_memory.cpp +++ b/src/axom/quest/examples/quest_shape_in_memory.cpp @@ -27,6 +27,10 @@ #error Shaping functionality requires Axom to be configured with Conduit #endif +#if !defined(AXOM_USE_RAJA) + #error Shaping test requires Axom to be configured with RAJA +#endif + #include "conduit_relay_io_blueprint.hpp" #if defined(AXOM_USE_MFEM) @@ -92,10 +96,10 @@ struct Input } // The shape to run. - std::string testShape {"tetmesh"}; + std::vector testShape; // The shapes this example is set up to run. const std::set - availableShapes {"tetmesh", "tri", "sphere", "cyl", "cone", "sor", "tet", "hex", "plane", "all"}; + availableShapes {"tetmesh", "tri", "sphere", "cyl", "cone", "sor", "tet", "hex", "plane"}; RuntimePolicy policy {RuntimePolicy::seq}; int outputOrder {2}; @@ -230,8 +234,12 @@ struct Input ->capture_default_str(); app.add_option("-s,--testShape", testShape) - ->description("The shape to run") - ->check(axom::CLI::IsMember(availableShapes)); + ->description( + "The shape(s) to run. Specifying multiple shapes will override scaling and translations " + "to shrink shapes and shift them to individual octants of the mesh.") + ->check(axom::CLI::IsMember(availableShapes)) + ->delimiter(',') + ->expected(1, 60); #ifdef AXOM_USE_CALIPER app.add_option("--caliper", annotationMode) @@ -333,6 +341,40 @@ struct Input }; // struct Input Input params; +/************************************************************ + * Shared variables. + ************************************************************/ + +const std::string topoName = "mesh"; +const std::string coordsetName = "coords"; +int cellCount = -1; +// Translation to individual octants (override) when running multiple shapes. +// Except that the plane is never moved. +std::vector> translations {{1, 1, -1}, + {-1, 1, -1}, + {-1, -1, -1}, + {1, -1, -1}, + {1, 1, 1}, + {-1, 1, 1}, + {-1, -1, 1}, + {1, -1, 1}}; +int translationIdx = 0; // To track what translations have been used. +std::map shapeReps; // Repetitions of the geometry. +std::map exactOverlapVols; +std::map errorToleranceRel; // Relative error tolerance. +std::map errorToleranceAbs; // Absolute error tolerance. +double vScale = 1.0; // Volume scale due to geometry scale. + +const auto hostAllocId = axom::execution_space::allocatorID(); +int arrayAllocId = axom::INVALID_ALLOCATOR_ID; + +// 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; + // Start property for all 3D shapes. axom::klee::TransformableGeometryProperties startProp {axom::klee::Dimensions::Three, axom::klee::LengthUnit::unspecified}; @@ -353,14 +395,23 @@ void addScaleOperator(axom::klee::CompositeOperator& compositeOp) } // Add translate operator. -void addTranslateOperator(axom::klee::CompositeOperator& compositeOp, - double shiftx, - double shifty, - double shiftz) +void addTranslateOperator(axom::klee::CompositeOperator& compositeOp) { - primal::Vector3D shift({shiftx, shifty, shiftz}); - auto translateOp = std::make_shared(shift, startProp); - compositeOp.addOperator(translateOp); + if(params.testShape.size() > 1) + { + const axom::NumericArray& shifts = + translations[(translationIdx++) % translations.size()]; + primal::Vector3D shift({shifts[0], shifts[1], shifts[2]}); + auto translateOp = std::make_shared(shift, startProp); + compositeOp.addOperator(translateOp); + } + else + { + // Use zero shift as a smoke test. + primal::Vector3D shift({0, 0, 0}); + auto translateOp = std::make_shared(shift, startProp); + compositeOp.addOperator(translateOp); + } } // Add operator to rotate x-axis to params.direction, if it is given. @@ -448,20 +499,6 @@ void printMfemMeshInfo(mfem::Mesh* mesh, const std::string& prefixMessage = "") * Shared variables. ************************************************************/ -const std::string topoName = "mesh"; -const std::string coordsetName = "coords"; -int cellCount = -1; - -const auto hostAllocId = axom::execution_space::allocatorID(); -int arrayAllocId = axom::INVALID_ALLOCATOR_ID; - -// 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) { switch(params.getBoxDim()) @@ -555,55 +592,6 @@ void finalizeLogger() } } -// Single triangle ShapeSet. -axom::klee::ShapeSet create2DShapeSet(sidre::DataStore& ds) -{ - sidre::Group* meshGroup = ds.getRoot()->createGroup("triangleMesh"); - const std::string topo = "mesh"; - const std::string coordset = "coords"; - axom::mint::UnstructuredMesh triangleMesh(2, - axom::mint::CellType::TRIANGLE, - meshGroup, - topo, - coordset); - - const double lll = 2.0; - - // Insert tet at origin. - triangleMesh.appendNode(0.0, 0.0); - triangleMesh.appendNode(lll, 0.0); - triangleMesh.appendNode(0.0, lll); - axom::IndexType conn[3] = {0, 1, 2}; - triangleMesh.appendCell(conn); - - SLIC_ERROR_ROOT_IF(!axom::mint::blueprint::isValidRootGroup(meshGroup), - "Triangle mesh blueprint is not valid"); - - axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Two, - axom::klee::LengthUnit::unspecified}; - axom::klee::Geometry triangleGeom(prop, triangleMesh.getSidreGroup(), topo, nullptr); - - std::vector shapes; - axom::klee::Shape triangleShape( - "triangle", - "AL", - {}, - {}, - axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}); - shapes.push_back( - axom::klee::Shape {"triangle", - "AL", - {}, - {}, - axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}}); - - axom::klee::ShapeSet shapeSet; - shapeSet.setShapes(shapes); - shapeSet.setDimensions(axom::klee::Dimensions::Two); - - return shapeSet; -} - axom::klee::Shape createShape_Sphere() { Point3D center = params.center.empty() ? Point3D {0, 0, 0} : Point3D {params.center.data()}; @@ -616,7 +604,7 @@ axom::klee::Shape createShape_Sphere() auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); addRotateOperator(*compositeOp); - addTranslateOperator(*compositeOp, 1, 1, 1); + addTranslateOperator(*compositeOp); const axom::IndexType levelOfRefinement = params.refinementLevel; axom::klee::Geometry sphereGeometry(prop, sphere, levelOfRefinement, compositeOp); @@ -662,7 +650,7 @@ axom::klee::Shape createShape_TetMesh(sidre::DataStore& ds) auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); addRotateOperator(*compositeOp); - addTranslateOperator(*compositeOp, -1, 1, 1); + addTranslateOperator(*compositeOp); axom::klee::Geometry tetMeshGeometry(prop, tetMesh.getSidreGroup(), topo, compositeOp); axom::klee::Shape tetShape("tetmesh", "TETMESH", {}, {}, tetMeshGeometry); @@ -672,7 +660,7 @@ axom::klee::Shape createShape_TetMesh(sidre::DataStore& ds) axom::klee::Geometry createGeometry_Sor(axom::primal::Point& sorBase, axom::primal::Vector& sorDirection, - axom::Array& discreteFunction, + axom::ArrayView discreteFunction, std::shared_ptr& compositeOp) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, @@ -690,7 +678,22 @@ axom::klee::Geometry createGeometry_Sor(axom::primal::Point& sorBase, return sorGeometry; } -axom::klee::Shape createShape_Sor() +double computeVolume_Sor(axom::Array& discreteFunction) +{ + using ConeType = axom::primal::Cone; + axom::IndexType segmentCount = discreteFunction.shape()[0]; + double vol = 0.0; + for(axom::IndexType s = 0; s < segmentCount - 1; ++s) + { + ConeType cone(discreteFunction(s, 1), + discreteFunction(s + 1, 1), + discreteFunction(s + 1, 0) - discreteFunction(s, 0)); + vol += cone.volume(); + } + return vol; +} + +axom::klee::Shape createShape_Sor(const std::string& shapeName) { Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; axom::primal::Vector sorDirection = params.direction.empty() @@ -700,36 +703,40 @@ axom::klee::Shape createShape_Sor() // discreteFunction are discrete z-r pairs describing the function // to be rotated around the z axis. axom::Array discreteFunction({numIntervals + 1, 2}, axom::ArrayStrideOrder::ROW); - double zLen = params.length < 0 ? 1.6 : params.length; + double zLen = params.length < 0 ? 2.40 : params.length; double zShift = -zLen / 2; - double maxR = params.radius < 0 ? 0.75 : params.radius; + double maxR = params.radius < 0 ? 1.10 : params.radius; double dz = zLen / numIntervals; - discreteFunction[0][0] = 0 * dz + zShift; - discreteFunction[0][1] = 0.0 * maxR; - discreteFunction[1][0] = 1 * dz + zShift; - discreteFunction[1][1] = 0.8 * maxR; - discreteFunction[2][0] = 2 * dz + zShift; - discreteFunction[2][1] = 0.4 * maxR; - discreteFunction[3][0] = 3 * dz + zShift; - discreteFunction[3][1] = 0.5 * maxR; - discreteFunction[4][0] = 4 * dz + zShift; - discreteFunction[4][1] = 1.0 * maxR; - discreteFunction[5][0] = 5 * dz + zShift; - discreteFunction[5][1] = 0.0; + discreteFunction(0, 0) = 0 * dz + zShift; + discreteFunction(0, 1) = 0.0 * maxR; + discreteFunction(1, 0) = 1 * dz + zShift; + discreteFunction(1, 1) = 0.8 * maxR; + discreteFunction(2, 0) = 2 * dz + zShift; + discreteFunction(2, 1) = 0.4 * maxR; + discreteFunction(3, 0) = 3 * dz + zShift; + discreteFunction(3, 1) = 0.5 * maxR; + discreteFunction(4, 0) = 4 * dz + zShift; + discreteFunction(4, 1) = 1.0 * maxR; + discreteFunction(5, 0) = 5 * dz + zShift; + discreteFunction(5, 1) = 1.0 * maxR; auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); - addTranslateOperator(*compositeOp, -1, -1, 1); + addTranslateOperator(*compositeOp); axom::klee::Geometry sorGeometry = createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape sorShape("sor", "SOR", {}, {}, sorGeometry); + axom::klee::Shape sorShape(shapeName, shapeName + ".mat", {}, {}, sorGeometry); + + exactOverlapVols[shapeName] = vScale * computeVolume_Sor(discreteFunction); + errorToleranceRel[shapeName] = 0.04; + errorToleranceAbs[shapeName] = 0.15; return sorShape; } -axom::klee::Shape createShape_Cylinder() +axom::klee::Shape createShape_Cylinder(const std::string& shapeName) { Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; axom::primal::Vector sorDirection = params.direction.empty() @@ -738,8 +745,8 @@ axom::klee::Shape createShape_Cylinder() // discreteFunction are discrete z-r pairs describing the function // to be rotated around the z axis. axom::Array discreteFunction({2, 2}, axom::ArrayStrideOrder::ROW); - double radius = params.radius < 0 ? 0.5 : params.radius; - double height = params.length < 0 ? 1.2 : params.length; + double radius = params.radius < 0 ? 0.695 : params.radius; + double height = params.length < 0 ? 2.78 : params.length; discreteFunction[0][0] = -height / 2; discreteFunction[0][1] = radius; discreteFunction[1][0] = height / 2; @@ -747,17 +754,22 @@ axom::klee::Shape createShape_Cylinder() auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); - addTranslateOperator(*compositeOp, 1, -1, 1); + addTranslateOperator(*compositeOp); axom::klee::Geometry sorGeometry = createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape sorShape("cyl", "CYL", {}, {}, sorGeometry); + axom::klee::Shape sorShape(shapeName, shapeName + ".mat", {}, {}, sorGeometry); + + exactOverlapVols[shapeName] = vScale * computeVolume_Sor(discreteFunction); + // error tolerance for 2 levels of refinement + errorToleranceRel[shapeName] = 0.05; + errorToleranceAbs[shapeName] = 0.2; return sorShape; } -axom::klee::Shape createShape_Cone() +axom::klee::Shape createShape_Cone(const std::string& shapeName) { Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; axom::primal::Vector sorDirection = params.direction.empty() @@ -766,9 +778,9 @@ axom::klee::Shape createShape_Cone() // discreteFunction are discrete z-r pairs describing the function // to be rotated around the z axis. axom::Array discreteFunction({2, 2}, axom::ArrayStrideOrder::ROW); - double baseRadius = params.radius < 0 ? 0.7 : params.radius; - double topRadius = params.radius2 < 0 ? 0.1 : params.radius2; - double height = params.length < 0 ? 1.3 : params.length; + double baseRadius = params.radius < 0 ? 1.23 : params.radius; + double topRadius = params.radius2 < 0 ? 0.176 : params.radius2; + double height = params.length < 0 ? 2.3 : params.length; discreteFunction[0][0] = -height / 2; discreteFunction[0][1] = baseRadius; discreteFunction[1][0] = height / 2; @@ -776,17 +788,21 @@ axom::klee::Shape createShape_Cone() auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); - addTranslateOperator(*compositeOp, 1, 1, -1); + addTranslateOperator(*compositeOp); axom::klee::Geometry sorGeometry = createGeometry_Sor(sorBase, sorDirection, discreteFunction, compositeOp); - axom::klee::Shape sorShape("cone", "CONE", {}, {}, sorGeometry); + axom::klee::Shape sorShape(shapeName, shapeName + ".mat", {}, {}, sorGeometry); + + exactOverlapVols[shapeName] = vScale * computeVolume_Sor(discreteFunction); + errorToleranceRel[shapeName] = 0.05; + errorToleranceAbs[shapeName] = 0.2; return sorShape; } -axom::klee::Shape createShape_Tet() +axom::klee::Shape createShape_Tet(const std::string& shapeName) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, axom::klee::LengthUnit::unspecified}; @@ -794,32 +810,35 @@ axom::klee::Shape createShape_Tet() SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); // Tetrahedron at origin. - const double len = params.length < 0 ? 0.8 : params.length; - const Point3D a {-len, -len, -len}; - const Point3D b {+len, -len, -len}; - const Point3D c {+len, +len, -len}; - const Point3D d {-len, +len, +len}; + const double len = params.length < 0 ? 1.55 : params.length; + const Point3D a {Point3D::NumericArray {1., 0., -1.} * len}; + const Point3D b {Point3D::NumericArray {-.8, 1, -1.} * len}; + const Point3D c {Point3D::NumericArray {-.8, -1, -1.} * len}; + const Point3D d {Point3D::NumericArray {0., 0., +1.} * len}; const primal::Tetrahedron tet {a, b, c, d}; auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); addRotateOperator(*compositeOp); - addTranslateOperator(*compositeOp, -1, 1, -1); + addTranslateOperator(*compositeOp); + exactOverlapVols[shapeName] = vScale * tet.volume(); + errorToleranceRel[shapeName] = 1e-6; + errorToleranceAbs[shapeName] = 1e-8; axom::klee::Geometry tetGeometry(prop, tet, compositeOp); - axom::klee::Shape tetShape("tet", "TET", {}, {}, tetGeometry); + axom::klee::Shape tetShape(shapeName, shapeName + ".mat", {}, {}, tetGeometry); return tetShape; } -axom::klee::Shape createShape_Hex() +axom::klee::Shape createShape_Hex(const std::string& shapeName) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, axom::klee::LengthUnit::unspecified}; SLIC_ASSERT(params.scaleFactors.empty() || params.scaleFactors.size() == 3); - const double md = params.length < 0 ? 0.6 : params.length / 2; + const double md = params.length < 0 ? 0.82 : params.length / 2; const double lg = 1.2 * md; const double sm = 0.8 * md; const Point3D p {-lg, -md, -sm}; @@ -835,15 +854,18 @@ axom::klee::Shape createShape_Hex() auto compositeOp = std::make_shared(startProp); addScaleOperator(*compositeOp); addRotateOperator(*compositeOp); - addTranslateOperator(*compositeOp, -1, -1, -1); + addTranslateOperator(*compositeOp); + exactOverlapVols[shapeName] = vScale * hex.volume(); + errorToleranceRel[shapeName] = 1e-6; + errorToleranceAbs[shapeName] = 1e-8; axom::klee::Geometry hexGeometry(prop, hex, compositeOp); - axom::klee::Shape hexShape("hex", "HEX", {}, {}, hexGeometry); + axom::klee::Shape hexShape(shapeName, shapeName + ".mat", {}, {}, hexGeometry); return hexShape; } -axom::klee::Shape createShape_Plane() +axom::klee::Shape createShape_Plane(const std::string& shapeName) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, axom::klee::LengthUnit::unspecified}; @@ -868,21 +890,21 @@ axom::klee::Shape createShape_Plane() const primal::Plane plane {normal, center, true}; axom::klee::Geometry planeGeometry(prop, plane, scaleOp); - axom::klee::Shape planeShape("plane", "PLANE", {}, {}, planeGeometry); + axom::klee::Shape planeShape(shapeName, shapeName + ".mat", {}, {}, planeGeometry); + + // Exact mesh overlap volume, assuming plane passes through center of box mesh. + using Pt3D = primal::Point; + Pt3D lower(params.boxMins.data()); + Pt3D upper(params.boxMaxs.data()); + auto diag = upper.array() - lower.array(); + double meshVolume = diag[0] * diag[1] * diag[2]; + exactOverlapVols[shapeName] = 0.5 * meshVolume; + errorToleranceRel[shapeName] = 1e-6; + errorToleranceAbs[shapeName] = 1e-8; return planeShape; } -//!@brief Create a ShapeSet with a single shape. -axom::klee::ShapeSet createShapeSet(const axom::klee::Shape& shape) -{ - axom::klee::ShapeSet shapeSet; - shapeSet.setShapes(std::vector {shape}); - shapeSet.setDimensions(axom::klee::Dimensions::Three); - - return shapeSet; -} - double volumeOfTetMesh(const axom::mint::UnstructuredMesh& tetMesh) { using TetType = axom::primal::Tetrahedron; @@ -933,6 +955,138 @@ double areaOfTriMesh(const axom::mint::UnstructuredMesh create2DShapeSet(sidre::DataStore& ds) +{ + sidre::Group* meshGroup = ds.getRoot()->createGroup("triangleMesh"); + AXOM_UNUSED_VAR(meshGroup); // variable is only referenced in debug configs + const std::string topo = "mesh"; + const std::string coordset = "coords"; + axom::mint::UnstructuredMesh triangleMesh(2, + axom::mint::CellType::TRIANGLE, + meshGroup, + topo, + coordset); + + double lll = 2.0; + + // Insert tet at origin. + triangleMesh.appendNode(0.0, 0.0); + triangleMesh.appendNode(lll, 0.0); + triangleMesh.appendNode(0.0, lll); + axom::IndexType conn[3] = {0, 1, 2}; + triangleMesh.appendCell(conn); + + SLIC_ASSERT(axom::mint::blueprint::isValidRootGroup(meshGroup)); + + axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Two, + axom::klee::LengthUnit::unspecified}; + axom::klee::Geometry triangleGeom(prop, triangleMesh.getSidreGroup(), topo, nullptr); + + std::vector shapes; + axom::klee::Shape triangleShape( + "triangle", + "AL", + {}, + {}, + axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}); + shapes.push_back( + axom::klee::Shape {"triangle", + "AL", + {}, + {}, + axom::klee::Geometry {prop, triangleMesh.getSidreGroup(), topo, nullptr}}); + + axom::klee::ShapeSet shapeSet; + shapeSet.setShapes(shapes); + shapeSet.setDimensions(axom::klee::Dimensions::Two); + + double shapeMeshVol = areaOfTriMesh(triangleMesh); + exactOverlapVols[triangleShape.getName()] = shapeMeshVol; + errorToleranceRel[triangleShape.getName()] = 1e-6; + errorToleranceAbs[triangleShape.getName()] = 1e-8; + + return shapes; +} + +axom::klee::Shape createShape_Sphere(const std::string& shapeName) +{ + Point3D center = params.center.empty() ? Point3D {0, 0, 0} : Point3D {params.center.data()}; + double radius = params.radius < 0 ? 1.0 : params.radius; + axom::primal::Sphere sphere {center, radius}; + + axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp); + + const axom::IndexType levelOfRefinement = params.refinementLevel; + axom::klee::Geometry sphereGeometry(prop, sphere, levelOfRefinement, compositeOp); + axom::klee::Shape sphereShape(shapeName, shapeName + ".mat", {}, {}, sphereGeometry); + exactOverlapVols[shapeName] = vScale * 4. / 3 * M_PI * radius * radius * radius; + errorToleranceRel[shapeName] = 0.1; + errorToleranceAbs[shapeName] = 0.38; + + return sphereShape; +} + +axom::klee::Shape createShape_TetMesh(sidre::DataStore& ds, const std::string& shapeName) +{ + // Shape a tetrahedal mesh. + sidre::Group* meshGroup = ds.getRoot()->createGroup(shapeName); + AXOM_UNUSED_VAR(meshGroup); // variable is only referenced in debug configs + const std::string topo = "mesh"; + const std::string coordset = "coords"; + axom::mint::UnstructuredMesh tetMesh(3, + axom::mint::CellType::TET, + meshGroup, + topo, + coordset); + + double lll = params.length < 0 ? 1.17 : params.length; + + // Insert tets around origin. + tetMesh.appendNode(-lll, -lll, -lll); + tetMesh.appendNode(+lll, -lll, -lll); + tetMesh.appendNode(-lll, +lll, -lll); + tetMesh.appendNode(-lll, -lll, +lll); + tetMesh.appendNode(+lll, +lll, +lll); + tetMesh.appendNode(-lll, +lll, +lll); + tetMesh.appendNode(+lll, +lll, -lll); + tetMesh.appendNode(+lll, -lll, +lll); + axom::IndexType conn0[4] = {0, 1, 2, 3}; + tetMesh.appendCell(conn0); + axom::IndexType conn1[4] = {4, 5, 6, 7}; + tetMesh.appendCell(conn1); + + SLIC_ASSERT(axom::mint::blueprint::isValidRootGroup(meshGroup)); + + axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addRotateOperator(*compositeOp); + addTranslateOperator(*compositeOp); + + axom::klee::Geometry tetMeshGeometry(prop, tetMesh.getSidreGroup(), topo, compositeOp); + axom::klee::Shape tetShape(shapeName, shapeName + ".mat", {}, {}, tetMeshGeometry); + + exactOverlapVols[shapeName] = vScale * volumeOfTetMesh(tetMesh); + errorToleranceRel[shapeName] = 1e-6; + errorToleranceAbs[shapeName] = 1e-8; + + return tetShape; +} + #if defined(AXOM_USE_MFEM) /*! * @brief Return the element volumes as a sidre::View. @@ -1086,7 +1240,7 @@ axom::sidre::View* getElementVolumes( * Most of this is lifted from IntersectionShaper::runShapeQueryImpl. */ template -axom::sidre::View* getElementVolumes( +axom::sidre::View* getElementVolumesImpl( sidre::Group* meshGrp, const std::string& volFieldName = std::string("elementVolumes")) { @@ -1296,6 +1450,34 @@ axom::sidre::View* getElementVolumes( return volSidreView; } +axom::sidre::View* getElementVolumes( + sidre::Group* meshGrp, + const std::string& volFieldName = std::string("elementVolumes")) +{ + switch(params.policy) + { +#if defined(AXOM_USE_CUDA) + case RuntimePolicy::cuda: + return getElementVolumesImpl>(meshGrp, volFieldName); + break; +#endif +#if defined(AXOM_USE_HIP) + case RuntimePolicy::hip: + return getElementVolumesImpl>(meshGrp, volFieldName); + break; +#endif +#if defined(AXOM_USE_OMP) + case RuntimePolicy::omp: + return getElementVolumesImpl(meshGrp, volFieldName); + break; +#endif + case RuntimePolicy::seq: + default: + return getElementVolumesImpl(meshGrp, volFieldName); + break; + } + return nullptr; +} #if defined(AXOM_USE_MFEM) /*! @@ -1349,7 +1531,7 @@ double sumMaterialVolumesImpl(sidre::Group* meshGrp, const std::string& material // Get cell volumes from meshGrp. const std::string volsName = "vol_" + material; - axom::sidre::View* elementVols = getElementVolumes(meshGrp, volsName); + axom::sidre::View* elementVols = getElementVolumesImpl(meshGrp, volsName); axom::ArrayView elementVolsView(elementVols->getData(), elementVols->getNumElements()); // Get material volume fractions @@ -1540,11 +1722,23 @@ int main(int argc, char** argv) exit(retval); } + if(params.testShape.size() > 1) + { + SLIC_WARNING( + "Multiple test configurations specified.\n" + "Scaling by half to shrink the geometries\n" + "and move them to individual octants so they don't overlap\n" + "with each other."); + params.scaleFactors.resize(3, 1.0); + for(auto& f : params.scaleFactors) f *= 0.5; + } + vScale = params.scaleFactors[0] * params.scaleFactors[1] * params.scaleFactors[2]; + axom::utilities::raii::AnnotationsWrapper annotations_raii_wrapper(params.annotationMode); const int arrayAllocId = axom::policyToDefaultAllocatorID(params.policy); - AXOM_ANNOTATE_SCOPE("quest shaping example"); + AXOM_ANNOTATE_BEGIN("quest shaping example"); AXOM_ANNOTATE_BEGIN("init"); // Storage for the shape geometry meshes. @@ -1553,58 +1747,58 @@ int main(int argc, char** argv) //--------------------------------------------------------------------------- // Create simple ShapeSet for the example. //--------------------------------------------------------------------------- - axom::klee::ShapeSet shapeSet; - - if(params.testShape == "tetmesh") - { - shapeSet = createShapeSet(createShape_TetMesh(ds)); - } - else if(params.testShape == "tet") - { - shapeSet = createShapeSet(createShape_Tet()); - } - else if(params.testShape == "tri") - { - SLIC_ERROR_IF(params.getBoxDim() != 2, "This example is only in 2D."); - shapeSet = create2DShapeSet(ds); - } - else if(params.testShape == "hex") + std::vector shapesVec; + for(const auto& tg : params.testShape) { - 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); + if(shapeReps.count(tg) == 0) + { + shapeReps[tg] = 0; + } + std::string name = axom::fmt::format("{}.{}", tg, shapeReps[tg]++); + + if(tg == "plane") + { + shapesVec.push_back(createShape_Plane(name)); + } + else if(tg == "hex") + { + shapesVec.push_back(createShape_Hex(name)); + } + else if(tg == "sphere") + { + shapesVec.push_back(createShape_Sphere(name)); + } + else if(tg == "tetmesh") + { + shapesVec.push_back(createShape_TetMesh(ds, name)); + } + else if(tg == "tet") + { + shapesVec.push_back(createShape_Tet(name)); + } + else if(tg == "sor") + { + shapesVec.push_back(createShape_Sor(name)); + } + else if(tg == "cyl") + { + shapesVec.push_back(createShape_Cylinder(name)); + } + else if(tg == "cone") + { + shapesVec.push_back(createShape_Cone(name)); + } + else if(tg == "tri") + { + SLIC_ERROR_IF(params.getBoxDim() != 2, "This example is only in 2D."); + shapesVec = create2DShapeSet(ds); + } } + axom::klee::ShapeSet shapeSet; + + shapeSet.setShapes(shapesVec); + shapeSet.setDimensions(params.getBoxDim() == 2 ? axom::klee::Dimensions::Two + : axom::klee::Dimensions::Three); // Save the discrete shapes for viz and testing. auto* shapeMeshGroup = ds.getRoot()->createGroup("shapeMeshGroup"); @@ -1794,16 +1988,19 @@ int main(int argc, char** argv) shape.getMaterial(), shapeFormat))); + const auto annotationName = "shaping:" + shape.getName(); + AXOM_ANNOTATE_BEGIN(annotationName); // Load the shape from file. This also applies any transformations. shaper->loadShape(shape); - slic::flushStreams(); + // slic::flushStreams(); // Generate a spatial index over the shape shaper->prepareShapeQuery(shapeSet.getDimensions(), shape); - slic::flushStreams(); + // slic::flushStreams(); // Query the mesh against this shape shaper->runShapeQuery(shape); + AXOM_ANNOTATE_END(annotationName); slic::flushStreams(); // Apply the replacement rules for this shape against the existing materials @@ -1945,43 +2142,49 @@ int main(int argc, char** argv) // shape mesh for closes shape. As long as the shapes don't overlap, this // should be a good correctness check. //--------------------------------------------------------------------------- - auto* meshVerificationGroup = ds.getRoot()->createGroup("meshVerification"); for(const auto& shape : shapeSet.getShapes()) { - axom::quest::DiscreteShape dShape(shape, meshVerificationGroup); - auto shapeMesh = - std::dynamic_pointer_cast>( - dShape.createMeshRepresentation()); - double shapeMeshVol = - params.getBoxDim() == 3 ? volumeOfTetMesh(*shapeMesh) : areaOfTriMesh(*shapeMesh); - SLIC_INFO(axom::fmt::format("{:-^80}", - axom::fmt::format("Shape '{}' discrete geometry has {} cells", - shape.getName(), - shapeMesh->getNumberOfCells()))); - - const std::string& materialName = shape.getMaterial(); - double shapeVol = -1; - if(params.useBlueprintSidre() || params.useBlueprintConduit()) - { - shapeVol = sumMaterialVolumes(compMeshGrp, materialName); - } + std::string fieldName = "shape_vol_frac_" + shape.getName(); + axom::ArrayView vfView = getFieldAsArrayView(fieldName); + axom::Array vfHostArray(vfView, axom::execution_space::allocatorID()); + vfView = vfHostArray.view(); + axom::sidre::View* elementVolsVu = nullptr; #if defined(AXOM_USE_MFEM) if(params.useMfem()) { - shapeVol = sumMaterialVolumes(shapingDC.get(), materialName); + elementVolsVu = getElementVolumes(shapingDC.get(), "elementVolumes"); } + else + { + elementVolsVu = getElementVolumes(compMeshGrp, "elementVolumes"); + } +#else + elementVolsVu = getElementVolumes(compMeshGrp, "elementVolumes"); #endif - double correctShapeVol = params.testShape == "plane" ? params.boxMeshVolume() / 2 : shapeMeshVol; + axom::ArrayView elementVolsView(elementVolsVu->getData(), + elementVolsVu->getNumElements()); + axom::Array elementVols(elementVolsView, hostAllocId); + elementVolsView = elementVols.view(); + using ReducePolicy = typename axom::execution_space::reduce_policy; + RAJA::ReduceSum shapedVolume(0); + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType i) { shapedVolume += vfView[i] * elementVolsView[i]; }); + double shapeVol = shapedVolume.get(); + double correctShapeVol = exactOverlapVols.at(shape.getName()); SLIC_ASSERT(correctShapeVol > 0.0); // Indicates error in the test setup. double diff = shapeVol - correctShapeVol; - bool err = !axom::utilities::isNearlyEqualRelative(shapeVol, correctShapeVol, 1e-6, 1e-8); + bool err = !axom::utilities::isNearlyEqualRelative(shapeVol, + correctShapeVol, + errorToleranceRel.at(shape.getName()), + errorToleranceAbs.at(shape.getName())); failCounts += err; SLIC_INFO(axom::fmt::format( "{:-^80}", axom::fmt::format("Material '{}' in shape '{}' has volume {} vs {}, diff of {}, {}.", - materialName, + shape.getMaterial(), shape.getName(), shapeVol, correctShapeVol, @@ -2022,6 +2225,10 @@ int main(int argc, char** argv) SLIC_INFO(axom::fmt::format("{:-^80}", "")); slic::flushStreams(); + AXOM_ANNOTATE_END("quest shaping example"); + + SLIC_INFO(axom::fmt::format("exiting with failure count {}", failCounts)); + finalizeLogger(); return failCounts; diff --git a/src/axom/quest/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index 14ad1e8188..ffa100f1e9 100644 --- a/src/axom/quest/tests/CMakeLists.txt +++ b/src/axom/quest/tests/CMakeLists.txt @@ -254,9 +254,11 @@ if(CONDUIT_FOUND AND RAJA_FOUND AND AXOM_ENABLE_SIDRE) blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) endif() - foreach(_meshType "bpSidre" "bpConduit") + set(_testgeoms "plane" "tet" "plane,tet") + set(_meshTypes "bpSidre" "bpConduit") + foreach(_meshType ${_meshTypes}) foreach(_policy ${_policies}) - foreach(_testgeom "tet") + foreach(_testgeom ${_testgeoms}) set(_testname "quest_mesh_clipper_${_meshType}_${_policy}_${_testgeom}") axom_add_test( NAME ${_testname} diff --git a/src/axom/quest/tests/quest_discretize.cpp b/src/axom/quest/tests/quest_discretize.cpp index f2d9ff1bfb..c232691abd 100644 --- a/src/axom/quest/tests/quest_discretize.cpp +++ b/src/axom/quest/tests/quest_discretize.cpp @@ -167,7 +167,7 @@ void degenerate_segment_test(const char* label, axom::Array& polyline, SCOPED_TRACE(label); - axom::Array generated; + axom::Array generated(0, 0, axom::execution_space::allocatorID()); if(expsuccess) { EXPECT_TRUE(axom::quest::discretize(polyline, len, gens, generated, octcount)); @@ -250,7 +250,7 @@ void segment_test(const char* label, axom::Array& polyline, int len) axom::Array handcut; discretized_segment(polyline[0], polyline[1], handcut); - axom::Array generatedDevice; + axom::Array generatedDevice(0, 0, axom::execution_space::allocatorID()); int octcount = 0; axom::quest::discretize(polyline, len, generations, generatedDevice, octcount); @@ -325,7 +325,7 @@ void multi_segment_test(const char* label, axom::Array& polyline, int l int generation = 0; - axom::Array generatedDevice; + axom::Array generatedDevice(0, 0, axom::execution_space::allocatorID()); int octcount = 0; axom::quest::discretize(polyline, len, generations, generatedDevice, octcount); diff --git a/src/axom/quest/tests/quest_mesh_clipper.cpp b/src/axom/quest/tests/quest_mesh_clipper.cpp index 1daaba635c..c1d0476f1f 100644 --- a/src/axom/quest/tests/quest_mesh_clipper.cpp +++ b/src/axom/quest/tests/quest_mesh_clipper.cpp @@ -27,6 +27,7 @@ #include "axom/sidre.hpp" #include "axom/klee.hpp" #include "axom/quest.hpp" +#include "axom/quest/detail/clipping/Plane3DClipper.hpp" #include "axom/quest/detail/clipping/TetClipper.hpp" #include "axom/fmt.hpp" @@ -101,7 +102,7 @@ struct Input // The shape to run. std::vector testGeom; // The shapes this example is set up to run. - const std::set availableShapes {"tet"}; // More geometries to come. + const std::set availableShapes {"tet", "plane"}; RuntimePolicy policy {RuntimePolicy::seq}; int refinementLevel {7}; @@ -110,6 +111,8 @@ struct Input std::string backgroundMaterial; + int screenLevel = -1; + // clang-format off enum class MeshType { bpSidre = 0, bpConduit = 1 }; const std::map meshTypeChoices @@ -137,10 +140,14 @@ struct Input void parse(int argc, char** argv, axom::CLI::App& app) { + app.add_option("--screenLevel", screenLevel) + ->description("Developer feature for MeshClipper.") + ->capture_default_str(); + app.add_option("-o,--outputFile", outputFile)->description("Path to output file(s)"); app.add_flag("-v,--verbose,!--no-verbose", m_verboseOutput) - ->description("Enable/disable verbose output") + ->description("Enable/disable verbose output, including SLIC_DEBUG") ->capture_default_str(); app.add_option("--meshType", meshType) @@ -270,7 +277,8 @@ const std::string matsetName = "matset"; const std::string coordsetName = "coords"; int cellCount = -1; // Translation to individual octants (override) when running multiple shapes. -// Except that the plane is never moved. +// Exception: the plane always placed at origin to facilitate finding its +// exact overlap volume. const double tDist = 0.9; // Bias toward origin to help keep shape inside domain. std::vector> translations {{tDist, tDist, -tDist}, {-tDist, tDist, -tDist}, @@ -387,6 +395,7 @@ void initializeLogger() #ifdef AXOM_USE_MPI int num_ranks = 1; MPI_Comm_size(MPI_COMM_WORLD, &num_ranks); + SLIC_ERROR_IF(num_ranks > 1, "Sorry, this test is serial."); if(num_ranks > 1) { std::string fmt = "[][]: \n"; @@ -417,6 +426,32 @@ void finalizeLogger() } } +double volumeOfTetMesh(const axom::mint::UnstructuredMesh& tetMesh) +{ + using TetType = axom::primal::Tetrahedron; + axom::StackArray nodesShape {tetMesh.getNumberOfNodes()}; + axom::ArrayView x(tetMesh.getCoordinateArray(0), nodesShape); + axom::ArrayView y(tetMesh.getCoordinateArray(1), nodesShape); + axom::ArrayView z(tetMesh.getCoordinateArray(2), nodesShape); + const axom::IndexType tetCount = tetMesh.getNumberOfCells(); + axom::Array tetVolumes(tetCount, tetCount); + double meshVolume = 0.0; + for(axom::IndexType ic = 0; ic < tetCount; ++ic) + { + const axom::IndexType* nodeIds = tetMesh.getCellNodeIDs(ic); + TetType tet; + for(int j = 0; j < 4; ++j) + { + auto cornerNodeId = nodeIds[j]; + tet[j][0] = x[cornerNodeId]; + tet[j][1] = y[cornerNodeId]; + tet[j][2] = z[cornerNodeId]; + } + meshVolume += tet.volume(); + } + return meshVolume; +} + axom::klee::Geometry createGeom_Tet(const std::string& geomName) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, @@ -443,6 +478,36 @@ axom::klee::Geometry createGeom_Tet(const std::string& geomName) return tetGeometry; } +axom::klee::Geometry createGeom_Plane(const std::string& geomName) +{ + axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + // Create a plane crossing center of mesh. No matter the normal, + // it cuts the mesh in half. + Point3D center {0.5 * + (axom::NumericArray(params.boxMins.data()) + + axom::NumericArray(params.boxMaxs.data()))}; + primal::Vector normal = params.direction.empty() + ? primal::Vector3D {1.0, 0.0, 0.0} + : primal::Vector3D {params.direction.data()}.unitVector(); + const primal::Plane plane {normal, center, true}; + + axom::klee::Geometry planeGeometry(prop, plane, {nullptr}); + + // Exact mesh overlap volume, assuming plane passes through center of box mesh. + using Pt3D = primal::Point; + Pt3D lower(params.boxMins.data()); + Pt3D upper(params.boxMaxs.data()); + auto diag = upper.array() - lower.array(); + double meshVolume = diag[0] * diag[1] * diag[2]; + exactGeomVols[geomName] = 0.5 * meshVolume; + errorToleranceRel[geomName] = 1e-6; + errorToleranceAbs[geomName] = 1e-8; + + return planeGeometry; +} + /*! @brief Return the element volumes as a sidre::View containing the volumes in an array. @@ -823,12 +888,16 @@ int main(int argc, char** argv) } std::string name = axom::fmt::format("{}.{}", tg, geomReps[tg]++); - if(tg == "tet") + if(tg == "plane") + { + geomStrategies.push_back( + std::make_shared(createGeom_Plane(name), name)); + } + else if(tg == "tet") { geomStrategies.push_back( std::make_shared(createGeom_Tet(name), name)); } - // More geometries to come. } { @@ -900,11 +969,19 @@ int main(int argc, char** argv) quest::experimental::MeshClipper clipper(sMesh, geomStrategies[i]); clipper.setVerbose(params.isVerbose()); + if(params.screenLevel >= 0) + { + clipper.setScreenLevel(params.screenLevel); + } + SLIC_INFO(axom::fmt::format("MeshClipper screen level: {}", clipper.getScreenLevel())); + axom::Array ovlap; AXOM_ANNOTATE_BEGIN(annotationName); clipper.clip(ovlap); AXOM_ANNOTATE_END(annotationName); + clipper.logClippingStats(); + // Save volume fractions in mesh, for plotting and checking. sMesh.setMatsetFromVolume(geomStrategies[i]->name(), ovlap.view(), false); diff --git a/src/axom/quest/util/mesh_helpers.cpp b/src/axom/quest/util/mesh_helpers.cpp index 8e4209341f..4b763f1d83 100644 --- a/src/axom/quest/util/mesh_helpers.cpp +++ b/src/axom/quest/util/mesh_helpers.cpp @@ -239,30 +239,41 @@ axom::sidre::Group* make_unstructured_blueprint_box_mesh_2d(axom::sidre::Group* void convert_blueprint_structured_explicit_to_unstructured_3d(axom::sidre::Group* meshGrp, const std::string& topoName, - axom::runtime_policy::Policy runtimePolicy) + axom::runtime_policy::Policy runtimePolicy, + const std::string& ugTopoName) { if(runtimePolicy == axom::runtime_policy::Policy::seq) { - convert_blueprint_structured_explicit_to_unstructured_impl_3d(meshGrp, topoName); + convert_blueprint_structured_explicit_to_unstructured_3d_impl( + meshGrp, + topoName, + ugTopoName.empty() ? topoName : ugTopoName); } #if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) if(runtimePolicy == axom::runtime_policy::Policy::omp) { - convert_blueprint_structured_explicit_to_unstructured_impl_3d(meshGrp, topoName); + convert_blueprint_structured_explicit_to_unstructured_3d_impl( + meshGrp, + topoName, + ugTopoName.empty() ? topoName : ugTopoName); } #endif #if defined(AXOM_RUNTIME_POLICY_USE_CUDA) if(runtimePolicy == axom::runtime_policy::Policy::cuda) { - convert_blueprint_structured_explicit_to_unstructured_impl_3d>(meshGrp, - topoName); + convert_blueprint_structured_explicit_to_unstructured_3d_impl>( + meshGrp, + topoName, + ugTopoName.empty() ? topoName : ugTopoName); } #endif #if defined(AXOM_RUNTIME_POLICY_USE_HIP) if(runtimePolicy == axom::runtime_policy::Policy::hip) { - convert_blueprint_structured_explicit_to_unstructured_impl_3d>(meshGrp, - topoName); + convert_blueprint_structured_explicit_to_unstructured_3d_impl>( + meshGrp, + topoName, + ugTopoName.empty() ? topoName : ugTopoName); } #endif } @@ -273,57 +284,73 @@ void convert_blueprint_structured_explicit_to_unstructured_2d(axom::sidre::Group { if(runtimePolicy == axom::runtime_policy::Policy::seq) { - convert_blueprint_structured_explicit_to_unstructured_impl_2d(meshGrp, topoName); + convert_blueprint_structured_explicit_to_unstructured_2d_impl(meshGrp, topoName); } #if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) if(runtimePolicy == axom::runtime_policy::Policy::omp) { - convert_blueprint_structured_explicit_to_unstructured_impl_2d(meshGrp, topoName); + convert_blueprint_structured_explicit_to_unstructured_2d_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_2d>(meshGrp, + convert_blueprint_structured_explicit_to_unstructured_2d_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_2d>(meshGrp, + convert_blueprint_structured_explicit_to_unstructured_2d_impl>(meshGrp, topoName); } #endif } template -void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre::Group* meshGrp, - const std::string& topoName) +void convert_blueprint_structured_explicit_to_unstructured_3d_impl(axom::sidre::Group* meshGrp, + const std::string& topoName, + const std::string& ugTopoName) { AXOM_ANNOTATE_SCOPE("convert_to_unstructured"); // Note: MSVC required `static` to pass DIM to the axom::for_all w/ C++14 // this restriction might be lifted w/ C++17 static constexpr int DIM = 3; - const std::string& coordsetName = - meshGrp->getView("topologies/" + topoName + "/coordset")->getString(); + sidre::Group* topologiesGrp = meshGrp->getGroup("topologies"); + axom::sidre::Group* topoGrp = topologiesGrp->getGroup(topoName); + axom::sidre::Group* ugTopoGrp = + ugTopoName == topoName ? topoGrp : topologiesGrp->createGroup(ugTopoName); + + const std::string& coordsetName = topoGrp->getView("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"); + if(!ugTopoGrp->hasView("coordset")) + { + ugTopoGrp->createViewString("coordset", coordsetName); + } + else + { + SLIC_ASSERT(ugTopoGrp->getView("coordset")->isString()); + SLIC_ASSERT(std::string(ugTopoGrp->getView("coordset")->getString()) == coordsetName); + } + + SLIC_ASSERT(std::string(topoGrp->getView("type")->getString()) == "structured"); + axom::sidre::View* ugTopoTypeView = + ugTopoGrp == topoGrp ? ugTopoGrp->getView("type") : ugTopoGrp->createView("type"); + ugTopoTypeView->setString("unstructured"); + axom::sidre::View* shapeView = ugTopoGrp->createView("elements/shape"); + shapeView->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")); + SLIC_ASSERT(!topoElemGrp->hasGroup("offsets")); + SLIC_ASSERT(!topoElemGrp->hasGroup("strides")); const axom::StackArray cShape { axom::IndexType(topoDimsGrp->getView("i")->getNode().to_value()), @@ -335,10 +362,10 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre:: const axom::StackArray connShape {cCount, 8}; axom::sidre::View* connView = - topoGrp->createViewWithShapeAndAllocate("elements/connectivity", - axom::sidre::detail::SidreTT::id, - 2, - connShape.begin()); + ugTopoGrp->createViewWithShapeAndAllocate("elements/connectivity", + axom::sidre::detail::SidreTT::id, + 2, + connShape.begin()); axom::ArrayView connArrayView( static_cast(connView->getVoidPtr()), connShape); @@ -370,15 +397,15 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre:: { 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. + * 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. + * 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; @@ -391,7 +418,7 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre:: valuesGrp->getView("z")->reshapeArray(2, vertsShape); // Make connectivity array 2D for the same reason. - auto* elementsGrp = topoGrp->getGroup("elements"); + auto* elementsGrp = ugTopoGrp->getGroup("elements"); auto* connView = elementsGrp->getView("connectivity"); curDim = connView->getShape(2, curShape); constexpr axom::IndexType NUM_VERTS_PER_HEX = 8; @@ -406,8 +433,6 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre:: // 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"); } @@ -423,7 +448,7 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre:: } template -void convert_blueprint_structured_explicit_to_unstructured_impl_2d(axom::sidre::Group* meshGrp, +void convert_blueprint_structured_explicit_to_unstructured_2d_impl(axom::sidre::Group* meshGrp, const std::string& topoName) { AXOM_ANNOTATE_SCOPE("convert_to_unstructured"); @@ -494,15 +519,15 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_2d(axom::sidre:: { AXOM_ANNOTATE_SCOPE("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. + * 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. + * 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] = {}; auto* valuesGrp = coordsetGrp->getGroup("values"); @@ -525,9 +550,6 @@ void convert_blueprint_structured_explicit_to_unstructured_impl_2d(axom::sidre:: // mint::Mesh requires connectivity strides, even though Blueprint doesn't. constexpr axom::IndexType BIT_SPECIFIC_NUM_VERTS_PER_QUAD = 4; elementsGrp->createViewScalar("stride", BIT_SPECIFIC_NUM_VERTS_PER_QUAD); - - // mint::Mesh requires field group, even though Blueprint doesn't. - meshGrp->createGroup("fields"); } return; diff --git a/src/axom/quest/util/mesh_helpers.hpp b/src/axom/quest/util/mesh_helpers.hpp index 6e367402d3..b9ec2ccd24 100644 --- a/src/axom/quest/util/mesh_helpers.hpp +++ b/src/axom/quest/util/mesh_helpers.hpp @@ -152,6 +152,8 @@ axom::sidre::Group* make_unstructured_blueprint_box_mesh_2d( * \param runtimePolicy Runtime policy, see axom::runtime_policy. * Memory in \c meshGrp must be compatible with the * specified policy. + * \param ugTopoName Name of the unstructured topology to create. + * Defaults to topoName. * * All input mesh data are expected to have the allocator id of * meshGrp->getDefaultAllocatorID(). On output, they will also have @@ -160,18 +162,20 @@ axom::sidre::Group* make_unstructured_blueprint_box_mesh_2d( */ void convert_blueprint_structured_explicit_to_unstructured_3d(axom::sidre::Group* meshGrp, const std::string& topoName, - axom::runtime_policy::Policy runtimePolicy); + axom::runtime_policy::Policy runtimePolicy, + const std::string& ugTopoName = ""); template -void convert_blueprint_structured_explicit_to_unstructured_impl_3d(axom::sidre::Group* meshGrp, - const std::string& topoName); +void convert_blueprint_structured_explicit_to_unstructured_3d_impl(axom::sidre::Group* meshGrp, + const std::string& topoName, + const std::string& ugTopoName); void convert_blueprint_structured_explicit_to_unstructured_2d(axom::sidre::Group* meshGrp, const std::string& topoName, axom::runtime_policy::Policy runtimePolicy); template -void convert_blueprint_structured_explicit_to_unstructured_impl_2d(axom::sidre::Group* meshGrp, +void convert_blueprint_structured_explicit_to_unstructured_2d_impl(axom::sidre::Group* meshGrp, const std::string& topoName); #if defined(AXOM_USE_CONDUIT) diff --git a/src/axom/sidre/core/Group.cpp b/src/axom/sidre/core/Group.cpp index 6fab301f48..0dc9bae329 100644 --- a/src/axom/sidre/core/Group.cpp +++ b/src/axom/sidre/core/Group.cpp @@ -1476,6 +1476,15 @@ bool Group::createNativeLayout(Node& n, const Attribute* attr) const */ bool Group::deepCopyToConduit(Node& dst, int tupleAllocId, int arrayAllocId, const Attribute* attr) const { + if(tupleAllocId == INVALID_ALLOCATOR_ID) + { + tupleAllocId = ConduitMemory::conduitAllocIdToAxom(dst.allocator()); + } + if(arrayAllocId == INVALID_ALLOCATOR_ID) + { + arrayAllocId = ConduitMemory::conduitAllocIdToAxom(dst.allocator()); + } + dst.set(DataType::object()); bool hasSavedViews = false; diff --git a/src/axom/sidre/core/View.cpp b/src/axom/sidre/core/View.cpp index 02306c7b6a..917a01573a 100644 --- a/src/axom/sidre/core/View.cpp +++ b/src/axom/sidre/core/View.cpp @@ -1036,7 +1036,7 @@ void View::createNativeLayout(Node& n) const * ************************************************************************* */ -void View::deepCopyToConduit(Node& dst) const +void View::deepCopyToConduit(Node& dst, int allocId) const { // see ATK-726 - Handle undescribed and unallocated views in Sidre's // createNativeLayout() @@ -1046,6 +1046,10 @@ void View::deepCopyToConduit(Node& dst) const const conduit::DataType& srcDtype = m_node.dtype(); dst.set(srcDtype); + if(allocId != INVALID_ALLOCATOR_ID) + { + dst.set_allocator(ConduitMemory::axomAllocIdToConduit(allocId)); + } if(isAllocated()) { // Using set_node to set dst: would reset dst's allocator id to diff --git a/src/axom/sidre/core/View.hpp b/src/axom/sidre/core/View.hpp index f32dc06c4f..cf7163d337 100644 --- a/src/axom/sidre/core/View.hpp +++ b/src/axom/sidre/core/View.hpp @@ -1000,7 +1000,7 @@ class View /*! * \brief Deep copy View into the given conduit::Node. */ - void deepCopyToConduit(Node& dst) const; + void deepCopyToConduit(Node& dst, int allocId = INVALID_ALLOCATOR_ID) const; /*! * \brief Copy metadata of the View to the given Conduit node