diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 611b404389..5aa467b238 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -48,6 +48,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ ### Changed +- Version of `quest::discretize` that approximates a surface-of-revolution from a polyline + now respects the allocator ID of the output `Array`. It no longer resets the ID to the + execution space default. - Evaluation methods for line integrals in `axom::primal` have been generalized, and `evaluate_scalar_line_integral` has been renamed to `evaluate_line_integral`. - Treatment of materials on strided-structured Blueprint meshes has changed in `axom::mir`. diff --git a/src/axom/primal/geometry/CoordinateTransformer.hpp b/src/axom/primal/geometry/CoordinateTransformer.hpp index d59304b2e8..bc0bd23e48 100644 --- a/src/axom/primal/geometry/CoordinateTransformer.hpp +++ b/src/axom/primal/geometry/CoordinateTransformer.hpp @@ -8,6 +8,7 @@ #include "axom/core/numerics/Matrix.hpp" #include "axom/core/utilities/Utilities.hpp" +#include "axom/core/NumericLimits.hpp" #include "axom/core/numerics/floating_point_limits.hpp" #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" @@ -157,7 +158,7 @@ class CoordinateTransformer * * Validity can be checked with isValid(). */ - AXOM_HOST_DEVICE void setInvalid() { m_P[0][0] = std::numeric_limits::quiet_NaN(); } + AXOM_HOST_DEVICE void setInvalid() { m_P[0][0] = axom::numeric_limits::quiet_NaN(); } //! @brief Whether transformer is valid. AXOM_HOST_DEVICE bool isValid() { return !std::isnan(m_P[0][0]); } @@ -377,7 +378,7 @@ class CoordinateTransformer * M = [ P v ] * [ 0 1 ] * - * Store m_P[0][0] = std::numeric_limits::quiet_NaN() to set this + * Store m_P[0][0] = axom::numeric_limits::quiet_NaN() to set this * CoordinateTransformer as invalid. See \a setInvalid() */ Matrx m_P; @@ -411,7 +412,7 @@ class CoordinateTransformer } else { - m[0][0] = std::numeric_limits::quiet_NaN(); + m[0][0] = axom::numeric_limits::quiet_NaN(); } } diff --git a/src/axom/quest/CMakeLists.txt b/src/axom/quest/CMakeLists.txt index a5420ffbd6..ba5fbf70ca 100644 --- a/src/axom/quest/CMakeLists.txt +++ b/src/axom/quest/CMakeLists.txt @@ -154,6 +154,8 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE) detail/clipping/TetClipper.hpp detail/clipping/HexClipper.hpp detail/clipping/SphereClipper.hpp + detail/clipping/MonotonicZSORClipper.hpp + detail/clipping/SORClipper.hpp detail/clipping/MeshClipperImpl.hpp) list(APPEND quest_sources MeshClipper.cpp MeshClipperStrategy.cpp @@ -161,6 +163,8 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE) detail/clipping/TetClipper.cpp detail/clipping/HexClipper.cpp detail/clipping/SphereClipper.cpp + detail/clipping/MonotonicZSORClipper.cpp + detail/clipping/SORClipper.cpp detail/clipping/MeshClipperImpl.hpp) endif() endif() diff --git a/src/axom/quest/IntersectionShaper.hpp b/src/axom/quest/IntersectionShaper.hpp index e3a626d7d0..0dca5dd4d8 100644 --- a/src/axom/quest/IntersectionShaper.hpp +++ b/src/axom/quest/IntersectionShaper.hpp @@ -824,7 +824,8 @@ class IntersectionShaper : public Shaper axom::fmt::format("{:-^80}", axom::fmt::format(" Refinement level set to {} ", m_level))); // Generate the Octahedra - // (octahedra m_octs will be on device) + // (Set m_octs's allocator id to where we want its data to live.) + 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); diff --git a/src/axom/quest/detail/Discretize_detail.hpp b/src/axom/quest/detail/Discretize_detail.hpp index a1ea0c6815..79f11263b5 100644 --- a/src/axom/quest/detail/Discretize_detail.hpp +++ b/src/axom/quest/detail/Discretize_detail.hpp @@ -7,6 +7,7 @@ #define AXOM_QUEST_DISCRETIZE_DETAIL_ #include "axom/primal/constants.hpp" +#include "math.h" namespace { @@ -152,12 +153,11 @@ int discrSeg(const Point2D &a, const Point2D &b, int levels, axom::ArrayView::allocatorID(); // Assert input assumptions - SLIC_ASSERT(b[0] - a[0] >= 0); SLIC_ASSERT(a[1] >= 0); SLIC_ASSERT(b[1] >= 0); // Deal with degenerate segments - if(b[0] - a[0] < axom::primal::PRIMAL_TINY) + if(fabs(b[0] - a[0]) < axom::primal::PRIMAL_TINY) { return 0; } @@ -256,7 +256,7 @@ namespace quest * less than the polyline's length). That is exponential growth. Use * appropriate caution. * - * This routine initializes an Array pointed to by \a out. + * This routine resizes and populates an Array pointed to by \a out. */ template bool discretize(const axom::ArrayView &polyline, @@ -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; @@ -273,11 +277,6 @@ bool discretize(const axom::ArrayView &polyline, { const Point2D &a = polyline[seg]; const Point2D &b = polyline[seg + 1]; - // invalid if a.x > b.x - if(a[0] > b[0]) - { - stillValid = false; - } if(a[1] < 0 || b[1] < 0) { stillValid = false; @@ -293,7 +292,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.clear(); + out.resize(axom::ArrayOptions::Uninitialized(), totaloctcount); axom::ArrayView out_view = out.view(); octcount = 0; @@ -303,6 +303,8 @@ bool discretize(const axom::ArrayView &polyline, discrSeg(polyline[seg], polyline[seg + 1], levels, out_view, octcount); octcount += segment_prism_count; } + // octcount may be < totaloctcount if there are degenerate segments. + out.resize(octcount); // TODO check for errors in each segment's computation return true; diff --git a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp index 5f1ccf8788..1e5782a59e 100644 --- a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp +++ b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp @@ -109,12 +109,19 @@ class MeshClipperImpl : public MeshClipper::Impl }); } - //! @brief Make a list of indices where labels have value LABEL_ON. + /*! + * @brief Make a list of indices where labels have value LABEL_ON, + * stored in the same allocator id as the labels. + */ void collectOnIndices(const axom::ArrayView& labels, axom::Array& onIndices) override { if(labels.empty()) { + if(onIndices.getAllocatorID() != labels.getAllocatorID()) + { + onIndices = axom::Array(0, 0, labels.getAllocatorID()); + } return; }; diff --git a/src/axom/quest/detail/clipping/MonotonicZSORClipper.cpp b/src/axom/quest/detail/clipping/MonotonicZSORClipper.cpp new file mode 100644 index 0000000000..5daac38a6b --- /dev/null +++ b/src/axom/quest/detail/clipping/MonotonicZSORClipper.cpp @@ -0,0 +1,736 @@ +// 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/core/numerics/matvecops.hpp" +#include "axom/core/utilities/Utilities.hpp" +#include "axom/primal/operators/squared_distance.hpp" +#include "axom/quest/Discretize.hpp" +#include "axom/quest/detail/clipping/MonotonicZSORClipper.hpp" +#include "axom/fmt.hpp" + +#include + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +MonotonicZSORClipper::MonotonicZSORClipper(const klee::Geometry& kGeom, const std::string& name) + : MeshClipperStrategy(kGeom) + , m_name(name.empty() ? std::string("FSor") : name) + , m_maxRadius(0.0) + , m_minRadius(numerics::floating_point_limits::max()) + , m_transformer() +{ + extractClipperInfo(); + + combineRadialSegments(m_sorCurve); + axom::Array turnIndices = findZSwitchbacks(m_sorCurve.view()); + if(turnIndices.size() > 2) + { + // The 2 "turns" allowed are the first and last points. Anything else is a switchback. + SLIC_ERROR( + "MonotonicZSORClipper does not work when a curve doubles back" + " in the axial direction. Use SORClipper instead."); + } + + for(auto& pt : m_sorCurve) + { + m_maxRadius = fmax(m_maxRadius, pt[1]); + m_minRadius = fmin(m_minRadius, pt[1]); + } + SLIC_ERROR_IF(m_minRadius < 0.0, + axom::fmt::format("MonotonicZSORClipper '{}' has a negative radius", m_name)); + + // Combine internal and external rotations into m_transformer. + m_transformer.applyRotation(Vector3DType({1, 0, 0}), m_sorDirection); + m_transformer.applyTranslation(m_sorOrigin.array()); + m_transformer.applyMatrix(m_extTrans); + m_invTransformer = m_transformer.getInverse(); + + for(const auto& pt : m_sorCurve) + { + m_curveBb.addPoint(pt); + } +} + +MonotonicZSORClipper::MonotonicZSORClipper(const klee::Geometry& kGeom, + const std::string& name, + axom::ArrayView discreteFunction, + const Point3DType& sorOrigin, + const Vector3DType& sorDirection, + axom::IndexType levelOfRefinement) + : MeshClipperStrategy(kGeom) + , m_name(name.empty() ? std::string("FSor") : name) + , m_sorCurve(discreteFunction, axom::execution_space::allocatorID()) + , m_maxRadius(0.0) + , m_minRadius(numerics::floating_point_limits::max()) + , m_sorOrigin(sorOrigin) + , m_sorDirection(sorDirection) + , m_levelOfRefinement(levelOfRefinement) + , m_transformer() +{ + combineRadialSegments(m_sorCurve); + axom::Array turnIndices = findZSwitchbacks(m_sorCurve.view()); + if(turnIndices.size() > 2) + { + // The 2 "turns" allowed are the first and last points. Anything else is a switchback. + SLIC_ERROR( + "MonotonicZSORClipper does not work when a curve doubles back" + " in the axial direction. Use SORClipper instead."); + } + + for(auto& pt : m_sorCurve) + { + m_maxRadius = fmax(m_maxRadius, pt[1]); + m_minRadius = fmin(m_minRadius, pt[1]); + } + SLIC_ERROR_IF(m_minRadius < 0.0, + axom::fmt::format("MonotonicZSORClipper '{}' has a negative radius", m_name)); + + // Combine internal and external rotations into m_transformer. + m_transformer.applyRotation(Vector3DType({1, 0, 0}), m_sorDirection); + m_transformer.applyTranslation(m_sorOrigin.array()); + m_transformer.applyMatrix(m_extTrans); + m_invTransformer = m_transformer.getInverse(); + + for(const auto& pt : m_sorCurve) + { + m_curveBb.addPoint(pt); + } +} + +bool MonotonicZSORClipper::labelCellsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& labels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "MonotonicZSORClipper requires a 3D mesh."); + + const int allocId = shapeMesh.getAllocatorID(); + const auto cellCount = shapeMesh.getCellCount(); + if(labels.size() < cellCount || labels.getAllocatorID() != allocId) + { + 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 MonotonicZSORClipper::labelTetsInOut(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::Array& tetLabels) +{ + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "MonotonicZSORClipper 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()) + { + 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; +} + +/* + * Implementation: (reverse) transform the mesh vertices to the r-z + * frame where the curve is defined as a r(z) function. It's easier to + * determine whether the point is in the sor that way. + */ +template +void MonotonicZSORClipper::labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView labels) +{ + axom::Array bbOn; + axom::Array bbUnder; + computeCurveBoxes(shapeMesh, bbOn, bbUnder); + const axom::ArrayView bbOnView = bbOn.view(); + const axom::ArrayView bbUnderView = bbUnder.view(); + + const auto cellCount = shapeMesh.getCellCount(); + auto meshHexes = shapeMesh.getCellsAsHexes(); + auto meshCellVolumes = shapeMesh.getCellVolumes(); + auto invTransformer = m_invTransformer; + constexpr double EPS = 1e-10; + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType cellId) { + if(axom::utilities::isNearlyEqual(meshCellVolumes[cellId], 0.0, EPS)) + { + labels[cellId] = LabelType::LABEL_OUT; + return; + } + auto cellHex = meshHexes[cellId]; + for(int vi = 0; vi < HexahedronType::NUM_HEX_VERTS; ++vi) + { + invTransformer.transform(cellHex[vi].array()); + } + BoundingBox2DType cellBbInRz = estimateBoundingBoxInRz(cellHex); + labels[cellId] = rzBbToLabel(cellBbInRz, bbOnView, bbUnderView); + }); +} + +template +void MonotonicZSORClipper::labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::ArrayView labels) +{ + axom::Array bbOn; + axom::Array bbUnder; + computeCurveBoxes(shapeMesh, bbOn, bbUnder); + const axom::ArrayView bbOnView = bbOn.view(); + const axom::ArrayView bbUnderView = bbUnder.view(); + + const auto cellCount = cellIds.size(); + auto meshHexes = shapeMesh.getCellsAsHexes(); + auto tetVolumes = shapeMesh.getTetVolumes(); + auto invTransformer = m_invTransformer; + constexpr double EPS = 1e-10; + + axom::for_all( + cellCount, + AXOM_LAMBDA(axom::IndexType ci) { + axom::IndexType cellId = cellIds[ci]; + + HexahedronType hex = meshHexes[cellId]; + for(int vi = 0; vi < HexahedronType::NUM_HEX_VERTS; ++vi) + { + invTransformer.transform(hex[vi].array()); + } + + TetrahedronType cellTets[NUM_TETS_PER_HEX]; + ShapeMesh::hexToTets(hex, cellTets); + + for(IndexType ti = 0; ti < NUM_TETS_PER_HEX; ++ti) + { + axom::IndexType tetId = cellId * NUM_TETS_PER_HEX + ti; + LabelType& tetLabel = labels[ci * NUM_TETS_PER_HEX + ti]; + if(axom::utilities::isNearlyEqual(tetVolumes[tetId], 0.0, EPS)) + { + tetLabel = LabelType::LABEL_OUT; + continue; + } + const TetrahedronType& tet = cellTets[ti]; + BoundingBox2DType bbInRz = estimateBoundingBoxInRz(tet); + tetLabel = rzBbToLabel(bbInRz, bbOnView, bbUnderView); + } + }); +} + +/* + Compute bounding box in rz space for a tet or hex geometry in + body frame (the 3D frame with the rotation along +x). + + 1. Rotate the tet or hex vertices into the rz plane. + 2. Compute bounding box for vertices. + 3. Expand 2D bounding box to contain edge that may + intersect SOR between vertices. +*/ +template +AXOM_HOST_DEVICE MonotonicZSORClipper::BoundingBox2DType MonotonicZSORClipper::estimateBoundingBoxInRz( + const PolyhedronType& vertices) +{ + MonotonicZSORClipper::BoundingBox2DType bbInRz; + + // Range of vertex angles in cylindrical coordinates. + primal::BoundingBox angleRange; + + for(IndexType vi = 0; vi < vertices.numVertices(); ++vi) + { + auto& vert = vertices[vi]; + Point2DType vertOnXPlane {vert[1], vert[2]}; + Point2DType vertOnRz { + vert[0], + std::sqrt(numerics::dot_product(vertOnXPlane.data(), vertOnXPlane.data(), 2))}; + bbInRz.addPoint(vertOnRz); + + double angle = atan2(vertOnXPlane[1], vertOnXPlane[0]); + angleRange.addPoint(primal::Point {angle}); + } + /* + The geometry can be closer to the rotation axis than its + individual vertices are, depending on the angle (about the axis) + between the vertices. Given the angle, extend the bottom of bbInRz + for the worst case. + */ + auto angleDiff = angleRange.range()[0]; + double factor = angleDiff > M_PI ? 0.0 : cos(angleDiff / 2); + auto newMin = bbInRz.getMin(); + newMin[1] *= factor; + bbInRz.addPoint(newMin); + return bbInRz; +} + +/* + Compute label based on a bounding box in rz space. + + - If bbInRz is close to any bbOn, label it ON. + - Else if bbInRz touches any bbUnder, label it IN. + It cannot possibly be partially outside, because it + doesn't cross the boundary or even touch any bbOn. + - Else, label bbInRz OUT. + + We expect bbOn and bbUnder to be small arrays, so we use + linear searches. If that's too slow, we can use a BVH. +*/ +AXOM_HOST_DEVICE inline MeshClipperStrategy::LabelType MonotonicZSORClipper::rzBbToLabel( + const BoundingBox2DType& bbInRz, + const axom::ArrayView& bbOn, + const axom::ArrayView& bbUnder) +{ + LabelType label = LabelType::LABEL_OUT; + + for(const auto& bbOn : bbOn) + { + double sqDist = axom::primal::squared_distance(bbInRz, bbOn); + if(sqDist <= 0.0) + { + label = LabelType::LABEL_ON; + } + } + + if(label == LabelType::LABEL_OUT) + { + for(const auto& bbUnder : bbUnder) + { + if(bbInRz.intersectsWith(bbUnder)) + { + label = LabelType::LABEL_IN; + } + } + } + + return label; +} + +/* +*/ +template +void MonotonicZSORClipper::computeCurveBoxes(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& bbOn, + axom::Array& bbUnder) +{ + /* + * Compute bounding boxes bbOn, which cover the curve segments, and + * bbUnder, which cover the space between bbOn and the z axis. bbOn + * includes end caps, the segments that join the curve to the + * z-axis. + */ + const int allocId = shapeMesh.getAllocatorID(); + const IndexType cellCount = shapeMesh.getCellCount(); + + axom::ArrayView cellLengths = shapeMesh.getCellLengths(); + + using ReducePolicy = typename axom::execution_space::reduce_policy; + using LoopPolicy = typename execution_space::loop_policy; + RAJA::ReduceSum sumCharLength(0.0); + RAJA::forall( + RAJA::RangeSegment(0, cellCount), + AXOM_LAMBDA(axom::IndexType cellId) { sumCharLength += cellLengths[cellId]; }); + double avgCharLength = sumCharLength.get() / cellCount; + + /* + Subdivide the SOR curve and place it with the correct allocator. + Create temporary sorCurve that is equivalent to m_sorCurve but + - with long segments subdivided into subsegments based on + characteristic length of mesh cells. + - with memory from allocId. + */ + axom::Array sorCurve = subdivideCurve(m_sorCurve, + 3 * avgCharLength /* maxMean */, + -1 /* maxDz, negative disables */, + -1 /* minDz, negative disables */); + sorCurve = axom::Array(sorCurve, allocId); + auto sorCurveView = sorCurve.view(); + + /* + Compute 2 sets of boxes. + - bbOn have boxes over each segment. + - bbUnder have boxes from the z axis to the bottom of bbOn. + Add add to bbOn boxes representing the vertical endcaps of the curve. + */ + auto segCount = sorCurve.size() - 1; + bbOn = axom::Array(segCount + 2, segCount + 2, allocId); + bbUnder = axom::Array(segCount, segCount, allocId); + auto bbOnView = bbOn.view(); + auto bbUnderView = bbUnder.view(); + + axom::for_all( + segCount, + AXOM_LAMBDA(axom::IndexType i) { + BoundingBox2DType& on = bbOnView[i]; + BoundingBox2DType& under = bbUnderView[i]; + on.addPoint(sorCurveView[i]); + on.addPoint(sorCurveView[i + 1]); + Point2DType underMin {on.getMin()[0], 0.0}; + Point2DType underMax {on.getMax()[0], on.getMin()[1]}; + under = BoundingBox2DType(underMin, underMax); + }); + + axom::Array endCaps(2, 2); + endCaps[0].addPoint(m_sorCurve.front()); + endCaps[0].addPoint(Point2DType {m_sorCurve.front()[0], 0.0}); + endCaps[1].addPoint(m_sorCurve.back()); + endCaps[1].addPoint(Point2DType {m_sorCurve.back()[0], 0.0}); + axom::copy(&bbOn[segCount], endCaps.data(), endCaps.size() * sizeof(BoundingBox2DType)); +} + +/* + * Replace SOR curve segments that have bounding boxes that overlap + * too much beyond what the segments actually overlap. + * + * Goal: Split up segments with excessively large bounding boxes, + * which reach too far beyond the SOR curve. These are long diagonal + * segments. But don't split up segments aligned close to z or r + * directions, because they don't have excessively large bounding + * boxes for their size. We do this by limiting the harmonic mean of + * the r and z sides of the bounding boxes. + */ +Array MonotonicZSORClipper::subdivideCurve( + const Array& sorCurveIn, + double maxMean, + double maxDz, + double minDz) +{ + Array sorCurveOut; + + if(sorCurveIn.empty()) + { + return sorCurveOut; + } + + // Reserve guessed total number of points needed + sorCurveOut.reserve(sorCurveIn.size() * 1.2 + 10); + sorCurveOut.push_back(sorCurveIn[0]); + + for(IndexType i = 1; i < sorCurveIn.size(); ++i) + { + const Point2DType& segStart = sorCurveIn[i - 1]; + const Point2DType& segEnd = sorCurveIn[i]; + + const auto delta = segEnd.array() - segStart.array(); + const auto absDelta = axom::abs(delta); + const double segDz = absDelta[0]; + const double segDr = absDelta[1]; + const double segMean = 2 * segDz * segDr / (segDz + segDr); + + int numSplitsByMean = + maxMean <= 0 && segMean > maxMean ? 0 : static_cast(std::ceil(segMean / maxMean)) - 1; + int numSplitsByDz = + maxDz <= 0 && segDz > maxDz ? 0 : static_cast(std::ceil(segDz / maxDz)) - 1; + + // Prevent dz from falling below minDz + int numSplitsByMinDz = minDz <= 0 && segDz > minDz ? 0 : static_cast(segDz / minDz) - 1; + + int numSplits = std::min(std::max(numSplitsByMean, numSplitsByDz), numSplitsByMinDz); + + for(int j = 1; j < numSplits; ++j) + { + double t = static_cast(j) / numSplits; + Point2DType newPt(segStart.array() + t * delta); + sorCurveOut.push_back(newPt); + } + sorCurveOut.push_back(segEnd); + } + + return sorCurveOut; +} + +bool MonotonicZSORClipper::getGeometryAsOcts(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& octs) +{ + AXOM_ANNOTATE_SCOPE("MonotonicZSORClipper::getGeometryAsOcts"); + switch(shapeMesh.getRuntimePolicy()) + { + case axom::runtime_policy::Policy::seq: + getGeometryAsOctsImpl(shapeMesh, octs); + break; +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) + case axom::runtime_policy::Policy::omp: + getGeometryAsOctsImpl(shapeMesh, octs); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) + case axom::runtime_policy::Policy::cuda: + getGeometryAsOctsImpl>(shapeMesh, octs); + break; +#endif +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) + case axom::runtime_policy::Policy::hip: + getGeometryAsOctsImpl>(shapeMesh, octs); + break; +#endif + default: + SLIC_ERROR("Axom Internal error: Unhandled execution policy."); + } + return true; +} + +/* + Compute octahedral geometry representation, with an execution policy. + + Side effect: m_sorCurve data is reallocated to the shapeMesh allocator, + if it's not there yet. +*/ +template +bool MonotonicZSORClipper::getGeometryAsOctsImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& octs) +{ + const int allocId = shapeMesh.getAllocatorID(); + octs = axom::Array(0, 0, allocId); + + const auto cellCount = shapeMesh.getCellCount(); + + // Compute an average characteristic length for the mesh cells. + using ReducePolicy = typename axom::execution_space::reduce_policy; + using LoopPolicy = typename execution_space::loop_policy; + axom::ArrayView cellVolumes = shapeMesh.getCellVolumes(); + RAJA::ReduceSum sumVolume(0.0); + RAJA::forall( + RAJA::RangeSegment(0, cellCount), + AXOM_LAMBDA(axom::IndexType cellId) { sumVolume += cellVolumes[cellId]; }); + double avgVolume = sumVolume.get() / cellCount; + double avgCharLength = pow(avgVolume, 1. / 3); + + axom::Array sorCurve = subdivideCurve(m_sorCurve, + 3 * avgCharLength /* maxMean */, + 3 * avgCharLength /* maxDz */, + 2 * avgCharLength /* minDz */); + + // Generate the Octahedra + int octCount = 0; + const bool good = axom::quest::discretize(sorCurve.view(), + int(sorCurve.size()), + m_levelOfRefinement, + octs, + octCount); + + AXOM_UNUSED_VAR(good); + SLIC_ASSERT(good); + SLIC_ASSERT(octCount == octs.size()); + + auto transformer = m_transformer; + auto octsView = octs.view(); + axom::for_all( + octCount, + AXOM_LAMBDA(axom::IndexType iOct) { + OctahedronType& oct = octsView[iOct]; + for(int iVert = 0; iVert < OctahedronType::NUM_VERTS; ++iVert) + { + transformer.transform(oct[iVert].array()); + } + }); + + SLIC_INFO(axom::fmt::format( + "MonotonicZSORClipper '{}' {}-level refinement got {} geometry octs from {} curve points.", + name(), + m_levelOfRefinement, + octs.size(), + sorCurve.size())); + + return true; +} + +/* + Combine consecutive radial segments in SOR curve. Change in place. +*/ +void MonotonicZSORClipper::combineRadialSegments(axom::Array& sorCurve) +{ + int ptCount = sorCurve.size(); + if(ptCount < 3) + { + return; + } + + constexpr double eps = 1e-14; + + // Set sorCurve[j] to sorCurve[i] where j <= i, skipping points + // joining consecutive radial segments. + + int j = 1; + bool prevIsRadial = axom::utilities::isNearlyEqual(sorCurve[j][0] - sorCurve[j - 1][0], eps); + bool curIsRadial = false; + for(int i = 2; i < ptCount; ++i) + { + curIsRadial = axom::utilities::isNearlyEqual(sorCurve[i][0] - sorCurve[i - 1][0], eps); + /* + Current and previous segments share point j. If both are + consecutive radial segments, discard point j by overwriting it + with point i. Else, copy point i to a new point j. + */ + if(!(curIsRadial && prevIsRadial)) + { + ++j; + } + sorCurve[j] = sorCurve[i]; + prevIsRadial = curIsRadial; + } + sorCurve.resize(j + 1); +} + +/* + Find points along the r(z) curve where the z-coordinate changes direction. + + Cases 1 and 2 below show direction changes at point o. Case 3 + shows a potential change at the radial segment, but not a real + change. (Radial segments have constant z and align with the radial + direction.) To decide between cases 2 and 3, defer until the + segment after the radial segment. (The next segment is not radial + because adjacent radials have been combined by combineRadialSegments.) + For case 2, prefer to split at the point closer to the axis of + rotation. + + r ^ + (or y) | (1) (2) (3) + | Single Radial Radial + | point segment segment w/o + | change change change + | + | \ \ \ + | \ \ \ + | o | | + | / o \ + | / / \ + +-------------------------------------> z (or x) +*/ +axom::Array MonotonicZSORClipper::findZSwitchbacks( + axom::ArrayView pts) +{ + const axom::IndexType segCount = pts.size() - 1; + SLIC_ASSERT(segCount > 0); + + // boundaryIdx is where curve's axial direction changes, plus end points. + axom::Array boundaryIdx(0, 2); + boundaryIdx.push_back(0); + + constexpr double eps = 1e-14; + + if(segCount > 1) + { + // Direction is whether z increases or decreases along the curve. + // curDir is the current direction, ignoring radial segments, + // which don't change z. + int curDir = axom::utilities::sign_of(pts[1][0] - pts[0][0], eps); + if(curDir == 0) + { + curDir = axom::utilities::sign_of(pts[2][0] - pts[1][0], eps); + } + + // Detect where z changes direction, and note those indices. + for(axom::IndexType i = 1; i < segCount; ++i) + { + int segDir = axom::utilities::sign_of(pts[i + 1][0] - pts[i][0], eps); + if(segDir == 0) + { + // Radial segment may or may not indicate change. Decide with next segment. + continue; + } + if(segDir != curDir) + { + // Direction change + int prevSegDir = axom::utilities::sign_of(pts[i][0] - pts[i - 1][0], eps); + if(prevSegDir != 0) + { + // Case 1, a clear turn not involving a radial segment. + boundaryIdx.push_back(i); + } + else + { + // Case 2, involving a radial segment. + // Use the radially-closer point of the segment. + int splitI = pts[i][1] < pts[i - 1][1] ? i : i - 1; + boundaryIdx.push_back(splitI); + } + curDir = segDir; + SLIC_ASSERT(curDir != 0); // curDir ignores radial segments. + } + } + } + boundaryIdx.push_back(pts.size() - 1); + return boundaryIdx; +} + +void MonotonicZSORClipper::extractClipperInfo() +{ + auto sorOriginArray = m_info.fetch_existing("sorOrigin").as_double_array(); + auto sorDirectionArray = m_info.fetch_existing("sorDirection").as_double_array(); + for(int d = 0; d < 3; ++d) + { + m_sorOrigin[d] = sorOriginArray[d]; + m_sorDirection[d] = sorDirectionArray[d]; + } + + auto discreteFunctionArray = m_info.fetch_existing("discreteFunction").as_double_array(); + auto n = discreteFunctionArray.number_of_elements(); + + SLIC_ERROR_IF( + n % 2 != 0, + axom::fmt::format( + "***MonotonicZSORClipper: Discrete function must have an even number of values. It has {}.", + n)); + + m_sorCurve.resize(axom::ArrayOptions::Uninitialized(), n / 2); + for(int i = 0; i < n / 2; ++i) + { + m_sorCurve[i] = Point2DType {discreteFunctionArray[i * 2], discreteFunctionArray[i * 2 + 1]}; + } + + m_levelOfRefinement = m_info.fetch_existing("levelOfRefinement").to_double(); +} + +} // namespace experimental +} // end namespace quest +} // end namespace axom diff --git a/src/axom/quest/detail/clipping/MonotonicZSORClipper.hpp b/src/axom/quest/detail/clipping/MonotonicZSORClipper.hpp new file mode 100644 index 0000000000..6d6108ffb4 --- /dev/null +++ b/src/axom/quest/detail/clipping/MonotonicZSORClipper.hpp @@ -0,0 +1,222 @@ +// 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_MONOTONICZSORCLIPPER_HPP +#define AXOM_QUEST_MONOTONICZSORCLIPPER_HPP + +#include "axom/klee/Geometry.hpp" +#include "axom/quest/MeshClipperStrategy.hpp" +#include "axom/primal/geometry/CoordinateTransformer.hpp" + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +/*! + * @brief Geometry clipping operations for simple 3D + * surface-of-revolution geometries. + * + * This implementation requires the SOR curve r(z) to be monotonic in z. + * The values of z that make up the polyline specifying the curve must + * be either non-decreasing or non-increasing. A curve may have segments + * with vertical slope, but it may not double back in z. + * For SOR curves that double back, use \a SORClipper. + * + * The SOR specification may include axis orientation and location + * in addition to any external transformation. + */ +class MonotonicZSORClipper : 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 + * + * \c kGeom.asHierarchy() must contain the following data: + * - sorOrigin: 3D coordinates of the point (z,r) = (0,0) + * - sorDirection: a vector in the direction of the SOR axis. + * - discreteFunction: The discretized r(z) curve as an array of (z,r) pairs. + * - levelOfRefinement: number of refinement levels used + * to approximate the sphere with octahedra. The number + * of octs grows exponentially with level + * (@see quest::discretize(const axom::ArrayView &polyline, int pointcount, int levels, axom::Array &out, int &octcount)). + * In practice, keep this number to 11 or less. + */ + MonotonicZSORClipper(const klee::Geometry& kGeom, const std::string& name = ""); + + /*! + * @brief Construct with parameters to override the specifications + * in the klee::Geometry. + */ + MonotonicZSORClipper(const klee::Geometry& kGeom, + const std::string& name, + axom::ArrayView discreteFunction, + const Point3DType& sorOrigin, + const Vector3DType& sorDirection, + axom::IndexType levelOfRefinement); + + virtual ~MonotonicZSORClipper() = 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 getGeometryAsOcts(quest::experimental::ShapeMesh& shappeMesh, + axom::Array& octs) override; + + axom::ArrayView getSorCurve() const { return m_sorCurve.view(); } + + //@{ + //! @name Utilities shared with SORClipper for handling SOR. + /*! + * @brief Find division points between curve sections where z (x) + * changes directions. + * + * @param sorCurve Set of at least 2 2D points describing a curve + * in r-z space (in host array). + * + * @return Indices of switchbacks, plus the first and last indices. + */ + static axom::Array findZSwitchbacks(axom::ArrayView pts); + + /* + * @brief Combine consecutive radial segments of the curve into a + * single segment. + * + * This step is necessary because some other steps assume there are + * no consecutive radial segments. + */ + static void combineRadialSegments(axom::Array& sorCurve); + //@} + +#if !defined(__CUDACC__) +private: +#endif + std::string m_name; + + /*! + * @brief The discrete r(z) curve as an array of y(x) points. + * + * This data is before internal or external transformations. + * It may include points on each end to connect the curve to + * the axis of rotation. + */ + axom::Array m_sorCurve; + + //! @brief Bounding box of points in m_sorCurve; + BoundingBox2DType m_curveBb; + + //! @brief Maximum radius of the SOR. + double m_maxRadius; + + //! @brief Minimum radius of the SOR. + double m_minRadius; + + //!@brief The point corresponding to z=0 on the SOR axis. + Point3DType m_sorOrigin; + + //!@brief SOR axis in 3D space, in the direction of increasing z. + Vector3DType m_sorDirection; + + //!@brief Level of refinement for discretizing curved + // analytical shapes and surfaces of revolutions. + axom::IndexType m_levelOfRefinement = 0; + + /*! + * @brief Boxes (in rz space) on the curve. + * + * The curve lies completely in these boxes and includes the planes + * of the base and top. Points in these boxes are require more + * computation to determine their signed distance. + */ + axom::Array m_bbOn; + + /*! + * @brief Boxes (in rz space) completely under the curve. + * + * These boxes lie completely under the curve. + */ + axom::Array m_bbUnder; + + //!@brief Internal and external transforms (includes m_sorDirection and m_sorOrigin). + axom::primal::experimental::CoordinateTransformer m_transformer; + + /*! + * @brief Inverse of m_transformer. + * + * Axom supports vector scaling. @see axom::klee::Scale. This means + * a SOR may be transformed into a shape that we cannot represent. + * Therefore, we don't transform the shape until after it's discretized. + * When needed, we will inverse-transform the mesh. + */ + axom::primal::experimental::CoordinateTransformer m_invTransformer; + + template + void labelCellsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView label); + + template + void labelTetsInOutImpl(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView cellIds, + axom::ArrayView tetLabels); + + template + void computeCurveBoxes(quest::experimental::ShapeMesh& shapeMesh, + axom::Array& bbOn, + axom::Array& bbUnder); + + /*! + * @brief Compute 2D bounding box of a polyhedron in the r-z plane. + * @tparam PolyhedronType Either TetrahedronType or HexahedronType. + */ + template + AXOM_HOST_DEVICE BoundingBox2DType estimateBoundingBoxInRz(const PolyhedronType& vertices); + + AXOM_HOST_DEVICE inline MeshClipperStrategy::LabelType rzBbToLabel( + const BoundingBox2DType& bbInRz, + const axom::ArrayView& bbOn, + const axom::ArrayView& bbUnder); + + // Extract clipper info from MeshClipperStrategy::m_info. + void extractClipperInfo(); + + /*! + * @brief Subdivide large segments of the SOR curve to make + * screening more precise. + * + * @param sorCurveIn [in] Un-divided SOR curve + * @param maxMean [in] Subdivide segment if the harmonic mean + * of its bounding box sides exceeds this value. + * @param maxDz [in] Subdivide segment if its length along the + * axis of symmetry exceeds this value. + * @param minDz [in] Don't subdivide segments below this dz. + */ + axom::Array subdivideCurve(const Array& sorCurveIn, + double maxMean, + double maxDz, + double minDz); + + //!@brief Compute geometry as octs, by policy. + template + bool getGeometryAsOctsImpl(quest::experimental::ShapeMesh& shappeMesh, + axom::Array& octs); +}; + +} // namespace experimental +} // namespace quest +} // namespace axom + +#endif // AXOM_QUEST_FSORCLIPPER_HPP diff --git a/src/axom/quest/detail/clipping/SORClipper.cpp b/src/axom/quest/detail/clipping/SORClipper.cpp new file mode 100644 index 0000000000..236e3aadad --- /dev/null +++ b/src/axom/quest/detail/clipping/SORClipper.cpp @@ -0,0 +1,205 @@ +// 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/Discretize.hpp" +#include "axom/quest/detail/clipping/SORClipper.hpp" +#include "axom/quest/MeshClipper.hpp" +#include "axom/spin/BVH.hpp" +#include "axom/fmt.hpp" + +#include + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +SORClipper::SORClipper(const klee::Geometry& kGeom, const std::string& name) + : MeshClipperStrategy(kGeom) + , m_name(name.empty() ? std::string("Sor") : name) + , m_maxRadius(0.0) + , m_minRadius(std::numeric_limits::max()) +{ + extractClipperInfo(); + + for(auto& pt : m_sorCurve) + { + m_maxRadius = fmax(m_maxRadius, pt[1]); + m_minRadius = fmin(m_minRadius, pt[1]); + } + SLIC_ERROR_IF(m_minRadius < 0.0, + axom::fmt::format("SORClipper '{}' has a negative radius", m_name)); + + for(const auto& pt : m_sorCurve) + { + m_curveBb.addPoint(pt); + } + + MonotonicZSORClipper::combineRadialSegments(m_sorCurve); + + axom::Array> sections; + splitIntoMonotonicSections(m_sorCurve.view(), sections); + for(int i = 0; i < sections.size(); ++i) + { + axom::ArrayView section = sections[i].view(); + std::string sectionName = axom::fmt::format("{}.section{:02d}", m_name, i); + m_fsorImpls.push_back(std::make_shared(kGeom, + sectionName, + section, + m_sorOrigin, + m_sorDirection, + m_levelOfRefinement)); + } +} + +bool SORClipper::specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics) +{ + /* + * The SOR curve has been split into SOR functions that do not double + * back on itself. We compute the overlaps for each section and + * accumulate them with the correct sign. (Functions going backward + * remove stuff from the functions above them, so they contribute a + * negative volume.) + * + * By convention, backward curves should generate negative volume, + * but the cone discretization functionality always generates + * positive volumes. We correct this by manually applying the + * correct sign. + */ + const axom::IndexType cellCount = ovlap.size(); + axom::Array tmpOvlap(cellCount, cellCount, ovlap.getAllocatorID()); + for(auto& fsorImpl : m_fsorImpls) + { + tmpOvlap.fill(0.0); + MeshClipper clipper(shapeMesh, fsorImpl); + clipper.setScreenLevel(m_screenLevel); + clipper.setVerbose(false); + clipper.clip(tmpOvlap); + auto sorCurve = fsorImpl->getSorCurve(); + const auto firstZ = sorCurve[0][0]; + const auto lastZ = sorCurve[sorCurve.size() - 1][0]; + int sign = axom::utilities::sign_of(lastZ - firstZ, 0.0); + accumulateData(ovlap, tmpOvlap.view(), double(sign), shapeMesh.getRuntimePolicy()); + MeshClipper::accumulateClippingStats(statistics, clipper.getClippingStats()); + } + return true; +} + +/* + * Split curve into sections that goes monotonically up or down the + * axis of symmetry. If x changes directions at a radial segment, + * split at the end with greater y value. A radial segment is one with + * constant x (or z), pointing in the y (or radial) direction. + * + * This method assumes there are no consecutive radial segments + * (combineRadialSegments has been called on pts). +*/ +void SORClipper::splitIntoMonotonicSections(axom::ArrayView pts, + axom::Array>& sections) +{ + AXOM_ANNOTATE_SCOPE("SORClipper::splitIntoMonotonicSections"); + axom::Array splitIdx = MonotonicZSORClipper::findZSwitchbacks(pts); + + const axom::IndexType sectionCount = splitIdx.size() - 1; + sections.clear(); + sections.resize(sectionCount); + for(axom::IndexType i = 0; i < sectionCount; ++i) + { + axom::IndexType firstInSection = splitIdx[i]; + axom::IndexType lastInSection = splitIdx[i + 1]; + auto& curSection = sections[i]; + curSection.reserve(lastInSection - firstInSection + 1); + for(axom::IndexType j = firstInSection; j <= lastInSection; ++j) + { + curSection.push_back(pts[j]); + } + } +} + +// Compute a += b. +template +void accumulateDataImpl(axom::ArrayView a, axom::ArrayView b, double scale) +{ + SLIC_ASSERT(a.size() == b.size()); + axom::for_all(a.size(), AXOM_LAMBDA(axom::IndexType i) { a[i] += scale * b[i]; }); +} + +void SORClipper::accumulateData(axom::ArrayView a, + axom::ArrayView b, + double scale, + axom::runtime_policy::Policy runtimePolicy) +{ + using RuntimePolicy = axom::runtime_policy::Policy; + if(runtimePolicy == RuntimePolicy::seq) + { + accumulateDataImpl(a, b, scale); + } +#ifdef AXOM_RUNTIME_POLICY_USE_OPENMP + else if(runtimePolicy == RuntimePolicy::omp) + { + accumulateDataImpl(a, b, scale); + } +#endif +#ifdef AXOM_RUNTIME_POLICY_USE_CUDA + else if(runtimePolicy == RuntimePolicy::cuda) + { + accumulateDataImpl>(a, b, scale); + } +#endif +#ifdef AXOM_RUNTIME_POLICY_USE_HIP + else if(runtimePolicy == RuntimePolicy::hip) + { + accumulateDataImpl>(a, b, scale); + } +#endif + else + { + SLIC_ERROR(axom::fmt::format("Unrecognized runtime policy {}", runtimePolicy)); + } + return; +} + +void SORClipper::extractClipperInfo() +{ + auto sorOriginArray = m_info.fetch_existing("sorOrigin").as_double_array(); + auto sorDirectionArray = m_info.fetch_existing("sorDirection").as_double_array(); + for(int d = 0; d < 3; ++d) + { + m_sorOrigin[d] = sorOriginArray[d]; + m_sorDirection[d] = sorDirectionArray[d]; + } + + auto discreteFunctionArray = m_info.fetch_existing("discreteFunction").as_double_array(); + auto n = discreteFunctionArray.number_of_elements(); + + SLIC_ERROR_IF( + n % 2 != 0, + axom::fmt::format( + "***SORClipper: Discrete function must have an even number of values. It has {}.", + n)); + + m_sorCurve.resize(axom::ArrayOptions::Uninitialized(), n / 2); + for(int i = 0; i < n / 2; ++i) + { + m_sorCurve[i] = Point2DType {discreteFunctionArray[i * 2], discreteFunctionArray[i * 2 + 1]}; + } + + m_levelOfRefinement = m_info.fetch_existing("levelOfRefinement").to_double(); + + if(m_info.has_child("screenLevel")) + { + m_screenLevel = m_info["screenLevel"].to_int(); + } +} + +} // namespace experimental +} // end namespace quest +} // end namespace axom diff --git a/src/axom/quest/detail/clipping/SORClipper.hpp b/src/axom/quest/detail/clipping/SORClipper.hpp new file mode 100644 index 0000000000..acab357afb --- /dev/null +++ b/src/axom/quest/detail/clipping/SORClipper.hpp @@ -0,0 +1,116 @@ +// 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_SORCLIPPER_HPP +#define AXOM_QUEST_SORCLIPPER_HPP + +#include "axom/klee/Geometry.hpp" +#include "axom/quest/MeshClipperStrategy.hpp" +#include "axom/quest/detail/clipping/MonotonicZSORClipper.hpp" +#include "axom/primal/geometry/CoordinateTransformer.hpp" + +namespace axom +{ +namespace quest +{ +namespace experimental +{ + +/*! + * @brief Geometry clipping operations for general 3D + * surface-of-revolution geometries. + * + * This implementation allows the SOR function to have non-monotonic + * z in the r(z) curve. For SOR curves that don't, the less complex + * MonotonicZSORClipper is sufficient. + * + * The SOR specification may include rotation and translation + * internally, in addition to any external transformation. +*/ +class SORClipper : 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 + * + * \c kGeom.asHierarchy() must contain the following data: + * - All the parameters used by MonotonicZSORClipper. + * The only difference is that the @a discreteFunction + * points may double back in z. + * - screenLevel: to override the default screen level of + * internal MeshClipper objects. @c MeshClipper::setScreenLevel(). + */ + SORClipper(const klee::Geometry& kGeom, const std::string& name = ""); + + virtual ~SORClipper() = default; + + const std::string& name() const override { return m_name; } + + bool specializedClipCells(quest::experimental::ShapeMesh& shapeMesh, + axom::ArrayView ovlap, + conduit::Node& statistics) override; + +#if !defined(__CUDACC__) +private: +#endif + std::string m_name; + + axom::Array> m_fsorImpls; + + /*! + * @brief The discrete r(z) curve as an array of (z,r) coordinates. + * + * This data is before internal or external transformations. + * It may include points on each end to connect the curve to + * the axis of rotation. + */ + axom::Array m_sorCurve; + + //! @brief Bounding box of points in m_sorCurve; + BoundingBox2DType m_curveBb; + + //! @brief Maximum radius of the SOR. + double m_maxRadius; + + //! @brief Minimum radius of the SOR. + double m_minRadius; + + //!@brief The point corresponding to z=0 on the SOR axis. + Point3DType m_sorOrigin; + + //!@brief SOR axis in 3D space, in the direction of increasing z. + Vector3DType m_sorDirection; + + //!@brief Level of refinement for discretizing curved + // analytical shapes and surfaces of revolutions. + axom::IndexType m_levelOfRefinement = 0; + + //!@brief Screen level setting to use in internal MeshClipper objects. + int m_screenLevel = 3; + + //!@brief Array implementation of a += scale*b. + void accumulateData(axom::ArrayView a, + axom::ArrayView b, + double scale, + axom::runtime_policy::Policy runtimePolicy); + + // Extract clipper info from MeshClipperStrategy::m_info. + void extractClipperInfo(); + + void splitIntoMonotonicSections(axom::ArrayView pts, + axom::Array>& sections); + + void initializeMonotonicZSORClippers(); +}; + +} // namespace experimental +} // namespace quest +} // namespace axom + +#endif // AXOM_QUEST_SORCLIPPER_HPP diff --git a/src/axom/quest/detail/clipping/SphereClipper.cpp b/src/axom/quest/detail/clipping/SphereClipper.cpp index 0640753da0..12edc4b933 100644 --- a/src/axom/quest/detail/clipping/SphereClipper.cpp +++ b/src/axom/quest/detail/clipping/SphereClipper.cpp @@ -187,10 +187,10 @@ AXOM_HOST_DEVICE inline MeshClipperStrategy::LabelType SphereClipper::polyhedron Otherwise, polyhedron is labeled either LABEL_ON or LABEL_IN. Sphere is convex, so polyhedron is IN only if all vertices are inside. - Some polyhedra may be LABEL_ON even though they are actually LABEL_OUT, - but this is a conservative error. The clip function will compute the - correct overlap volume. The purpose of labeling is bypass the - clip function where we can do it efficiently. + Some polyhedra may be labeled ON even though they are actually outside. + This is a conservative error that is fixed ty the clip function later. + The purpose of labeling is bypass the clip function where we can do + it efficiently. */ BoundingBox3DType bb(verts[0]); auto vertCount = Polyhedron::numVertices(); diff --git a/src/axom/quest/detail/clipping/TetClipper.cpp b/src/axom/quest/detail/clipping/TetClipper.cpp index 13e993db54..3ec97b3b46 100644 --- a/src/axom/quest/detail/clipping/TetClipper.cpp +++ b/src/axom/quest/detail/clipping/TetClipper.cpp @@ -44,7 +44,7 @@ TetClipper::TetClipper(const klee::Geometry& kGeom, const std::string& name) bool TetClipper::labelCellsInOut(quest::experimental::ShapeMesh& shapeMesh, axom::Array& cellLabels) { - SLIC_ERROR_IF(shapeMesh.dimension() != 3, "FSorClipper requires a 3D mesh."); + SLIC_ERROR_IF(shapeMesh.dimension() != 3, "TetClipper requires a 3D mesh."); const int allocId = shapeMesh.getAllocatorID(); const auto cellCount = shapeMesh.getCellCount(); diff --git a/src/axom/quest/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index 9d6f9779ba..4539875b78 100644 --- a/src/axom/quest/tests/CMakeLists.txt +++ b/src/axom/quest/tests/CMakeLists.txt @@ -258,7 +258,7 @@ if(CONDUIT_FOUND AND RAJA_FOUND AND AXOM_ENABLE_SIDRE) blt_list_append(TO _policies ELEMENTS "hip" IF AXOM_ENABLE_HIP) endif() - set(_testgeoms "plane" "tet" "hex" "sphere" "plane,hex,tet") + set(_testgeoms "plane" "tet" "hex" "sphere" "sor" "cyl" "cone" "plane,hex,tet") set(_meshTypes "bpSidre" "bpConduit") foreach(_meshType ${_meshTypes}) foreach(_policy ${_policies}) diff --git a/src/axom/quest/tests/quest_discretize.cpp b/src/axom/quest/tests/quest_discretize.cpp index f2d9ff1bfb..b2ae04937b 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)); @@ -219,10 +219,6 @@ void run_degen_segment_tests() polyline[0] = {1., 1.}; polyline[1] = {1.5, -.1}; degenerate_segment_test("b.y < 0", polyline, 2, false); - - polyline[0] = {.5, 1.}; - polyline[1] = {0., 1.}; - degenerate_segment_test("a.x > b.x", polyline, 2, false); } template @@ -250,7 +246,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 +321,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 4b68bbc8ab..b8f0ec0c41 100644 --- a/src/axom/quest/tests/quest_mesh_clipper.cpp +++ b/src/axom/quest/tests/quest_mesh_clipper.cpp @@ -29,6 +29,8 @@ #include "axom/quest.hpp" #include "axom/quest/detail/clipping/HexClipper.hpp" #include "axom/quest/detail/clipping/Plane3DClipper.hpp" +#include "axom/quest/detail/clipping/MonotonicZSORClipper.hpp" +#include "axom/quest/detail/clipping/SORClipper.hpp" #include "axom/quest/detail/clipping/SphereClipper.hpp" #include "axom/quest/detail/clipping/TetClipper.hpp" @@ -106,7 +108,8 @@ struct Input // The shape to run. std::vector testGeom; // The shapes this example is set up to run. - const std::set availableShapes {"tet", "hex", "sphere", "plane"}; // More geometries to come. + const std::set + availableShapes {"sphere", "cyl", "cone", "sor", "tet", "hex", "plane"}; RuntimePolicy policy {RuntimePolicy::seq}; int refinementLevel {7}; @@ -460,6 +463,149 @@ axom::klee::Geometry createGeom_Sphere(const std::string& geomName) return sphereGeometry; } +/* + * Utility function to make a SOR geometry from the specifications in the arguments. + */ +axom::klee::Geometry makeSorGeometry(axom::primal::Point& sorBase, + axom::primal::Vector& sorDirection, + axom::ArrayView discreteFunction, + std::shared_ptr& compositeOp) +{ + axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, + axom::klee::LengthUnit::unspecified}; + + const axom::IndexType levelOfRefinement = params.refinementLevel; + axom::klee::Geometry sorGeometry(prop, + discreteFunction, + sorBase, + sorDirection, + levelOfRefinement, + compositeOp); + return sorGeometry; +} + +double computeVolume_Sor(axom::ArrayView 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::Geometry createGeom_Sor(const std::string& geomName) +{ + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; + axom::primal::Vector sorDirection = params.direction.empty() + ? primal::Vector3D {1.0, 0.0, 0.0} + : primal::Vector3D {params.direction.data()}; + // discreteFunction is discrete (z,r) pairs describing the r(z) function + // to be rotated around the z axis. + using Point2DType = axom::primal::Point; + double zLen = 0.5 * (params.length < 0 ? 2.40 : params.length); + double maxR = params.radius < 0 ? 1.10 : params.radius; + axom::Array discretePts(0, 10); + discretePts.push_back(Point2DType({-1.0 * zLen, 1.0 * maxR})); + discretePts.push_back(Point2DType({0.4 * zLen, 1.0 * maxR})); + discretePts.push_back(Point2DType({0.4 * zLen, 0.7 * maxR})); + discretePts.push_back(Point2DType({1.0 * zLen, 0.7 * maxR})); + discretePts.push_back(Point2DType({1.0 * zLen, 0.4 * maxR})); + discretePts.push_back(Point2DType({0.5 * zLen, 0.4 * maxR})); + discretePts.push_back(Point2DType({0.5 * zLen, 0.3 * maxR})); + discretePts.push_back(Point2DType({0.0 * zLen, 0.3 * maxR})); + discretePts.push_back(Point2DType({0.0 * zLen, 0.5 * maxR})); + discretePts.push_back(Point2DType({0.2 * zLen, 0.5 * maxR})); + discretePts.push_back(Point2DType({0.2 * zLen, 0.7 * maxR})); + discretePts.push_back(Point2DType({-1.0 * zLen, 0.7 * maxR})); + axom::ArrayView discreteFunction((const double*)discretePts.data(), + discretePts.size(), + 2); + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp); + + axom::klee::Geometry sorGeometry = + makeSorGeometry(sorBase, sorDirection, discreteFunction, compositeOp); + sorGeometry.asHierarchy()["screenLevel"] = params.screenLevel; + + exactGeomVols[geomName] = vScale * computeVolume_Sor(discreteFunction); + // Tolerance should account for discretization errors. + errorToleranceRel[geomName] = params.refinementLevel <= 5 ? 0.007 : 0.0063; + errorToleranceAbs[geomName] = errorToleranceRel[geomName] * exactGeomVols[geomName]; + + return sorGeometry; +} + +axom::klee::Geometry createGeom_Cylinder(const std::string& geomName) +{ + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; + axom::primal::Vector sorDirection = params.direction.empty() + ? primal::Vector3D {1.0, 0.0, 0.0} + : primal::Vector3D {params.direction.data()}; + // 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.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; + discreteFunction(1, 1) = radius; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp); + + axom::klee::Geometry sorGeometry = + makeSorGeometry(sorBase, sorDirection, discreteFunction, compositeOp); + + exactGeomVols[geomName] = vScale * computeVolume_Sor(discreteFunction); + // Tolerance should account for discretization errors. + errorToleranceRel[geomName] = params.refinementLevel <= 5 ? 0.00075 : 0.00005; + errorToleranceAbs[geomName] = errorToleranceRel[geomName] * exactGeomVols[geomName]; + + return sorGeometry; +} + +axom::klee::Geometry createGeom_Cone(const std::string& geomName) +{ + Point3D sorBase = params.center.empty() ? Point3D {0.0, 0.0, 0.0} : Point3D {params.center.data()}; + axom::primal::Vector sorDirection = params.direction.empty() + ? primal::Vector3D {1.0, 0.0, 0.0} + : primal::Vector3D {params.direction.data()}; + // 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 ? 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; + discreteFunction(1, 1) = topRadius; + + auto compositeOp = std::make_shared(startProp); + addScaleOperator(*compositeOp); + addTranslateOperator(*compositeOp); + + axom::klee::Geometry sorGeometry = + makeSorGeometry(sorBase, sorDirection, discreteFunction, compositeOp); + + exactGeomVols[geomName] = vScale * computeVolume_Sor(discreteFunction); + // Tolerance should account for discretization errors. + errorToleranceRel[geomName] = params.refinementLevel <= 5 ? 0.00075 : 0.00005; + errorToleranceAbs[geomName] = errorToleranceRel[geomName] * exactGeomVols[geomName]; + + return sorGeometry; +} + axom::klee::Geometry createGeom_Tet(const std::string& geomName) { axom::klee::TransformableGeometryProperties prop {axom::klee::Dimensions::Three, @@ -946,6 +1092,23 @@ int main(int argc, char** argv) geomStrategies.push_back( std::make_shared(createGeom_Tet(name), name)); } + else if(tg == "cyl") + { + geomStrategies.push_back( + std::make_shared(createGeom_Cylinder(name), + name)); + } + else if(tg == "cone") + { + geomStrategies.push_back( + std::make_shared(createGeom_Cone(name), + name)); + } + else if(tg == "sor") + { + geomStrategies.push_back( + std::make_shared(createGeom_Sor(name), name)); + } // More geometries to come. }