Commit b11575fb authored by Nghia Truong's avatar Nghia Truong
Browse files

PERF: Parallel implementation for mapping classes

parent 386cbccc
......@@ -52,7 +52,7 @@ public:
///
/// \brief Destructor
///
~GeometryMap() = default;
virtual ~GeometryMap() = default;
///
/// \brief Compute the map
......@@ -118,7 +118,7 @@ public:
/// \param idx
/// \return index of Master corresponding to the idx of Slave
///
virtual size_t getMapIdx(const size_t& idx) { return 0; }
virtual size_t getMapIdx(const size_t&) { return 0; }
protected:
///
......
......@@ -20,6 +20,11 @@
=========================================================================*/
#include "imstkOneToOneMap.h"
#include "imstkParallelUtils.h"
#undef min
#undef max
#include <climits>
namespace imstk
{
......@@ -33,36 +38,53 @@ OneToOneMap::compute()
}
// returns the first matching vertex
auto findMatchingVertex = [](std::shared_ptr<PointSet> masterMesh, const Vec3d& p) -> size_t
auto findMatchingVertex = [](const std::shared_ptr<PointSet>& masterMesh, const Vec3d& p, size_t& nodeId)
{
for (size_t nodeId = 0; nodeId < masterMesh->getNumVertices(); ++nodeId)
for (size_t idx = 0; idx < masterMesh->getNumVertices(); ++idx)
{
if (masterMesh->getInitialVertexPosition(nodeId) == p)
if (masterMesh->getInitialVertexPosition(idx) == p)
{
return nodeId;
nodeId = idx;
return true;
}
}
return -1;
return false;
};
auto meshMaster = std::dynamic_pointer_cast<PointSet>(m_master);
auto meshSlave = std::dynamic_pointer_cast<PointSet>(m_slave);
m_oneToOneMap.clear();
for (size_t nodeId = 0; nodeId < meshSlave->getNumVertices(); ++nodeId)
{
// Find the enclosing or closest tetrahedron
size_t matchingNodeId = findMatchingVertex(meshMaster, meshSlave->getVertexPosition(nodeId));
bool bValid = true;
if (matchingNodeId < 0)
ParallelUtils::parallelFor(meshSlave->getNumVertices(),
[&](const size_t nodeId)
{
LOG(WARNING) << "Could not find matching node for the node " << nodeId;
continue;
}
// Find the enclosing or closest tetrahedron
size_t matchingNodeId;
if (!findMatchingVertex(meshMaster, meshSlave->getVertexPosition(nodeId), matchingNodeId))
{
LOG(WARNING) << "Could not find matching node for the node " << nodeId;
bValid = false;
return;
}
// add to the map
// Note: This replaces the map if one with <nodeId> already exists
m_oneToOneMap[nodeId] = matchingNodeId;
});
if (!bValid)
{
m_oneToOneMap.clear();
return;
}
// add to the map
// Note: This replaces the map if one with <nodeId> already exists
m_oneToOneMap[nodeId] = matchingNodeId;
// Copy data from map to vector for parallel processing
m_oneToOneMapVector.resize(0);
for (auto kv: m_oneToOneMap)
{
m_oneToOneMapVector.push_back({kv.first, kv.second});
}
}
......@@ -77,8 +99,8 @@ OneToOneMap::isValid() const
for (auto const& mapValue : m_oneToOneMap)
{
if (mapValue.first >= 0 && mapValue.first < numVertSlave &&
mapValue.second >= 0 && mapValue.second < numVertMaster)
if (mapValue.first < numVertSlave &&
mapValue.second < numVertMaster)
{
continue;
}
......@@ -94,7 +116,14 @@ OneToOneMap::isValid() const
void
OneToOneMap::setMap(const std::map<size_t, size_t>& sourceMap)
{
this->m_oneToOneMap = sourceMap;
m_oneToOneMap = sourceMap;
// Copy data from map to vector for parallel processing
m_oneToOneMapVector.resize(0);
for (auto kv: m_oneToOneMap)
{
m_oneToOneMapVector.push_back({kv.first, kv.second});
}
}
void
......@@ -114,13 +143,16 @@ OneToOneMap::apply()
return;
}
LOG_IF(FATAL, (m_oneToOneMap.size() != m_oneToOneMapVector.size())) << "Internal data is corrupted";
auto meshMaster = std::dynamic_pointer_cast<PointSet>(m_master);
auto meshSlave = std::dynamic_pointer_cast<PointSet>(m_slave);
for (auto const& mapValue : m_oneToOneMap)
{
meshSlave->setVertexPosition(mapValue.first, meshMaster->getVertexPosition(mapValue.second));
}
ParallelUtils::parallelFor(m_oneToOneMapVector.size(),
[&](const size_t idx){
const auto& mapValue = m_oneToOneMapVector[idx];
meshSlave->setVertexPosition(mapValue.first, meshMaster->getVertexPosition(mapValue.second));
});
}
void
......@@ -131,7 +163,7 @@ OneToOneMap::print() const
// Print the one-to-one map
LOG(INFO) << "[slaveVertId, masterVertexId]\n";
for (auto const& mapValue : this->m_oneToOneMap)
for (auto const& mapValue : m_oneToOneMap)
{
LOG(INFO) << "[" << mapValue.first << ", " << mapValue.second << "]\n";
}
......
......@@ -48,7 +48,7 @@ public:
///
/// \brief Default destructor
///
~OneToOneMap() = default;
virtual ~OneToOneMap() override = default;
///
/// \brief Compute the tetra-triangle mesh map
......@@ -88,13 +88,13 @@ public:
///
/// \brief
///
inline size_t getMapIdx(const size_t& idx) override
{
return m_oneToOneMap[idx];
}
size_t getMapIdx(const size_t& idx) override { return m_oneToOneMap[idx]; }
protected:
std::map<size_t, size_t> m_oneToOneMap; ///> One to one mapping data
// This vector is for parallel processing, it should contain identical data as m_oneToOneMap
std::vector<std::pair<size_t, size_t>> m_oneToOneMapVector; ///> One to one mapping data
};
} // imstk
......
......@@ -20,6 +20,7 @@
=========================================================================*/
#include "imstkTetraTriangleMap.h"
#include "imstkParallelUtils.h"
namespace imstk
{
......@@ -37,26 +38,46 @@ TetraTriangleMap::compute()
m_verticesEnclosingTetraId.clear();
m_verticesWeights.clear();
for (const Vec3d& surfVertPos : triMesh->getVertexPositions())
{
// Find the enclosing or closest tetrahedron
size_t closestTetId = findEnclosingTetrahedron(tetMesh, surfVertPos);
if (closestTetId == std::numeric_limits<size_t>::max())
{
closestTetId = findClosestTetrahedron(tetMesh, surfVertPos);
}
if (closestTetId == std::numeric_limits<size_t>::max())
{
LOG(WARNING) << "Could not find closest tetrahedron";
return;
}
// Compute the weights
TetrahedralMesh::WeightsArray weights = { 0.0, 0.0, 0.0, 0.0 };
tetMesh->computeBarycentricWeights(closestTetId, surfVertPos, weights);
ParallelUtils::SpinLock lock;
bool bValid = true;
m_verticesEnclosingTetraId.push_back(closestTetId); // store nearest tetrahedron
m_verticesWeights.push_back(weights); // store weights
ParallelUtils::parallelFor(triMesh->getNumVertices(),
[&](const size_t vertexIdx)
{
if (!bValid)
{
return;
}
const Vec3d& surfVertPos = triMesh->getVertexPositions()[vertexIdx];
// Find the enclosing or closest tetrahedron
size_t closestTetId = findEnclosingTetrahedron(tetMesh, surfVertPos);
if (closestTetId == std::numeric_limits<size_t>::max())
{
closestTetId = findClosestTetrahedron(tetMesh, surfVertPos);
}
if (closestTetId == std::numeric_limits<size_t>::max())
{
LOG(WARNING) << "Could not find closest tetrahedron";
bValid = false;
return;
}
// Compute the weights
TetrahedralMesh::WeightsArray weights = { 0.0, 0.0, 0.0, 0.0 };
tetMesh->computeBarycentricWeights(closestTetId, surfVertPos, weights);
lock.lock();
m_verticesEnclosingTetraId.push_back(closestTetId); // store nearest tetrahedron
m_verticesWeights.push_back(weights); // store weights
lock.unlock();
});
// Clear result if could not find closest tet
if (!bValid)
{
m_verticesEnclosingTetraId.clear();
m_verticesWeights.clear();
}
}
......@@ -80,19 +101,21 @@ TetraTriangleMap::apply()
auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh> (m_master);
auto triMesh = std::dynamic_pointer_cast<SurfaceMesh> (m_slave);
Vec3d newPos;
for (size_t vertexId = 0; vertexId < triMesh->getNumVertices(); ++vertexId)
{
newPos.setZero();
size_t tetId = m_verticesEnclosingTetraId.at(vertexId);
auto tetVerts = tetMesh->getTetrahedronVertices(tetId);
auto weights = m_verticesWeights.at(vertexId);
for (size_t i = 0; i < 4; ++i)
ParallelUtils::parallelFor(triMesh->getNumVertices(),
[&](const size_t vertexId)
{
newPos += tetMesh->getVertexPosition(tetVerts[i]) * weights[i];
}
triMesh->setVertexPosition(vertexId, newPos);
}
Vec3d newPos(0, 0, 0);
size_t tetId = m_verticesEnclosingTetraId[vertexId];
auto tetVerts = tetMesh->getTetrahedronVertices(tetId);
auto weights = m_verticesWeights[vertexId];
for (size_t i = 0; i < 4; ++i)
{
newPos += tetMesh->getVertexPosition(tetVerts[i]) * weights[i];
}
// This writes newPos to position array at individual vertexId, thus should not cause race condition
triMesh->setVertexPosition(vertexId, newPos);
});
}
void
......@@ -123,8 +146,7 @@ TetraTriangleMap::isValid() const
for (size_t tetId = 0; tetId < m_verticesEnclosingTetraId.size(); ++tetId)
{
if (!(m_verticesEnclosingTetraId.at(tetId) < totalElementsMaster &&
m_verticesEnclosingTetraId.at(tetId) >= 0))
if (!(m_verticesEnclosingTetraId[tetId] < totalElementsMaster))
{
return false;
}
......@@ -155,10 +177,9 @@ TetraTriangleMap::setSlave(std::shared_ptr<Geometry> slave)
}
size_t
TetraTriangleMap::findClosestTetrahedron(std::shared_ptr<TetrahedralMesh> tetraMesh,
const Vec3d& pos)
TetraTriangleMap::findClosestTetrahedron(std::shared_ptr<TetrahedralMesh> tetraMesh, const Vec3d& pos)
{
double closestDistance = MAX_D;
double closestDistanceSqr = MAX_D;
size_t closestTetrahedron = std::numeric_limits<size_t>::max();
for (size_t tetId = 0; tetId < tetraMesh->getNumTetrahedra(); ++tetId)
{
......@@ -168,11 +189,11 @@ TetraTriangleMap::findClosestTetrahedron(std::shared_ptr<TetrahedralMesh> tetraM
{
center += tetraMesh->getInitialVertexPosition(vert[i]);
}
center = center / 4;
double dist = (pos - center).norm();
if (dist < closestDistance)
center = center / 4.0;
double distSqr = (pos - center).squaredNorm();
if (distSqr < closestDistanceSqr)
{
closestDistance = dist;
closestDistanceSqr = distSqr;
closestTetrahedron = tetId;
}
}
......@@ -180,8 +201,7 @@ TetraTriangleMap::findClosestTetrahedron(std::shared_ptr<TetrahedralMesh> tetraM
}
size_t
TetraTriangleMap::findEnclosingTetrahedron(std::shared_ptr<TetrahedralMesh> tetraMesh,
const Vec3d& pos)
TetraTriangleMap::findEnclosingTetrahedron(std::shared_ptr<TetrahedralMesh> tetraMesh, const Vec3d& pos)
{
Vec3d boundingBoxMin;
Vec3d boundingBoxMax;
......
......@@ -48,7 +48,7 @@ public:
///
/// \brief Destructor
///
virtual ~TetraTriangleMap() = default;
virtual ~TetraTriangleMap() override = default;
///
/// \brief Compute the tetra-triangle mesh map
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment