Skip to content

Commit

Permalink
use unordered set in fast edge counter, and move to utility file for …
Browse files Browse the repository at this point in the history
…shared use
  • Loading branch information
nyoungbq committed Jan 8, 2025
1 parent e8a3457 commit 33f62d4
Show file tree
Hide file tree
Showing 2 changed files with 148 additions and 121 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,13 @@
#include "simplnx/DataStructure/DataGroup.hpp"
#include "simplnx/DataStructure/Geometry/IGeometry.hpp"
#include "simplnx/DataStructure/Geometry/TriangleGeom.hpp"
#include "simplnx/Utilities/GeometryHelpers.hpp"

#include "EbsdLib/Core/Orientation.hpp"
#include "EbsdLib/Core/OrientationTransformation.hpp"

#include <unordered_set>

using namespace nx::core;

namespace
Expand All @@ -19,134 +22,16 @@ using VertsStore = AbstractDataStore<INodeGeometry0D::SharedVertexList::value_ty
constexpr double k_Multiplier = 1.0 / (4.0 * Constants::k_PiD);
constexpr float64 k_ScaleFactor = 1.0;

constexpr uint64 k_MaxOptimizedValue = static_cast<uint64>(std::numeric_limits<uint32>::max());

template<typename T>
usize SafeEdgeCount(const AbstractDataStore<T>& faceStore)
{
const usize numFaces = faceStore.getNumberOfTuples();
const usize numComp = faceStore.getNumberOfComponents();
T v0 = 0;
T v1 = 0;

std::set<std::pair<T, T>> edgeSet;

for(usize i = 0; i < numFaces; i++)
{
const usize offset = i * numComp;

for(usize j = 0; j < numComp; j++)
{
if(j == (numComp - 1))
{
if(faceStore[offset + j] > faceStore[offset + 0])
{
v0 = faceStore[offset + 0];
v1 = faceStore[offset + j];
}
else
{
v0 = faceStore[offset + j];
v1 = faceStore[offset + 0];
}
}
else
{
if(faceStore[offset + j] > faceStore[offset + j + 1])
{
v0 = faceStore[offset + j + 1];
v1 = faceStore[offset + j];
}
else
{
v0 = faceStore[offset + j];
v1 = faceStore[offset + j + 1];
}
}
std::pair<T, T> edge = std::make_pair(v0, v1);
edgeSet.insert(edge);
}
}

return edgeSet.size();
}

template<typename T>
usize FastEdgeCount(const AbstractDataStore<T>& faceStore)
{
const usize numFaces = faceStore.getNumberOfTuples();
const usize numComp = faceStore.getNumberOfComponents();
uint32 v0 = 0;
uint32 v1 = 0;

std::set<uint64> edgeSet;

for(usize i = 0; i < numFaces; i++)
{
const usize offset = i * numComp;

for(usize j = 0; j < numComp; j++)
{
if(j == (numComp - 1))
{
if(faceStore[offset + j] > faceStore[offset + 0])
{
v0 = static_cast<uint32>(faceStore[offset + 0]);
v1 = static_cast<uint32>(faceStore[offset + j]);
}
else
{
v0 = static_cast<uint32>(faceStore[offset + j]);
v1 = static_cast<uint32>(faceStore[offset + 0]);
}
}
else
{
if(faceStore[offset + j] > faceStore[offset + j + 1])
{
v0 = static_cast<uint32>(faceStore[offset + j + 1]);
v1 = static_cast<uint32>(faceStore[offset + j]);
}
else
{
v0 = static_cast<uint32>(faceStore[offset + j]);
v1 = static_cast<uint32>(faceStore[offset + j + 1]);
}
}

edgeSet.insert(static_cast<uint64>(v0) << 32 | v1);
}
}

return edgeSet.size();
}

template<typename T>
usize FindNumEdges(const AbstractDataStore<T>& faceStore, usize numVertices = (k_MaxOptimizedValue + 1))
{
// This case may seem niche, but it is designed with Indexing types in mind specifically IGeometry::MeshIndexType
if constexpr(!std::is_signed_v<T>)
{
if(numVertices < k_MaxOptimizedValue)
{
// speedier method because max vertices value fits into uint32
return FastEdgeCount(faceStore);
}
}
// Slower method to avoid overflow
return SafeEdgeCount(faceStore);
}

usize FindEulerCharacteristic(usize numVertices, usize numFaces, usize numRegions)
{
return numVertices + numFaces - (2 * numRegions);
}

template<typename T>
template <typename T>
bool ValidateMesh(const AbstractDataStore<T>& faceStore, usize numVertices, usize numRegions)
{
// Expensive call
usize numEdges = FindNumEdges(faceStore, numVertices);
usize numEdges = GeometryHelpers::Connectivity::FindNumEdges(faceStore, numVertices);

return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions);
}
Expand Down Expand Up @@ -274,7 +159,7 @@ Result<> ComputeTriangleGeomShapes::operator()()
// the assumption here is face labels contains information on region ids, that it is contiguous in the values, and that 0 is an invalid id
// (ie the max function means that if the values in array are [1,2,4,5] it will assume there are 5 regions)
usize numRegions = *std::max_element(faceLabels.begin(), faceLabels.end());

