Skip to content

Commit

Permalink
Merge pull request #1178 from LLNL/feature/gunney/marching-cubes-ghosts
Browse files Browse the repository at this point in the history
Ghost layer support for `MarchingCubes`  I've tested this branch on Cuda and HIP.
  • Loading branch information
gunney1 committed Nov 17, 2023
2 parents e4acf4a + 12efcf5 commit 35ae725
Show file tree
Hide file tree
Showing 15 changed files with 3,020 additions and 821 deletions.
17 changes: 17 additions & 0 deletions src/axom/core/Array.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -344,6 +344,23 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE>>

/// @}

/*!
@brief Convert 1D Array into a StackArray.
*/
template <int LENGTH1D, typename TT = T, int TDIM = DIM>
AXOM_HOST_DEVICE inline
typename std::enable_if<TDIM == 1, axom::StackArray<TT, LENGTH1D>>::type
to_stack_array() const
{
axom::StackArray<TT, LENGTH1D> rval;
IndexType copyCount = LENGTH1D <= m_num_elements ? LENGTH1D : m_num_elements;
for(IndexType i = 0; i < copyCount; ++i)
{
rval[i] = m_data[i];
}
return rval;
}

/// @}

/// \name Array methods to modify the data.
Expand Down
116 changes: 116 additions & 0 deletions src/axom/quest/ArrayIndexer.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
// Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and
// other Axom Project Developers. See the top-level LICENSE file for details.
//
// SPDX-License-Identifier: (BSD-3-Clause)

#ifndef QUEST_ARRAYINDEXER_HPP_
#define QUEST_ARRAYINDEXER_HPP_

#include "axom/slic.hpp"
#include "axom/core/StackArray.hpp"
#include "axom/core/numerics/matvecops.hpp"

namespace axom
{
/*!
@brief Indexing into a multidimensional structured array.
Supports row-major and column-major ordering and arbitrary
permutations of the ordering.
*/
template <typename T, int DIM>
class ArrayIndexer
{
axom::StackArray<T, DIM> m_strides;
axom::StackArray<std::uint16_t, DIM> m_slowestDirs;

public:
/*!
@brief Constructor for row- or column major indexing.
@param [in] lengths Lengths of the array
@param [in] order: c is column major; r is row major.
*/
ArrayIndexer(const axom::StackArray<T, DIM>& lengths, char order)
{
SLIC_ASSERT(order == 'c' || order == 'r');
if(order == 'r')
{
for(int d = 0; d < DIM; ++d)
{
m_slowestDirs[d] = DIM - 1 - d;
}
m_strides[0] = 1;
for(int d = 1; d < DIM; ++d)
{
m_strides[d] = m_strides[d - 1] * lengths[d - 1];
}
}
else
{
for(int d = 0; d < DIM; ++d)
{
m_slowestDirs[d] = d;
}
m_strides[DIM - 1] = 1;
for(int d = DIM - 2; d >= 0; --d)
{
m_strides[d] = m_strides[d + 1] * lengths[d + 1];
}
}
}

//!@brief Constructor for arbitrary-stride indexing.
ArrayIndexer(const axom::StackArray<T, DIM>& strides) : m_strides(strides)
{
for(int d = 0; d < DIM; ++d)
{
m_slowestDirs[d] = d;
}
for(int s = 0; s < DIM; ++s)
{
for(int d = s; d < DIM; ++d)
{
if(m_strides[m_slowestDirs[s]] < m_strides[m_slowestDirs[d]])
{
std::swap(m_slowestDirs[s], m_slowestDirs[d]);
}
}
}
}

//!@brief Index directions, ordered from slowest to fastest.
inline AXOM_HOST_DEVICE const axom::StackArray<std::uint16_t, DIM>& slowestDirs() const
{
return m_slowestDirs;
}

//!@brief Strides.
inline AXOM_HOST_DEVICE const axom::StackArray<axom::IndexType, DIM>& strides() const
{
return m_strides;
}

//!@brief Convert multidimensional index to flat index.
inline AXOM_HOST_DEVICE T toFlatIndex(const axom::StackArray<T, DIM>& multiIndex) const
{
T rval = numerics::dot_product(multiIndex.begin(), m_strides.begin(), DIM);
return rval;
}

//!@brief Convert flat index to multidimensional index.
inline AXOM_HOST_DEVICE axom::StackArray<T, DIM> toMultiIndex(T flatIndex) const
{
axom::StackArray<T, DIM> multiIndex;
for(int d = 0; d < DIM; ++d)
{
int dir = m_slowestDirs[d];
multiIndex[dir] = flatIndex / m_strides[dir];
flatIndex -= multiIndex[dir] * m_strides[dir];
}
return multiIndex;
}
};

} // end namespace axom

