Commit dda1f518 authored by Nghia Truong's avatar Nghia Truong

BUG: Fix bug of Eigen library that cause slow computation due to lazy evaluation

parent e93c27fc
......@@ -314,12 +314,12 @@ void
pointToSphere(const Vec3r& point, uint32_t pointIdx, Sphere* const sphere,
const std::shared_ptr<CollisionData>& colData)
{
const auto sphereCenter = sphere->getPosition();
const auto sphereRadius = sphere->getRadius();
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const Vec3r sphereCenter = sphere->getPosition();
const auto sphereRadius = sphere->getRadius();
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto pc = sphereCenter - point;
const auto distSqr = pc.squaredNorm();
const Vec3r pc = sphereCenter - point;
const auto distSqr = pc.squaredNorm();
if (distSqr < sphereRadiusSqr)
{
const Vec3r direction = distSqr > Real(1e-12) ? pc / std::sqrt(distSqr) : Vec3r(0, 0, 0);
......@@ -344,12 +344,12 @@ void
pointToSpherePicking(const Vec3r& point, uint32_t pointIdx, Sphere* const sphere,
const std::shared_ptr<CollisionData>& colData)
{
const auto sphereCenter = sphere->getPosition();
const auto sphereRadius = sphere->getRadius();
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const Vec3r sphereCenter = sphere->getPosition();
const auto sphereRadius = sphere->getRadius();
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto pc = sphereCenter - point;
const auto distSqr = pc.squaredNorm();
const Vec3r pc = sphereCenter - point;
const auto distSqr = pc.squaredNorm();
if (distSqr < sphereRadiusSqr)
{
colData->NodePickData.safeAppend({ pc, pointIdx, 0 });
......
......@@ -75,9 +75,9 @@ PbdAreaConstraint::solvePositionConstraint(PbdModel& model)
n /= 2 * A;
const auto grad0 = e2.cross(n);
const auto grad1 = e3.cross(n);
const auto grad2 = e1.cross(n);
const Vec3d grad0 = e2.cross(n);
const Vec3d grad1 = e3.cross(n);
const Vec3d grad2 = e1.cross(n);
auto lambda = im0 * grad0.squaredNorm() + im1 * grad1.squaredNorm() + im2 * grad2.squaredNorm();
......
......@@ -77,7 +77,7 @@ PbdConstantDensityConstraint::solvePositionConstraint(PbdModel& model)
double
PbdConstantDensityConstraint::wPoly6(const Vec3d& pi, const Vec3d& pj)
{
double rLengthSqr = (pi - pj).squaredNorm();
double rLengthSqr = (Vec3d(pi - pj)).squaredNorm();
return (rLengthSqr > m_maxDistSqr || rLengthSqr < 1e-20) ?
0 :
......
......@@ -93,10 +93,10 @@ PbdDihedralConstraint::solvePositionConstraint(PbdModel& model)
return false;
}
const auto grad0 = -(l / A1) * n1;
const auto grad1 = -(l / A2) * n2;
const auto grad2 = (e.dot(e1) / (A1 * l)) * n1 + (e.dot(e3) / (A2 * l)) * n2;
const auto grad3 = (e.dot(e2) / (A1 * l)) * n1 + (e.dot(e4) / (A2 * l)) * n2;
const Vec3d grad0 = -(l / A1) * n1;
const Vec3d grad1 = -(l / A2) * n2;
const Vec3d grad2 = (e.dot(e1) / (A1 * l)) * n1 + (e.dot(e3) / (A2 * l)) * n2;
const Vec3d grad3 = (e.dot(e2) / (A1 * l)) * n1 + (e.dot(e4) / (A2 * l)) * n2;
auto lambda = im0 * grad0.squaredNorm() +
im1 * grad1.squaredNorm() +
......
......@@ -72,7 +72,7 @@ neighborSearchBruteForce(StdVectorOfVec3r& particles, std::vector<std::vector<si
}
const auto qpos = particles[q];
const auto d2 = (ppos - qpos).squaredNorm();
const auto d2 = (Vec3r(ppos - qpos)).squaredNorm();
if (d2 < radiusSqr)
{
pneighbors.push_back(q);
......@@ -139,7 +139,7 @@ neighborSearchBruteForce(StdVectorOfVec3r& setA, StdVectorOfVec3r& setB, std::ve
for (size_t q = 0; q < setB.size(); ++q)
{
const auto qpos = setB[q];
const auto d2 = (ppos - qpos).squaredNorm();
const auto d2 = (Vec3r(ppos - qpos)).squaredNorm();
if (d2 < radiusSqr)
{
pneighbors.push_back(q);
......
......@@ -118,7 +118,8 @@ GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>& result,
for (auto q : m_Grid.getCellData(cellX, cellY, cellZ).particleIndices)
{
const auto qpos = setB[q];
const auto d2 = (ppos - qpos).squaredNorm();
const Vec3r diff = ppos - qpos;
const auto d2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
if (d2 < m_SearchRadiusSqr)
{
pneighbors.push_back(q);
......
......@@ -171,7 +171,8 @@ SpatialHashTableSeparateChaining::getPointsInSphere(std::vector<size_t>& result,
for (auto it = first; it != last; ++it)
{
const Vec3d& qpos = it->point;
const auto d2 = (ppos - qpos).squaredNorm();
const Vec3d diff = ppos - qpos;
const auto d2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
if (d2 < radiusSqr)
{
result.push_back(it->ID);
......
......@@ -405,7 +405,7 @@ SPHModel::computeSurfaceTension()
const auto K_ij = Real(2) * m_Parameters->m_RestDensity / (pdensity + qdensity);
// Cohesion acc
auto r = qInfo.xpq;
const auto r = qInfo.xpq;
const auto d2 = r.squaredNorm();
if (d2 > Real(1e-20))
{
......
......@@ -205,7 +205,7 @@ TetraTriangleMap::findClosestTetrahedron(std::shared_ptr<TetrahedralMesh> tetraM
center += tetraMesh->getInitialVertexPosition(vert[i]);
}
center = center / 4.0;
double distSqr = (pos - center).squaredNorm();
double distSqr = (Vec3d(pos - center)).squaredNorm();
if (distSqr < closestDistanceSqr)
{
closestDistanceSqr = distSqr;
......
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