diff --git a/src/axom/quest/CMakeLists.txt b/src/axom/quest/CMakeLists.txt index 9ffccdc970..1cdc5997cc 100644 --- a/src/axom/quest/CMakeLists.txt +++ b/src/axom/quest/CMakeLists.txt @@ -150,13 +150,16 @@ 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/SphereClipper.hpp detail/clipping/MeshClipperImpl.hpp) list(APPEND quest_sources MeshClipper.cpp + MeshClipperStrategy.cpp + detail/clipping/Plane3DClipper.cpp detail/clipping/TetClipper.cpp detail/clipping/SphereClipper.cpp - MeshClipperStrategy.cpp) + detail/clipping/MeshClipperImpl.hpp) endif() endif() diff --git a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp index 3acee0c513..5f1ccf8788 100644 --- a/src/axom/quest/detail/clipping/MeshClipperImpl.hpp +++ b/src/axom/quest/detail/clipping/MeshClipperImpl.hpp @@ -1047,17 +1047,16 @@ class MeshClipperImpl : public MeshClipper::Impl */ AXOM_HOST_DEVICE static 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])}; + // Oct vertex indices of the four tets that the oct decomposes into. + IndexType tIds[4][4] = {{0, 3, 1, 2}, {0, 3, 2, 4}, {0, 3, 4, 5}, {0, 3, 5, 1}}; double octVol = 0.0; for(int i = 0; i < 4; ++i) { - double tetVol = tets[i].signedVolume(); - SLIC_ASSERT(tetVol >= -EPS); // Tet may be degenerate but not inverted. - octVol += axom::utilities::abs(tetVol); + TetrahedronType tet(oct[tIds[i][0]], oct[tIds[i][1]], oct[tIds[i][2]], oct[tIds[i][3]]); + double tetVol = tet.signedVolume(); + octVol -= tetVol; // Octs from the discretized geometries are inverted w.r.t. tIDs. } + SLIC_ASSERT(octVol > 0.0); return octVol; } diff --git a/src/axom/quest/detail/clipping/Plane3DClipper.cpp b/src/axom/quest/detail/clipping/Plane3DClipper.cpp new file mode 100644 index 0000000000..d3b87c26fb --- /dev/null +++ b/src/axom/quest/detail/clipping/Plane3DClipper.cpp @@ -0,0 +1,412 @@ +// 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; + + axom::ReduceSum missSum {0}; + + 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); + if(overlap.numVertices() >= 4) + { + auto volume = overlap.volume(); + vol += volume; + } + else + { + missSum += 1; + } + } + ovlap[cellId] = vol; + }); + + statistics["clipsOn"].set_int64(cellIds.size() * NUM_TETS_PER_HEX); + statistics["clipsSum"].set_int64(cellIds.size() * NUM_TETS_PER_HEX); + statistics["missSum"].set_int64(missSum.get()); +} + +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; + + 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); + }); + + // Because the tet screening is perfect, all tets in tetIds are on the plane. + statistics["onSum"].set_int64(tetCount); + statistics["clipsSum"].set_int64(tetCount); +} + +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..b7406f33c2 --- /dev/null +++ b/src/axom/quest/detail/clipping/Plane3DClipper.hpp @@ -0,0 +1,111 @@ +// 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. + * + * @internal Because this class provides screening via the + * labelCellsInOut and labelTetsInOut methods, the + * specializedClipCells methods below are not essential. They are + * implemented only to avoid crashing when the MeshClipper's screen + * level is overridden (for performance comparisons). The + * specializedClipTets method is much faster and is the one used for + * the default screen level. If the screen level override is remove, + * the non-essential methods can also be removed. + */ + 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/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index ba1329069a..1e85dbd181 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 "tet" "sphere" "tet,sphere") + set(_testgeoms "plane" "tet" "sphere" "plane,tet,sphere") set(_meshTypes "bpSidre" "bpConduit") foreach(_meshType ${_meshTypes}) foreach(_policy ${_policies}) @@ -280,7 +280,6 @@ if(CONDUIT_FOUND AND RAJA_FOUND AND AXOM_ENABLE_SIDRE) --refinements 5 --scale .99 .99 .99 --dir 8 4 2 - --screenLevel 2 --meshType ${_meshType} inline_mesh --min -2 -2 -2 --max 2 2 2 --res 16 16 16 NUM_MPI_TASKS ${_nranks} diff --git a/src/axom/quest/tests/quest_mesh_clipper.cpp b/src/axom/quest/tests/quest_mesh_clipper.cpp index 98a91a9f5f..daa3cc5b79 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/SphereClipper.hpp" #include "axom/quest/detail/clipping/TetClipper.hpp" @@ -102,7 +103,7 @@ struct Input // The shape to run. std::vector testGeom; // The shapes this example is set up to run. - const std::set availableShapes {"tet", "sphere"}; // More geometries to come. + const std::set availableShapes {"tet", "sphere", "plane"}; // More geometries to come. RuntimePolicy policy {RuntimePolicy::seq}; int refinementLevel {7}; @@ -277,8 +278,8 @@ const std::string matsetName = "matset"; const std::string coordsetName = "coords"; int cellCount = -1; // Translation to individual octants (override) when running multiple shapes. -// Exception: the plane always placed at origin to facilitate finding its -// exact overlap volume. +// Exception: the plane always placed at the center of the box mesh +// 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}, @@ -482,6 +483,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. @@ -862,7 +893,12 @@ int main(int argc, char** argv) } std::string name = axom::fmt::format("{}.{}", tg, geomReps[tg]++); - if(tg == "sphere") + if(tg == "plane") + { + geomStrategies.push_back( + std::make_shared(createGeom_Plane(name), name)); + } + else if(tg == "sphere") { geomStrategies.push_back( std::make_shared(createGeom_Sphere(name), name));