Skip to content
179 changes: 123 additions & 56 deletions src/axom/quest/MeshClipper.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,16 @@ MeshClipper::MeshClipper(quest::experimental::ShapeMesh& shapeMesh,
, m_strategy(strategy)
, m_impl(newImpl())
, m_verbose(false)
{ }
, m_screenLevel(3)
{
// Initialize statistics used by this class.
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);
}

void MeshClipper::clip(axom::Array<double>& ovlap)
{
Expand Down Expand Up @@ -58,9 +67,22 @@ void MeshClipper::clip(axom::ArrayView<double> 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<LabelType> 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;

Expand All @@ -79,7 +101,8 @@ void MeshClipper::clip(axom::ArrayView<double> ovlap)

if(m_verbose)
{
logLabelStats(cellLabels, "cells");
getLabelCounts(cellLabels, cellsInCount, cellsOnCount, cellsOutCount);
logClippingStats();
}

AXOM_ANNOTATE_BEGIN("MeshClipper::processInOut");
Expand All @@ -90,16 +113,24 @@ void MeshClipper::clip(axom::ArrayView<double> ovlap)
m_impl->collectOnIndices(cellLabels.view(), cellsOnBdry);

axom::Array<LabelType> 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<axom::IndexType> tetsOnBdry;

if(withTetInOut)
{
if(m_verbose)
{
logLabelStats(tetLabels, "tets");
getLabelCounts(tetLabels, tetsInCount, tetsOnCount, tetsOutCount);
logClippingStats();
}

m_impl->collectOnIndices(tetLabels.view(), tetsOnBdry);
m_impl->remapTetIndices(cellsOnBdry, tetsOnBdry);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gut check: Was this change intentional? Or perhaps a bug fix? I don't see an API change in this PR for Impl::remapTetIndices()

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It fixes an error I made when I split up the huge PR into smaller ones. I didn't catch the error immediately because this code wasn't exercised until the PRs I submitted later. So this is a correction to a code that wasn't being exercised, but I thought it made sense to keep the correction as close as possible to the place where the error was made. In the original huge PR, it was correct.


Expand All @@ -116,27 +147,25 @@ void MeshClipper::clip(axom::ArrayView<double> ovlap)
//
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);
}

if(!done)
{
AXOM_ANNOTATE_SCOPE("MeshClipper::clip3D_limited");
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
{
Expand All @@ -145,15 +174,15 @@ void MeshClipper::clip(axom::ArrayView<double> ovlap)
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;
}
}

Expand Down Expand Up @@ -197,52 +226,90 @@ std::unique_ptr<MeshClipper::Impl> MeshClipper::newImpl()
return impl;
}

void MeshClipper::logLabelStats(axom::ArrayView<const LabelType> labels, const std::string& labelType)
#if defined(AXOM_USE_MPI)
template <typename T>
void globalReduce(axom::Array<T>& 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<axom::IndexType>::type, MPI_SUM, 0, MPI_COMM_WORLD);
axom::Array<T> localValues(values);
MPI_Allreduce(localValues.data(),
values.data(),
values.size(),
axom::mpi_traits<T>::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<axom::IndexType>::type,
MPI_SUM,
0,
MPI_COMM_WORLD);
MPI_Reduce(&localCellInCount,
&maxLocalCellInCount,
1,
axom::mpi_traits<axom::IndexType>::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<std::int64_t> 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<std::int64_t> 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
Expand Down
Loading