#endif // QUEST_ARRAYINDEXER_HPP_
2 changes: 2 additions & 0 deletions src/axom/quest/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@ set( quest_headers
interface/signed_distance.hpp

util/mesh_helpers.hpp
MeshViewUtil.hpp
ArrayIndexer.hpp
)

set( quest_sources
Expand Down
63 changes: 43 additions & 20 deletions src/axom/quest/MarchingCubes.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,16 @@
//
// SPDX-License-Identifier: (BSD-3-Clause)

#include "axom/core/execution/execution_space.hpp"
#include "axom/quest/MarchingCubes.hpp"
#include "axom/quest/detail/MarchingCubesImpl.hpp"
#include "conduit_blueprint.hpp"
#include "axom/fmt.hpp"
#include "axom/config.hpp"

// Implementation requires Conduit.
#ifdef AXOM_USE_CONDUIT
#include "conduit_blueprint.hpp"

#include "axom/core/execution/execution_space.hpp"
#include "axom/quest/MarchingCubes.hpp"
#include "axom/quest/detail/MarchingCubesImpl.hpp"
#include "axom/fmt.hpp"

namespace axom
{
Expand All @@ -20,7 +25,9 @@ MarchingCubes::MarchingCubes(RuntimePolicy runtimePolicy,
: m_runtimePolicy(runtimePolicy)
, m_singles()
, m_topologyName(topologyName)
, m_fcnFieldName()
, m_fcnPath()
, m_maskFieldName(maskField)
, m_maskPath(maskField.empty() ? std::string() : "fields/" + maskField)
{
const bool isMultidomain = conduit::blueprint::mesh::is_multi_domain(bpMesh);
Expand All @@ -40,6 +47,7 @@ MarchingCubes::MarchingCubes(RuntimePolicy runtimePolicy,

void MarchingCubes::setFunctionField(const std::string& fcnField)
{
m_fcnFieldName = fcnField;
m_fcnPath = "fields/" + fcnField;
for(auto& s : m_singles)
{
Expand All @@ -49,7 +57,7 @@ void MarchingCubes::setFunctionField(const std::string& fcnField)

void MarchingCubes::computeIsocontour(double contourVal)
{
SLIC_ASSERT_MSG(!m_fcnPath.empty(),
SLIC_ASSERT_MSG(!m_fcnFieldName.empty(),
"You must call setFunctionField before computeIsocontour.");

for(int dId = 0; dId < m_singles.size(); ++dId)
Expand Down Expand Up @@ -118,13 +126,15 @@ void MarchingCubes::populateContourMesh(
auto* domainIdPtr =
mesh.getFieldPtr<axom::IndexType>(domainIdField,
axom::mint::CELL_CENTERED);
// TODO: Verify that UnstructuredMesh only supports host memory.

int userDomainId = single->getDomainId(dId);

axom::detail::ArrayOps<axom::IndexType, MemorySpace::Dynamic>::fill(
domainIdPtr,
nPrev,
nNew - nPrev,
execution_space<axom::SEQ_EXEC>::allocatorID(),
dId);
userDomainId);
}
}
SLIC_ASSERT(mesh.getNumberOfNodes() == contourNodeCount);
Expand All @@ -139,7 +149,9 @@ MarchingCubesSingleDomain::MarchingCubesSingleDomain(RuntimePolicy runtimePolicy
, m_dom(nullptr)
, m_ndim(0)
, m_topologyName(topologyName)
, m_fcnFieldName()
, m_fcnPath()
, m_maskFieldName(maskField)
, m_maskPath(maskField.empty() ? std::string() : "fields/" + maskField)
{
SLIC_ASSERT_MSG(
Expand Down Expand Up @@ -170,8 +182,8 @@ void MarchingCubesSingleDomain::setDomain(const conduit::Node& dom)

m_dom = &dom;

m_ndim = conduit::blueprint::mesh::coordset::dims(
dom.fetch_existing("coordsets/coords"));
m_ndim = conduit::blueprint::mesh::topology::dims(
dom.fetch_existing(axom::fmt::format("topologies/{}", m_topologyName)));
SLIC_ASSERT(m_ndim >= 2 && m_ndim <= 3);

const conduit::Node& coordsValues =
Expand All @@ -184,22 +196,31 @@ void MarchingCubesSingleDomain::setDomain(const conduit::Node& dom)

void MarchingCubesSingleDomain::setFunctionField(const std::string& fcnField)
{
m_fcnFieldName = fcnField;
m_fcnPath = "fields/" + fcnField;
SLIC_ASSERT(m_dom->has_path(m_fcnPath));
SLIC_ASSERT(m_dom->fetch_existing(m_fcnPath + "/association").as_string() ==
"vertex");
SLIC_ASSERT(m_dom->has_path(m_fcnPath + "/values"));
}

int MarchingCubesSingleDomain::getDomainId(int defaultId) const
{
int rval = defaultId;
if(m_dom->has_path("state/domain_id"))
{
rval = m_dom->fetch_existing("state/domain_id").as_int();
}
return rval;
}

void MarchingCubesSingleDomain::computeIsocontour(double contourVal)
{
SLIC_ASSERT_MSG(!m_fcnPath.empty(),
SLIC_ASSERT_MSG(!m_fcnFieldName.empty(),
"You must call setFunctionField before computeIsocontour.");

allocateImpl();
const std::string coordsetPath = "coordsets/" +
m_dom->fetch_existing("topologies/" + m_topologyName + "/coordset").as_string();
m_impl->initialize(*m_dom, coordsetPath, m_fcnPath, m_maskPath);
m_impl->initialize(*m_dom, m_topologyName, m_fcnFieldName, m_maskFieldName);
m_impl->setContourValue(contourVal);
m_impl->markCrossings();
m_impl->scanCrossings();
Expand All @@ -217,7 +238,7 @@ void MarchingCubesSingleDomain::allocateImpl()
: std::unique_ptr<ImplBase>(
new MarchingCubesImpl<3, axom::SEQ_EXEC, axom::SEQ_EXEC>);
}
#ifdef _AXOM_MC_USE_OPENMP
#ifdef _AXOM_MARCHINGCUBES_USE_OPENMP
else if(m_runtimePolicy == RuntimePolicy::omp)
{
m_impl = m_ndim == 2
Expand All @@ -226,8 +247,8 @@ void MarchingCubesSingleDomain::allocateImpl()
: std::unique_ptr<ImplBase>(
new MarchingCubesImpl<3, axom::OMP_EXEC, axom::SEQ_EXEC>);
}
#endif
#ifdef _AXOM_MC_USE_CUDA
#endif
#ifdef _AXOM_MARCHINGCUBES_USE_CUDA
else if(m_runtimePolicy == RuntimePolicy::cuda)
{
m_impl = m_ndim == 2
Expand All @@ -236,8 +257,8 @@ void MarchingCubesSingleDomain::allocateImpl()
: std::unique_ptr<ImplBase>(
new MarchingCubesImpl<3, axom::CUDA_EXEC<256>, axom::CUDA_EXEC<1>>);
}
#endif
#ifdef _AXOM_MC_USE_HIP
#endif
#ifdef _AXOM_MARCHINGCUBES_USE_HIP
else if(m_runtimePolicy == RuntimePolicy::hip)
{
m_impl = m_ndim == 2
Expand All @@ -246,7 +267,7 @@ void MarchingCubesSingleDomain::allocateImpl()
: std::unique_ptr<ImplBase>(
new MarchingCubesImpl<3, axom::HIP_EXEC<256>, axom::HIP_EXEC<1>>);
}
#endif
#endif
else
{
SLIC_ERROR(axom::fmt::format(
Expand All @@ -255,5 +276,7 @@ void MarchingCubesSingleDomain::allocateImpl()
}
}

#endif // AXOM_USE_CONDUIT

} // end namespace quest
} // end namespace axom

0 comments on commit 35ae725

Please sign in to comment.