From adffcf743c1b7def25e6fe97828f87bf5994f471 Mon Sep 17 00:00:00 2001
From: Sreekanth Arikatla <sreekanth.arikatla@kitware.com>
Date: Tue, 25 Apr 2017 21:33:19 -0400
Subject: [PATCH] ENH: Move the LPC from the nonlinear system to CG solver
Move the Linear Projection Constraints (LPC) from the nonlinear system to CG solver for simplicity.
LPC are only used in CG.
---
Base/Core/imstkAPIUtilities.h | 1 -
Base/Solvers/imstkConjugateGradient.cpp | 15 +++++++-----
Base/Solvers/imstkConjugateGradient.h | 31 +++++++++++++++++++++++++
Base/Solvers/imstkNewtonSolver.cpp | 2 +-
Base/Solvers/imstkNonlinearSystem.h | 21 ++++++++++++++---
Examples/Sandbox/main.cpp | 10 ++++----
6 files changed, 63 insertions(+), 17 deletions(-)
diff --git a/Base/Core/imstkAPIUtilities.h b/Base/Core/imstkAPIUtilities.h
index 41adcc731..3daa959cb 100644
--- a/Base/Core/imstkAPIUtilities.h
+++ b/Base/Core/imstkAPIUtilities.h
@@ -209,7 +209,6 @@ createNonLinearSystem(std::shared_ptr<imstk::FEMDeformableBodyModel> dynaModel)
{
linProj.push_back(LinearProjectionConstraint(i, true));
}
- nlSystem->setLinearProjectors(linProj);
nlSystem->setUnknownVector(dynaModel->getUnknownVec());
nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction());
diff --git a/Base/Solvers/imstkConjugateGradient.cpp b/Base/Solvers/imstkConjugateGradient.cpp
index 793f453b5..2ba8bd966 100644
--- a/Base/Solvers/imstkConjugateGradient.cpp
+++ b/Base/Solvers/imstkConjugateGradient.cpp
@@ -68,7 +68,7 @@ ConjugateGradient::solve(Vectord& x)
return;
}
- if (m_linearSystem->getLinearProjectors().size() == 0)
+ if (m_FixedLinearProjConstraints->size() == 0)
{
x = m_cgSolver.solve(m_linearSystem->getRHSVector());
}
@@ -81,16 +81,17 @@ ConjugateGradient::solve(Vectord& x)
void
ConjugateGradient::modifiedCGSolve(Vectord& x)
{
- const auto &linearProjectors = m_linearSystem->getLinearProjectors();
const auto &b = m_linearSystem->getRHSVector();
const auto &A = m_linearSystem->getMatrix();
// Set the initial guess to zero
x.setZero();
- applyLinearProjectionFilter(x, linearProjectors, true);
+ applyLinearProjectionFilter(x, *m_DynamicLinearProjConstraints, true);
+ applyLinearProjectionFilter(x, *m_FixedLinearProjConstraints, true);
auto res = b;
- applyLinearProjectionFilter(res, linearProjectors, false);
+ applyLinearProjectionFilter(res, *m_DynamicLinearProjConstraints, false);
+ applyLinearProjectionFilter(res, *m_FixedLinearProjConstraints, false);
auto c = res;
auto delta = res.dot(res);
auto deltaPrev = delta;
@@ -103,7 +104,8 @@ ConjugateGradient::modifiedCGSolve(Vectord& x)
while (delta > eps)
{
q = A * c;
- applyLinearProjectionFilter(q, linearProjectors, false);
+ applyLinearProjectionFilter(q, *m_DynamicLinearProjConstraints, false);
+ applyLinearProjectionFilter(q, *m_FixedLinearProjConstraints, false);
dotval = c.dot(q);
if (dotval != 0.0)
{
@@ -120,7 +122,8 @@ ConjugateGradient::modifiedCGSolve(Vectord& x)
delta = res.dot(res);
c *= delta / deltaPrev;
c += res;
- applyLinearProjectionFilter(c, linearProjectors, false);
+ applyLinearProjectionFilter(c, *m_DynamicLinearProjConstraints, false);
+ applyLinearProjectionFilter(c, *m_FixedLinearProjConstraints, false);
if (++iterNum >= m_maxIterations)
{
diff --git a/Base/Solvers/imstkConjugateGradient.h b/Base/Solvers/imstkConjugateGradient.h
index d9dd6a035..6a27c2a0c 100644
--- a/Base/Solvers/imstkConjugateGradient.h
+++ b/Base/Solvers/imstkConjugateGradient.h
@@ -98,6 +98,34 @@ public:
/// \brief Apply a filter to the vector supplied
///
void applyLinearProjectionFilter(Vectord& x, const std::vector<LinearProjectionConstraint>& linProj, const bool setVal);
+
+ /// \brief Get the vector denoting the filter
+ ///
+ void setLinearProjectors(std::vector<LinearProjectionConstraint>* f)
+ {
+ m_FixedLinearProjConstraints = f;
+ }
+
+ /// \brief Get the vector denoting the filter
+ ///
+ std::vector<LinearProjectionConstraint>& getLinearProjectors()
+ {
+ return *m_FixedLinearProjConstraints;
+ }
+
+ /// \brief Get the vector denoting the filter
+ ///
+ void setDynamicLinearProjectors(std::vector<LinearProjectionConstraint>* f)
+ {
+ m_DynamicLinearProjConstraints = f;
+ }
+
+ /// \brief Get the vector denoting the filter
+ ///
+ std::vector<LinearProjectionConstraint>& getDynamicLinearProjectors()
+ {
+ return *m_DynamicLinearProjConstraints;
+ }
private:
///
@@ -107,6 +135,9 @@ private:
///> Pointer to the Eigen's Conjugate gradient solver
Eigen::ConjugateGradient<SparseMatrixd> m_cgSolver;
+
+ std::vector<LinearProjectionConstraint> *m_FixedLinearProjConstraints;
+ std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints;
};
} // imstk
diff --git a/Base/Solvers/imstkNewtonSolver.cpp b/Base/Solvers/imstkNewtonSolver.cpp
index 8938b3917..a01e56688 100644
--- a/Base/Solvers/imstkNewtonSolver.cpp
+++ b/Base/Solvers/imstkNewtonSolver.cpp
@@ -144,7 +144,7 @@ NewtonSolver::updateJacobian(const Vectord& x)
auto &b = m_nonLinearSystem->m_F(x, m_isSemiImplicit);
auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(A, b);
- linearSystem->setLinearProjectors(m_nonLinearSystem->getLinearProjectors());
+ //linearSystem->setLinearProjectors(m_nonLinearSystem->getLinearProjectors()); // TODO: Left for near future reference. Clear in future.
m_linearSolver->setSystem(linearSystem);
return b.dot(b);
diff --git a/Base/Solvers/imstkNonlinearSystem.h b/Base/Solvers/imstkNonlinearSystem.h
index b6a5d5480..95d04c3e1 100644
--- a/Base/Solvers/imstkNonlinearSystem.h
+++ b/Base/Solvers/imstkNonlinearSystem.h
@@ -92,7 +92,7 @@ namespace imstk
/// \brief Get the vector denoting the filter
///
- void setLinearProjectors(std::vector<LinearProjectionConstraint>& f)
+ /*void setLinearProjectors(std::vector<LinearProjectionConstraint>& f)
{
m_LinearProjConstraints = &f;
}
@@ -102,7 +102,7 @@ namespace imstk
std::vector<LinearProjectionConstraint>& getLinearProjectors()
{
return *m_LinearProjConstraints;
- }
+ }*/
///
/// \brief Set the update function
@@ -120,6 +120,20 @@ namespace imstk
m_FUpdatePrevState = updateFunc;
}
+ /// \brief Get the vector denoting the filter
+ ///
+ /*void setDynamicLinearProjectors(std::vector<LinearProjectionConstraint>* f)
+ {
+ m_DynamicLinearProjConstraints = f;
+ }
+
+ /// \brief Get the vector denoting the filter
+ ///
+ std::vector<LinearProjectionConstraint>& getDynamicLinearProjectors()
+ {
+ return *m_DynamicLinearProjConstraints;
+ }*/
+
public:
VectorFunctionType m_F; ///> Nonlinear function
MatrixFunctionType m_dF; ///> Gradient of the Nonlinear function with respect to the unknown vector
@@ -127,7 +141,8 @@ public:
UpdateFunctionType m_FUpdate;
UpdatePrevStateFunctionType m_FUpdatePrevState;
- std::vector<LinearProjectionConstraint> *m_LinearProjConstraints;
+ /*std::vector<LinearProjectionConstraint> *m_LinearProjConstraints;
+ std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints;*/
};
} // imstk
diff --git a/Examples/Sandbox/main.cpp b/Examples/Sandbox/main.cpp
index f224938bf..150b8bbee 100644
--- a/Examples/Sandbox/main.cpp
+++ b/Examples/Sandbox/main.cpp
@@ -1256,11 +1256,9 @@ void testDeformableBody()
for (auto i : dynaModel->getFixNodeIds())
{
auto s = LinearProjectionConstraint(i, false);
- s.setProjectorToDirichlet(i);
- s.setValue(Vec3d(0.001, 0, 0));
+ s.setProjectorToDirichlet(i, Vec3d(0.001, 0, 0));
projList.push_back(s);
}
- nlSystem->setLinearProjectors(projList);
nlSystem->setUnknownVector(dynaModel->getUnknownVec());
nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
@@ -1271,9 +1269,9 @@ void testDeformableBody()
// create a non-linear solver and add to the scene
auto nlSolver = std::make_shared<NewtonSolver>();
+ cgLinSolver->setLinearProjectors(&projList);
nlSolver->setLinearSolver(cgLinSolver);
nlSolver->setSystem(nlSystem);
- //nlSolver->setToFullyImplicit();
scene->addNonlinearSolver(nlSolver);
// Display UPS
@@ -2286,7 +2284,6 @@ void testDeformableBodyCollision()
{
linProj.push_back(LinearProjectionConstraint(id, true));
}
- nlSystem->setLinearProjectors(linProj);
nlSystem->setUnknownVector(dynaModel->getUnknownVec());
nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction());
@@ -2294,6 +2291,7 @@ void testDeformableBodyCollision()
// create a non-linear solver and add to the scene
auto nlSolver = std::make_shared<NewtonSolver>();
auto cgLinSolver = std::make_shared<ConjugateGradient>();// create a linear solver to be used in the NL solver
+ cgLinSolver->setLinearProjectors(&linProj);
nlSolver->setLinearSolver(cgLinSolver);
nlSolver->setSystem(nlSystem);
scene->addNonlinearSolver(nlSolver);
@@ -2386,7 +2384,6 @@ void liverToolInteraction()
{
linProj.push_back(LinearProjectionConstraint(id, true));
}
- nlSystem->setLinearProjectors(linProj);
nlSystem->setUnknownVector(dynaModel->getUnknownVec());
nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction());
@@ -2394,6 +2391,7 @@ void liverToolInteraction()
// create a non-linear solver and add to the scene
auto nlSolver = std::make_shared<NewtonSolver>();
auto cgLinSolver = std::make_shared<ConjugateGradient>();// create a linear solver to be used in the NL solver
+ cgLinSolver->setLinearProjectors(&linProj);
nlSolver->setLinearSolver(cgLinSolver);
nlSolver->setSystem(nlSystem);
//nlSolver->setToFullyImplicit();
--
GitLab