diff --git a/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/LinearSolverConstraintCorrection.inl b/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/LinearSolverConstraintCorrection.inl index 4943e17ab4c..305166be82e 100644 --- a/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/LinearSolverConstraintCorrection.inl +++ b/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/LinearSolverConstraintCorrection.inl @@ -203,7 +203,7 @@ void LinearSolverConstraintCorrection::addComplianceInConstraintSpace } // use the Linear solver to compute J*A^-1*J^T, where A is the mechanical linear system matrix - l_linearSolver->setSystemLHVector(sofa::core::MultiVecDerivId::null()); + l_linearSolver->getLinearSystem()->setSystemSolution(sofa::core::MultiVecDerivId::null()); l_linearSolver->addJMInvJt(W, &m_constraintJacobian, factor); addRegularization(W); @@ -247,9 +247,10 @@ void LinearSolverConstraintCorrection< DataTypes >::computeMotionCorrection(cons { if (mstate && l_linearSolver.get()) { - l_linearSolver->setSystemRHVector(f); - l_linearSolver->setSystemLHVector(dx); + l_linearSolver->getLinearSystem()->setRHS(f); + l_linearSolver->getLinearSystem()->setSystemSolution(dx); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(dx); } } @@ -373,9 +374,10 @@ void LinearSolverConstraintCorrection::applyContactForce(const linear } } } - l_linearSolver->setSystemRHVector(forceID); - l_linearSolver->setSystemLHVector(dxID); + l_linearSolver->getLinearSystem()->setRHS(forceID); + l_linearSolver->getLinearSystem()->setSystemSolution(dxID); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(dxID); //TODO: tell the solver not to recompute the matrix @@ -541,13 +543,13 @@ void LinearSolverConstraintCorrection::resetForUnbuiltResolution(SRea core::VecDerivId forceID(core::VecDerivId::V_FIRST_DYNAMIC_INDEX); core::VecDerivId dxID = core::vec_id::write_access::dx; - l_linearSolver->setSystemRHVector(forceID); - l_linearSolver->setSystemLHVector(dxID); + l_linearSolver->getLinearSystem()->setRHS(forceID); + l_linearSolver->getLinearSystem()->setSystemSolution(dxID); - systemMatrix_buf = l_linearSolver->getSystemBaseMatrix(); - systemRHVector_buf = l_linearSolver->getSystemRHBaseVector(); - systemLHVector_buf = l_linearSolver->getSystemLHBaseVector(); + systemMatrix_buf = l_linearSolver->getLinearSystem()->getSystemBaseMatrix(); + systemRHVector_buf = l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(); + systemLHVector_buf = l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector(); systemLHVector_buf_fullvector = dynamic_cast*>(systemLHVector_buf); // Cast checking whether the LH vector is a FullVector to improve performances constexpr const auto derivDim = Deriv::total_size; diff --git a/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/PrecomputedConstraintCorrection.inl b/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/PrecomputedConstraintCorrection.inl index 6a81c4d3955..b50aac185e2 100644 --- a/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/PrecomputedConstraintCorrection.inl +++ b/Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/PrecomputedConstraintCorrection.inl @@ -355,9 +355,6 @@ void PrecomputedConstraintCorrection::bwdInit() fact *= eulerSolver->getPositionIntegrationFactor(); // here, we compute a compliance eulerSolver->solve(core::execparams::defaultInstance(), dt, core::vec_id::write_access::position, core::vec_id::write_access::velocity); - - if (linearSolver) - linearSolver->freezeSystemMatrix(); // do not recompute the matrix for the rest of the precomputation } for (unsigned int v = 0; v < nbNodes; v++) @@ -374,10 +371,6 @@ void PrecomputedConstraintCorrection::bwdInit() msg_info() << tmpStr.str(); } - // Do not recompute the matrix for the rest of the precomputation - if (linearSolver) - linearSolver->freezeSystemMatrix(); - saveCompliance(invName); // Restore gravity diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.h index d64eef84488..a3bd4624c0d 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.h @@ -62,21 +62,19 @@ class AsyncSparseLDLSolver : public SparseLDLSolver newInvertDataReady { false }; - bool m_hasUpdatedMatrix { false }; + Data< bool > d_enableAssembly; }; #if !defined(SOFA_COMPONENT_LINEARSOLVER_ASYNCSPARSELDLSOLVER_CPP) diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl index 2965a02d60f..21ada3b9aa6 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl @@ -29,44 +29,58 @@ namespace sofa::component::linearsolver::direct { template -void AsyncSparseLDLSolver::init() +AsyncSparseLDLSolver::AsyncSparseLDLSolver() + : d_enableAssembly(initData(&d_enableAssembly, true, "enableAssembly", "Allow assembly of the linear system")) { - Inherit1::init(); - waitForAsyncTask = true; - m_asyncThreadInvertData = &m_secondInvertData; - m_mainThreadInvertData = static_cast(this->invertData.get()); } template -void AsyncSparseLDLSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) +void AsyncSparseLDLSolver::init() { - if (isAsyncFactorizationFinished() || !m_asyncResult.valid()) + Inherit1::init(); + + if (!this->isComponentStateInvalid()) { - SCOPED_TIMER_VARNAME(setSystemMBKMatrixTimer, "setSystemMBKMatrix"); - Inherit1::setSystemMBKMatrix(mparams); - m_hasUpdatedMatrix = true; + if (this->l_linearSystem) + { + this->l_linearSystem->d_enableAssembly.setParent(&d_enableAssembly); + } + + waitForAsyncTask = true; + m_asyncThreadInvertData = &m_secondInvertData; + m_mainThreadInvertData = static_cast(this->invertData.get()); } } +template +void AsyncSparseLDLSolver::reset() +{ + d_enableAssembly.setValue(true); + waitForAsyncTask = true; +} + template void AsyncSparseLDLSolver::solveSystem() { - SCOPED_TIMER_VARNAME(invertDataCopyTimer, "AsyncSolve"); + SCOPED_TIMER_VARNAME(invertDataCopyTimer, "solveSystem"); if (newInvertDataReady) { swapInvertData(); } - if (this->linearSystem.needInvert) + if (this->d_factorizationInvalidation.getValue()) { if (this->invertData == nullptr) { - this->getMatrixInvertData(this->getSystemMatrix()); + this->getMatrixInvertData(this->l_linearSystem->getSystemMatrix()); m_mainThreadInvertData = static_cast(this->invertData.get()); } launchAsyncFactorization(); - this->linearSystem.needInvert = false; + this->d_factorizationInvalidation.setValue(false); + + //matrix assembly is temporarily stopped until the next factorization + d_enableAssembly.setValue(false); } if (waitForAsyncTask) @@ -81,15 +95,11 @@ void AsyncSparseLDLSolver::solveSystem() swapInvertData(); } - this->solve(*this->getSystemMatrix(), *this->getSystemLHVector(), *this->getSystemRHVector()); + auto* A = this->l_linearSystem->getSystemMatrix(); + auto* b = this->l_linearSystem->getRHSVector(); + auto* x = this->l_linearSystem->getSolutionVector(); - if (!this->linearSystem.solutionVecId.isNull()) - { - if (this->l_linearSystem) - { - this->l_linearSystem->dispatchSystemSolution(this->linearSystem.solutionVecId); - } - } + this->solve(*A, *x, *b); } template @@ -107,8 +117,8 @@ void AsyncSparseLDLSolver::invert(TMatrix& M) } template -bool AsyncSparseLDLSolver::addJMInvJtLocal(TMatrix* M, ResMatrixType* result, - const JMatrixType* J, SReal fact) +bool AsyncSparseLDLSolver::addJMInvJtLocal( + TMatrix* M, ResMatrixType* result, const JMatrixType* J, SReal fact) { SOFA_UNUSED(M); @@ -119,18 +129,6 @@ bool AsyncSparseLDLSolver::addJMInvJtLocal(TMa return Inherit1::doAddJMInvJtLocal(result, J, fact, m_mainThreadInvertData); } -template -bool AsyncSparseLDLSolver::hasUpdatedMatrix() -{ - return m_hasUpdatedMatrix; -} - -template -void AsyncSparseLDLSolver::updateSystemMatrix() -{ - m_hasUpdatedMatrix = false; -} - template AsyncSparseLDLSolver::~AsyncSparseLDLSolver() { @@ -154,9 +152,14 @@ void AsyncSparseLDLSolver::launchAsyncFactoriz template void AsyncSparseLDLSolver::asyncFactorization() { + SCOPED_TIMER_TR("asyncFactorization"); + newInvertDataReady = false; - this->invert(*this->getSystemMatrix()); + this->invert(*this->l_linearSystem->getSystemMatrix()); newInvertDataReady = true; + + //factorization is finished: matrix assembly is authorized once again + d_enableAssembly.setValue(true); } template diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/BTDLinearSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/BTDLinearSolver.inl index 5628a90ab71..ee609309bcf 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/BTDLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/BTDLinearSolver.inl @@ -394,7 +394,7 @@ void BTDLinearSolver::init_partial_solve() { constexpr Index bsize = Matrix::getSubMatrixDim(); - const Index nb = this->getSystemRHVector()->size() / bsize; + const Index nb = this->l_linearSystem->getRHSVector()->size() / bsize; //TODO => optimisation ?? bwdContributionOnLH.clear(); @@ -419,7 +419,7 @@ void BTDLinearSolver::init_partial_solve() { Vec_dRH[i]=0; Vec_dRH[i].resize(bsize); - _rh_buf.asub(i,bsize) = this->getSystemRHVector()->asub(i,bsize) ; + _rh_buf.asub(i,bsize) = this->l_linearSystem->getRHSVector()->asub(i,bsize) ; } @@ -456,7 +456,7 @@ void BTDLinearSolver::bwdAccumulateRHinBloc(Index indMaxBloc) { // evaluate the Right Hand Term for the block b - RHbloc = this->getSystemRHVector()->asub(b,bsize) ; + RHbloc = this->l_linearSystem->getRHSVector()->asub(b,bsize) ; // compute the contribution on LH created by RH _acc_lh_bloc += Minv.asub(b,b,bsize,bsize) * RHbloc; @@ -474,7 +474,7 @@ void BTDLinearSolver::bwdAccumulateRHinBloc(Index indMaxBloc) b = current_bloc; // compute the block which indice is current_bloc - this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) + + this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) + bwdContributionOnLH.asub(b, bsize); dmsg_info_when(showProblem) << "LH[" << b << "] = Minv[" << b << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << ")) + bwdLH(" << b << ")"; @@ -500,7 +500,7 @@ void BTDLinearSolver::bwdAccumulateLHGlobal( ) dmsg_info_when(showProblem) << "bwdLH[" << current_bloc - 1 << "] = H[" << current_bloc - 1 << "][" << current_bloc << "] *( bwdLH[" << current_bloc << "] + Minv[" << current_bloc << "][" << current_bloc << "] * RH[" << current_bloc << "])"; // BwdLH += Minv*RH - _acc_lh_bloc += Minv.asub(current_bloc,current_bloc,bsize,bsize) * this->getSystemRHVector()->asub(current_bloc,bsize) ; + _acc_lh_bloc += Minv.asub(current_bloc,current_bloc,bsize,bsize) * this->l_linearSystem->getRHSVector()->asub(current_bloc,bsize) ; current_bloc--; // BwdLH(n-1) = H(n-1)(n)*BwdLH(n) @@ -540,7 +540,7 @@ void BTDLinearSolver::fwdAccumulateRHGlobal(Index indMinBloc) { // fwdRH(n) += RH(n) - _acc_rh_bloc += this->getSystemRHVector()->asub(current_bloc,bsize); + _acc_rh_bloc += this->l_linearSystem->getRHSVector()->asub(current_bloc,bsize); // fwdRH(n+1) = H(n+1)(n) * fwdRH(n) _acc_rh_bloc = -(lambda[current_bloc].t() * _acc_rh_bloc); @@ -556,7 +556,7 @@ void BTDLinearSolver::fwdAccumulateRHGlobal(Index indMinBloc) Index b = current_bloc; // compute the block which indice is _indMaxFwdLHComputed - this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) + + this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) + bwdContributionOnLH.asub(b, bsize); dmsg_info_when(showProblem) << "LH[" << b << "] = Minv[" << b << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << ")) + bwdLH(" << b << ")"; @@ -584,13 +584,13 @@ void BTDLinearSolver::fwdComputeLHinBloc(Index indMaxBloc) { dmsg_info_when(showProblem) << " fwdRH[" << b + 1 << "] = H[" << b + 1 << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << "))"; // fwdRH(n+1) = H(n+1)(n) * (fwdRH(n) + RH(n)) - fwdContributionOnRH.asub(b+1, bsize) = (-lambda[b].t())* ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) ; + fwdContributionOnRH.asub(b+1, bsize) = (-lambda[b].t())* ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) ; } _indMaxFwdLHComputed++; b++; // compute the block which indice is _indMaxFwdLHComputed - this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) + + this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) + bwdContributionOnLH.asub(b, bsize); dmsg_info_when(showProblem) << "LH["<::partial_solve(ListIndex& Iout, ListIndex& { constexpr Index bsize = Matrix::getSubMatrixDim(); Vector *Result_partial_Solve = new Vector(); - (*Result_partial_Solve) = (*this->getSystemLHVector()); + (*Result_partial_Solve) = (*this->l_linearSystem->getSolutionVector()); - solve(*this->getSystemMatrix(),*this->getSystemLHVector(), *this->getSystemRHVector()); + solve(*this->l_linearSystem->getSystemMatrix(),*this->l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector()); Vector *Result = new Vector(); - (*Result) = (*this->getSystemLHVector()); + (*Result) = (*this->l_linearSystem->getSolutionVector()); Vector *DR = new Vector(); (*DR) = (*Result); diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h index ad0bfbdb5d2..afb07d31872 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h @@ -82,12 +82,11 @@ class PrecomputedLinearSolver : public sofa::component::linearsolver::MatrixLine Data d_jmjt_twostep; ///< Use two step algorithm to compute JMinvJt Data d_use_file; ///< Dump system matrix in a file - Data init_Tolerance; + Data init_Tolerance; PrecomputedLinearSolver(); - void solve (TMatrix& M, TVector& x, TVector& b) override; + void solve (TMatrix& M, TVector& solution, TVector& rh) override; void invert(TMatrix& M) override; - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; void loadMatrix(TMatrix& M); void loadMatrixWithCholeskyDecomposition(TMatrix& M); bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; @@ -120,8 +119,8 @@ protected : private : bool first; unsigned systemSize; - double dt; - double factInt; + Real dt; + Real factInt; std::vector isActiveDofs; }; diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl index 149aca006eb..bf5100e7363 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl @@ -53,28 +53,16 @@ PrecomputedLinearSolver::PrecomputedLinearSolver() } template -void PrecomputedLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) +void PrecomputedLinearSolver::solve (TMatrix& M, TVector& solution, TVector& rh) { - // Update the matrix only the first time - if (first) - { - first = false; - Inherit::setSystemMBKMatrix(mparams); - loadMatrix(*this->getSystemMatrix()); - } -} - -//Solve x = R * M^-1 * R^t * b -template -void PrecomputedLinearSolver::solve (TMatrix& , TVector& z, TVector& r) -{ - z = internalData.Minv * r; + SOFA_UNUSED(M); + solution = internalData.Minv * rh; } template void PrecomputedLinearSolver::loadMatrix(TMatrix& M) { - systemSize = this->getSystemMatrix()->rowSize(); + systemSize = this->l_linearSystem->getSystemMatrix()->rowSize(); internalData.Minv.resize(systemSize,systemSize); dt = this->getContext()->getDt(); @@ -151,7 +139,18 @@ void PrecomputedLinearSolver::loadMatrixWithCholeskyDecompositi } template -void PrecomputedLinearSolver::invert(TMatrix& /*M*/) {} +void PrecomputedLinearSolver::invert(TMatrix& /*M*/) +{ + if (first) + { + if (this->l_linearSystem) + { + loadMatrix(*this->l_linearSystem->getSystemMatrix()); + this->l_linearSystem->d_enableAssembly.setValue(false); + } + first = false; + } +} template template void PrecomputedLinearSolver::computeActiveDofs(JMatrix& J) @@ -194,7 +193,6 @@ bool PrecomputedLinearSolver::addJMInvJt(linearalgebra::BaseMat msg_error() << "The construction of the matrix when the solver is used only as cvonstraint " "correction is not implemented. You first need to save the matrix into a file. " ; - setSystemMBKMatrix(&mparams); } if (SparseMatrix* j = dynamic_cast*>(J)) diff --git a/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt b/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt index 50514787466..300a11f361a 100644 --- a/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt +++ b/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt @@ -20,6 +20,8 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/ShewchukPCGLinearSolver.inl ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.h ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.inl + ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PreconditionedMatrixFreeSystem.h + ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PreconditionedMatrixFreeSystem.inl ) set(SOURCE_FILES @@ -32,6 +34,7 @@ set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MatrixLinearSystem[GraphScattered].cpp ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MinResLinearSolver.cpp ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.cpp + ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PreconditionedMatrixFreeSystem.cpp ) sofa_find_package(Sofa.Simulation.Core REQUIRED) diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.h index 46783dd8559..77df33d5274 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.h @@ -64,10 +64,6 @@ class CGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver< void init() override; void reinit() override {}; - void resetSystem() override; - - void setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) override; - /// Solve iteratively the linear system Ax=b following a conjugate gradient descent void solve (Matrix& A, Vector& x, Vector& b) override; }; diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.inl index 51a042632df..725aff86e01 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/CGLinearSolver.inl @@ -68,27 +68,6 @@ void CGLinearSolver::init() equilibriumReached = false; } -/// Clear graph and clean the RHS / LHS vectors -template -void CGLinearSolver::resetSystem() -{ - d_graph.beginEdit()->clear(); - d_graph.endEdit(); - - Inherit::resetSystem(); -} - -/// For unbuilt approach (e.g. with GraphScattered types), -/// it passes the coefficients multiplying the matrices M, B and K from the ODE to the LinearSolver (MechanicalOperations::setKFactor) and includes a resetSystem -/// In other cases -/// the global system matrix is setup (pass coefficients with MechanicalOperations::setKFactor) and built it by calling the addMBKToMatrix visitor -template -void CGLinearSolver::setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) -{ - SCOPED_TIMER("CG-setSystemMBKMatrix"); - Inherit::setSystemMBKMatrix(mparams); -} - /// Solve iteratively the linear system Ax=b following a conjugate gradient descent template void CGLinearSolver::solve(Matrix& A, Vector& x, Vector& b) diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].cpp index 1da21eb1bbd..b5912ebc54f 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].cpp @@ -21,6 +21,7 @@ ******************************************************************************/ #define SOFA_COMPONENT_LINEARSOLVER_MATRIXFREESYSTEM_GRAPHSCATTERED_CPP #include +#include #include namespace sofa::component::linearsystem @@ -30,11 +31,14 @@ using sofa::component::linearsolver::GraphScatteredMatrix; using sofa::component::linearsolver::GraphScatteredVector; template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API MatrixFreeSystem; +template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API CompositeLinearSystem; void registerMatrixFreeSystemGraphScattered(sofa::core::ObjectFactory* factory) { factory->registerObjects(core::ObjectRegistrationData("Matrix-free (unbuilt) linear system.") - .add< MatrixFreeSystem >()); + .add< MatrixFreeSystem >(true)); + factory->registerObjects(core::ObjectRegistrationData("Matrix-free (unbuilt) linear system.") + .add< CompositeLinearSystem >()); } } //namespace sofa::component::linearsystem diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].h index 68d0ba15516..a7232d6705b 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixFreeSystem[GraphScattered].h @@ -23,6 +23,7 @@ #include #include +#include #include namespace sofa::component::linearsystem @@ -32,7 +33,8 @@ using sofa::component::linearsolver::GraphScatteredMatrix; using sofa::component::linearsolver::GraphScatteredVector; #if !defined(SOFA_COMPONENT_LINEARSOLVER_MATRIXFREESYSTEM_GRAPHSCATTERED_CPP) - extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API MatrixFreeSystem; +extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API MatrixFreeSystem; +extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API CompositeLinearSystem; #endif } //namespace sofa::component::linearsystem diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.cpp index 0102405a53b..20ae5d57cbb 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.cpp @@ -34,56 +34,6 @@ namespace sofa::component::linearsolver using sofa::core::behavior::LinearSolver; using sofa::core::objectmodel::BaseContext; - -template<> -void MatrixLinearSolver::resetSystem() -{ - if (l_linearSystem) - { - l_linearSystem->clearSystem(); - } - linearSystem.solutionVecId = core::MultiVecDerivId::null(); - linearSystem.needInvert = true; - -} - -template<> -void MatrixLinearSolver::resizeSystem(Size n) -{ - if (l_linearSystem) - { - l_linearSystem->resizeSystem(n); - } - linearSystem.solutionVecId = core::MultiVecDerivId::null(); - linearSystem.needInvert = true; -} - -template<> -void MatrixLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - resetSystem(); - if (auto* matrix = getSystemMatrix()) - { - matrix->setMBKFacts(mparams); - } -} - -template<> -void MatrixLinearSolver::rebuildSystem(SReal /*massFactor*/, SReal /*forceFactor*/) -{ -} - -template<> -void MatrixLinearSolver::setSystemLHVector(core::MultiVecDerivId v) -{ - linearSystem.solutionVecId = v; - getSystemLHVector()->set(v); -} - -template<> -void MatrixLinearSolver::applySystemSolution() -{} - template<> void MatrixLinearSolver::applyConstraintForce(const sofa::core::ConstraintParams* /*cparams*/, sofa::core::MultiVecDerivId /*dx*/, const linearalgebra::BaseVector* /*f*/) { @@ -100,15 +50,6 @@ GraphScatteredVector* MatrixLinearSolver -linearalgebra::BaseMatrix* MatrixLinearSolver::getSystemBaseMatrix() { return nullptr; } - -template<> -linearalgebra::BaseVector* MatrixLinearSolver::getSystemRHBaseVector() { return nullptr; } - -template<> -linearalgebra::BaseVector* MatrixLinearSolver::getSystemLHBaseVector() { return nullptr; } - template<> void MatrixLinearSolver::checkLinearSystem() { diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.h index a48dd2c2c89..f0f134841bc 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.h @@ -63,7 +63,8 @@ class BaseMatrixLinearSolver : public sofa::core::behavior::LinearSolver virtual void solve(Matrix& M, Vector& solution, Vector& rh) = 0; - virtual Matrix * getSystemMatrix() = 0; + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual Matrix * getSystemMatrix() final = delete; }; @@ -191,82 +192,30 @@ class MatrixLinearSolver : public BaseMatrixLinea void init() override; /// Reset the current linear system. - void resetSystem() override; - - /// Reset the current linear system. - void resizeSystem(Size n); - - /// Set the linear system matrix, combining the mechanical M,B,K matrices using the given coefficients - /// - /// Note that this automatically resizes the linear system to the number of active degrees of freedoms - /// - /// @todo Should we put this method in a specialized class for mechanical systems, or express it using more general terms (i.e. coefficients of the second order ODE to solve) - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; - - /// Rebuild the system using a mass and force factor - void rebuildSystem(SReal massFactor, SReal forceFactor) override; - - /// Set the linear system right-hand term vector, from the values contained in the (Mechanical/Physical)State objects - void setSystemRHVector(core::MultiVecDerivId v) override; - - /// Set the initial estimate of the linear system left-hand term vector, from the values contained in the (Mechanical/Physical)State objects - /// This vector will be replaced by the solution of the system once solveSystem is called - void setSystemLHVector(core::MultiVecDerivId v) override; - - /// Get the linear system matrix, or nullptr if this solver does not build it - Matrix* getSystemMatrix() override; + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + void resizeSystem(Size n) = delete; /// Get the linear system right-hand term vector, or nullptr if this solver does not build it + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_DEPRECATED_ASSEMBLY_API() Vector* getSystemRHVector() { return l_linearSystem ? l_linearSystem->getRHSVector() : nullptr; } /// Get the linear system left-hand term vector, or nullptr if this solver does not build it + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_DEPRECATED_ASSEMBLY_API() Vector* getSystemLHVector() { return l_linearSystem ? l_linearSystem->getSolutionVector() : nullptr; } - /// Get the linear system matrix, or nullptr if this solver does not build it - linearalgebra::BaseMatrix* getSystemBaseMatrix() override; - - /// Get the linear system right-hand term vector, or nullptr if this solver does not build it - linearalgebra::BaseVector* getSystemRHBaseVector() override { return l_linearSystem ? l_linearSystem->getRHSVector() : nullptr; } - - /// Get the linear system left-hand term vector, or nullptr if this solver does not build it - linearalgebra::BaseVector* getSystemLHBaseVector() override { return l_linearSystem ? l_linearSystem->getSolutionVector() : nullptr; } - /// Returns the linear system component associated to the linear solver - sofa::component::linearsystem::TypedMatrixLinearSystem* getLinearSystem() const { return l_linearSystem.get(); } + sofa::component::linearsystem::TypedMatrixLinearSystem* getLinearSystem() const override { return l_linearSystem.get(); } /// Solve the system as constructed using the previous methods void solveSystem() override; /// Apply the solution of the system to all the objects - void applySystemSolution(); + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + void applySystemSolution() = delete; /// Invert the system, this method is optional because it's call when solveSystem() is called for the first time void invertSystem() override; - void prepareVisitor(simulation::Visitor* v) - { - v->setTags(this->getTags()); - } - - void prepareVisitor(simulation::BaseMechanicalVisitor* v) - { - prepareVisitor(static_cast(v)); - } - - template - void executeVisitor(T v) - { - prepareVisitor(&v); - v.execute( this->getContext() ); - } - - template - void executeVisitor(T* v) - { - prepareVisitor(v); - v->execute( this->getContext() ); - } - /// Implementing the GetCustomTemplateName is mandatory to have a custom template name parameters /// instead of the default one generated automatically by the SOFA_CLASS() macro. static std::string GetCustomTemplateName() @@ -279,6 +228,8 @@ class MatrixLinearSolver : public BaseMatrixLinea return ThreadManager::isAsyncSolver(); } + virtual void invertIfInvalidated(Matrix& M) final; + void invert(Matrix& /*M*/) override {} void solve(Matrix& M, Vector& solution, Vector& rh) override = 0; @@ -301,12 +252,41 @@ class MatrixLinearSolver : public BaseMatrixLinea ///< matrix of the linear system and J is any matrix with compatible dimensions Data d_parallelInverseProduct; -public: + SingleLink< + MatrixLinearSolver, + sofa::component::linearsystem::TypedMatrixLinearSystem, + BaseLink::FLAG_STOREPATH|BaseLink::FLAG_STRONGLINK + > l_linearSystem; MatrixInvertData * getMatrixInvertData(linearalgebra::BaseMatrix * m); protected: + void prepareVisitor(simulation::Visitor* v) + { + v->setTags(this->getTags()); + } + + void prepareVisitor(simulation::BaseMechanicalVisitor* v) + { + prepareVisitor(static_cast(v)); + } + + template + void executeVisitor(T v) + { + prepareVisitor(&v); + v.execute( this->getContext() ); + } + + template + void executeVisitor(T* v) + { + prepareVisitor(v); + v->execute( this->getContext() ); + } + + virtual void checkLinearSystem(); /** @@ -351,45 +331,15 @@ class MatrixLinearSolver : public BaseMatrixLinea virtual MatrixInvertData * createInvertData(); - struct LinearSystemData - { - bool needInvert; - Matrix* systemMatrix; - Vector* systemRHVector; - Vector* systemLHVector; - core::MultiVecDerivId solutionVecId; - -#if SOFA_CORE_ENABLE_CRSMULTIMATRIXACCESSOR - core::behavior::CRSMultiMatrixAccessor matrixAccessor; -#else - core::behavior::DefaultMultiMatrixAccessor matrixAccessor; -#endif // SOFA_CORE_ENABLE_CRSMULTIMATRIXACCESSOR - - LinearSystemData() - : needInvert(true), systemMatrix(nullptr), systemRHVector(nullptr), systemLHVector(nullptr), - solutionVecId(core::MultiVecDerivId::null()) - {} - ~LinearSystemData() - { - if (systemMatrix) deleteMatrix(systemMatrix); - if (systemRHVector) deletePersistentVector(systemRHVector); - if (systemLHVector) deletePersistentVector(systemLHVector); - } - }; - - LinearSystemData linearSystem; + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved linearSystem; - SReal currentMFactor, currentBFactor, currentKFactor; + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved currentMFactor, currentBFactor, currentKFactor; bool singleThreadAddJMInvJtLocal(Matrix * /*M*/,ResMatrixType * result,const JMatrixType * J, SReal fact); -protected: - SingleLink< - MatrixLinearSolver, - sofa::component::linearsystem::TypedMatrixLinearSystem, - BaseLink::FLAG_STOREPATH|BaseLink::FLAG_STRONGLINK - > l_linearSystem; - + Data d_factorizationInvalidation; }; ////////////////////////////////////////////////////////////// @@ -414,36 +364,9 @@ class MatrixLinearSolver SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::resetSystem(); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::resizeSystem(Size); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::rebuildSystem(SReal massFactor, SReal forceFactor); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::setSystemLHVector(core::MultiVecDerivId v); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -void MatrixLinearSolver::applySystemSolution(); - template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API GraphScatteredVector* MatrixLinearSolver::createPersistentVector(); -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -linearalgebra::BaseMatrix* MatrixLinearSolver::getSystemBaseMatrix(); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -linearalgebra::BaseVector* MatrixLinearSolver::getSystemRHBaseVector(); - -template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API -linearalgebra::BaseVector* MatrixLinearSolver::getSystemLHBaseVector(); - template<> SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void MatrixLinearSolver::applyConstraintForce(const sofa::core::ConstraintParams* /*cparams*/, sofa::core::MultiVecDerivId /*dx*/, const linearalgebra::BaseVector* /*f*/); diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.inl index b931d0d17d6..495dfd83a70 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MatrixLinearSolver.inl @@ -50,10 +50,12 @@ MatrixLinearSolver::MatrixLinearSolver() "where M is the matrix of the linear system and J is any " "matrix with compatible dimensions")) , invertData() - , linearSystem() - , currentMFactor(), currentBFactor(), currentKFactor() , l_linearSystem(initLink("linearSystem", "The linear system to solve")) + , d_factorizationInvalidation(initData(&d_factorizationInvalidation, false, "factorizationInvalidation", "Internal data for the detection of cache invalidation of the matrix factorization")) { + d_factorizationInvalidation.setReadOnly(true); + d_factorizationInvalidation.setDisplayed(false); + this->addUpdateCallback("parallelInverseProduct", {&d_parallelInverseProduct}, [this](const core::DataTracker& tracker) -> sofa::core::objectmodel::ComponentState { @@ -125,7 +127,7 @@ void MatrixLinearSolver::doCheckLinearSystem() } else { - msg_warning() << "A linear system has been found, but not the expected type." + msg_warning() << "A linear system has been found, but not the expected type. " << "Add a linear system with a compatible type to your scene to remove this warning.\n" << "A component of type " << TLinearSystemType::GetClass()->className << " (template " << TLinearSystemType::GetClass()->templateName << ") will be automatically added for you in Node " @@ -179,17 +181,24 @@ void MatrixLinearSolver::doCheckLinearSystem() } } } + + // this serves as an observer on the matrix to detect when the matrix is re-initialized or cleared + if (l_linearSystem) + { + d_factorizationInvalidation.setParent(&l_linearSystem->d_matrixChanged); + } } -template -template -void MatrixLinearSolver::createDefaultLinearSystem() +template +template +void MatrixLinearSolver::createDefaultLinearSystem() { if (auto system = sofa::core::objectmodel::New()) { - this->addSlave(system); - system->setName( this->getContext()->getNameHelper().resolveName(system->getClassName(), {})); + this->getContext()->addObject(system); + system->setName( + this->getContext()->getNameHelper().resolveName(system->getClassName(), {})); system->f_printLog.setValue(this->f_printLog.getValue()); l_linearSystem.set(system); } @@ -217,125 +226,10 @@ MatrixInvertData * MatrixLinearSolver::createInvertData() return nullptr; } -template -void MatrixLinearSolver::resetSystem() -{ - if (!this->frozen) - { - if (auto* systemMatrix = this->getSystemMatrix()) - { - systemMatrix->clear(); - } - linearSystem.needInvert = true; - } - if (auto* rhs = this->getSystemRHVector()) - { - rhs->clear(); - } - if (auto* solution = this->getSystemLHVector()) - { - solution->clear(); - } - linearSystem.solutionVecId = core::MultiVecDerivId::null(); -} - -template -void MatrixLinearSolver::resizeSystem(Size n) -{ - if (this->getSystemRHVector()) - { - this->getSystemRHVector()->resize(n); - } - - if (this->getSystemLHVector()) - { - this->getSystemLHVector()->resize(n); - } - - linearSystem.needInvert = true; -} - -template -void MatrixLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - this->currentMFactor = mparams->mFactor(); - this->currentBFactor = sofa::core::mechanicalparams::bFactor(mparams); - this->currentKFactor = mparams->kFactor(); - - if (!this->frozen) - { - if (l_linearSystem) - { - l_linearSystem->buildSystemMatrix(mparams); - resizeSystem(l_linearSystem->getMatrixSize()[0]); - } - } - -} - -template -void MatrixLinearSolver::rebuildSystem(SReal massFactor, SReal forceFactor) -{ - sofa::core::MechanicalParams mparams; - mparams.setMFactor(this->currentMFactor*massFactor); - mparams.setBFactor(this->currentBFactor*forceFactor); - mparams.setKFactor(this->currentKFactor*forceFactor); - if (!this->frozen) - { - simulation::common::MechanicalOperations mops(&mparams, this->getContext()); - if (!linearSystem.systemMatrix) linearSystem.systemMatrix = createMatrix(); - linearSystem.matrixAccessor.setGlobalMatrix(linearSystem.systemMatrix); - linearSystem.matrixAccessor.clear(); - mops.getMatrixDimension(&(linearSystem.matrixAccessor)); - linearSystem.matrixAccessor.setupMatrices(); - resizeSystem(linearSystem.matrixAccessor.getGlobalDimension()); - linearSystem.systemMatrix->clear(); - mops.addMBK_ToMatrix(&(linearSystem.matrixAccessor), mparams.mFactor(), mparams.bFactor(), mparams.kFactor()); - linearSystem.matrixAccessor.computeGlobalMatrix(); - } - - this->invertSystem(); -} - -template -void MatrixLinearSolver::setSystemRHVector(core::MultiVecDerivId v) -{ - if (l_linearSystem) - { - l_linearSystem->setRHS(v); - } -} - -template -void MatrixLinearSolver::setSystemLHVector(core::MultiVecDerivId v) -{ - linearSystem.solutionVecId = v; - if (l_linearSystem) - { - l_linearSystem->setSystemSolution(v); - } -} - -template -Matrix* MatrixLinearSolver::getSystemMatrix() -{ - return l_linearSystem ? l_linearSystem->getSystemMatrix() : nullptr; -} - -template -linearalgebra::BaseMatrix* MatrixLinearSolver::getSystemBaseMatrix() -{ - if (!l_linearSystem) - { - return nullptr; - } - return l_linearSystem->getSystemMatrix(); -} - template void MatrixLinearSolver::solveSystem() { - auto* systemMatrix = getSystemMatrix(); + auto* systemMatrix = l_linearSystem->getSystemMatrix(); if (!systemMatrix) { msg_error() << "System matrix is not setup properly"; @@ -343,29 +237,10 @@ void MatrixLinearSolver::solveSystem() } // Step 1: Invert the system, e.g. factorization of the matrix - if (linearSystem.needInvert) - { - this->invert(*systemMatrix); - linearSystem.needInvert = false; - } + this->invertIfInvalidated(*systemMatrix); // Step 2: Solve the system based on the system inversion - this->solve(*systemMatrix, *this->getSystemLHVector(), *this->getSystemRHVector()); - - // Step 3: Apply the solution - applySystemSolution(); -} - -template -void MatrixLinearSolver::applySystemSolution() -{ - if (!linearSystem.solutionVecId.isNull()) - { - if (l_linearSystem) - { - l_linearSystem->dispatchSystemSolution(linearSystem.solutionVecId); - } - } + this->solve(*systemMatrix, *l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector()); } template @@ -392,17 +267,27 @@ void MatrixLinearSolver::deleteMatrix(Matrix* v) delete v; } -template -void MatrixLinearSolver::invertSystem() +template +void MatrixLinearSolver::invertSystem() { - if (linearSystem.needInvert && l_linearSystem) + if (l_linearSystem) { - this->invert(*l_linearSystem->getSystemMatrix()); - linearSystem.needInvert = false; + this->invertIfInvalidated(*l_linearSystem->getSystemMatrix()); } } -template +template +void MatrixLinearSolver::invertIfInvalidated(Matrix& M) +{ + // since this Data is linked to the linear system, the linear system may have modified this value + if (d_factorizationInvalidation.getValue()) + { + invert(M); + d_factorizationInvalidation.setValue(false); + } +} + +template bool MatrixLinearSolver::addJMInvJtLocal(Matrix* M, ResMatrixType* result, const JMatrixType* J, const SReal fact) { if (!this->isComponentStateValid()) @@ -417,18 +302,14 @@ bool MatrixLinearSolver::addJMInvJtLocal(Matrix* M, ResMatrixTyp static_assert(std::is_same_v>, "This function supposes a SparseMatrix"); - auto* systemMatrix = getSystemMatrix(); + auto* systemMatrix = l_linearSystem->getSystemMatrix(); if (!systemMatrix) { msg_error() << "System matrix is not setup properly"; return false; } - if (linearSystem.needInvert) - { - this->invert(*systemMatrix); - linearSystem.needInvert = false; - } + this->invertIfInvalidated(*systemMatrix); simulation::TaskScheduler* taskScheduler = simulation::MainTaskSchedulerFactory::createInRegistry(); assert(taskScheduler); @@ -486,29 +367,28 @@ bool MatrixLinearSolver::singleThreadAddJMInvJtLocal(Matrix* M, SOFA_UNUSED(M); static_assert(std::is_same_v>, "This function supposes a SparseMatrix"); - auto* systemMatrix = getSystemMatrix(); + auto* systemMatrix = l_linearSystem->getSystemMatrix(); if (!systemMatrix) { msg_error() << "System matrix is not setup properly"; return false; } - if (linearSystem.needInvert) - { - this->invert(*systemMatrix); - linearSystem.needInvert = false; - } + auto* rhsVector = l_linearSystem->getRHSVector(); + auto* lhsVector = l_linearSystem->getSolutionVector(); + + this->invertIfInvalidated(*systemMatrix); for (typename JMatrixType::Index row = 0; row < J->rowSize(); ++row) { // STEP 1 : put each line of matrix Jt in the right hand term of the system for (typename JMatrixType::Index i = 0; i < J->colSize(); ++i) { - this->getSystemRHVector()->set(i, J->element(row, i)); // linearSystem.systemMatrix->rowSize() + rhsVector->set(i, J->element(row, i)); // linearSystem.systemMatrix->rowSize() } // STEP 2 : solve the system : - this->solve(*systemMatrix, *this->getSystemLHVector(), *this->getSystemRHVector()); + this->solve(*systemMatrix, *lhsVector, *rhsVector); // STEP 3 : project the result using matrix J for (const auto& [row2, line] : *J) @@ -516,7 +396,7 @@ bool MatrixLinearSolver::singleThreadAddJMInvJtLocal(Matrix* M, Real acc = 0; for (const auto& [col2, val2] : line) { - acc += val2 * getSystemLHVector()->element(col2); + acc += val2 * lhsVector->element(col2); } acc *= fact; result->add(row2, row, acc); @@ -529,12 +409,13 @@ bool MatrixLinearSolver::singleThreadAddJMInvJtLocal(Matrix* M, template bool MatrixLinearSolver::addMInvJtLocal(Matrix * /*M*/,ResMatrixType * result,const JMatrixType * J, SReal fact) { + auto* rhsVector = l_linearSystem->getRHSVector(); for (typename JMatrixType::Index row = 0; row < J->rowSize(); ++row) { // STEP 1 : put each line of matrix Jt in the right hand term of the system for (typename JMatrixType::Index i = 0; i < J->colSize(); ++i) { - getSystemRHVector()->set(i, J->element(row, i)); // linearSystem.systemMatrix->rowSize() + rhsVector->set(i, J->element(row, i)); // linearSystem.systemMatrix->rowSize() } // STEP 2 : solve the system : @@ -543,7 +424,7 @@ bool MatrixLinearSolver::addMInvJtLocal(Matrix * /*M*/,ResMatrixT // STEP 3 : project the result using matrix J for (typename JMatrixType::Index i = 0; i < J->colSize(); ++i) { - result->add(row, i, getSystemRHVector()->element(i) * fact); + result->add(row, i, rhsVector->element(i) * fact); } } @@ -560,7 +441,7 @@ bool MatrixLinearSolver::addJMInvJt(linearalgebra::BaseMatrix* re const JMatrixType * j_local = internalData.getLocalJ(J); ResMatrixType * res_local = internalData.getLocalRes(result); - const bool res = addJMInvJtLocal(getSystemMatrix(), res_local, j_local, fact); + const bool res = addJMInvJtLocal(l_linearSystem->getSystemMatrix(), res_local, j_local, fact); internalData.addLocalRes(result); return res; } @@ -572,7 +453,7 @@ bool MatrixLinearSolver::addMInvJt(linearalgebra::BaseMatrix* res const JMatrixType * j_local = internalData.getLocalJ(J); ResMatrixType * res_local = internalData.getLocalRes(result); - const bool res = addMInvJtLocal(getSystemMatrix(), res_local, j_local, fact); + const bool res = addMInvJtLocal(l_linearSystem->getSystemMatrix(), res_local, j_local, fact); internalData.addLocalRes(result); return res; } @@ -582,7 +463,7 @@ bool MatrixLinearSolver::buildComplianceMatrix(const sofa::core:: { JMatrixType * j_local = internalData.getLocalJ(); j_local->clear(); - j_local->resize(result->rowSize(), getSystemMatrix()->colSize()); + j_local->resize(result->rowSize(), l_linearSystem->getSystemMatrix()->colSize()); if (result->rowSize() == 0) { @@ -610,12 +491,15 @@ bool MatrixLinearSolver::buildComplianceMatrix(const sofa::core:: template void MatrixLinearSolver::applyConstraintForce(const sofa::core::ConstraintParams* cparams, sofa::core::MultiVecDerivId dx, const linearalgebra::BaseVector* f) { - getSystemRHVector()->clear(); - getSystemRHVector()->resize(getSystemMatrix()->colSize()); + auto* systemMatrix = l_linearSystem->getSystemMatrix(); + auto* rhsVector = l_linearSystem->getRHSVector(); + + rhsVector->clear(); + rhsVector->resize(systemMatrix->colSize()); /// rhs = J^t * f - internalData.projectForceInConstraintSpace(getSystemRHVector(), f); + internalData.projectForceInConstraintSpace(rhsVector, f); /// lhs = M^-1 * rhs - this->solve(*getSystemMatrix(), *getSystemLHVector(), *getSystemRHVector()); + this->solve(*systemMatrix, *rhsVector, *rhsVector); l_linearSystem->dispatchSystemSolution(dx); l_linearSystem->dispatchSystemRHS(cparams->lambda()); @@ -624,17 +508,18 @@ void MatrixLinearSolver::applyConstraintForce(const sofa::core::C template void MatrixLinearSolver::computeResidual(const core::ExecParams* params,linearalgebra::BaseVector* f) { - getSystemRHVector()->clear(); - getSystemRHVector()->resize(getSystemMatrix()->colSize()); + auto* rhsVector = l_linearSystem->getRHSVector(); + rhsVector->clear(); + rhsVector->resize(l_linearSystem->getSystemBaseMatrix()->colSize()); /// rhs = J^t * f - internalData.projectForceInConstraintSpace(getSystemRHVector(), f); + internalData.projectForceInConstraintSpace(rhsVector, f); sofa::simulation::common::VectorOperations vop( params, this->getContext() ); sofa::core::behavior::MultiVecDeriv force(&vop, core::vec_id::write_access::force ); // force += rhs - executeVisitor( MechanicalMultiVectorPeqBaseVectorVisitor(core::execparams::defaultInstance(), force, getSystemRHVector(), &(linearSystem.matrixAccessor)) ); + executeVisitor( MechanicalMultiVectorPeqBaseVectorVisitor(core::execparams::defaultInstance(), force, rhsVector, nullptr) ); } diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.h index 2f072621a82..ae8dc6032cc 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.h @@ -52,13 +52,10 @@ class MinResLinearSolver : public sofa::component::linearsolver::MatrixLinearSol MinResLinearSolver(); public: - void resetSystem() override; - void setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) override; - /// Solve Mx=b void solve (Matrix& M, Vector& x, Vector& b) override; - void parse(core::objectmodel::BaseObjectDescription *arg) override; + void parse(core::objectmodel::BaseObjectDescription* arg) override; }; #if !defined(SOFA_COMPONENT_LINEARSOLVER_MINRESLINEARSOLVER_CPP) diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.inl index d4d87e13f6f..7ba7b3fc744 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/MinResLinearSolver.inl @@ -49,23 +49,6 @@ MinResLinearSolver::MinResLinearSolver() f_tolerance.setRequired(true); } -template -void MinResLinearSolver::resetSystem() -{ - f_graph.beginEdit()->clear(); - f_graph.endEdit(); - - Inherit::resetSystem(); -} - -template -void MinResLinearSolver::setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) -{ - Inherit::setSystemMBKMatrix(mparams); -} - - - /// Solve Ax=b /// code issued (and modified) from tminres (https://code.google.com/p/tminres/) // - Umberto Villa, Emory University - uvilla@emory.edu diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp index 72f4df45bed..f76c68e98f8 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp @@ -30,7 +30,7 @@ namespace sofa::component::linearsolver::iterative void registerPCGLinearSolver(sofa::core::ObjectFactory* factory) { - factory->registerObjects(core::ObjectRegistrationData("Linear system solver using the Shewchuk conjugate gradient iterative algorithm.") + factory->registerObjects(core::ObjectRegistrationData("Linear solver using the preconditioned conjugate gradient iterative algorithm.") .add< PCGLinearSolver >()); } diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h index f54a08d8c00..38ddcd331b3 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h @@ -31,69 +31,74 @@ namespace sofa::component::linearsolver::iterative { -/// Linear system solver using the conjugate gradient iterative algorithm +/// Linear system solver using the preconditioned conjugate gradient iterative algorithm template class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver { + public: - SOFA_CLASS(SOFA_TEMPLATE2(PCGLinearSolver,TMatrix,TVector),SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector)); - typedef TMatrix Matrix; - typedef TVector Vector; - typedef sofa::component::linearsolver::MatrixLinearSolver Inherit; + SOFA_CLASS( + SOFA_TEMPLATE2(PCGLinearSolver,TMatrix,TVector), + SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, TMatrix, TVector)); + + using Matrix = TMatrix; + using Vector = TVector; + using Real = typename Matrix::Real; + using Inherit = sofa::component::linearsolver::MatrixLinearSolver; Data d_maxIter; ///< Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop - Data d_tolerance; ///< Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm) + Data d_tolerance; ///< Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm) Data d_use_precond; ///< Use a preconditioner SingleLink l_preconditioner; ///< Link towards the linear solver used to precondition the conjugate gradient - Data d_update_step; ///< Number of steps before the next refresh of preconditioners - Data d_build_precond; ///< Build the preconditioners, if false build the preconditioner only at the initial step - Data > > d_graph; ///< Graph of residuals at each iteration - -protected: - PCGLinearSolver(); + core::objectmodel::lifecycle::DeprecatedData d_update_step; ///< Number of steps before the next refresh of preconditioners + Data > > d_graph; ///< Graph of residuals at each iteration -public: void solve (Matrix& M, Vector& x, Vector& b) override; void init() override; - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; + void bwdInit() override; -private : +private: unsigned next_refresh_step; bool first; int newton_iter; protected: - /// This method is separated from the rest to be able to use custom/optimized versions depending on the types of vectors. - /// It computes: p = p*beta + r - inline void cgstep_beta(Vector& p, Vector& r, double beta); + PCGLinearSolver(); + + void ensureRequiredLinearSystemType(); + + /// This method is separated from the rest to be able to use custom/optimized versions depending on + /// the types of vectors. It computes: p = p*beta + r + inline void cgstep_beta(Vector& p, Vector& r, Real beta); + /// This method is separated from the rest to be able to use custom/optimized versions depending on the types of vectors. /// It computes: x += p*alpha, r -= q*alpha - inline void cgstep_alpha(Vector& x,Vector& p,double alpha); + inline void cgstep_alpha(Vector& x,Vector& p, Real alpha); void handleEvent(sofa::core::objectmodel::Event* event) override; - + void checkLinearSystem() override; }; template -inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta) +inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, Real beta) { p *= beta; p += r; //z; } template -inline void PCGLinearSolver::cgstep_alpha(Vector& x,Vector& p,double alpha) +inline void PCGLinearSolver::cgstep_alpha(Vector& x,Vector& p, Real alpha) { x.peq(p,alpha); // x = x + alpha p } template<> -inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta); +inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, Real beta); template<> -inline void PCGLinearSolver::cgstep_alpha(Vector& x,Vector& p,double alpha); +inline void PCGLinearSolver::cgstep_alpha(Vector& x,Vector& p, Real alpha); #if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP) template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver; diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl index 4335dc9b165..e87042d3d74 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl @@ -20,15 +20,15 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include #include - +#include #include -#include -#include -#include #include #include - +#include +#include + #include namespace sofa::component::linearsolver::iterative @@ -36,12 +36,11 @@ namespace sofa::component::linearsolver::iterative template PCGLinearSolver::PCGLinearSolver() - : d_maxIter(initData(&d_maxIter, (unsigned)25, "iterations", "Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop") ) + : d_maxIter(initData(&d_maxIter, 25u, "iterations", "Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop") ) , d_tolerance(initData(&d_tolerance, 1e-5, "tolerance", "Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm)") ) , d_use_precond(initData(&d_use_precond, true, "use_precond", "Use a preconditioner") ) , l_preconditioner(initLink("preconditioner", "Link towards the linear solver used to precondition the conjugate gradient")) - , d_update_step(initData(&d_update_step, (unsigned)1, "update_step", "Number of steps before the next refresh of preconditioners") ) - , d_build_precond(initData(&d_build_precond, true, "build_precond", "Build the preconditioners, if false build the preconditioner only at the initial step") ) + , d_update_step(this, "v25.12", "v26.06", "update_step", "Instead, use the Data 'assemblingRate' in the associated PreconditionedMatrixFreeSystem") , d_graph(initData(&d_graph, "graph", "Graph of residuals at each iteration") ) , next_refresh_step(0) , newton_iter(0) @@ -51,30 +50,33 @@ PCGLinearSolver::PCGLinearSolver() this->f_listening.setValue(true); } -template -void PCGLinearSolver::init() +template +void PCGLinearSolver::init() { Inherit1::init(); // Find linear solvers if (l_preconditioner.empty()) { - msg_info() << "Link \"preconditioner\" to the desired linear solver should be set to precondition the conjugate gradient."; + msg_info() << "Link '" << l_preconditioner.getName() << "' to the desired linear solver " + "should be set to precondition the conjugate gradient. Without preconditioner, the " + "solver will act as a regular conjugate gradient solver."; } else { if (l_preconditioner.get() == nullptr) { msg_error() << "No preconditioner found at path: " << l_preconditioner.getLinkedPath(); - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Invalid); return; } else { - if (l_preconditioner.get()->getTemplateName() == "GraphScattered") + if (l_preconditioner->getTemplateName() == "GraphScattered") { - msg_error() << "Can not use the preconditioner " << l_preconditioner.get()->getName() << " because it is templated on GraphScatteredType"; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + msg_error() << "Cannot use the preconditioner " << l_preconditioner->getName() + << " because it is templated on GraphScatteredType"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; } else @@ -84,86 +86,106 @@ void PCGLinearSolver::init() } } + ensureRequiredLinearSystemType(); + if (this->isComponentStateInvalid()) + return; + first = true; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Valid); + } } -template -void PCGLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) +template +void PCGLinearSolver::ensureRequiredLinearSystemType() { - sofa::helper::AdvancedTimer::valSet("PCG::buildMBK", 1); - + if (this->l_linearSystem) { - SCOPED_TIMER("PCG::setSystemMBKMatrix"); - Inherit::setSystemMBKMatrix(mparams); + auto* preconditionedMatrix = + dynamic_cast*>(this->l_linearSystem.get()); + if (!preconditionedMatrix) + { + msg_error() << "This linear solver is designed to work with a " + << PreconditionedMatrixFreeSystem::GetClass()->className + << " linear system, but a " << this->l_linearSystem->getClassName() + << " was found"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } } +} - if (l_preconditioner.get()==nullptr) return; +template +void PCGLinearSolver::bwdInit() +{ + if (this->isComponentStateInvalid()) + return; - if (first) //We initialize all the preconditioners for the first step + //link the linear systems of both the preconditioner and the PCG + if (l_preconditioner && this->l_linearSystem) { - l_preconditioner.get()->setSystemMBKMatrix(mparams); - first = false; - next_refresh_step = 1; - } - else if (d_build_precond.getValue()) - { - sofa::helper::AdvancedTimer::valSet("PCG::PrecondBuildMBK", 1); - SCOPED_TIMER_VARNAME(mbkTimer, "PCG::PrecondSetSystemMBKMatrix"); - - if (d_update_step.getValue() > 0) + if (auto* preconditionerLinearSystem = l_preconditioner->getLinearSystem()) { - if (next_refresh_step >= d_update_step.getValue()) + if (auto* preconditionedMatrix = dynamic_cast*>(this->l_linearSystem.get())) { - l_preconditioner.get()->setSystemMBKMatrix(mparams); - next_refresh_step=1; + msg_info() << "Linking the preconditioner linear system (" << preconditionerLinearSystem->getPathName() << ") to the PCG linear system (" << preconditionedMatrix->getPathName() << ")"; + //this link is essential to ensure that the preconditioner matrix is assembled + preconditionedMatrix->l_preconditionerSystem.set(preconditionerLinearSystem); } else { - next_refresh_step++; + msg_error() << "The preconditioned linear system (" << preconditionedMatrix->getPathName() << ") is not a PreconditionedMatrixFreeSystem"; + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Invalid); + return; } } } - - l_preconditioner.get()->updateSystemMatrix(); } -template<> -inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta) +template <> +inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, Real beta) { p.eq(r,p,beta); // p = p*beta + r } template<> -inline void PCGLinearSolver::cgstep_alpha(Vector& x, Vector& p, double alpha) +inline void PCGLinearSolver::cgstep_alpha(Vector& x, Vector& p, Real alpha) { x.peq(p,alpha); // x = x + alpha p } -template -void PCGLinearSolver::handleEvent(sofa::core::objectmodel::Event* event) { +template +void PCGLinearSolver::handleEvent(sofa::core::objectmodel::Event* event) +{ /// this event shoul be launch before the addKToMatrix if (sofa::simulation::AnimateBeginEvent::checkEventType(event)) { newton_iter = 0; - std::map < std::string, sofa::type::vector >& graph = * d_graph.beginEdit(); + std::map >& graph = *d_graph.beginEdit(); graph.clear(); } } +template +void PCGLinearSolver::checkLinearSystem() +{ + // a PreconditionedMatrixFreeSystem component is created in the absence of a linear system + this->template doCheckLinearSystem >(); +} -template +template void PCGLinearSolver::solve (Matrix& M, Vector& x, Vector& b) { SCOPED_TIMER_VARNAME(solveTimer, "PCGLinearSolver::solve"); - std::map < std::string, sofa::type::vector >& graph = * d_graph.beginEdit(); -// sofa::type::vector& graph_error = graph["Error"]; + std::map < std::string, sofa::type::vector >& graph = * d_graph.beginEdit(); +// sofa::type::vector& graph_error = graph["Error"]; newton_iter++; char name[256]; sprintf(name,"Error %d",newton_iter); - sofa::type::vector& graph_error = graph[std::string(name)]; + sofa::type::vector& graph_error = graph[std::string(name)]; const core::ExecParams* params = core::execparams::defaultInstance(); typename Inherit::TempVectorContainer vtmp(this, params, M, x, b); @@ -173,8 +195,8 @@ void PCGLinearSolver::solve (Matrix& M, Vector& x, Vector& b) const bool apply_precond = l_preconditioner.get()!=nullptr && d_use_precond.getValue(); - const double b_norm = b.dot(b); - const double tol = d_tolerance.getValue() * b_norm; + const Real b_norm = b.dot(b); + const Real tol = d_tolerance.getValue() * b_norm; r = M * x; cgstep_beta(r,b,-1);// r = -1 * r + b = b - (M * x) @@ -182,24 +204,25 @@ void PCGLinearSolver::solve (Matrix& M, Vector& x, Vector& b) if (apply_precond) { SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond"); - l_preconditioner.get()->setSystemLHVector(w); - l_preconditioner.get()->setSystemRHVector(r); - l_preconditioner.get()->solveSystem(); + l_preconditioner->getLinearSystem()->setSystemSolution(w); + l_preconditioner->getLinearSystem()->setRHS(r); + l_preconditioner->solveSystem(); + l_preconditioner->getLinearSystem()->dispatchSystemSolution(w); } else { w = r; } - double r_norm = r.dot(w); + Real r_norm = r.dot(w); graph_error.push_back(r_norm/b_norm); unsigned iter=1; while ((iter <= d_maxIter.getValue()) && (r_norm > tol)) { s = M * w; - const double dtq = w.dot(s); - double alpha = r_norm / dtq; + const Real dtq = w.dot(s); + Real alpha = r_norm / dtq; cgstep_alpha(x,w,alpha);//for(int i=0; i::solve (Matrix& M, Vector& x, Vector& b) if (apply_precond) { SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond"); - l_preconditioner.get()->setSystemLHVector(s); - l_preconditioner.get()->setSystemRHVector(r); - l_preconditioner.get()->solveSystem(); + l_preconditioner->getLinearSystem()->setSystemSolution(s); + l_preconditioner->getLinearSystem()->setRHS(r); + l_preconditioner->solveSystem(); + l_preconditioner->getLinearSystem()->dispatchSystemSolution(s); } else { s = r; } - const double deltaOld = r_norm; + const Real deltaOld = r_norm; r_norm = r.dot(s); graph_error.push_back(r_norm/b_norm); - double beta = r_norm / deltaOld; + Real beta = r_norm / deltaOld; cgstep_beta(w,s,beta);//for (int i=0; i. * + ******************************************************************************* + * Authors: The SOFA Team and external contributors (see Authors.txt) * + * * + * Contact information: contact@sofa-framework.org * + ******************************************************************************/ +#define SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PRECONDITIONEDMATRIXFREESYSTEM_CPP +#include +#include + +namespace sofa::component::linearsolver::iterative +{ + +void registerPreconditionedMatrixFreeSystem(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(core::ObjectRegistrationData("A matrix-free linear system that must be used with a preconditioned matrix-free solver") + .add< PreconditionedMatrixFreeSystem >()); +} + +template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PreconditionedMatrixFreeSystem; + +} diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h new file mode 100644 index 00000000000..cf67bdddd4e --- /dev/null +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h @@ -0,0 +1,68 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +namespace sofa::component::linearsolver::iterative +{ + +/** + * A matrix-free linear system that must be used with a preconditioned matrix-free solver + * + * This component is like a @MatrixFreeSystem (its base class), but also has a link to another + * linear system that assembles a matrix. This other linear system is used by a preconditioner + * in the context of a preconditioned solver. + */ +template +class PreconditionedMatrixFreeSystem + : public sofa::component::linearsystem::MatrixFreeSystem +{ +public: + SOFA_CLASS(SOFA_TEMPLATE2(PreconditionedMatrixFreeSystem, TMatrix, TVector), + SOFA_TEMPLATE2(sofa::component::linearsystem::MatrixFreeSystem, TMatrix, TVector)); + + void init() override; + void reset() override; + void buildSystemMatrix(const core::MechanicalParams* mparams) override; + void resizeSystem(sofa::Size n) override; + void clearSystem() override; + + ///< The matrix system of the preconditioner + SingleLink l_preconditionerSystem; + + Data d_assemblingRate; + + void reinitAssemblyCounter(); + +protected: + PreconditionedMatrixFreeSystem(); + + unsigned int m_assemblyCounter {}; +}; + + +#if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PRECONDITIONEDMATRIXFREESYSTEM_CPP) + extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PreconditionedMatrixFreeSystem; +#endif + +} diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl new file mode 100644 index 00000000000..1e5defbb51a --- /dev/null +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl @@ -0,0 +1,99 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::linearsolver::iterative +{ + +template +PreconditionedMatrixFreeSystem::PreconditionedMatrixFreeSystem() + : l_preconditionerSystem(initLink("preconditionerSystem", "Link toward the linear system of the preconditioner")) + , d_assemblingRate(initData(&d_assemblingRate, 1u, "assemblingRate", + "Rate of update of the preconditioner matrix, in number of time steps or Newton iterations")) +{ + this->addAlias(&d_assemblingRate, "update_step"); +} + +template +void PreconditionedMatrixFreeSystem::init() +{ + linearsystem::MatrixFreeSystem::init(); + + reinitAssemblyCounter(); +} + +template +void PreconditionedMatrixFreeSystem::reinitAssemblyCounter() +{ + m_assemblyCounter = d_assemblingRate.getValue(); // to assemble the first time +} + +template +void PreconditionedMatrixFreeSystem::reset() +{ + reinitAssemblyCounter(); +} + +template +void PreconditionedMatrixFreeSystem::buildSystemMatrix( + const core::MechanicalParams* mparams) +{ + //this component builds its own matrix... + linearsystem::MatrixFreeSystem::buildSystemMatrix(mparams); + + //... and also the one from the preconditioner + if (l_preconditionerSystem) + { + if (++m_assemblyCounter >= d_assemblingRate.getValue()) + { + l_preconditionerSystem->buildSystemMatrix(mparams); + m_assemblyCounter = 0; + } + } +} + +template +void PreconditionedMatrixFreeSystem::resizeSystem(sofa::Size n) +{ + linearsystem::MatrixFreeSystem::resizeSystem(n); + + if (l_preconditionerSystem) + { + l_preconditionerSystem->resizeSystem(n); + m_assemblyCounter = 0; + } +} + +template +void PreconditionedMatrixFreeSystem::clearSystem() +{ + linearsystem::MatrixFreeSystem::clearSystem(); + + if (l_preconditionerSystem) + { + l_preconditionerSystem->clearSystem(); + m_assemblyCounter = 0; + } +} + +} // namespace sofa::component::linearsolver::iterative diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/config.h.in b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/config.h.in index 13e67323c21..7448cce2c43 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/config.h.in +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/config.h.in @@ -37,3 +37,16 @@ namespace sofa::component::linearsolver::iterative constexpr const char* MODULE_VERSION = "@PROJECT_VERSION@"; } // namespace sofa::component::linearsolver::iterative +#ifdef SOFA_BUILD_SOFA_COMPONENT_LINEARSOLVER_ITERATIVE +#define SOFA_ITERATIVE_SOLVER_ATTRIBUTE_DEPRECATED_ASSEMBLY_API() +#else +#define SOFA_ITERATIVE_SOLVER_ATTRIBUTE_DEPRECATED_ASSEMBLY_API() \ + SOFA_ATTRIBUTE_DEPRECATED("v25.12", "v26.06", "The assembly of the linear system is no longer the responsibility of the solver. Instead, a linear system component lives along with the linear solver. This component is in charge of the assembly.") +#endif + +#ifdef SOFA_BUILD_SOFA_COMPONENT_LINEARSOLVER_ITERATIVE +#define SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() +#else +#define SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() \ + SOFA_ATTRIBUTE_DISABLED("v25.12", "v25.12", "The assembly of the linear system is no longer the responsibility of the solver. Instead, a linear system component lives along with the linear solver. This component is in charge of the assembly.") +#endif diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/init.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/init.cpp index f6fa6c78938..461a8f6437c 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/init.cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/init.cpp @@ -34,6 +34,7 @@ namespace sofa::component::linearsolver::iterative extern void registerCGLinearSolver(sofa::core::ObjectFactory* factory); extern void registerMinResLinearSolver(sofa::core::ObjectFactory* factory); extern void registerPCGLinearSolver(sofa::core::ObjectFactory* factory); +extern void registerPreconditionedMatrixFreeSystem(sofa::core::ObjectFactory* factory); extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -63,6 +64,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) linearsystem::registerMatrixFreeSystemGraphScattered(factory); registerMinResLinearSolver(factory); registerPCGLinearSolver(factory); + registerPreconditionedMatrixFreeSystem(factory); } void init() diff --git a/Sofa/Component/LinearSolver/Preconditioner/CMakeLists.txt b/Sofa/Component/LinearSolver/Preconditioner/CMakeLists.txt index 3ffd9601673..1a8cf860609 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/CMakeLists.txt +++ b/Sofa/Component/LinearSolver/Preconditioner/CMakeLists.txt @@ -17,6 +17,7 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERPRECONDITIONER_SOURCE_DIR}/WarpPreconditioner.h ${SOFACOMPONENTLINEARSOLVERPRECONDITIONER_SOURCE_DIR}/WarpPreconditioner.inl ${SOFACOMPONENTLINEARSOLVERPRECONDITIONER_SOURCE_DIR}/RotationMatrixSystem.h + ${SOFACOMPONENTLINEARSOLVERPRECONDITIONER_SOURCE_DIR}/RotationMatrixSystem.inl ${SOFACOMPONENTLINEARSOLVERPRECONDITIONER_SOURCE_DIR}/PrecomputedMatrixSystem.h ) diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.h index f3b21e2d882..1d6f02e61c8 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.h +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.h @@ -46,7 +46,6 @@ class JacobiPreconditioner : public sofa::component::linearsolver::MatrixLinearS protected: JacobiPreconditioner(); public: - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; void solve (Matrix& M, Vector& x, Vector& b) override; void invert(Matrix& M) override; diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.inl b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.inl index 646a2f11b61..912a004396c 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.inl +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/JacobiPreconditioner.inl @@ -39,12 +39,6 @@ JacobiPreconditioner::JacobiPreconditioner() { } -template -void JacobiPreconditioner::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - Inherit::setSystemMBKMatrix(mparams); -} - /// Solve P^-1 Mx= P^-1 b // P[i][j] = M[i][j] ssi i=j //P^-1[i][j] = 1/M[i][j] diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h index 81760f03e5b..74ae12683ac 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h @@ -122,14 +122,11 @@ class PrecomputedWarpPreconditioner : public sofa::component::linearsolver::Matr public: void solve (TMatrix& M, TVector& x, TVector& b) override; void invert(TMatrix& M) override; - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; void draw(const core::visual::VisualParams* vparams) override; void init() override; void loadMatrix(TMatrix& M); - bool hasUpdatedMatrix() override {return false;} - TBaseMatrix * getSystemMatrixInv() { return internalData.MinvPtr; diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl index 985e20688fc..5d7475a50ea 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl @@ -67,31 +67,7 @@ PrecomputedWarpPreconditioner::PrecomputedWarpPreconditioner() template void PrecomputedWarpPreconditioner::checkLinearSystem() { - if (!this->l_linearSystem) - { - auto* matrixLinearSystem=this->getContext()->template get >(); - if(!matrixLinearSystem) - { - this->template createDefaultLinearSystem >(); - } - } -} - -template -void PrecomputedWarpPreconditioner::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - // Update the matrix only the first time - if (first) - { - first = false; - init_mFact = mparams->mFactor(); - init_bFact = sofa::core::mechanicalparams::bFactor(mparams); - init_kFact = mparams->kFactor(); - Inherit::setSystemMBKMatrix(mparams); - loadMatrix(*this->getSystemMatrix()); - } - - this->linearSystem.needInvert = usePrecond; + this->template doCheckLinearSystem>(); } //Solve x = R * M^-1 * R^t * b @@ -337,7 +313,6 @@ void PrecomputedWarpPreconditioner::loadMatrixWithSolver() mparams.setMFactor(init_mFact); mparams.setBFactor(init_bFact); mparams.setKFactor(init_kFact); - linearSolver->setSystemMBKMatrix(&mparams); } helper::WriteAccessor > dataForce = *mstate->write(core::vec_id::write_access::externalForce); @@ -397,13 +372,12 @@ void PrecomputedWarpPreconditioner::loadMatrixWithSolver() if (linearSolver) { - linearSolver->setSystemRHVector(rhId); - linearSolver->setSystemLHVector(lhId); + linearSolver->getLinearSystem()->setRHS(rhId); + linearSolver->getLinearSystem()->setSystemSolution(lhId); linearSolver->solveSystem(); + linearSolver->getLinearSystem()->dispatchSystemSolution(lhId); } - if (linearSolver && pid_j*dof_on_node+d == 0) linearSolver->freezeSystemMatrix(); // do not recompute the matrix for the rest of the precomputation - if(pid_j*dof_on_node+d < 2) { EulerSolver->f_printLog.setValue(false); @@ -429,8 +403,6 @@ void PrecomputedWarpPreconditioner::loadMatrixWithSolver() msg_info() << "Precomputing constraint correction : " << std::fixed << 100.0f << " %" ; /////////////////////////////////////////////////////////////////////////////////////////////// - if (linearSolver) linearSolver->freezeSystemMatrix(); // do not recompute the matrix for the rest of the precomputation - ///////////////////////// RESET PARAMETERS AT THEIR PREVIOUS VALUE ///////////////////////////////// // gravity is reset at its previous value this->getContext()->setGravity(gravity); diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.cpp b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.cpp index ffa1961911f..c99f627fd00 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.cpp +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.cpp @@ -19,7 +19,8 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include +#define SOFA_COMPONENT_LINEARSOLVER_PRECONDITIONER_ROTATIONMATRIXSYSTEM_CPP +#include #include #include #include @@ -33,7 +34,7 @@ template class SOFA_COMPONENT_LINEARSOLVER_PRECONDITIONER_API RotationMatrixSyst void registerRotationMatrixSystem(sofa::core::ObjectFactory* factory) { - factory->registerObjects(core::ObjectRegistrationData("Rotation matrix warpping the main linear system.") + factory->registerObjects(core::ObjectRegistrationData("Rotation matrix warping the main linear system.") .add, FullVector > >()); } diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.h index abac06f291f..6b1805ff953 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.h +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.h @@ -1,6 +1,29 @@ -#pragma once +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once #include #include +#include + namespace sofa::component::linearsolver::preconditioner { @@ -9,5 +32,31 @@ class RotationMatrixSystem : public linearsystem::TypedMatrixLinearSystem d_assemblingRate; + + SingleLink l_mainAssembledSystem; + SingleLink l_rotationFinder; + + void init() override; + void reset() override; + + void buildSystemMatrix(const core::MechanicalParams* mparams) override; + + + +protected: + RotationMatrixSystem(); + + // count the number of steps since the last assembly + unsigned int m_assemblyCounter {}; + + void reinitAssemblyCounter(); + void updateMatrixWithRotations(); }; -} + +#if !defined(SOFA_COMPONENT_LINEARSOLVER_PRECONDITIONER_ROTATIONMATRIXSYSTEM_CPP) +extern template class SOFA_COMPONENT_LINEARSOLVER_PRECONDITIONER_API RotationMatrixSystem< RotationMatrix, FullVector >; +#endif + +} // namespace sofa::component::linearsolver::preconditioner diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.inl b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.inl index dbfcf7235aa..0411aaa8c46 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.inl +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/RotationMatrixSystem.inl @@ -1,5 +1,148 @@ -#pragma once +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once #include +#include namespace sofa::component::linearsolver::preconditioner { + +template +RotationMatrixSystem::RotationMatrixSystem() + : d_assemblingRate(initData(&d_assemblingRate, 1u, "assemblingRate", + "Rate of update of the preconditioner matrix")) + , l_mainAssembledSystem(initLink("mainSystem", "Main assembled linear system that will be warped")) + , l_rotationFinder(initLink("rotationFinder", "Link toward the rotation finder used to compute the rotation matrix")) +{} + +template +void RotationMatrixSystem::init() +{ + linearsystem::TypedMatrixLinearSystem::init(); + + if (!l_rotationFinder) + { + if (auto* rotationFinder = this->getContext()->template get()) + { + l_rotationFinder.set(rotationFinder); + msg_info() << "Rotation finder found: '" << l_rotationFinder->getPathName() << "'"; + } + else + { + msg_error() << "A RotationFinder is required by " << this->getClassName() << " but has not " + "been found. The list of available rotation finders is: " + << core::ObjectFactory::getInstance()->listClassesDerivedFrom(); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + } + + if (l_mainAssembledSystem.empty()) + { + msg_info() << "Link \"" << l_mainAssembledSystem.getName() << "\" to the desired linear system should be set to ensure right behavior." << msgendl + << "First assembled linear system found in current context will be used (if any)."; + + const auto listSystems = this->getContext()->template getObjects(sofa::core::objectmodel::BaseContext::Local); + for (const auto& system : listSystems) + { + if (system->getTemplateName() != "GraphScattered" && system != this) + { + l_mainAssembledSystem.set(system); + break; + } + } + } + + if (l_mainAssembledSystem.get() == nullptr) + { + msg_error() << "No linear system component found at path: " << l_mainAssembledSystem.getLinkedPath() << ", nor in current context: " << this->getContext()->name; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + if (l_mainAssembledSystem->getTemplateName() == "GraphScattered") + { + msg_error() << "Cannot use the solver " << l_mainAssembledSystem->getName() + << " because it is templated on GraphScatteredType"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + msg_info() << "Linear system path used: '" << l_mainAssembledSystem->getPathName() << "'"; + + reinitAssemblyCounter(); + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void RotationMatrixSystem::reset() +{ + reinitAssemblyCounter(); +} + +template +void RotationMatrixSystem::buildSystemMatrix( + const core::MechanicalParams* mparams) +{ + bool isMainSystemAssembled = false; + + if (l_mainAssembledSystem) + { + if (++m_assemblyCounter >= d_assemblingRate.getValue()) + { + l_mainAssembledSystem->buildSystemMatrix(mparams); + m_assemblyCounter = 0; + isMainSystemAssembled = true; + } + } + + if (isMainSystemAssembled) + { + const auto matrixSize = l_mainAssembledSystem->getMatrixSize(); + this->resizeSystem(matrixSize[0]); + this->getSystemMatrix()->setIdentity(); + } + else + { + updateMatrixWithRotations(); + } +} + +template +void RotationMatrixSystem::reinitAssemblyCounter() +{ + m_assemblyCounter = d_assemblingRate.getValue(); // to assemble the first time +} + +template +void RotationMatrixSystem::updateMatrixWithRotations() +{ + if (l_rotationFinder) + { + l_rotationFinder->getRotations(this->getSystemBaseMatrix()); + } } + +} // namespace sofa::component::linearsolver::preconditioner diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h index d1cacfed720..94958d3e855 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h @@ -53,50 +53,38 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol using Index = typename TMatrix::Index; SingleLink l_linearSolver; ///< Link towards the linear solver used to build the warp conditioner - - Data d_useRotationFinder; ///< Which rotation Finder to use - Data d_updateStep; ///< Number of steps before the next refresh of the system matrix in the main solver + core::objectmodel::lifecycle::DeprecatedData d_useRotationFinder; + core::objectmodel::lifecycle::DeprecatedData d_updateStep; protected: WarpPreconditioner(); public: - ~WarpPreconditioner(); void init() override; void bwdInit() override; - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; - void invert(Matrix& M) override; - void solve(Matrix& M, Vector& solution, Vector& rh) override; + void solve(Matrix& R, Vector& solution, Vector& rhs) override; bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; bool addMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; - Index getSystemDimention(const sofa::core::MechanicalParams* mparams); - void computeResidual(const core::ExecParams* params, linearalgebra::BaseVector* /*f*/) override; - void updateSystemMatrix() override; - protected: void checkLinearSystem() override; + void ensureRequiredLinearSystemType(); private : int updateSystemSize,currentSystemSize; - int indexwork; - bool first; unsigned nextRefreshStep {}; - TRotationMatrix * rotationWork[2]; - std::vector rotationFinders; - JMatrixType j_local; }; diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.inl b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.inl index d4f761f15ef..19c2553567c 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.inl +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.inl @@ -44,39 +44,43 @@ namespace sofa::component::linearsolver::preconditioner template WarpPreconditioner::WarpPreconditioner() : l_linearSolver(initLink("linearSolver", "Link towards the linear solver used to build the warp conditioner")) - , d_useRotationFinder(initData(&d_useRotationFinder, (unsigned)0, "useRotationFinder", "Which rotation Finder to use" ) ) - , d_updateStep(initData(&d_updateStep, 1u, "update_step", "Number of steps before the next refresh of the system matrix in the main solver" ) ) + , d_useRotationFinder(this, "v25.12", "v26.06", "useRotationFinder", "This Data has been replaced with the Link 'rotationFinder' in RotationMatrixSystem") + , d_updateStep(this, "v25.12", "v26.06", "update_step", "Instead, use the Data 'assemblingRate' in the RotationMatrixSystem" ) { - rotationWork[0] = nullptr; - rotationWork[1] = nullptr; - - first = true; - indexwork = 0; } -template -WarpPreconditioner::~WarpPreconditioner() +template +void WarpPreconditioner::init() { - if (rotationWork[0]) delete rotationWork[0]; - if (rotationWork[1]) delete rotationWork[1]; + Inherit1::init(); - rotationWork[0] = nullptr; - rotationWork[1] = nullptr; + ensureRequiredLinearSystemType(); } template -void WarpPreconditioner::init() +void WarpPreconditioner::ensureRequiredLinearSystemType() { - Inherit1::init(); - first = true; + if (this->l_linearSystem) + { + auto* preconditionedMatrix = + dynamic_cast*>(this->l_linearSystem.get()); + if (!preconditionedMatrix) + { + msg_error() << "This linear solver is designed to work with a " + << RotationMatrixSystem::GetClass()->className + << " linear system, but a " << this->l_linearSystem->getClassName() + << " was found"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } + } } -template +template void WarpPreconditioner::bwdInit() { if (l_linearSolver.empty()) { - msg_info() << "Link \"linearSolver\" to the desired linear solver should be set to ensure right behavior." << msgendl + msg_info() << "Link \"" << l_linearSolver.getName() << "\" to the desired linear solver should be set to ensure right behavior." << msgendl << "First LinearSolver found in current context will be used."; l_linearSolver.set( this->getContext()->template get(sofa::core::objectmodel::BaseContext::Local) ); } @@ -84,174 +88,80 @@ void WarpPreconditioner::bwdInit() if (l_linearSolver.get() == nullptr) { msg_error() << "No LinearSolver component found at path: " << l_linearSolver.getLinkedPath() << ", nor in current context: " << this->getContext()->name; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; } - else - { - if (l_linearSolver.get()->getTemplateName() == "GraphScattered") - { - msg_error() << "Can not use the solver " << l_linearSolver.get()->getName() << " because it is templated on GraphScatteredType"; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } - else - { - msg_info() << "LinearSolver path used: '" << l_linearSolver.getLinkedPath() << "'"; - } - } - - const sofa::core::objectmodel::BaseContext * c = this->getContext(); - c->get(&rotationFinders, sofa::core::objectmodel::BaseContext::Local); - - std::stringstream tmpStr; - tmpStr << "Found " << rotationFinders.size() << " Rotation finders" << msgendl; - for (unsigned i=0; igetName() << msgendl; - } - msg_info() << tmpStr.str(); - - first = true; - indexwork = 0; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); -} - -template -typename WarpPreconditioner::Index -WarpPreconditioner::getSystemDimention(const sofa::core::MechanicalParams* mparams) { - simulation::common::MechanicalOperations mops(mparams, this->getContext()); - - this->linearSystem.matrixAccessor.setGlobalMatrix(this->getSystemMatrix()); - this->linearSystem.matrixAccessor.clear(); - mops.getMatrixDimension(&(this->linearSystem.matrixAccessor)); - this->linearSystem.matrixAccessor.setupMatrices(); - return this->linearSystem.matrixAccessor.getGlobalDimension(); -} - -template -void WarpPreconditioner::setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) -{ - this->currentMFactor = mparams->mFactor(); - this->currentBFactor = sofa::core::mechanicalparams::bFactor(mparams); - this->currentKFactor = mparams->kFactor(); - if (!this->frozen) - { - simulation::common::MechanicalOperations mops(mparams, this->getContext()); - if (!this->l_linearSystem) - { - msg_error() << "No linear system associated to this component. Cannot assemble the matrix."; - return; - } - } - if (first || ( d_updateStep.getValue() > 0 && nextRefreshStep >= d_updateStep.getValue()) || (d_updateStep.getValue() == 0)) + if (l_linearSolver->getTemplateName() == "GraphScattered") { - l_linearSolver.get()->setSystemMBKMatrix(mparams); - nextRefreshStep = 1; - } - - if (first) - { - updateSystemSize = getSystemDimention(mparams); - this->resizeSystem(updateSystemSize); - - first = false; - - if (!rotationWork[indexwork]) rotationWork[indexwork] = new TRotationMatrix(); - - rotationWork[indexwork]->resize(updateSystemSize,updateSystemSize); - rotationFinders[d_useRotationFinder.getValue()]->getRotations(rotationWork[indexwork]); - - if (l_linearSolver.get()->isAsyncSolver()) indexwork = (indexwork==0) ? 1 : 0; - - if (!rotationWork[indexwork]) rotationWork[indexwork] = new TRotationMatrix(); - - rotationWork[indexwork]->resize(updateSystemSize,updateSystemSize); - rotationFinders[d_useRotationFinder.getValue()]->getRotations(rotationWork[indexwork]); - - this->l_linearSystem->resizeSystem(updateSystemSize); - this->l_linearSystem->getSystemMatrix()->setIdentity(); // identity rotationa after update - + msg_error() << "Cannot use the solver " << l_linearSolver->getName() + << " because it is templated on GraphScatteredType"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; } - else if (l_linearSolver.get()->hasUpdatedMatrix()) - { - updateSystemSize = getSystemDimention(mparams); - this->resizeSystem(updateSystemSize); - if (!rotationWork[indexwork]) rotationWork[indexwork] = new TRotationMatrix(); + msg_info() << "LinearSolver path used: '" << l_linearSolver.getLinkedPath() << "'"; - rotationWork[indexwork]->resize(updateSystemSize,updateSystemSize); - rotationFinders[d_useRotationFinder.getValue()]->getRotations(rotationWork[indexwork]); - - if (l_linearSolver.get()->isAsyncSolver()) indexwork = (indexwork==0) ? 1 : 0; - - this->l_linearSystem->getSystemMatrix()->setIdentity(); // identity rotationa after update - } - else + if (!this->isComponentStateInvalid()) { - currentSystemSize = getSystemDimention(sofa::core::mechanicalparams::defaultInstance()); - - this->l_linearSystem->resizeSystem(currentSystemSize); - this->l_linearSystem->getSystemMatrix()->clear(); - - rotationFinders[d_useRotationFinder.getValue()]->getRotations(this->l_linearSystem->getSystemMatrix()); - - this->l_linearSystem->getSystemMatrix()->opMulTM(this->l_linearSystem->getSystemMatrix(),rotationWork[indexwork]); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); } } template void WarpPreconditioner::invert(Matrix& /*Rcur*/) {} -template -void WarpPreconditioner::updateSystemMatrix() { - ++nextRefreshStep; - l_linearSolver.get()->updateSystemMatrix(); -} - template void WarpPreconditioner::checkLinearSystem() { - if (!this->l_linearSystem) - { - auto* matrixLinearSystem = this->getContext()->template get >(); - if (!matrixLinearSystem) - { - this->template createDefaultLinearSystem >(); - } - } + // a RotationMatrixSystem component is created in the absence of a linear system + this->template doCheckLinearSystem>(); } - -/// Solve the system as constructed using the previous methods template -void WarpPreconditioner::solve(Matrix& Rcur, Vector& solution, Vector& rh) { - Rcur.opMulTV(l_linearSolver.get()->getSystemRHBaseVector(),&rh); - - l_linearSolver.get()->solveSystem(); - - Rcur.opMulV(&solution,l_linearSolver.get()->getSystemLHBaseVector()); +void WarpPreconditioner::solve(Matrix& R, Vector& solution, Vector& rhs) +{ + // The matrix A in l_linearSolver is rotated such as R * A * R^T + // The new linear system to solve is then R * A * R^T * x = b + // This is solved in 3 steps: + + // Step 1: + // R * A * R^T * x = b <=> A * R^T * x = R^T * b + // R^T * b is computed in this step: + R.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rhs); + + // Step 2: + // Let's define y = R^T * x, then the linear system is A * y = R^T * b. + // This step solves A * y = R^T * b using the linear solver where y is the unknown. + l_linearSolver->solveSystem(); + + // Step 3: + // Since y = R^T * x, x is deduced: x = R * y + R.opMulV(&solution, l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector()); } /// Solve the system as constructed using the previous methods template -bool WarpPreconditioner::addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) { +bool WarpPreconditioner::addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) +{ if (J->rowSize()==0 || !l_linearSolver.get()) return true; this->l_linearSystem->getSystemMatrix()->rotateMatrix(&j_local,J); - return l_linearSolver.get()->addJMInvJt(result,&j_local,fact); + return l_linearSolver->addJMInvJt(result,&j_local,fact); } template -bool WarpPreconditioner::addMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) { +bool WarpPreconditioner::addMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) +{ this->l_linearSystem->getSystemMatrix()->rotateMatrix(&j_local,J); - return l_linearSolver.get()->addMInvJt(result,&j_local,fact); + return l_linearSolver->addMInvJt(result,&j_local,fact); } template -void WarpPreconditioner::computeResidual(const core::ExecParams* params, linearalgebra::BaseVector* f) { - l_linearSolver.get()->computeResidual(params,f); +void WarpPreconditioner::computeResidual(const core::ExecParams* params, linearalgebra::BaseVector* f) +{ + l_linearSolver->computeResidual(params,f); } } // namespace sofa::component::linearsolver::preconditioner diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h index a0821e3a770..71f14b14ea7 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h @@ -28,9 +28,9 @@ namespace sofa::component::linearsystem { /** - * Component acting like a linear system, but delegates the linear system functionalities to a list of real linear systems. + * Component acting like a linear system, but delegating the linear system functionalities to a list of real linear systems. * - * Using this component allows to assemble more than one global matrix. It is useful when only a partial assembly is + * Using this component allows assembling more than one global matrix. It is useful when only a partial assembly is * necessary. For example, the first linear system is the global one, and the second one could be only the stiffness * matrix. */ @@ -50,6 +50,7 @@ class CompositeLinearSystem : public TypedMatrixLinearSystem TVector* getRHSVector() const override; TVector* getSolutionVector() const override; [[nodiscard]] sofa::linearalgebra::BaseMatrix* getSystemBaseMatrix() const override; + void buildSystemMatrix(const core::MechanicalParams* mparams) override; void resizeSystem(sofa::Size n) override; void clearSystem() override; void setRHS(core::MultiVecDerivId v) override; @@ -58,15 +59,8 @@ class CompositeLinearSystem : public TypedMatrixLinearSystem void dispatchSystemRHS(core::MultiVecDerivId v) override; protected: - void allocateSystem() override; - void resizeVectors(sofa::Size n) override; - - void preAssembleSystem(const core::MechanicalParams* /*mparams*/) override; - void assembleSystem(const core::MechanicalParams* /*mparams*/) override; - void postAssembleSystem(const core::MechanicalParams* /*mparams*/) override; - ///< List of linear systems to assemble - MultiLink < MyType, TypedMatrixLinearSystem, BaseLink::FLAG_DUPLICATE > l_linearSystems; + MultiLink < MyType, sofa::core::behavior::BaseMatrixLinearSystem, BaseLink::FLAG_DUPLICATE > l_linearSystems; ///< Among the list of linear systems, which one is to be used by the linear solver SingleLink < MyType, TypedMatrixLinearSystem, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK > l_solverLinearSystem; diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.inl index 245ac5e3af7..21bf157fee0 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.inl +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.inl @@ -55,25 +55,31 @@ void CompositeLinearSystem::init() else { if (l_solverLinearSystem) + { + if (std::none_of(l_linearSystems.begin(), l_linearSystems.end(), [&](auto linearSystem){return linearSystem == l_solverLinearSystem;})) + { + l_linearSystems.add(l_solverLinearSystem); + } + } + else { bool found = false; - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) + for (const auto& linearSystem : l_linearSystems) { - if (l_solverLinearSystem == l_linearSystems[i]) + if (auto* typedLinearSystem = dynamic_cast*>(linearSystem)) { + l_solverLinearSystem.set(typedLinearSystem); found = true; + break; } } if (!found) { - l_linearSystems.add(l_solverLinearSystem); + msg_error() << "At least one linear system with a compatible type must be provided"; + sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); } } - else - { - l_solverLinearSystem.set(l_linearSystems[0]); - } } } @@ -101,6 +107,19 @@ linearalgebra::BaseMatrix* CompositeLinearSystem::getSystemBas return l_solverLinearSystem ? l_solverLinearSystem->getSystemBaseMatrix() : nullptr; } +template +void CompositeLinearSystem::buildSystemMatrix( + const core::MechanicalParams* mparams) +{ + for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) + { + if (l_linearSystems[i]) + { + l_linearSystems[i]->buildSystemMatrix(mparams); + } + } +} + template void CompositeLinearSystem::resizeSystem(sofa::Size n) { @@ -167,64 +186,4 @@ void CompositeLinearSystem::dispatchSystemRHS(core::MultiVecDe } } -template -void CompositeLinearSystem::allocateSystem() -{ - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) - { - if (l_linearSystems[i]) - { - l_linearSystems[i]->allocateSystem(); - } - } -} - -template -void CompositeLinearSystem::resizeVectors(sofa::Size n) -{ - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) - { - if (l_linearSystems[i]) - { - l_linearSystems[i]->resizeVectors(n); - } - } -} - -template -void CompositeLinearSystem::preAssembleSystem(const core::MechanicalParams* mechanical_params) -{ - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) - { - if (l_linearSystems[i]) - { - l_linearSystems[i]->preAssembleSystem(mechanical_params); - } - } -} - -template -void CompositeLinearSystem::assembleSystem(const core::MechanicalParams* mechanical_params) -{ - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) - { - if (l_linearSystems[i]) - { - l_linearSystems[i]->assembleSystem(mechanical_params); - } - } -} - -template -void CompositeLinearSystem::postAssembleSystem(const core::MechanicalParams* mechanical_params) -{ - for (unsigned int i = 0 ; i < l_linearSystems.size(); ++i) - { - if (l_linearSystems[i]) - { - l_linearSystems[i]->postAssembleSystem(mechanical_params); - } - } -} - } diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixFreeSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixFreeSystem.h index 2ac3f125c27..3c907cfbff1 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixFreeSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixFreeSystem.h @@ -33,7 +33,10 @@ class SOFA_COMPONENT_LINEARSYSTEM_API MatrixFreeSystem : public TypedMatrixLinea public: SOFA_CLASS(SOFA_TEMPLATE2(MatrixFreeSystem, TMatrix, TVector), SOFA_TEMPLATE2(TypedMatrixLinearSystem, TMatrix, TVector)); - void assembleSystem(const core::MechanicalParams* /* mparams */) override {} + void assembleSystem(const core::MechanicalParams* mparams) override + { + this->getSystemMatrix()->setMBKFacts(mparams); + } void associateLocalMatrixToComponents(const core::MechanicalParams* /* mparams */) override {} }; diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.cpp b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.cpp index d36535b0ac1..be6131f7605 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.cpp +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.cpp @@ -59,7 +59,11 @@ void registerMatrixLinearSystem(sofa::core::ObjectFactory* factory) .add >, FullVector > >() .add >, FullVector > >() .add >, FullVector > >() - .add >, FullVector > >()); + .add >, FullVector > >() + .add, FullVector > >() + .add, FullVector > >() + .add, FullVector > >() + ); } } //namespace sofa::component::linearsystem diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h index 509bb5a5a90..87e20f62e1d 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h @@ -71,28 +71,38 @@ class TypedMatrixLinearSystem : public sofa::core::behavior::BaseMatrixLinearSys virtual TVector* getSolutionVector() const; linearalgebra::BaseMatrix* getSystemBaseMatrix() const override; + linearalgebra::BaseVector* getSystemRHSBaseVector() const override; + linearalgebra::BaseVector* getSystemSolutionBaseVector() const override; /// Set the size of the matrix to n x n, and the size of RHS and solution to n - virtual void resizeSystem(sofa::Size n); + void resizeSystem(sofa::Size n) override; - virtual void clearSystem(); + void clearSystem() override; /// Assemble the right-hand side of the linear system, from the values contained in the (Mechanical/Physical)State objects /// Warning: it assumes m_mappingGraph is already built - virtual void setRHS(core::MultiVecDerivId v); + void setRHS(core::MultiVecDerivId v) override; /// Set the initial estimate of the linear system solution vector, from the values contained in the (Mechanical/Physical)State objects /// This vector will be replaced by the solution of the system once the system is solved /// Warning: it assumes m_mappingGraph is already built - virtual void setSystemSolution(core::MultiVecDerivId v); + void setSystemSolution(core::MultiVecDerivId v) override; - virtual void dispatchSystemSolution(core::MultiVecDerivId v); - virtual void dispatchSystemRHS(core::MultiVecDerivId v); + void dispatchSystemSolution(core::MultiVecDerivId v) override; + void dispatchSystemRHS(core::MultiVecDerivId v) override; core::objectmodel::BaseContext* getSolveContext(); + /** + * This Data is used only to notify other components that the system matrix changed (resize, + * clear) + */ + Data d_matrixChanged; + protected: + TypedMatrixLinearSystem(); + LinearSystemData m_linearSystem; /// Relationships between the mechanical states and their associated components @@ -121,20 +131,20 @@ class TypedMatrixLinearSystem : public sofa::core::behavior::BaseMatrixLinearSys /** * Returns the factor to apply to the contributions depending on the contribution type */ - template + template static SReal getContributionFactor( const sofa::core::MechanicalParams* mparams, const sofa::core::matrixaccumulator::get_component_type* object) { if constexpr (c == core::matrixaccumulator::Contribution::STIFFNESS) { - return sofa::core::mechanicalparams::kFactorIncludingRayleighDamping(mparams, - object->rayleighStiffness.getValue()); + return sofa::core::mechanicalparams::kFactorIncludingRayleighDamping( + mparams, object->rayleighStiffness.getValue()); } else if constexpr (c == core::matrixaccumulator::Contribution::MASS) { - return sofa::core::mechanicalparams::mFactorIncludingRayleighDamping(mparams, - object->rayleighMass.getValue()); + return sofa::core::mechanicalparams::mFactorIncludingRayleighDamping( + mparams, object->rayleighMass.getValue()); } else if constexpr (c == core::matrixaccumulator::Contribution::DAMPING) { diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl index 58fcf531f8a..21e5cef1dfa 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl @@ -34,6 +34,14 @@ namespace sofa::component::linearsystem { +template +TypedMatrixLinearSystem::TypedMatrixLinearSystem() + : d_matrixChanged(initData(&d_matrixChanged, false, "factorizationInvalidation", "Internal Data indicating a change in the matrix")) +{ + d_matrixChanged.setReadOnly(true); + d_matrixChanged.setDisplayed(false); +} + template void TypedMatrixLinearSystem::preAssembleSystem(const core::MechanicalParams* mparams) { @@ -69,6 +77,8 @@ void TypedMatrixLinearSystem::preAssembleSystem(const core::Me } associateLocalMatrixToComponents(mparams); + + d_matrixChanged.setValue(true); } template @@ -136,17 +146,43 @@ linearalgebra::BaseMatrix* TypedMatrixLinearSystem::getSystemB return nullptr; } } +template +linearalgebra::BaseVector* TypedMatrixLinearSystem::getSystemRHSBaseVector() const +{ + if constexpr (std::is_base_of_v) + { + return getRHSVector(); + } + else + { + return nullptr; + } +} +template +linearalgebra::BaseVector* TypedMatrixLinearSystem::getSystemSolutionBaseVector() const +{ + if constexpr (std::is_base_of_v) + { + return getSolutionVector(); + } + else + { + return nullptr; + } +} template void TypedMatrixLinearSystem::resizeSystem(sofa::Size n) { m_linearSystem.resizeSystem(n); + d_matrixChanged.setValue(true); } template void TypedMatrixLinearSystem::clearSystem() { m_linearSystem.clearSystem(); + d_matrixChanged.setValue(true); } template @@ -197,12 +233,14 @@ void TypedMatrixLinearSystem::dispatchSystemRHS(core::MultiVec template core::objectmodel::BaseContext* TypedMatrixLinearSystem::getSolveContext() { - auto* linearSolver = this->getContext()->template get(core::objectmodel::BaseContext::Local); + auto* linearSolver = this->getContext()->template get( + core::objectmodel::BaseContext::Local); if (linearSolver) { return linearSolver->getContext(); } - linearSolver = this->getContext()->template get(core::objectmodel::BaseContext::SearchUp); + linearSolver = this->getContext()->template get( + core::objectmodel::BaseContext::SearchUp); if (linearSolver) { return linearSolver->getContext(); @@ -211,4 +249,4 @@ core::objectmodel::BaseContext* TypedMatrixLinearSystem::getSo return this->getContext(); } -} +} // namespace sofa::component::linearsystem diff --git a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/BaseLinearMultiStepMethod.cpp b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/BaseLinearMultiStepMethod.cpp index a437453d8d0..5b1d7acf378 100644 --- a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/BaseLinearMultiStepMethod.cpp +++ b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/BaseLinearMultiStepMethod.cpp @@ -194,9 +194,10 @@ struct ResidualFunction : newton_raphson::BaseNonLinearFunction { vop.v_clear(dv); - linearSolver->setSystemLHVector(dv); - linearSolver->setSystemRHVector(rhs); + linearSolver->getLinearSystem()->setSystemSolution(dv); + linearSolver->getLinearSystem()->setRHS(rhs); linearSolver->solveSystem(); + linearSolver->getLinearSystem()->dispatchSystemSolution(dv); } void computeDxFromDv() diff --git a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/EulerImplicitSolver.cpp b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/EulerImplicitSolver.cpp index 9d829a9020f..212e2dd7fb8 100644 --- a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/EulerImplicitSolver.cpp +++ b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/EulerImplicitSolver.cpp @@ -178,9 +178,10 @@ void EulerImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa:: { SCOPED_TIMER("MBKSolve"); - l_linearSolver->setSystemLHVector(x); - l_linearSolver->setSystemRHVector(b); + l_linearSolver->getLinearSystem()->setSystemSolution(x); + l_linearSolver->getLinearSystem()->setRHS(b); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(x); } #ifdef SOFA_DUMP_VISITOR_INFO simulation::Visitor::printCloseNode("SystemSolution"); diff --git a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/NewmarkImplicitSolver.cpp b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/NewmarkImplicitSolver.cpp index e69abebe1e6..01b572682b2 100644 --- a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/NewmarkImplicitSolver.cpp +++ b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/NewmarkImplicitSolver.cpp @@ -135,9 +135,10 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa const core::MatricesFactors::K kFact ( -h * h * beta - h * rK * gamma); mop.setSystemMBKMatrix(mFact, bFact, kFact, l_linearSolver.get()); - l_linearSolver->setSystemLHVector(aResult); - l_linearSolver->setSystemRHVector(b); + l_linearSolver->getLinearSystem()->setSystemSolution(aResult); + l_linearSolver->getLinearSystem()->setRHS(b); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(aResult); msg_info() << "a1 = " << aResult; diff --git a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/StaticSolver.cpp b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/StaticSolver.cpp index 85de14bca27..8b279100cb9 100644 --- a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/StaticSolver.cpp +++ b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/StaticSolver.cpp @@ -150,9 +150,10 @@ struct StaticResidualFunction : newton_raphson::BaseNonLinearFunction { SCOPED_TIMER("MBKSolve"); - linearSolver->setSystemLHVector(dx); - linearSolver->setSystemRHVector(force); + linearSolver->getLinearSystem()->setSystemSolution(dx); + linearSolver->getLinearSystem()->setRHS(force); linearSolver->solveSystem(); + linearSolver->getLinearSystem()->dispatchSystemSolution(dx); } SReal squaredNormDx() override diff --git a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/VariationalSymplecticSolver.cpp b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/VariationalSymplecticSolver.cpp index 2211fba1849..bb7cfc7c926 100644 --- a/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/VariationalSymplecticSolver.cpp +++ b/Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/VariationalSymplecticSolver.cpp @@ -224,9 +224,10 @@ void VariationalSymplecticSolver::solve(const core::ExecParams* params, SReal dt { SCOPED_TIMER("MBKSolve"); // resolution of matrix*res=b - l_linearSolver->setSystemLHVector(res); - l_linearSolver->setSystemRHVector(b); + l_linearSolver->getLinearSystem()->setSystemSolution(res); + l_linearSolver->getLinearSystem()->setRHS(b); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(res); } /// Updates of q(k,i) /// @@ -307,9 +308,10 @@ void VariationalSymplecticSolver::solve(const core::ExecParams* params, SReal dt mop.setSystemMBKMatrix(core::MatricesFactors::M(1), core::MatricesFactors::B(0), core::MatricesFactors::K(0), l_linearSolver.get()); // resolution of matrix*b=newp - l_linearSolver->setSystemLHVector(b); - l_linearSolver->setSystemRHVector(newp); + l_linearSolver->getLinearSystem()->setSystemSolution(b); + l_linearSolver->getLinearSystem()->setRHS(newp); l_linearSolver->solveSystem(); // b = inv(matrix)*newp = Minv*newp + l_linearSolver->getLinearSystem()->dispatchSystemSolution(b); const auto hamiltonianKineticEnergy = 0.5*(newp.dot(b)); diff --git a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp index ebf86e5ed79..9f33d9df4f0 100644 --- a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp +++ b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp @@ -337,9 +337,10 @@ void EulerExplicitSolver::assembleSystemMatrix(sofa::simulation::common::Mechani void EulerExplicitSolver::solveSystem(core::MultiVecDerivId solution, core::MultiVecDerivId rhs) const { SCOPED_TIMER("MBKSolve"); - l_linearSolver->setSystemLHVector(solution); - l_linearSolver->setSystemRHVector(rhs); + l_linearSolver->getLinearSystem()->setSystemSolution(solution); + l_linearSolver->getLinearSystem()->setRHS(rhs); l_linearSolver->solveSystem(); + l_linearSolver->getLinearSystem()->dispatchSystemSolution(solution); } } // namespace sofa::component::odesolver::forward diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedronFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedronFEMForceField.h index 9bed7b8db99..81171ce7709 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedronFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedronFEMForceField.h @@ -79,10 +79,15 @@ class HexahedronFEMForceFieldInternalData * 0---------1-->X */ template -class HexahedronFEMForceField : virtual public BaseLinearElasticityFEMForceField, public sofa::core::behavior::BaseRotationFinder +class HexahedronFEMForceField : + virtual public BaseLinearElasticityFEMForceField, + public sofa::core::behavior::BaseRotationFinder { public: - SOFA_CLASS(SOFA_TEMPLATE(HexahedronFEMForceField, DataTypes), SOFA_TEMPLATE(BaseLinearElasticityFEMForceField, DataTypes)); + SOFA_CLASS2( + SOFA_TEMPLATE(HexahedronFEMForceField, DataTypes), + SOFA_TEMPLATE(BaseLinearElasticityFEMForceField, DataTypes), + sofa::core::behavior::BaseRotationFinder); typedef BaseLinearElasticityFEMForceField InheritForceField; typedef typename DataTypes::VecCoord VecCoord; diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseForceField.cpp b/Sofa/framework/Core/src/sofa/core/behavior/BaseForceField.cpp index 03f74235b4e..d4cc43b94c2 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseForceField.cpp +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseForceField.cpp @@ -38,8 +38,13 @@ BaseForceField::BaseForceField() void BaseForceField::addMBKdx(const MechanicalParams* mparams, MultiVecDerivId dfId) { - if (sofa::core::mechanicalparams::kFactorIncludingRayleighDamping(mparams,rayleighStiffness.getValue()) != 0.0 || sofa::core::mechanicalparams::bFactor(mparams) != 0.0) + const auto kFactor = sofa::core::mechanicalparams::kFactorIncludingRayleighDamping(mparams,rayleighStiffness.getValue()); + const auto bFactor = sofa::core::mechanicalparams::bFactor(mparams); + + if (kFactor != 0.0 || bFactor != 0.0) + { addDForce(mparams, dfId); + } } void BaseForceField::addBToMatrix(const MechanicalParams* /*mparams*/, const sofa::core::behavior::MultiMatrixAccessor* /*matrix*/) diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp index ae24e93ce84..5a341c1277a 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp @@ -27,15 +27,22 @@ namespace sofa::core::behavior BaseMatrixLinearSystem::BaseMatrixLinearSystem() : Inherit1() , d_matrixSize(initData(&d_matrixSize, "matrixSize", "Size of the global matrix")) +, d_enableAssembly(initData(&d_enableAssembly, true, "enableAssembly", "Allows to assemble the system matrix")) { d_matrixSize.setReadOnly(true); + + d_enableAssembly.setReadOnly(true); + d_enableAssembly.setDisplayed(false); } void BaseMatrixLinearSystem::buildSystemMatrix(const core::MechanicalParams* mparams) { - preAssembleSystem(mparams); - assembleSystem(mparams); - postAssembleSystem(mparams); + if (d_enableAssembly.getValue()) + { + preAssembleSystem(mparams); + assembleSystem(mparams); + postAssembleSystem(mparams); + } } void BaseMatrixLinearSystem::preAssembleSystem(const core::MechanicalParams* mparams) diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h index 62d5742ca6c..47462077153 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h @@ -24,6 +24,7 @@ #include #include +#include namespace sofa::core::behavior { @@ -40,19 +41,39 @@ class SOFA_CORE_API BaseMatrixLinearSystem : public virtual core::objectmodel::B protected: BaseMatrixLinearSystem(); +public: + /// Size of the linear system Data< sofa::type::Vec2u > d_matrixSize; -public: + Data< bool > d_enableAssembly; - /// Returns the system matrix as a linearalgebra::BaseMatrix* + /// Returns the system matrix as a sofa::linearalgebra::BaseMatrix* virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() const { return nullptr; } + virtual linearalgebra::BaseVector* getSystemRHSBaseVector() const { return nullptr; } + virtual linearalgebra::BaseVector* getSystemSolutionBaseVector() const { return nullptr; } + /// Construct and assemble the linear system matrix - void buildSystemMatrix(const core::MechanicalParams* mparams); + virtual void buildSystemMatrix(const core::MechanicalParams* mparams); sofa::type::Vec2u getMatrixSize() const { return d_matrixSize.getValue(); } + /// Set the size of the matrix to n x n, and the size of RHS and solution to n + virtual void resizeSystem(sofa::Size n) = 0; + + virtual void clearSystem() = 0; + + /// Assemble the right-hand side of the linear system from the values contained in the (Mechanical/Physical)State objects + virtual void setRHS(core::MultiVecDerivId v) = 0; + + /// Set the initial estimate of the linear system solution vector, from the values contained in the (Mechanical/Physical)State objects + /// This vector will be replaced by the solution of the system once the system is solved + virtual void setSystemSolution(core::MultiVecDerivId v) = 0; + + virtual void dispatchSystemSolution(core::MultiVecDerivId v) = 0; + virtual void dispatchSystemRHS(core::MultiVecDerivId v) = 0; + protected: virtual void preAssembleSystem(const core::MechanicalParams* /*mparams*/); virtual void assembleSystem(const core::MechanicalParams* /*mparams*/); diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseRotationFinder.h b/Sofa/framework/Core/src/sofa/core/behavior/BaseRotationFinder.h index 3b0c62482a0..9f91a70fab8 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseRotationFinder.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseRotationFinder.h @@ -30,6 +30,8 @@ namespace sofa::core::behavior class BaseRotationFinder : public virtual sofa::core::objectmodel::BaseObject { public: + SOFA_ABSTRACT_CLASS(BaseRotationFinder, sofa::core::objectmodel::BaseObject); + virtual void getRotations(linearalgebra::BaseMatrix * m, int offset = 0) = 0; }; diff --git a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.cpp b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.cpp index 45511304aa9..5a47af3681d 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.cpp +++ b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.cpp @@ -24,12 +24,28 @@ namespace sofa::core::behavior { -LinearSolver::LinearSolver() - : frozen(false) -{} +void LinearSolver::init_partial_solve() +{ + msg_warning() << "partial_solve is not implemented yet."; +} + +void LinearSolver::partial_solve(std::list&, + std::list&, bool) +{ + msg_warning() << "partial_solve is not implemented yet"; +} -LinearSolver::~LinearSolver() -{} +void LinearSolver::applyConstraintForce(const sofa::core::ConstraintParams*, + sofa::core::MultiVecDerivId, + const linearalgebra::BaseVector*) +{ + msg_error() << "applyConstraintForce has not been implemented."; +} + +void LinearSolver::computeResidual(const core::ExecParams*, linearalgebra::BaseVector*) +{ + msg_error() << "computeResidual has not been implemented."; +} } // namespace sofa::core::behavior diff --git a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h index a45ca12b695..7bc3a8402c9 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h @@ -21,9 +21,11 @@ ******************************************************************************/ #pragma once +#include +#include #include #include -#include +#include namespace sofa::core::behavior { @@ -35,54 +37,62 @@ namespace sofa::core::behavior class SOFA_CORE_API LinearSolver : public BaseLinearSolver { public: - SOFA_ABSTRACT_CLASS(LinearSolver, BaseLinearSolver); + SOFA_ABSTRACT_CLASS(LinearSolver, BaseLinearSolver) SOFA_BASE_CAST_IMPLEMENTATION(LinearSolver) -protected: - LinearSolver(); - ~LinearSolver() override; -public: /// Reset the current linear system. - virtual void resetSystem() = 0; + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void resetSystem() final = delete; /// Set the linear system matrix, combining the mechanical M,B,K matrices using the given coefficients /// /// @todo Should we put this method in a specialized class for mechanical systems, or express it using more general terms (i.e. coefficients of the second order ODE to solve) - virtual void setSystemMBKMatrix(const MechanicalParams* mparams) = 0; + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void setSystemMBKMatrix(const MechanicalParams* mparams) final = delete; /// Rebuild the system using a mass and force factor /// Experimental API used to investigate convergence issues. - SOFA_ATTRIBUTE_DEPRECATED__REBUILDSYSTEM() virtual void rebuildSystem(SReal /*massFactor*/, SReal /*forceFactor*/){} + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void rebuildSystem(SReal /*massFactor*/, SReal /*forceFactor*/) final = delete; + + virtual sofa::core::behavior::BaseMatrixLinearSystem* getLinearSystem() const = 0; - /// Indicate if the solver update the system in parallel + /// Indicate if the solver updates the system in parallel virtual bool isAsyncSolver() { return false; } /// Returns true if the solver supports non-symmetric systems virtual bool supportNonSymmetricSystem() const { return false; } /// Indicate if the solver updated the system after the last call of setSystemMBKMatrix (should return true if isParallelSolver return false) - virtual bool hasUpdatedMatrix() { return true; } + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual bool hasUpdatedMatrix() final = delete; /// This function is use for the preconditioner it must be called at each time step event if setSystemMBKMatrix is not called - virtual void updateSystemMatrix() {} + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void updateSystemMatrix() final = delete; /// Set the linear system right-hand term vector, from the values contained in the (Mechanical/Physical)State objects - virtual void setSystemRHVector(core::MultiVecDerivId v) = 0; + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void setSystemRHVector(core::MultiVecDerivId v) final = delete; /// Set the initial estimate of the linear system left-hand term vector, from the values contained in the (Mechanical/Physical)State objects /// This vector will be replaced by the solution of the system once solveSystem is called - virtual void setSystemLHVector(core::MultiVecDerivId v) = 0; + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void setSystemLHVector(core::MultiVecDerivId v) final = delete; /// Solve the system as constructed using the previous methods virtual void solveSystem() = 0; /// - virtual void init_partial_solve() { msg_warning() << "partial_solve is not implemented yet."; } + virtual void init_partial_solve(); /// - virtual void partial_solve(std::list& /*I_last_Disp*/, std::list& /*I_last_Dforce*/, bool /*NewIn*/) { msg_warning() << "partial_solve is not implemented yet"; } + virtual void partial_solve(std::list& /*I_last_Disp*/, + std::list& /*I_last_Dforce*/, + bool /*NewIn*/); - /// Invert the system, this method is optional because it's called when solveSystem() is called for the first time + /// Invert the system, this method is optional because it's called when solveSystem() is called + /// for the first time virtual void invertSystem() {} /// Multiply the inverse of the system matrix by the transpose of the given matrix J @@ -118,20 +128,19 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver } /// Apply the contactforce dx = Minv * J^t * f and store the result in dx VecId - virtual void applyConstraintForce(const sofa::core::ConstraintParams* /*cparams*/,sofa::core::MultiVecDerivId /*dx*/, const linearalgebra::BaseVector* /*f*/) { - msg_error() << "applyConstraintForce has not been implemented."; - } + virtual void applyConstraintForce(const sofa::core::ConstraintParams* /*cparams*/, + sofa::core::MultiVecDerivId /*dx*/, + const linearalgebra::BaseVector* /*f*/); /// Compute the residual in the newton iterations due to the constraints forces /// i.e. compute mparams->dF() = J^t lambda /// the result is written in mparams->dF() SOFA_ATTRIBUTE_DEPRECATED__COMPUTERESIDUAL() - virtual void computeResidual(const core::ExecParams* /*params*/, linearalgebra::BaseVector* /*f*/) { - msg_error() << "computeResidual has not been implemented."; - } - + virtual void computeResidual(const core::ExecParams* /*params*/, + linearalgebra::BaseVector* /*f*/); - /// Multiply the inverse of the system matrix by the transpose of the given matrix, and multiply the result with the given matrix J + /// Multiply the inverse of the system matrix by the transpose of the given matrix, and multiply + /// the result with the given matrix J /// /// This method can compute the Schur complement of the constrained system: /// W = H A^{-1} H^T, where: @@ -142,7 +151,8 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver /// @param result the variable where the result will be added /// @param J the matrix J to use /// @param fact integrator parameter - /// @return false if the solver does not support this operation, of it the system matrix is not invertible + /// @return false if the solver does not support this operation, of it the system matrix is not + /// invertible virtual bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) { SOFA_UNUSED(result); @@ -152,16 +162,20 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver } /// Get the linear system matrix, or nullptr if this solver does not build it - virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() { return nullptr; } + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() final = delete; /// Get the linear system right-hand term vector, or nullptr if this solver does not build it - virtual linearalgebra::BaseVector* getSystemRHBaseVector() { return nullptr; } + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual linearalgebra::BaseVector* getSystemRHBaseVector() final = delete; /// Get the linear system left-hand term vector, or nullptr if this solver does not build it - virtual linearalgebra::BaseVector* getSystemLHBaseVector() { return nullptr; } + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual linearalgebra::BaseVector* getSystemLHBaseVector() final = delete; /// Get the linear system inverse matrix, or nullptr if this solver does not build it - virtual linearalgebra::BaseMatrix* getSystemInverseBaseMatrix() { return nullptr; } + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual linearalgebra::BaseMatrix* getSystemInverseBaseMatrix() final = delete; /// Read the Matrix solver from a file virtual bool readFile(std::istream& /*in*/) { return false;} @@ -170,13 +184,13 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver virtual bool writeFile(std::ostream& /*out*/) {return false;} /// Ask the solver to no longer update the system matrix - virtual void freezeSystemMatrix() { frozen = true; } - - + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void freezeSystemMatrix() = delete; protected: - bool frozen; + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved frozen; }; } // namespace sofa::core::behavior diff --git a/Sofa/framework/Core/src/sofa/core/behavior/RotationFinder.h b/Sofa/framework/Core/src/sofa/core/behavior/RotationFinder.h index 3a3db40d205..7adde6580c2 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/RotationFinder.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/RotationFinder.h @@ -33,6 +33,8 @@ template class RotationFinder : public BaseRotationFinder { public: + SOFA_CLASS(SOFA_TEMPLATE(RotationFinder, DataTypes), BaseRotationFinder); + typedef typename DataTypes::VecCoord VecCoord; typedef typename DataTypes::VecDeriv VecDeriv; typedef typename DataTypes::Coord Coord; diff --git a/Sofa/framework/Core/src/sofa/core/config.h.in b/Sofa/framework/Core/src/sofa/core/config.h.in index a728ab0e371..44582e0de78 100644 --- a/Sofa/framework/Core/src/sofa/core/config.h.in +++ b/Sofa/framework/Core/src/sofa/core/config.h.in @@ -161,3 +161,11 @@ SOFA_ATTRIBUTE_DEPRECATED("v25.12", "v26.06", "Use getContactDistance or setCont #define SOFA_ATTRIBUTE_DEPRECATED__COMPUTERESIDUAL() \ SOFA_ATTRIBUTE_DEPRECATED("v25.12", "v26.06", "The method computeResidual is not used") #endif + + +#ifdef SOFA_BUILD_SOFA_CORE +#define SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() +#else +#define SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() \ + SOFA_ATTRIBUTE_DISABLED("v25.12", "v25.12", "The assembly of the linear system is no longer the responsibility of the solver. Instead, a linear system component lives along with the linear solver. This component is in charge of the assembly.") +#endif diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp index e99915bcfc0..d6c0e0b020d 100644 --- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp @@ -400,7 +400,10 @@ void MechanicalOperations::resetSystem(core::behavior::LinearSolver* linearSolve { if (linearSolver) { - linearSolver->resetSystem(); + if (auto* linearSystem = linearSolver->getLinearSystem()) + { + linearSystem->clearSystem(); + } } } @@ -414,25 +417,34 @@ void MechanicalOperations::setSystemMBKMatrix( mparams.setBFactor(b.get()); mparams.setKFactor(k.get()); mparams.setSupportOnlySymmetricMatrix(!linearSolver->supportNonSymmetricSystem()); - linearSolver->setSystemMBKMatrix(&mparams); + if (auto* linearSystem = linearSolver->getLinearSystem()) + { + linearSystem->buildSystemMatrix(&mparams); + } } } -void MechanicalOperations::setSystemRHVector(core::MultiVecDerivId v, - core::behavior::LinearSolver* linearSolver) +void MechanicalOperations::setSystemRHVector( + core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver) { if (linearSolver) { - linearSolver->setSystemRHVector(v); + if (auto* linearSystem = linearSolver->getLinearSystem()) + { + linearSystem->setRHS(v); + } } } -void MechanicalOperations::setSystemLHVector(core::MultiVecDerivId v, - core::behavior::LinearSolver* linearSolver) +void MechanicalOperations::setSystemLHVector( + core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver) { if (linearSolver) { - linearSolver->setSystemLHVector(v); + if (auto* linearSystem = linearSolver->getLinearSystem()) + { + linearSystem->setSystemSolution(v); + } } } @@ -444,12 +456,25 @@ void MechanicalOperations::solveSystem(core::behavior::LinearSolver* linearSolve } } +void MechanicalOperations::solveSystem( + core::behavior::LinearSolver* linearSolver, core::MultiVecDerivId v) +{ + if (linearSolver) + { + linearSolver->solveSystem(); + if (auto* linearSystem = linearSolver->getLinearSystem()) + { + linearSystem->dispatchSystemSolution(v); + } + } +} + void MechanicalOperations::print(std::ostream& out, core::behavior::LinearSolver* linearSolver) { if (linearSolver) { - const linearalgebra::BaseMatrix* m = linearSolver->getSystemBaseMatrix(); + const linearalgebra::BaseMatrix* m = linearSolver->getLinearSystem()->getSystemBaseMatrix(); if (!m) { return; diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h index 1406df545a7..883607d20c2 100644 --- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h @@ -107,6 +107,7 @@ class SOFA_SIMULATION_CORE_API MechanicalOperations void setSystemRHVector(core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver); void setSystemLHVector(core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver); void solveSystem(core::behavior::LinearSolver* linearSolver); + void solveSystem(core::behavior::LinearSolver* linearSolver, core::MultiVecDerivId v); void print( std::ostream& out, core::behavior::LinearSolver* linearSolver); /// @} diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/mechanicalvisitor/MechanicalMultiVectorPeqBaseVectorVisitor.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/mechanicalvisitor/MechanicalMultiVectorPeqBaseVectorVisitor.cpp index 4d697728c91..24fa5be55ff 100644 --- a/Sofa/framework/Simulation/Core/src/sofa/simulation/mechanicalvisitor/MechanicalMultiVectorPeqBaseVectorVisitor.cpp +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/mechanicalvisitor/MechanicalMultiVectorPeqBaseVectorVisitor.cpp @@ -42,7 +42,7 @@ MechanicalMultiVectorPeqBaseVectorVisitor::Result MechanicalMultiVectorPeqBaseVe mm->addFromBaseVectorSameSize(dest.getId(mm), src, o); offset = (int)o; } - //if (!matrix) offset += mm->getMatrixSize(); + if (!matrix) offset += mm->getMatrixSize(); return RESULT_CONTINUE; } diff --git a/examples/Component/LinearSolver/FEMBAR-common.xml b/examples/Component/LinearSolver/FEMBAR-common.xml index 26e3896d17f..5ff5b29ba12 100644 --- a/examples/Component/LinearSolver/FEMBAR-common.xml +++ b/examples/Component/LinearSolver/FEMBAR-common.xml @@ -13,6 +13,7 @@ + diff --git a/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn b/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn index 89b5871bb32..ffd5be46917 100644 --- a/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn +++ b/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn @@ -2,7 +2,7 @@ - + - \ No newline at end of file + diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn index b8a29eb500e..2f7a003c9cf 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn @@ -2,8 +2,10 @@ + - + + diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn index fcb6eea9c34..ce794336778 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn @@ -2,9 +2,19 @@ - - - + + + + + + + + + + + + + diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn index aa4c61c823e..fe9a78ee749 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn @@ -2,9 +2,18 @@ - - - + + + + + + + + + + + + diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn index 5198f1f8e9e..4ce13948640 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn @@ -2,9 +2,14 @@ - - + + + + + + + diff --git a/examples/Component/ODESolver/Backward/BDFOdeSolver.scn b/examples/Component/ODESolver/Backward/BDFOdeSolver.scn index 208f56a87d0..4be6d5e0e46 100644 --- a/examples/Component/ODESolver/Backward/BDFOdeSolver.scn +++ b/examples/Component/ODESolver/Backward/BDFOdeSolver.scn @@ -30,8 +30,8 @@ absoluteEstimateDifferenceThreshold="1e-5" relativeEstimateDifferenceThreshold="1e-5"/> - - + + diff --git a/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn b/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn index 967751f69f0..23cae4c230c 100644 --- a/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn +++ b/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn @@ -39,7 +39,8 @@ - + +