if(!ValidateMesh(triangleList, verts.getNumberOfTuples(), numRegions))
{
return MakeErrorResult(-64720, fmt::format("The Euler Characteristic of the shape was found to be unequal to 2, this implies the shape may not be watertight or is malformed."));
Expand Down
142 changes: 142 additions & 0 deletions src/simplnx/Utilities/GeometryHelpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@

#include <Eigen/Dense>

#include <unordered_set>

namespace nx::core
{
namespace GeometryHelpers
Expand All @@ -32,6 +34,146 @@ SIMPLNX_EXPORT std::string GenerateGeometryInfo(const nx::core::SizeVec3& dims,

namespace Connectivity
{
namespace detail
{
constexpr uint64 k_MaxOptimizedValue = static_cast<uint64>(std::numeric_limits<uint32>::max());

/**
* @brief !!! DO NOT CALL DIRECTLY !!! Prefer FindNumEdges()
* @tparam T indexing type of the face store
* @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList
* @return usize This is the number of edges
*/
template <typename T>
usize SafeEdgeCount(const AbstractDataStore<T>& faceStore)
{
const usize numFaces = faceStore.getNumberOfTuples();
const usize numComp = faceStore.getNumberOfComponents();
T v0 = 0;
T v1 = 0;

std::set<std::pair<T, T>> edgeSet;

for(usize i = 0; i < numFaces; i++)
{
const usize offset = i * numComp;

for(usize j = 0; j < numComp; j++)
{
if(j == (numComp - 1))
{
if(faceStore[offset + j] > faceStore[offset + 0])
{
v0 = faceStore[offset + 0];
v1 = faceStore[offset + j];
}
else
{
v0 = faceStore[offset + j];
v1 = faceStore[offset + 0];
}
}
else
{
if(faceStore[offset + j] > faceStore[offset + j + 1])
{
v0 = faceStore[offset + j + 1];
v1 = faceStore[offset + j];
}
else
{
v0 = faceStore[offset + j];
v1 = faceStore[offset + j + 1];
}
}
std::pair<T, T> edge = std::make_pair(v0, v1);
edgeSet.insert(edge);
}
}

return edgeSet.size();
}

/**
* @brief !!! DO NOT CALL DIRECTLY !!! Prefer FindNumEdges()
* @tparam T indexing type of the face store
* @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList
* @return usize This is the number of edges
*/
template <typename T>
usize FastEdgeCount(const AbstractDataStore<T>& faceStore)
{
const usize numFaces = faceStore.getNumberOfTuples();
const usize numComp = faceStore.getNumberOfComponents();
uint32 v0 = 0;
uint32 v1 = 0;

std::unordered_set<uint64> edgeSet;

for(usize i = 0; i < numFaces; i++)
{
const usize offset = i * numComp;

for(usize j = 0; j < numComp; j++)
{
if(j == (numComp - 1))
{
if(faceStore[offset + j] > faceStore[offset + 0])
{
v0 = static_cast<uint32>(faceStore[offset + 0]);
v1 = static_cast<uint32>(faceStore[offset + j]);
}
else
{
v0 = static_cast<uint32>(faceStore[offset + j]);
v1 = static_cast<uint32>(faceStore[offset + 0]);
}
}
else
{
if(faceStore[offset + j] > faceStore[offset + j + 1])
{
v0 = static_cast<uint32>(faceStore[offset + j + 1]);
v1 = static_cast<uint32>(faceStore[offset + j]);
}
else
{
v0 = static_cast<uint32>(faceStore[offset + j]);
v1 = static_cast<uint32>(faceStore[offset + j + 1]);
}
}

edgeSet.insert(static_cast<uint64>(v0) << 32 | v1);
}
}

return edgeSet.size();
}
} // namespace detail

/**
* @brief !!! EXPENSIVE !!! This function is a wrapper method for implicitly determining the correct edge calculation/counting algorithm
* @tparam T indexing type of the face store
* @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList
* @param numVertices This value is optional, since it is exclusively used for determining if faster algorithm is viable
* @return usize This is the number of edges
*/
template <typename T>
usize FindNumEdges(const AbstractDataStore<T>& faceStore, usize numVertices = (detail::k_MaxOptimizedValue + 1))
{
// This case may seem niche, but it is designed with Indexing types in mind specifically IGeometry::MeshIndexType
if constexpr(!std::is_signed_v<T>)
{
if(numVertices < detail::k_MaxOptimizedValue)
{
// speedier method because max vertices value fits into uint32
return detail::FastEdgeCount(faceStore);
}
}
// Slower method to avoid overflow
return detail::SafeEdgeCount(faceStore);
}

/**
* @brief
* @tparam T
Expand Down

0 comments on commit 33f62d4

Please sign in to comment.