diff --git a/Base/Core/imstkAPIUtilities.h b/Base/Core/imstkAPIUtilities.h index 41adcc7311390d94f2821cb497cae0e564f6aa86..3daa959cb104c54e5c772b45756b21ea3ccbf529 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 793f453b529352aa064c94900a8bc939f8937d37..2ba8bd96627588efce34e875006a5728b4cbae37 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 d9dd6a035ebf559f8a017c862b798955a91909a5..6a27c2a0c1d3e58fa2b3239a00a0b6aae24f4e08 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 8938b39178ca6c249c3884017cafd56a7256954c..a01e566881651116eef52f22ad0724f976a88993 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 b6a5d54808b50a4e609a5e273f9d6c79a58559aa..95d04c3e1b68f52bfa85252ec642651249e7dc46 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 f224938bf0b363aa68d2db5e984fb44335f103b8..150b8bbee217004a82038925c3c8f191511b5bea 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();