From 63514a4b4e6db478fea1abb1e33e72dad055238a Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 6 Aug 2025 16:28:58 +0200 Subject: [PATCH 01/48] Start removing assembly function from linear solvers --- .../LinearSolverConstraintCorrection.inl | 22 +- .../PrecomputedConstraintCorrection.inl | 7 - .../direct/AsyncSparseLDLSolver.h | 1 - .../direct/AsyncSparseLDLSolver.inl | 36 +--- .../linearsolver/direct/BTDLinearSolver.inl | 24 +-- .../direct/PrecomputedLinearSolver.h | 1 - .../direct/PrecomputedLinearSolver.inl | 15 +- .../linearsolver/iterative/CGLinearSolver.h | 4 - .../linearsolver/iterative/CGLinearSolver.inl | 21 -- .../iterative/MatrixLinearSolver.cpp | 59 ----- .../iterative/MatrixLinearSolver.h | 120 +---------- .../iterative/MatrixLinearSolver.inl | 201 +++--------------- .../iterative/MinResLinearSolver.h | 5 +- .../iterative/MinResLinearSolver.inl | 17 -- .../linearsolver/iterative/PCGLinearSolver.h | 1 - .../iterative/PCGLinearSolver.inl | 54 +---- .../preconditioner/JacobiPreconditioner.h | 1 - .../preconditioner/JacobiPreconditioner.inl | 6 - .../PrecomputedWarpPreconditioner.h | 1 - .../PrecomputedWarpPreconditioner.inl | 27 +-- .../preconditioner/WarpPreconditioner.h | 2 - .../preconditioner/WarpPreconditioner.inl | 113 +++------- .../linearsystem/CompositeLinearSystem.h | 4 +- .../component/linearsystem/MatrixFreeSystem.h | 5 +- .../linearsystem/TypedMatrixLinearSystem.h | 24 ++- .../linearsystem/TypedMatrixLinearSystem.inl | 24 +++ .../backward/BaseLinearMultiStepMethod.cpp | 5 +- .../backward/EulerImplicitSolver.cpp | 5 +- .../backward/NewmarkImplicitSolver.cpp | 5 +- .../odesolver/backward/StaticSolver.cpp | 5 +- .../backward/VariationalSymplecticSolver.cpp | 10 +- .../odesolver/forward/EulerExplicitSolver.cpp | 5 +- .../src/sofa/core/behavior/BaseForceField.cpp | 7 +- .../core/behavior/BaseMatrixLinearSystem.h | 21 +- .../src/sofa/core/behavior/LinearSolver.cpp | 26 ++- .../src/sofa/core/behavior/LinearSolver.h | 75 ++----- .../sofa/simulation/MechanicalOperations.cpp | 43 +++- .../sofa/simulation/MechanicalOperations.h | 1 + 38 files changed, 271 insertions(+), 732 deletions(-) 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..e719d452ecb 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 @@ -63,7 +63,6 @@ class AsyncSparseLDLSolver : public SparseLDLSolver::init() m_mainThreadInvertData = static_cast(this->invertData.get()); } -template -void AsyncSparseLDLSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - if (isAsyncFactorizationFinished() || !m_asyncResult.valid()) - { - SCOPED_TIMER_VARNAME(setSystemMBKMatrixTimer, "setSystemMBKMatrix"); - Inherit1::setSystemMBKMatrix(mparams); - m_hasUpdatedMatrix = true; - } -} - template void AsyncSparseLDLSolver::solveSystem() { @@ -58,16 +47,12 @@ void AsyncSparseLDLSolver::solveSystem() swapInvertData(); } - if (this->linearSystem.needInvert) + if (this->invertData == nullptr) { - if (this->invertData == nullptr) - { - this->getMatrixInvertData(this->getSystemMatrix()); - m_mainThreadInvertData = static_cast(this->invertData.get()); - } - launchAsyncFactorization(); - this->linearSystem.needInvert = false; + this->getMatrixInvertData(this->l_linearSystem->getSystemMatrix()); + m_mainThreadInvertData = static_cast(this->invertData.get()); } + launchAsyncFactorization(); if (waitForAsyncTask) { @@ -81,15 +66,8 @@ void AsyncSparseLDLSolver::solveSystem() swapInvertData(); } - this->solve(*this->getSystemMatrix(), *this->getSystemLHVector(), *this->getSystemRHVector()); - - if (!this->linearSystem.solutionVecId.isNull()) - { - if (this->l_linearSystem) - { - this->l_linearSystem->dispatchSystemSolution(this->linearSystem.solutionVecId); - } - } + this->solve(*this->l_linearSystem->getSystemMatrix(), + *this->l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector()); } template @@ -155,7 +133,7 @@ template void AsyncSparseLDLSolver::asyncFactorization() { newInvertDataReady = false; - this->invert(*this->getSystemMatrix()); + this->invert(*this->l_linearSystem->getSystemMatrix()); newInvertDataReady = true; } 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..8a72beba63b 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 @@ -87,7 +87,6 @@ class PrecomputedLinearSolver : public sofa::component::linearsolver::MatrixLine PrecomputedLinearSolver(); void solve (TMatrix& M, TVector& x, TVector& b) 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; 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..fdf602322ff 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 @@ -52,18 +52,6 @@ PrecomputedLinearSolver::PrecomputedLinearSolver() first = true; } -template -void PrecomputedLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - // 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) @@ -74,7 +62,7 @@ void PrecomputedLinearSolver::solve (TMatrix& , TVector& z, TVe 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(); @@ -194,7 +182,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/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/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..1a582396eb9 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 @@ -62,9 +62,6 @@ class BaseMatrixLinearSolver : public sofa::core::behavior::LinearSolver virtual void invert(Matrix& M) = 0; virtual void solve(Matrix& M, Vector& solution, Vector& rh) = 0; - - virtual Matrix * getSystemMatrix() = 0; - }; /// Empty class used for default solver implementation without multi-threading support @@ -190,56 +187,12 @@ 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; - - /// Get the linear system right-hand term vector, or nullptr if this solver does not build it - 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 - 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(); - /// Invert the system, this method is optional because it's call when solveSystem() is called for the first time void invertSystem() override; @@ -301,7 +254,11 @@ 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); @@ -351,45 +308,7 @@ 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; - - SReal 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; - }; ////////////////////////////////////////////////////////////// @@ -414,36 +333,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..49cab16112b 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,8 +50,6 @@ 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")) { this->addUpdateCallback("parallelInverseProduct", {&d_parallelInverseProduct}, @@ -217,125 +215,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 +226,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->invert(*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 @@ -395,10 +259,9 @@ void MatrixLinearSolver::deleteMatrix(Matrix* v) template void MatrixLinearSolver::invertSystem() { - if (linearSystem.needInvert && l_linearSystem) + if (l_linearSystem) { this->invert(*l_linearSystem->getSystemMatrix()); - linearSystem.needInvert = false; } } @@ -417,18 +280,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->invert(*systemMatrix); simulation::TaskScheduler* taskScheduler = simulation::MainTaskSchedulerFactory::createInRegistry(); assert(taskScheduler); @@ -486,29 +345,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->invert(*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 +374,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 +387,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 +402,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 +419,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 +431,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 +441,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 +469,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 +486,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, &(linearSystem.matrixAccessor)) ); } 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.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h index f54a08d8c00..46b81730b42 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 @@ -56,7 +56,6 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver public: void solve (Matrix& M, Vector& x, Vector& b) override; void init() override; - void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; private : unsigned next_refresh_step; 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..1fd9152fe44 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 @@ -88,46 +88,6 @@ void PCGLinearSolver::init() sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); } -template -void PCGLinearSolver::setSystemMBKMatrix(const core::MechanicalParams* mparams) -{ - sofa::helper::AdvancedTimer::valSet("PCG::buildMBK", 1); - - { - SCOPED_TIMER("PCG::setSystemMBKMatrix"); - Inherit::setSystemMBKMatrix(mparams); - } - - if (l_preconditioner.get()==nullptr) return; - - if (first) //We initialize all the preconditioners for the first step - { - 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 (next_refresh_step >= d_update_step.getValue()) - { - l_preconditioner.get()->setSystemMBKMatrix(mparams); - next_refresh_step=1; - } - else - { - next_refresh_step++; - } - } - } - - l_preconditioner.get()->updateSystemMatrix(); -} - template<> inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta) { @@ -182,9 +142,10 @@ 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 { @@ -207,9 +168,10 @@ void PCGLinearSolver::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 { 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..fe1886839c5 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,7 +122,6 @@ 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; 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..26bc4a3d1d7 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 @@ -77,23 +77,6 @@ void PrecomputedWarpPreconditioner::checkLinearSystem() } } -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; -} - //Solve x = R * M^-1 * R^t * b template void PrecomputedWarpPreconditioner::solve (TMatrix& /*M*/, TVector& z, TVector& r) @@ -337,7 +320,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 +379,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 +410,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/WarpPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h index d1cacfed720..f080f8e955e 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 @@ -66,8 +66,6 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol 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; 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..c30fc282af1 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 @@ -118,93 +118,24 @@ void WarpPreconditioner::bwdInit() 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) +WarpPreconditioner::getSystemDimention(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)) - { - 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 - - } - else if (l_linearSolver.get()->hasUpdatedMatrix()) - { - updateSystemSize = getSystemDimention(mparams); - this->resizeSystem(updateSystemSize); - - 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; - - this->l_linearSystem->getSystemMatrix()->setIdentity(); // identity rotationa after update - } - else - { - 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()); + simulation::common::MechanicalOperations mops(mparams, this->getContext()); - this->l_linearSystem->getSystemMatrix()->opMulTM(this->l_linearSystem->getSystemMatrix(),rotationWork[indexwork]); - } + // this->linearSystem.matrixAccessor.setGlobalMatrix(this->l_linearSystem->getSystemMatrix()); + // this->linearSystem.matrixAccessor.clear(); + // mops.getMatrixDimension(&(this->linearSystem.matrixAccessor)); + // this->linearSystem.matrixAccessor.setupMatrices(); + // return this->linearSystem.matrixAccessor.getGlobalDimension(); + return 0; } template void WarpPreconditioner::invert(Matrix& /*Rcur*/) {} template -void WarpPreconditioner::updateSystemMatrix() { +void WarpPreconditioner::updateSystemMatrix() +{ ++nextRefreshStep; l_linearSolver.get()->updateSystemMatrix(); } @@ -225,33 +156,37 @@ void WarpPreconditioner::checkLinearSystem() /// 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); +void WarpPreconditioner::solve(Matrix& Rcur, Vector& solution, Vector& rh) +{ + Rcur.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(),&rh); - l_linearSolver.get()->solveSystem(); + l_linearSolver->solveSystem(); - Rcur.opMulV(&solution,l_linearSolver.get()->getSystemLHBaseVector()); + Rcur.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..32f42d78f6f 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. */ 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/TypedMatrixLinearSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h index 509bb5a5a90..254cb55835c 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h @@ -71,23 +71,25 @@ 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(); @@ -121,20 +123,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..e475943e464 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl @@ -136,6 +136,30 @@ 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) 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/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.h b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h index 62d5742ca6c..d69ba52d3a5 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 { @@ -45,14 +46,32 @@ class SOFA_CORE_API BaseMatrixLinearSystem : public virtual core::objectmodel::B public: - /// 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); 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/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..2609ac77c60 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h @@ -21,9 +21,10 @@ ******************************************************************************/ #pragma once +#include #include #include -#include +#include namespace sofa::core::behavior { @@ -37,24 +38,11 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver public: SOFA_ABSTRACT_CLASS(LinearSolver, BaseLinearSolver); SOFA_BASE_CAST_IMPLEMENTATION(LinearSolver) -protected: - LinearSolver(); - - ~LinearSolver() override; public: - /// Reset the current linear system. - virtual void resetSystem() = 0; - - /// 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; - /// 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*/){} + 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 @@ -66,23 +54,19 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver /// This function is use for the preconditioner it must be called at each time step event if setSystemMBKMatrix is not called virtual void updateSystemMatrix() {} - /// 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; - - /// 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; - /// 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 +102,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 +125,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); @@ -151,32 +135,11 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver return false; } - /// Get the linear system matrix, or nullptr if this solver does not build it - virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() { return nullptr; } - - /// Get the linear system right-hand term vector, or nullptr if this solver does not build it - virtual linearalgebra::BaseVector* getSystemRHBaseVector() { return nullptr; } - - /// Get the linear system left-hand term vector, or nullptr if this solver does not build it - virtual linearalgebra::BaseVector* getSystemLHBaseVector() { return nullptr; } - - /// Get the linear system inverse matrix, or nullptr if this solver does not build it - virtual linearalgebra::BaseMatrix* getSystemInverseBaseMatrix() { return nullptr; } - /// Read the Matrix solver from a file virtual bool readFile(std::istream& /*in*/) { return false;} /// Read the Matrix solver from a file virtual bool writeFile(std::ostream& /*out*/) {return false;} - - /// Ask the solver to no longer update the system matrix - virtual void freezeSystemMatrix() { frozen = true; } - - - -protected: - - bool frozen; }; } // namespace sofa::core::behavior 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); /// @} From 3f0c14fa4d1c4467bf4cde74721eddf1308f5f5e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 6 Aug 2025 16:57:45 +0200 Subject: [PATCH 02/48] move internal functions to protected --- .../iterative/MatrixLinearSolver.h | 49 ++++++++++--------- 1 file changed, 25 insertions(+), 24 deletions(-) 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 1a582396eb9..e411219b8c1 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 @@ -196,30 +196,6 @@ class MatrixLinearSolver : public BaseMatrixLinea /// 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() @@ -264,6 +240,31 @@ class MatrixLinearSolver : public BaseMatrixLinea 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(); /** From 5f4c3a75428d9eafbadf9b3e95e9af5ef3c50bfc Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 11:41:36 +0200 Subject: [PATCH 03/48] cache factorization + cache invalidation --- .../iterative/MatrixLinearSolver.h | 4 ++ .../iterative/MatrixLinearSolver.inl | 44 ++++++++++++++----- .../linearsystem/TypedMatrixLinearSystem.h | 8 ++++ .../linearsystem/TypedMatrixLinearSystem.inl | 20 +++++++-- 4 files changed, 62 insertions(+), 14 deletions(-) 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 e411219b8c1..84c00ee0e75 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 @@ -208,6 +208,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; @@ -310,6 +312,8 @@ class MatrixLinearSolver : public BaseMatrixLinea virtual MatrixInvertData * createInvertData(); bool singleThreadAddJMInvJtLocal(Matrix * /*M*/,ResMatrixType * result,const JMatrixType * J, SReal fact); + + Data d_factorizationInvalidation; }; ////////////////////////////////////////////////////////////// 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 49cab16112b..d3b0e14eef1 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 @@ -51,7 +51,11 @@ MatrixLinearSolver::MatrixLinearSolver() "matrix with compatible dimensions")) , invertData() , 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 { @@ -177,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(), {})); + system->setName( + this->getContext()->getNameHelper().resolveName(system->getClassName(), {})); system->f_printLog.setValue(this->f_printLog.getValue()); l_linearSystem.set(system); } @@ -226,7 +237,7 @@ void MatrixLinearSolver::solveSystem() } // Step 1: Invert the system, e.g. factorization of the matrix - this->invert(*systemMatrix); + this->invertIfInvalidated(*systemMatrix); // Step 2: Solve the system based on the system inversion this->solve(*systemMatrix, *l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector()); @@ -256,16 +267,27 @@ void MatrixLinearSolver::deleteMatrix(Matrix* v) delete v; } -template -void MatrixLinearSolver::invertSystem() +template +void MatrixLinearSolver::invertSystem() { if (l_linearSystem) { - this->invert(*l_linearSystem->getSystemMatrix()); + 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()) @@ -287,7 +309,7 @@ bool MatrixLinearSolver::addJMInvJtLocal(Matrix* M, ResMatrixTyp return false; } - this->invert(*systemMatrix); + this->invertIfInvalidated(*systemMatrix); simulation::TaskScheduler* taskScheduler = simulation::MainTaskSchedulerFactory::createInRegistry(); assert(taskScheduler); @@ -355,7 +377,7 @@ bool MatrixLinearSolver::singleThreadAddJMInvJtLocal(Matrix* M, auto* rhsVector = l_linearSystem->getRHSVector(); auto* lhsVector = l_linearSystem->getSolutionVector(); - this->invert(*systemMatrix); + this->invertIfInvalidated(*systemMatrix); for (typename JMatrixType::Index row = 0; row < J->rowSize(); ++row) { diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h index 254cb55835c..87e20f62e1d 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.h @@ -93,8 +93,16 @@ class TypedMatrixLinearSystem : public sofa::core::behavior::BaseMatrixLinearSys 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 diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/TypedMatrixLinearSystem.inl index e475943e464..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 @@ -165,12 +175,14 @@ 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 @@ -221,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(); @@ -235,4 +249,4 @@ core::objectmodel::BaseContext* TypedMatrixLinearSystem::getSo return this->getContext(); } -} +} // namespace sofa::component::linearsystem From 181bd27f3d5d1fc0aac415eef1267dabe88d9c3a Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 13:29:13 +0200 Subject: [PATCH 04/48] remove more "assembly" functions --- .../direct/AsyncSparseLDLSolver.h | 3 --- .../direct/AsyncSparseLDLSolver.inl | 24 +++++++------------ .../PrecomputedWarpPreconditioner.h | 2 -- .../preconditioner/WarpPreconditioner.h | 2 -- .../preconditioner/WarpPreconditioner.inl | 7 ------ .../src/sofa/core/behavior/LinearSolver.h | 6 ----- 6 files changed, 8 insertions(+), 36 deletions(-) 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 e719d452ecb..2e1fc4c48ee 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 @@ -68,9 +68,6 @@ class AsyncSparseLDLSolver : public SparseLDLSolver::solveSystem() swapInvertData(); } - if (this->invertData == nullptr) + if (this->d_factorizationInvalidation.getValue()) { - this->getMatrixInvertData(this->l_linearSystem->getSystemMatrix()); - m_mainThreadInvertData = static_cast(this->invertData.get()); + if (this->invertData == nullptr) + { + this->getMatrixInvertData(this->l_linearSystem->getSystemMatrix()); + m_mainThreadInvertData = static_cast(this->invertData.get()); + } + launchAsyncFactorization(); + this->d_factorizationInvalidation.setValue(false); } - launchAsyncFactorization(); if (waitForAsyncTask) { @@ -97,18 +101,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() { 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 fe1886839c5..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 @@ -127,8 +127,6 @@ class PrecomputedWarpPreconditioner : public sofa::component::linearsolver::Matr 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/WarpPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/WarpPreconditioner.h index f080f8e955e..2d9d91fec21 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 @@ -78,8 +78,6 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol void computeResidual(const core::ExecParams* params, linearalgebra::BaseVector* /*f*/) override; - void updateSystemMatrix() override; - protected: void checkLinearSystem() override; 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 c30fc282af1..2378f48b902 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 @@ -133,13 +133,6 @@ WarpPreconditioner::getSystemDimention(const sof template void WarpPreconditioner::invert(Matrix& /*Rcur*/) {} -template -void WarpPreconditioner::updateSystemMatrix() -{ - ++nextRefreshStep; - l_linearSolver.get()->updateSystemMatrix(); -} - template void WarpPreconditioner::checkLinearSystem() { diff --git a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h index 2609ac77c60..7f044dd4e38 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h @@ -48,12 +48,6 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver /// 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; } - - /// This function is use for the preconditioner it must be called at each time step event if setSystemMBKMatrix is not called - virtual void updateSystemMatrix() {} - /// Solve the system as constructed using the previous methods virtual void solveSystem() = 0; From 68ae25a92c53380b134df0ea648670e2e7d87951 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 14:32:21 +0200 Subject: [PATCH 05/48] fix AsyncSparseLDLSolver --- .../direct/AsyncSparseLDLSolver.h | 6 ++- .../direct/AsyncSparseLDLSolver.inl | 47 ++++++++++++++++--- .../MatrixFreeSystem[GraphScattered].cpp | 2 +- .../core/behavior/BaseMatrixLinearSystem.cpp | 13 +++-- .../core/behavior/BaseMatrixLinearSystem.h | 4 +- 5 files changed, 58 insertions(+), 14 deletions(-) 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 2e1fc4c48ee..efa7629b642 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,17 +62,19 @@ class AsyncSparseLDLSolver : public SparseLDLSolver newInvertDataReady { false }; - bool m_hasUpdatedMatrix { false }; + Data< bool > d_authorizeAssembly; }; #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 10b2c986493..ef0649739a9 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 @@ -28,19 +28,41 @@ namespace sofa::component::linearsolver::direct { +template +AsyncSparseLDLSolver::AsyncSparseLDLSolver() + : d_authorizeAssembly(initData(&d_authorizeAssembly, true, "authorizeAssembly", "Allow assembly of the linear system")) +{ +} + template void AsyncSparseLDLSolver::init() { Inherit1::init(); + + if (!this->isComponentStateInvalid()) + { + if (this->l_linearSystem) + { + this->l_linearSystem->d_authorizeAssembly.setParent(&d_authorizeAssembly); + } + + waitForAsyncTask = true; + m_asyncThreadInvertData = &m_secondInvertData; + m_mainThreadInvertData = static_cast(this->invertData.get()); + } +} + +template +void AsyncSparseLDLSolver::reset() +{ + d_authorizeAssembly.setValue(true); waitForAsyncTask = true; - m_asyncThreadInvertData = &m_secondInvertData; - m_mainThreadInvertData = static_cast(this->invertData.get()); } template void AsyncSparseLDLSolver::solveSystem() { - SCOPED_TIMER_VARNAME(invertDataCopyTimer, "AsyncSolve"); + SCOPED_TIMER_VARNAME(invertDataCopyTimer, "solveSystem"); if (newInvertDataReady) { @@ -56,6 +78,9 @@ void AsyncSparseLDLSolver::solveSystem() } launchAsyncFactorization(); this->d_factorizationInvalidation.setValue(false); + + //matrix assembly is temporarily stopped until the next factorization + d_authorizeAssembly.setValue(false); } if (waitForAsyncTask) @@ -70,8 +95,11 @@ void AsyncSparseLDLSolver::solveSystem() swapInvertData(); } - this->solve(*this->l_linearSystem->getSystemMatrix(), - *this->l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector()); + auto* A = this->l_linearSystem->getSystemMatrix(); + auto* b = this->l_linearSystem->getRHSVector(); + auto* x = this->l_linearSystem->getSolutionVector(); + + this->solve(*A, *x, *b); } template @@ -89,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); @@ -124,9 +152,14 @@ void AsyncSparseLDLSolver::launchAsyncFactoriz template void AsyncSparseLDLSolver::asyncFactorization() { + SCOPED_TIMER_TR("asyncFactorization"); + newInvertDataReady = false; this->invert(*this->l_linearSystem->getSystemMatrix()); newInvertDataReady = true; + + //factorization is finished: matrix assembly is authorized once again + d_authorizeAssembly.setValue(true); } template 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..c369962134c 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 @@ -34,7 +34,7 @@ template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API MatrixFreeSystemregisterObjects(core::ObjectRegistrationData("Matrix-free (unbuilt) linear system.") - .add< MatrixFreeSystem >()); + .add< MatrixFreeSystem >(true)); } } //namespace sofa::component::linearsystem diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp index ae24e93ce84..f42dacd8b6f 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_authorizeAssembly(initData(&d_authorizeAssembly, true, "authorizeAssembly", "Allows to assemble the system matrix")) { d_matrixSize.setReadOnly(true); + + d_authorizeAssembly.setReadOnly(true); + d_authorizeAssembly.setDisplayed(false); } void BaseMatrixLinearSystem::buildSystemMatrix(const core::MechanicalParams* mparams) { - preAssembleSystem(mparams); - assembleSystem(mparams); - postAssembleSystem(mparams); + if (d_authorizeAssembly.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 d69ba52d3a5..45abb642054 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h @@ -41,10 +41,12 @@ 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_authorizeAssembly; /// Returns the system matrix as a sofa::linearalgebra::BaseMatrix* virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() const { return nullptr; } From 79c49e89901ea4858d6ed9a58474f5102599c34a Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 15:09:29 +0200 Subject: [PATCH 06/48] add missing templates to the factory --- .../src/sofa/component/linearsystem/MatrixLinearSystem.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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 From b2beca4c36a85d870c52d2ef86797155e6fa3031 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 16:25:10 +0200 Subject: [PATCH 07/48] CompositeLinearSystem can support multiple types of linear systems --- .../MatrixFreeSystem[GraphScattered].cpp | 4 + .../MatrixFreeSystem[GraphScattered].h | 4 +- .../linearsystem/CompositeLinearSystem.h | 10 +- .../linearsystem/CompositeLinearSystem.inl | 93 ++++++------------- .../core/behavior/BaseMatrixLinearSystem.h | 2 +- .../FEMBAR_PCG_JacobiPreconditioner.scn | 9 +- 6 files changed, 43 insertions(+), 79 deletions(-) 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 c369962134c..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 >(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/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h index 32f42d78f6f..71f14b14ea7 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/CompositeLinearSystem.h @@ -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/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h index 45abb642054..2b9ab50cd50 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h @@ -55,7 +55,7 @@ class SOFA_CORE_API BaseMatrixLinearSystem : public virtual core::objectmodel::B 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(); } diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn index b8a29eb500e..c8739f7482e 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn @@ -2,8 +2,13 @@ - - + + + + + + + From dcc559fb4ede2528223b1e2a19e4aca9ed42c050 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 7 Aug 2025 16:56:23 +0200 Subject: [PATCH 08/48] introduce preconditioned systems --- .../LinearSolver/Iterative/CMakeLists.txt | 3 + .../linearsolver/iterative/PCGLinearSolver.h | 5 +- .../iterative/PCGLinearSolver.inl | 63 +++++++++++------ .../PreconditionedMatrixFreeSystem.cpp | 37 ++++++++++ .../PreconditionedMatrixFreeSystem.h | 53 +++++++++++++++ .../PreconditionedMatrixFreeSystem.inl | 68 +++++++++++++++++++ .../component/linearsolver/iterative/init.cpp | 2 + .../FEMBAR_PCG_JacobiPreconditioner.scn | 8 +-- 8 files changed, 211 insertions(+), 28 deletions(-) create mode 100644 Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp create mode 100644 Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h create mode 100644 Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl 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/PCGLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h index 46b81730b42..7c733d0450f 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 @@ -46,8 +46,6 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver 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: @@ -56,6 +54,7 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver public: void solve (Matrix& M, Vector& x, Vector& b) override; void init() override; + void bwdInit() override; private : unsigned next_refresh_step; @@ -72,7 +71,7 @@ private : void handleEvent(sofa::core::objectmodel::Event* event) override; - + void checkLinearSystem(); }; template 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 1fd9152fe44..aaab72a3b2b 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,17 +20,18 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include #include - #include -#include -#include -#include #include #include - +#include +#include + #include +#include "PreconditionedMatrixFreeSystem.h" + namespace sofa::component::linearsolver::iterative { @@ -40,8 +41,6 @@ PCGLinearSolver::PCGLinearSolver() , 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_graph(initData(&d_graph, "graph", "Graph of residuals at each iteration") ) , next_refresh_step(0) , newton_iter(0) @@ -51,30 +50,35 @@ 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_error() << "Link \"preconditioner\" to the desired linear solver should be set to " + "precondition the conjugate gradient."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; } 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") { - 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() << "Can not use the preconditioner " + << l_preconditioner.get()->getName() + << " because it is templated on GraphScatteredType"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; } else @@ -85,10 +89,25 @@ void PCGLinearSolver::init() } first = true; - sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Valid); } -template<> +template +void PCGLinearSolver::bwdInit() +{ + if (l_preconditioner && this->l_linearSystem) + { + if (auto* preconditionerLinearSystem = l_preconditioner->getLinearSystem()) + { + if (auto* preconditionedMatrix = dynamic_cast*>(this->l_linearSystem.get())) + { + preconditionedMatrix->l_preconditionerSystem.set(preconditionerLinearSystem); + } + } + } +} + +template <> inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta) { p.eq(r,p,beta); // p = p*beta + r @@ -100,19 +119,25 @@ inline void PCGLinearSolver -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() +{ + this->template doCheckLinearSystem >(); +} -template +template void PCGLinearSolver::solve (Matrix& M, Vector& x, Vector& b) { SCOPED_TIMER_VARNAME(solveTimer, "PCGLinearSolver::solve"); diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp new file mode 100644 index 00000000000..2fe52650324 --- /dev/null +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp @@ -0,0 +1,37 @@ +/****************************************************************************** + * 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 * + ******************************************************************************/ +#define SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PRECONDITIONEDMATRIXFREESYSTEM_CPP +#include +#include + +namespace sofa::component::linearsolver::iterative +{ + +void registerPreconditionedMatrixFreeSystem(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(core::ObjectRegistrationData("Matrix-free (unbuilt) linear system used in conjunction with built linear system serving in a preconditioner.") + .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..aba4f2d58d7 --- /dev/null +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h @@ -0,0 +1,53 @@ +/****************************************************************************** +* 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 +{ + +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 buildSystemMatrix(const core::MechanicalParams* mparams) override; + void resizeSystem(sofa::Size n) override; + void clearSystem() override; + + ///< The matrix system of the preconditioner + SingleLink l_preconditionerSystem; + +protected: + PreconditionedMatrixFreeSystem(); +}; + + +#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..c5930529211 --- /dev/null +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl @@ -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 + +namespace sofa::component::linearsolver::iterative +{ + +template +PreconditionedMatrixFreeSystem::PreconditionedMatrixFreeSystem() + : l_preconditionerSystem(initLink("preconditionerSystem", "Link toward the linear system of the preconditioner")) +{ +} + +template +void PreconditionedMatrixFreeSystem::buildSystemMatrix( + const core::MechanicalParams* mparams) +{ + linearsystem::MatrixFreeSystem::buildSystemMatrix(mparams); + + if (l_preconditionerSystem) + { + l_preconditionerSystem->buildSystemMatrix(mparams); + } +} + +template +void PreconditionedMatrixFreeSystem::resizeSystem(sofa::Size n) +{ + linearsystem::MatrixFreeSystem::resizeSystem(n); + + if (l_preconditionerSystem) + { + l_preconditionerSystem->resizeSystem(n); + } +} + +template +void PreconditionedMatrixFreeSystem::clearSystem() +{ + linearsystem::MatrixFreeSystem::clearSystem(); + + if (l_preconditionerSystem) + { + l_preconditionerSystem->clearSystem(); + } +} + +} // namespace sofa::component::linearsolver::iterative 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/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn index c8739f7482e..626d6ba18ff 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn @@ -2,12 +2,8 @@ - - - - - - + + From 505c882517273ec4d380a11f60383236acc4159e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 10:38:54 +0200 Subject: [PATCH 09/48] a bit of cleaning --- .../linearsolver/iterative/PCGLinearSolver.inl | 14 +++++++------- .../iterative/PreconditionedMatrixFreeSystem.cpp | 2 +- .../iterative/PreconditionedMatrixFreeSystem.h | 7 +++++++ .../Iterative/FEMBAR_ShewchukPCGLinearSolver.scn | 4 ++-- 4 files changed, 17 insertions(+), 10 deletions(-) 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 aaab72a3b2b..42e2fdf7177 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 @@ -58,10 +58,9 @@ void PCGLinearSolver::init() // Find linear solvers if (l_preconditioner.empty()) { - msg_error() << "Link \"preconditioner\" to the desired linear solver should be set to " - "precondition the conjugate gradient."; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; + 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 { @@ -73,10 +72,9 @@ void PCGLinearSolver::init() } else { - if (l_preconditioner.get()->getTemplateName() == "GraphScattered") + if (l_preconditioner->getTemplateName() == "GraphScattered") { - msg_error() << "Can not use the preconditioner " - << l_preconditioner.get()->getName() + msg_error() << "Can not use the preconditioner " << l_preconditioner->getName() << " because it is templated on GraphScatteredType"; this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; @@ -95,6 +93,7 @@ void PCGLinearSolver::init() template void PCGLinearSolver::bwdInit() { + //link the linear systems of both the preconditioner and the PCG if (l_preconditioner && this->l_linearSystem) { if (auto* preconditionerLinearSystem = l_preconditioner->getLinearSystem()) @@ -134,6 +133,7 @@ void PCGLinearSolver::handleEvent(sofa::core::objectmodel::Event template void PCGLinearSolver::checkLinearSystem() { + // a PreconditionedMatrixFreeSystem component is created in the absence of a linear system this->template doCheckLinearSystem >(); } diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp index 2fe52650324..6d637f9804f 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp @@ -28,7 +28,7 @@ namespace sofa::component::linearsolver::iterative void registerPreconditionedMatrixFreeSystem(sofa::core::ObjectFactory* factory) { - factory->registerObjects(core::ObjectRegistrationData("Matrix-free (unbuilt) linear system used in conjunction with built linear system serving in a preconditioner.") + factory->registerObjects(core::ObjectRegistrationData("A matrix-free linear system that must be used with a preconditioned matrix-free solver") .add< 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 index aba4f2d58d7..3f24920b21f 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h @@ -26,6 +26,13 @@ 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 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 + From 6daf4f3ad1d32955b6c1b4f325bee55179aa4f9d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 10:45:14 +0200 Subject: [PATCH 10/48] convert double to Real --- .../linearsolver/iterative/PCGLinearSolver.h | 18 +++++++------ .../iterative/PCGLinearSolver.inl | 26 +++++++++---------- 2 files changed, 23 insertions(+), 21 deletions(-) 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 7c733d0450f..cfe2d1b3f52 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 @@ -40,13 +40,14 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver typedef TMatrix Matrix; typedef TVector Vector; + using Real = typename Matrix::Real; typedef sofa::component::linearsolver::MatrixLinearSolver Inherit; 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_graph; ///< Graph of residuals at each iteration + Data > > d_graph; ///< Graph of residuals at each iteration protected: PCGLinearSolver(); @@ -64,10 +65,11 @@ private : 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); + 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; @@ -75,23 +77,23 @@ private : }; 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 42e2fdf7177..0e45efc7f6a 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 @@ -107,13 +107,13 @@ void PCGLinearSolver::bwdInit() } template <> -inline void PCGLinearSolver::cgstep_beta(Vector& p, Vector& r, double beta) +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 } @@ -125,7 +125,7 @@ void PCGLinearSolver::handleEvent(sofa::core::objectmodel::Event if (sofa::simulation::AnimateBeginEvent::checkEventType(event)) { newton_iter = 0; - std::map >& graph = *d_graph.beginEdit(); + std::map >& graph = *d_graph.beginEdit(); graph.clear(); } } @@ -142,13 +142,13 @@ 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); @@ -158,8 +158,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) @@ -177,15 +177,15 @@ void PCGLinearSolver::solve (Matrix& M, Vector& x, Vector& b) 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) 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 Date: Fri, 8 Aug 2025 10:54:17 +0200 Subject: [PATCH 11/48] add an error if wrong type --- .../linearsolver/iterative/PCGLinearSolver.inl | 11 +++++++++++ .../Component/ODESolver/Backward/BDFOdeSolver.scn | 4 ++-- 2 files changed, 13 insertions(+), 2 deletions(-) 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 0e45efc7f6a..61472e77ac8 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 @@ -93,6 +93,9 @@ void PCGLinearSolver::init() template void PCGLinearSolver::bwdInit() { + if (this->isComponentStateInvalid()) + return; + //link the linear systems of both the preconditioner and the PCG if (l_preconditioner && this->l_linearSystem) { @@ -100,8 +103,16 @@ void PCGLinearSolver::bwdInit() { if (auto* preconditionedMatrix = dynamic_cast*>(this->l_linearSystem.get())) { + 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 + { + msg_error() << "The preconditioned linear system (" << preconditionedMatrix->getPathName() << ") is not a PreconditionedMatrixFreeSystem"; + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Invalid); + return; + } } } } 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"/> - - + + From 294cc4fdfcdb2301f899f073ac9e7e8f83ad0305 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 11:23:56 +0200 Subject: [PATCH 12/48] Restoration of the parameter controlling the assembly rate of the preconditioned matrix. It was the Data 'update_step' in PCGLinearSolver, but now it's the Data 'assemblingRate' in PreconditionedMatrixFreeSystem. Benchmark of examples\Component\LinearSolver\Preconditioner\FEMBAR_PCG_JacobiPreconditioner.scn on 1000 time steps with variation of `assemblingRate`: 1: 700.128 FPS 2: 968.381 FPS 3: 1074.87 FPS 4: 1156.56 FPS 10: 1383.11 FPS 20: 1478.05 FPS 30: 1561.64 FPS 60: 1563.85 FPS --- .../PreconditionedMatrixFreeSystem.h | 6 +++++ .../PreconditionedMatrixFreeSystem.inl | 25 ++++++++++++++++++- .../FEMBAR_PCG_JacobiPreconditioner.scn | 1 + 3 files changed, 31 insertions(+), 1 deletion(-) 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 index 3f24920b21f..fdec51f2df8 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h @@ -41,6 +41,8 @@ class PreconditionedMatrixFreeSystem 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; @@ -48,8 +50,12 @@ class PreconditionedMatrixFreeSystem ///< The matrix system of the preconditioner SingleLink l_preconditionerSystem; + Data d_assemblingRate; + protected: PreconditionedMatrixFreeSystem(); + + unsigned int m_assemblyCounter {}; }; 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 index c5930529211..5557240f88f 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl @@ -28,7 +28,24 @@ 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(); + + m_assemblyCounter = d_assemblingRate.getValue(); // to assemble the first time +} + +template +void PreconditionedMatrixFreeSystem::reset() +{ + m_assemblyCounter = d_assemblingRate.getValue(); // to assemble the first time } template @@ -39,7 +56,11 @@ void PreconditionedMatrixFreeSystem::buildSystemMatrix( if (l_preconditionerSystem) { - l_preconditionerSystem->buildSystemMatrix(mparams); + if (++m_assemblyCounter >= d_assemblingRate.getValue()) + { + l_preconditionerSystem->buildSystemMatrix(mparams); + m_assemblyCounter = 0; + } } } @@ -51,6 +72,7 @@ void PreconditionedMatrixFreeSystem::resizeSystem(sofa::Size n if (l_preconditionerSystem) { l_preconditionerSystem->resizeSystem(n); + m_assemblyCounter = 0; } } @@ -62,6 +84,7 @@ void PreconditionedMatrixFreeSystem::clearSystem() if (l_preconditionerSystem) { l_preconditionerSystem->clearSystem(); + m_assemblyCounter = 0; } } diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn index 626d6ba18ff..d8896c44a1e 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn @@ -2,6 +2,7 @@ + From 74783e307670bad4669f52c255f6e6299ea25e59 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 11:37:42 +0200 Subject: [PATCH 13/48] restore removed Data with a deprecation layer --- .../linearsolver/iterative/PCGLinearSolver.h | 1 + .../linearsolver/iterative/PCGLinearSolver.inl | 11 +++++++++++ .../iterative/PreconditionedMatrixFreeSystem.h | 2 ++ .../iterative/PreconditionedMatrixFreeSystem.inl | 8 +++++++- 4 files changed, 21 insertions(+), 1 deletion(-) 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 cfe2d1b3f52..962bd28efd0 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 @@ -47,6 +47,7 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver 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_graph; ///< Graph of residuals at each iteration protected: 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 61472e77ac8..199ff1e0fed 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 @@ -41,6 +41,7 @@ PCGLinearSolver::PCGLinearSolver() , 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_graph(initData(&d_graph, "graph", "Graph of residuals at each iteration") ) , next_refresh_step(0) , newton_iter(0) @@ -55,6 +56,16 @@ void PCGLinearSolver::init() { Inherit1::init(); + if (d_update_step.isSet()) + { + if (auto* preconditionedMatrix = dynamic_cast*>(this->l_linearSystem.get())) + { + msg_deprecated() << "The Data '" << d_update_step.getName() << "' is deprecated. Instead use the Data '" << preconditionedMatrix->d_assemblingRate.getName() << "' in " << preconditionedMatrix->getPathName() << "."; + preconditionedMatrix->d_assemblingRate.setValue(d_update_step.getValue()); + preconditionedMatrix->reinitAssemblyCounter(); + } + } + // Find linear solvers if (l_preconditioner.empty()) { 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 index fdec51f2df8..cf67bdddd4e 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.h @@ -52,6 +52,8 @@ class PreconditionedMatrixFreeSystem Data d_assemblingRate; + void reinitAssemblyCounter(); + protected: PreconditionedMatrixFreeSystem(); 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 index 5557240f88f..b4a710b7c50 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl @@ -39,13 +39,19 @@ 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() { - m_assemblyCounter = d_assemblingRate.getValue(); // to assemble the first time + reinitAssemblyCounter(); } template From 4bcf3ddba991709dc59ff3af359b00bd2ef4a0d7 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 13:52:55 +0200 Subject: [PATCH 14/48] restore initial parameter for regression tests --- .../linearsolver/iterative/PreconditionedMatrixFreeSystem.inl | 2 ++ .../Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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 index b4a710b7c50..1e5defbb51a 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.inl @@ -58,8 +58,10 @@ 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()) diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn index d8896c44a1e..2f7a003c9cf 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn @@ -2,7 +2,7 @@ - + From ae0577c11c129cb3676d626068a0875bbfe58f7b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 8 Aug 2025 15:23:21 +0200 Subject: [PATCH 15/48] fix PrecomputedLinearSolver --- .../direct/PrecomputedLinearSolver.h | 8 ++++---- .../direct/PrecomputedLinearSolver.inl | 19 +++++++++++++++---- 2 files changed, 19 insertions(+), 8 deletions(-) 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 8a72beba63b..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,10 +82,10 @@ 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 loadMatrix(TMatrix& M); void loadMatrixWithCholeskyDecomposition(TMatrix& M); @@ -119,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 fdf602322ff..ada8a553931 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 @@ -52,11 +52,11 @@ PrecomputedLinearSolver::PrecomputedLinearSolver() first = true; } -//Solve x = R * M^-1 * R^t * b template -void PrecomputedLinearSolver::solve (TMatrix& , TVector& z, TVector& r) +void PrecomputedLinearSolver::solve (TMatrix& M, TVector& solution, TVector& rh) { - z = internalData.Minv * r; + SOFA_UNUSED(M); + solution = internalData.Minv * rh; } template @@ -139,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_authorizeAssembly.setValue(false); + } + first = false; + } +} template template void PrecomputedLinearSolver::computeActiveDofs(JMatrix& J) From 3536bbc212b86cac8ff8a80daf8d2d65b4c463d4 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 15 Sep 2025 16:12:09 +0200 Subject: [PATCH 16/48] mark WarpPreconditioner destructor as `override` --- .../component/linearsolver/preconditioner/WarpPreconditioner.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 2d9d91fec21..72404e637bb 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 @@ -61,7 +61,7 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol WarpPreconditioner(); public: - ~WarpPreconditioner(); + ~WarpPreconditioner() override; void init() override; void bwdInit() override; From bd2dcccedace11b08f8c05eda08da4a09bcf427c Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 15 Sep 2025 17:06:32 +0200 Subject: [PATCH 17/48] minor cleaning --- .../linearsolver/preconditioner/WarpPreconditioner.inl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 2378f48b902..13a130c8973 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 @@ -84,15 +84,15 @@ 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") + if (l_linearSolver->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); + 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 From 949bedac617386d9e94b7939c524a66a0332e148 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 15 Sep 2025 17:06:47 +0200 Subject: [PATCH 18/48] simplify `checkLinearSystem` in WarpPreconditioner implementation --- .../linearsolver/preconditioner/WarpPreconditioner.inl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) 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 13a130c8973..4dcac228808 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 @@ -136,14 +136,8 @@ void WarpPreconditioner::invert(Matrix& /*Rcur*/ 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>(); } From bea6494d49d859e48ae5c9538c2852997dfc1723 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 15 Sep 2025 17:07:50 +0200 Subject: [PATCH 19/48] add license header to RotationMatrixSystem header file --- .../preconditioner/RotationMatrixSystem.h | 23 ++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) 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..78e08699cd9 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,4 +1,25 @@ -#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 From 310b266bb922e0db8a675abf5365f6c68d32a67b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 15 Sep 2025 17:08:18 +0200 Subject: [PATCH 20/48] fix typo in registration description of RotationMatrixSystem --- .../linearsolver/preconditioner/RotationMatrixSystem.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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..0691f08151a 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 @@ -33,7 +33,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 > >()); } From 8f8e03061d3b4967525863d62a86af05c10472d1 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 09:39:17 +0200 Subject: [PATCH 21/48] modernize code style --- .../iterative/PCGLinearSolver.cpp | 2 +- .../linearsolver/iterative/PCGLinearSolver.h | 24 ++++++++++--------- .../iterative/PCGLinearSolver.inl | 7 +++--- 3 files changed, 17 insertions(+), 16 deletions(-) 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 962bd28efd0..c4d6cb01c20 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,17 +31,21 @@ 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; + 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; - typedef sofa::component::linearsolver::MatrixLinearSolver Inherit; + 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) @@ -50,20 +54,18 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver Data d_update_step; ///< Number of steps before the next refresh of preconditioners Data > > d_graph; ///< Graph of residuals at each iteration -protected: - PCGLinearSolver(); - -public: void solve (Matrix& M, Vector& x, Vector& b) override; void init() override; void bwdInit() override; -private : +private: unsigned next_refresh_step; bool first; int newton_iter; protected: + PCGLinearSolver(); + /// 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); @@ -74,7 +76,7 @@ private : void handleEvent(sofa::core::objectmodel::Event* event) override; - void checkLinearSystem(); + void checkLinearSystem() override; }; template 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 199ff1e0fed..07d5648d0e6 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 @@ -22,6 +22,7 @@ #pragma once #include #include +#include #include #include #include @@ -30,18 +31,16 @@ #include -#include "PreconditionedMatrixFreeSystem.h" - 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_update_step(initData(&d_update_step, 1u, "update_step", "Number of steps before the next refresh of preconditioners") ) , d_graph(initData(&d_graph, "graph", "Graph of residuals at each iteration") ) , next_refresh_step(0) , newton_iter(0) From d00e57bfa4c04bda6bd815b3c8d8218b72b6e434 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 15:19:34 +0200 Subject: [PATCH 22/48] add missing `RotationMatrixSystem.inl` to CMakeLists.txt --- Sofa/Component/LinearSolver/Preconditioner/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) 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 ) From 978df3939a4ebcb0ca95ade6dd171d5361b18a28 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 15:23:31 +0200 Subject: [PATCH 23/48] replace `addSlave` with `addObject` in MatrixLinearSolver This is to allow the added object to be initialized. Slaves are not initialized because visitors don't visit them. They are owned by their master, then it's up to their master to initialize them or not. There is no reason for a linear system to be a slave. --- .../component/linearsolver/iterative/MatrixLinearSolver.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 d3b0e14eef1..944ec2c8487 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 @@ -196,7 +196,7 @@ void MatrixLinearSolver::createDefaultLinearSystem() { if (auto system = sofa::core::objectmodel::New()) { - this->addSlave(system); + this->getContext()->addObject(system); system->setName( this->getContext()->getNameHelper().resolveName(system->getClassName(), {})); system->f_printLog.setValue(this->f_printLog.getValue()); From 29dc74758fa237791416b8bdf5b2f1cfcd48fcfc Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 15:38:30 +0200 Subject: [PATCH 24/48] add SOFA class macros for proper class hierarchy --- .../solidmechanics/fem/elastic/HexahedronFEMForceField.h | 9 +++++++-- .../Core/src/sofa/core/behavior/BaseRotationFinder.h | 2 ++ .../Core/src/sofa/core/behavior/RotationFinder.h | 2 ++ 3 files changed, 11 insertions(+), 2 deletions(-) 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/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/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; From db723247ffdb34b10fca45193a415a7497e48d0e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 15:38:51 +0200 Subject: [PATCH 25/48] Missing space --- .../component/linearsolver/iterative/MatrixLinearSolver.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 944ec2c8487..e1502d3413f 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 @@ -127,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 " From cd5ef4172ea5051f910071de0b81da05fdd33dae Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 16 Sep 2025 15:45:18 +0200 Subject: [PATCH 26/48] proper initialization of systems and solvers --- .../preconditioner/RotationMatrixSystem.cpp | 3 +- .../preconditioner/RotationMatrixSystem.h | 33 ++++- .../preconditioner/RotationMatrixSystem.inl | 124 +++++++++++++++++- .../preconditioner/WarpPreconditioner.h | 11 -- .../preconditioner/WarpPreconditioner.inl | 60 +-------- 5 files changed, 164 insertions(+), 67 deletions(-) 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 0691f08151a..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 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 78e08699cd9..d8a60de438b 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 @@ -22,6 +22,8 @@ #pragma once #include #include +#include + namespace sofa::component::linearsolver::preconditioner { @@ -30,5 +32,34 @@ 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(); + + unsigned int m_assemblyCounter {}; + + std::array, 2> rotationWork; + std::size_t indexwork {}; + + void reinitAssemblyCounter(); + void updateMatrixWithRotations(); + void ensureValidRotationWork(); }; -} + +#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..d475a269be6 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,127 @@ -#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()->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; + } + } + + 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 isAssembled = false; + + if (l_mainAssembledSystem) + { + if (++m_assemblyCounter >= d_assemblingRate.getValue()) + { + l_mainAssembledSystem->buildSystemMatrix(mparams); + m_assemblyCounter = 0; + isAssembled = true; + } + } + + if (isAssembled) + { + 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) + { + ensureValidRotationWork(); + l_rotationFinder->getRotations(rotationWork[indexwork].get()); + } +} + +template +void RotationMatrixSystem::ensureValidRotationWork() +{ + if (indexwork < rotationWork.size()) + { + rotationWork[indexwork] = std::make_unique(); + } + else + { + msg_error() << "Wrong index"; + } } + +} // 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 72404e637bb..1b2d81c16f4 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 @@ -54,16 +54,11 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol 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 - protected: WarpPreconditioner(); public: - ~WarpPreconditioner() override; - void init() override; void bwdInit() override; void invert(Matrix& M) override; @@ -74,8 +69,6 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol 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; protected: @@ -86,12 +79,8 @@ 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 4dcac228808..208ced75ecf 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,32 +44,7 @@ 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" ) ) -{ - rotationWork[0] = nullptr; - rotationWork[1] = nullptr; - - first = true; - indexwork = 0; -} - -template -WarpPreconditioner::~WarpPreconditioner() -{ - if (rotationWork[0]) delete rotationWork[0]; - if (rotationWork[1]) delete rotationWork[1]; - - rotationWork[0] = nullptr; - rotationWork[1] = nullptr; -} - -template -void WarpPreconditioner::init() -{ - Inherit1::init(); - first = true; -} +{} template void WarpPreconditioner::bwdInit() @@ -101,33 +76,10 @@ void WarpPreconditioner::bwdInit() } } - 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; + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); } - 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->l_linearSystem->getSystemMatrix()); - // this->linearSystem.matrixAccessor.clear(); - // mops.getMatrixDimension(&(this->linearSystem.matrixAccessor)); - // this->linearSystem.matrixAccessor.setupMatrices(); - // return this->linearSystem.matrixAccessor.getGlobalDimension(); - return 0; } template @@ -145,10 +97,12 @@ void WarpPreconditioner::checkLinearSystem() template void WarpPreconditioner::solve(Matrix& Rcur, Vector& solution, Vector& rh) { - Rcur.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(),&rh); + // rh = Rcur^T * rhs + Rcur.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rh); l_linearSolver->solveSystem(); + // solution = Rcur * solution Rcur.opMulV(&solution, l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector()); } From 48f707b9f50c2e949582e5740ca4f9429478d39a Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 11:06:27 +0200 Subject: [PATCH 27/48] validate `l_linearSystem` type in `PCGLinearSolver` --- .../linearsolver/iterative/PCGLinearSolver.h | 6 ++-- .../iterative/PCGLinearSolver.inl | 28 ++++++++++++++++++- 2 files changed, 31 insertions(+), 3 deletions(-) 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 c4d6cb01c20..8aa89d788c3 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 @@ -66,8 +66,10 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver protected: PCGLinearSolver(); - /// 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 + 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. 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 07d5648d0e6..90264b8524c 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 @@ -96,8 +96,34 @@ void PCGLinearSolver::init() } } + ensureRequiredLinearSystemType(); + if (this->isComponentStateInvalid()) + return; + first = true; - this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Valid); + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue( sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void PCGLinearSolver::ensureRequiredLinearSystemType() +{ + 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 " + << PreconditionedMatrixFreeSystem::GetClass()->className + << " linear system, but a " << this->l_linearSystem->getClassName() + << " was found"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } + } } template From cc4cd1bffc261f7681d0192c1cb06492fcbb6563 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 11:22:45 +0200 Subject: [PATCH 28/48] fix typo in error message of `PCGLinearSolver` --- .../sofa/component/linearsolver/iterative/PCGLinearSolver.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 90264b8524c..fad121e3742 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 @@ -84,7 +84,7 @@ void PCGLinearSolver::init() { if (l_preconditioner->getTemplateName() == "GraphScattered") { - msg_error() << "Can not use the preconditioner " << l_preconditioner->getName() + 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; From 8c248740517e75cc9d3043357bb670682709e3dd Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 11:27:32 +0200 Subject: [PATCH 29/48] validate `l_mainAssembledSystem` and handle missing or invalid linear systems in `RotationMatrixSystem` --- .../preconditioner/RotationMatrixSystem.inl | 40 ++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) 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 d475a269be6..245912e0515 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 @@ -55,6 +55,39 @@ void RotationMatrixSystem::init() } } + 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()->getObjects(sofa::core::objectmodel::BaseContext::Local); + for (const auto& system : listSystems) + { + if (system->getTemplateName() != "GraphScattered") + { + 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()) @@ -107,6 +140,8 @@ void RotationMatrixSystem::updateMatrixWithRotations() if (l_rotationFinder) { ensureValidRotationWork(); + const auto matrixSize = l_mainAssembledSystem->getMatrixSize(); + rotationWork[indexwork]->resize(matrixSize[0], matrixSize[1]); l_rotationFinder->getRotations(rotationWork[indexwork].get()); } } @@ -116,7 +151,10 @@ void RotationMatrixSystem::ensureValidRotationWork() { if (indexwork < rotationWork.size()) { - rotationWork[indexwork] = std::make_unique(); + if (auto& matrix = rotationWork[indexwork]; matrix == nullptr) + { + matrix = std::make_unique(); + } } else { From 1774746757c8d2562a484a767d388f5f867bea06 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 11:33:07 +0200 Subject: [PATCH 30/48] validate `l_linearSystem` type and ensure compatibility with `RotationMatrixSystem` in `WarpPreconditioner` --- .../preconditioner/WarpPreconditioner.h | 3 +- .../preconditioner/WarpPreconditioner.inl | 31 +++++++++++++++++-- 2 files changed, 31 insertions(+), 3 deletions(-) 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 1b2d81c16f4..027195cd7af 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 @@ -59,6 +59,7 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol public: + void init() override; void bwdInit() override; void invert(Matrix& M) override; @@ -74,6 +75,7 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol protected: void checkLinearSystem() override; + void ensureRequiredLinearSystemType(); private : @@ -81,7 +83,6 @@ private : unsigned nextRefreshStep {}; - 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 208ced75ecf..4fea2e277ee 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,9 +44,36 @@ namespace sofa::component::linearsolver::preconditioner template WarpPreconditioner::WarpPreconditioner() : l_linearSolver(initLink("linearSolver", "Link towards the linear solver used to build the warp conditioner")) -{} +{ +} -template +template +void WarpPreconditioner::init() +{ + Inherit1::init(); + + ensureRequiredLinearSystemType(); +} + +template +void WarpPreconditioner::ensureRequiredLinearSystemType() +{ + 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 void WarpPreconditioner::bwdInit() { if (l_linearSolver.empty()) From 51d5f2ac92212185f50123c0a3685ab32a8af443 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 11:34:55 +0200 Subject: [PATCH 31/48] don't link to itself --- .../linearsolver/preconditioner/RotationMatrixSystem.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 245912e0515..34bdec7a73f 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 @@ -63,7 +63,7 @@ void RotationMatrixSystem::init() const auto listSystems = this->getContext()->getObjects(sofa::core::objectmodel::BaseContext::Local); for (const auto& system : listSystems) { - if (system->getTemplateName() != "GraphScattered") + if (system->getTemplateName() != "GraphScattered" && system != this) { l_mainAssembledSystem.set(system); break; From 733afcfe66b131757579e8cb693133af138f0464 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 17 Sep 2025 17:27:53 +0200 Subject: [PATCH 32/48] refactor `RotationMatrixSystem` and `WarpPreconditioner` for clearer logic and enhanced matrix handling --- .../preconditioner/RotationMatrixSystem.h | 4 -- .../preconditioner/RotationMatrixSystem.inl | 29 +++--------- .../preconditioner/WarpPreconditioner.h | 2 +- .../preconditioner/WarpPreconditioner.inl | 45 ++++++++++--------- .../FEMBAR_PCG_WarpPreconditioner.scn | 12 +++-- 5 files changed, 41 insertions(+), 51 deletions(-) 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 d8a60de438b..c8958f8d143 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 @@ -50,12 +50,8 @@ class RotationMatrixSystem : public linearsystem::TypedMatrixLinearSystem, 2> rotationWork; - std::size_t indexwork {}; - void reinitAssemblyCounter(); void updateMatrixWithRotations(); - void ensureValidRotationWork(); }; #if !defined(SOFA_COMPONENT_LINEARSOLVER_PRECONDITIONER_ROTATIONMATRIXSYSTEM_CPP) 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 34bdec7a73f..352699390ff 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 @@ -106,7 +106,7 @@ template void RotationMatrixSystem::buildSystemMatrix( const core::MechanicalParams* mparams) { - bool isAssembled = false; + bool isMainSystemAssembled = false; if (l_mainAssembledSystem) { @@ -114,12 +114,14 @@ void RotationMatrixSystem::buildSystemMatrix( { l_mainAssembledSystem->buildSystemMatrix(mparams); m_assemblyCounter = 0; - isAssembled = true; + isMainSystemAssembled = true; } } - if (isAssembled) + if (isMainSystemAssembled) { + const auto matrixSize = l_mainAssembledSystem->getMatrixSize(); + this->resizeSystem(matrixSize[0]); this->getSystemMatrix()->setIdentity(); } else @@ -139,26 +141,7 @@ void RotationMatrixSystem::updateMatrixWithRotations() { if (l_rotationFinder) { - ensureValidRotationWork(); - const auto matrixSize = l_mainAssembledSystem->getMatrixSize(); - rotationWork[indexwork]->resize(matrixSize[0], matrixSize[1]); - l_rotationFinder->getRotations(rotationWork[indexwork].get()); - } -} - -template -void RotationMatrixSystem::ensureValidRotationWork() -{ - if (indexwork < rotationWork.size()) - { - if (auto& matrix = rotationWork[indexwork]; matrix == nullptr) - { - matrix = std::make_unique(); - } - } - else - { - msg_error() << "Wrong index"; + l_rotationFinder->getRotations(this->getSystemBaseMatrix()); } } 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 027195cd7af..dd95af7cdfa 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 @@ -64,7 +64,7 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol void invert(Matrix& M) override; - void solve(Matrix& M, Vector& solution, Vector& rh) override; + void solve(Matrix& M, Vector& solution, Vector& rhs) override; bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; 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 4fea2e277ee..9ecf54883e3 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 @@ -78,7 +78,7 @@ 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) ); } @@ -89,20 +89,17 @@ void WarpPreconditioner::bwdInit() this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; } - else + + if (l_linearSolver->getTemplateName() == "GraphScattered") { - if (l_linearSolver->getTemplateName() == "GraphScattered") - { - 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 - { - msg_info() << "LinearSolver path used: '" << l_linearSolver.getLinkedPath() << "'"; - } + 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; } + msg_info() << "LinearSolver path used: '" << l_linearSolver.getLinkedPath() << "'"; + if (!this->isComponentStateInvalid()) { this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); @@ -119,18 +116,26 @@ void WarpPreconditioner::checkLinearSystem() this->template doCheckLinearSystem>(); } - -/// Solve the system as constructed using the previous methods template -void WarpPreconditioner::solve(Matrix& Rcur, Vector& solution, Vector& rh) +void WarpPreconditioner::solve(Matrix& M, Vector& solution, Vector& rhs) { - // rh = Rcur^T * rhs - Rcur.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rh); - + // 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: + M.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rhs); + + // Step 2: + // Solve A * R^T * x = R^T * b using the linear solver. The solution stored in the linear + // solver is then y = R^T * x l_linearSolver->solveSystem(); - // solution = Rcur * solution - Rcur.opMulV(&solution, l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector()); + // Step 3: + // y = R^T * x <=> x = R * y + M.opMulV(&solution, l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector()); } /// Solve the system as constructed using the previous methods diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn index aa4c61c823e..70cf11b7de0 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn @@ -2,9 +2,15 @@ - - - + + + + + + + + + From d31ed9241e00e5a467663e4b9fad922a0627c496 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 18 Sep 2025 09:29:59 +0200 Subject: [PATCH 33/48] some comments --- .../Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn index 70cf11b7de0..51b2a59147e 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn @@ -2,12 +2,15 @@ + + + From cdfbd4677b3831e5ae961c04af19fdf1572061d0 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 18 Sep 2025 10:29:20 +0200 Subject: [PATCH 34/48] proper deprecation --- .../linearsolver/iterative/PCGLinearSolver.h | 2 +- .../linearsolver/iterative/PCGLinearSolver.inl | 12 +----------- .../linearsolver/preconditioner/WarpPreconditioner.h | 2 ++ .../preconditioner/WarpPreconditioner.inl | 2 ++ 4 files changed, 6 insertions(+), 12 deletions(-) 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 8aa89d788c3..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 @@ -51,7 +51,7 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver 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 + 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 void solve (Matrix& M, Vector& x, Vector& b) override; 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 fad121e3742..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 @@ -40,7 +40,7 @@ PCGLinearSolver::PCGLinearSolver() , 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, 1u, "update_step", "Number of steps before the next refresh of preconditioners") ) + , 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) @@ -55,16 +55,6 @@ void PCGLinearSolver::init() { Inherit1::init(); - if (d_update_step.isSet()) - { - if (auto* preconditionedMatrix = dynamic_cast*>(this->l_linearSystem.get())) - { - msg_deprecated() << "The Data '" << d_update_step.getName() << "' is deprecated. Instead use the Data '" << preconditionedMatrix->d_assemblingRate.getName() << "' in " << preconditionedMatrix->getPathName() << "."; - preconditionedMatrix->d_assemblingRate.setValue(d_update_step.getValue()); - preconditionedMatrix->reinitAssemblyCounter(); - } - } - // Find linear solvers if (l_preconditioner.empty()) { 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 dd95af7cdfa..059b1c1525e 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,6 +53,8 @@ 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 + core::objectmodel::lifecycle::DeprecatedData d_useRotationFinder; + core::objectmodel::lifecycle::DeprecatedData d_updateStep; protected: WarpPreconditioner(); 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 9ecf54883e3..642ace8efbf 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,6 +44,8 @@ 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(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" ) { } From d828970a704a5a07dec21413dab7939e01509903 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 18 Sep 2025 10:36:44 +0200 Subject: [PATCH 35/48] update scene --- .../FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) 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 @@ - - + + + + + + + From 581d3455553be24cb424f5ebb5571c6fe5283aaa Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 18 Sep 2025 13:22:26 +0200 Subject: [PATCH 36/48] minor changes --- .../preconditioner/WarpPreconditioner.h | 2 +- .../preconditioner/WarpPreconditioner.inl | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) 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 059b1c1525e..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 @@ -66,7 +66,7 @@ class WarpPreconditioner : public sofa::component::linearsolver::MatrixLinearSol void invert(Matrix& M) override; - void solve(Matrix& M, Vector& solution, Vector& rhs) override; + void solve(Matrix& R, Vector& solution, Vector& rhs) override; bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; 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 642ace8efbf..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 @@ -119,7 +119,7 @@ void WarpPreconditioner::checkLinearSystem() } template -void WarpPreconditioner::solve(Matrix& M, Vector& solution, Vector& rhs) +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 @@ -127,17 +127,17 @@ void WarpPreconditioner::solve(Matrix& M, Vector // Step 1: // R * A * R^T * x = b <=> A * R^T * x = R^T * b - // R^T * b is computed: - M.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rhs); + // R^T * b is computed in this step: + R.opMulTV(l_linearSolver->getLinearSystem()->getSystemRHSBaseVector(), &rhs); // Step 2: - // Solve A * R^T * x = R^T * b using the linear solver. The solution stored in the linear - // solver is then y = R^T * x + // 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: - // y = R^T * x <=> x = R * y - M.opMulV(&solution, l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector()); + // 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 From 223262a0639ff7ff82dd7b6f4c74008551cc906b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 18 Sep 2025 13:48:16 +0200 Subject: [PATCH 37/48] add `template` keyword to `get` and `getObjects` calls in `RotationMatrixSystem` for proper type resolution --- .../linearsolver/preconditioner/RotationMatrixSystem.inl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 352699390ff..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 @@ -40,7 +40,7 @@ void RotationMatrixSystem::init() if (!l_rotationFinder) { - if (auto* rotationFinder = this->getContext()->get()) + if (auto* rotationFinder = this->getContext()->template get()) { l_rotationFinder.set(rotationFinder); msg_info() << "Rotation finder found: '" << l_rotationFinder->getPathName() << "'"; @@ -60,7 +60,7 @@ void RotationMatrixSystem::init() 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()->getObjects(sofa::core::objectmodel::BaseContext::Local); + const auto listSystems = this->getContext()->template getObjects(sofa::core::objectmodel::BaseContext::Local); for (const auto& system : listSystems) { if (system->getTemplateName() != "GraphScattered" && system != this) From 539c5ec281aa8817eab838b543c2e0a890470cbd Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 19 Sep 2025 10:08:15 +0200 Subject: [PATCH 38/48] invert includes order to allow object factory to find the component module --- .../linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp index 6d637f9804f..1fc4dc6073e 100644 --- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp +++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PreconditionedMatrixFreeSystem.cpp @@ -20,8 +20,8 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #define SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PRECONDITIONEDMATRIXFREESYSTEM_CPP -#include #include +#include namespace sofa::component::linearsolver::iterative { From fa2e52bec8d49a509620badf8c3663030c314acf Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 19 Sep 2025 10:11:08 +0200 Subject: [PATCH 39/48] apply changes --- .../FEM/Heterogeneous-TetrahedronFEMForceField.scn | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 @@ - + + From e0aee69a3d2b6ae940b34a279ac1b945eea174cc Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 9 Oct 2025 17:44:26 +0200 Subject: [PATCH 40/48] deprecation --- .../iterative/MatrixLinearSolver.h | 26 ++++++++ .../linearsolver/iterative/config.h.in | 6 ++ .../src/sofa/core/behavior/LinearSolver.h | 61 ++++++++++++++++++- Sofa/framework/Core/src/sofa/core/config.h.in | 8 +++ 4 files changed, 99 insertions(+), 2 deletions(-) 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 84c00ee0e75..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 @@ -62,6 +62,10 @@ class BaseMatrixLinearSolver : public sofa::core::behavior::LinearSolver virtual void invert(Matrix& M) = 0; virtual void solve(Matrix& M, Vector& solution, Vector& rh) = 0; + + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual Matrix * getSystemMatrix() final = delete; + }; /// Empty class used for default solver implementation without multi-threading support @@ -187,12 +191,28 @@ class MatrixLinearSolver : public BaseMatrixLinea void init() override; + /// Reset the current linear system. + 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; } + /// Returns the linear system component associated to the linear solver 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 + 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; @@ -311,6 +331,12 @@ class MatrixLinearSolver : public BaseMatrixLinea virtual MatrixInvertData * createInvertData(); + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved linearSystem; + + SOFA_ITERATIVE_SOLVER_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved currentMFactor, currentBFactor, currentKFactor; + bool singleThreadAddJMInvJtLocal(Matrix * /*M*/,ResMatrixType * result,const JMatrixType * J, SReal fact); Data d_factorizationInvalidation; 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..c70503f0ffe 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,9 @@ 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_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/framework/Core/src/sofa/core/behavior/LinearSolver.h b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h index 7f044dd4e38..7bc3a8402c9 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/LinearSolver.h @@ -21,6 +21,7 @@ ******************************************************************************/ #pragma once +#include #include #include #include @@ -36,9 +37,23 @@ 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) -public: + + /// Reset the current linear system. + 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) + 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_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void rebuildSystem(SReal /*massFactor*/, SReal /*forceFactor*/) final = delete; virtual sofa::core::behavior::BaseMatrixLinearSystem* getLinearSystem() const = 0; @@ -48,6 +63,23 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver /// 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) + 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 + 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 + 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 + 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; @@ -129,11 +161,36 @@ class SOFA_CORE_API LinearSolver : public BaseLinearSolver return false; } + /// Get the linear system matrix, or nullptr if this solver does not build it + 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 + 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 + 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 + 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;} /// Read the Matrix solver from a file virtual bool writeFile(std::ostream& /*out*/) {return false;} + + /// Ask the solver to no longer update the system matrix + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + virtual void freezeSystemMatrix() = delete; + +protected: + + SOFA_CORE_ATTRIBUTE_REMOVE_ASSEMBLY_API() + DeprecatedAndRemoved frozen; }; } // namespace sofa::core::behavior 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 From dae0bbc59a222c8566d306b8ee245e96c476f39b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 10 Oct 2025 08:57:52 +0200 Subject: [PATCH 41/48] forgot a file --- .../src/sofa/component/linearsolver/iterative/config.h.in | 7 +++++++ 1 file changed, 7 insertions(+) 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 c70503f0ffe..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,6 +37,13 @@ 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 From e4d0d9bebe23f56fa5d0b082205a38b73b61de1c Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 13 Oct 2025 15:15:31 +0200 Subject: [PATCH 42/48] simplify `checkLinearSystem` logic using `doCheckLinearSystem` --- .../preconditioner/PrecomputedWarpPreconditioner.inl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) 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 26bc4a3d1d7..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,14 +67,7 @@ PrecomputedWarpPreconditioner::PrecomputedWarpPreconditioner() template void PrecomputedWarpPreconditioner::checkLinearSystem() { - if (!this->l_linearSystem) - { - auto* matrixLinearSystem=this->getContext()->template get >(); - if(!matrixLinearSystem) - { - this->template createDefaultLinearSystem >(); - } - } + this->template doCheckLinearSystem>(); } //Solve x = R * M^-1 * R^t * b From 43c5df158db9ce777bf51900bd6268adb9717a7f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 13 Oct 2025 15:16:16 +0200 Subject: [PATCH 43/48] add components that are created implicitly in the scene --- .../Component/LinearSolver/FEMBAR-common.xml | 1 + .../FEMBAR_PCG_PrecomputedWarpPreconditioner.scn | 16 +++++++++++++--- 2 files changed, 14 insertions(+), 3 deletions(-) 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/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 @@ - - - + + + + + + + + + + + + + From 9f2932cf5bc06f3abb41e3c18dee54580e1454ac Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 13 Oct 2025 15:47:32 +0200 Subject: [PATCH 44/48] add a description to a class member --- .../component/linearsolver/preconditioner/RotationMatrixSystem.h | 1 + 1 file changed, 1 insertion(+) 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 c8958f8d143..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 @@ -48,6 +48,7 @@ class RotationMatrixSystem : public linearsystem::TypedMatrixLinearSystem Date: Mon, 13 Oct 2025 15:47:43 +0200 Subject: [PATCH 45/48] see if regression test is happy --- .../Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn index 51b2a59147e..fe9a78ee749 100644 --- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn +++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn @@ -7,7 +7,7 @@ - + From e81b8ae529708e0483cd6c83c228e69510351fd8 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 16 Oct 2025 07:58:50 +0200 Subject: [PATCH 46/48] rename `d_authorizeAssembly` to `d_enableAssembly` --- .../linearsolver/direct/AsyncSparseLDLSolver.h | 2 +- .../linearsolver/direct/AsyncSparseLDLSolver.inl | 10 +++++----- .../linearsolver/direct/PrecomputedLinearSolver.inl | 2 +- .../src/sofa/core/behavior/BaseMatrixLinearSystem.cpp | 8 ++++---- .../src/sofa/core/behavior/BaseMatrixLinearSystem.h | 2 +- 5 files changed, 12 insertions(+), 12 deletions(-) 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 efa7629b642..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 @@ -97,7 +97,7 @@ class AsyncSparseLDLSolver : public SparseLDLSolver newInvertDataReady { false }; - Data< bool > d_authorizeAssembly; + 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 ef0649739a9..7865cd68b01 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 @@ -30,7 +30,7 @@ namespace sofa::component::linearsolver::direct template AsyncSparseLDLSolver::AsyncSparseLDLSolver() - : d_authorizeAssembly(initData(&d_authorizeAssembly, true, "authorizeAssembly", "Allow assembly of the linear system")) + : d_enableAssembly(initData(&d_enableAssembly, true, "authorizeAssembly", "Allow assembly of the linear system")) { } @@ -43,7 +43,7 @@ void AsyncSparseLDLSolver::init() { if (this->l_linearSystem) { - this->l_linearSystem->d_authorizeAssembly.setParent(&d_authorizeAssembly); + this->l_linearSystem->d_enableAssembly.setParent(&d_enableAssembly); } waitForAsyncTask = true; @@ -55,7 +55,7 @@ void AsyncSparseLDLSolver::init() template void AsyncSparseLDLSolver::reset() { - d_authorizeAssembly.setValue(true); + d_enableAssembly.setValue(true); waitForAsyncTask = true; } @@ -80,7 +80,7 @@ void AsyncSparseLDLSolver::solveSystem() this->d_factorizationInvalidation.setValue(false); //matrix assembly is temporarily stopped until the next factorization - d_authorizeAssembly.setValue(false); + d_enableAssembly.setValue(false); } if (waitForAsyncTask) @@ -159,7 +159,7 @@ void AsyncSparseLDLSolver::asyncFactorization( newInvertDataReady = true; //factorization is finished: matrix assembly is authorized once again - d_authorizeAssembly.setValue(true); + d_enableAssembly.setValue(true); } template 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 ada8a553931..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 @@ -146,7 +146,7 @@ void PrecomputedLinearSolver::invert(TMatrix& /*M*/) if (this->l_linearSystem) { loadMatrix(*this->l_linearSystem->getSystemMatrix()); - this->l_linearSystem->d_authorizeAssembly.setValue(false); + this->l_linearSystem->d_enableAssembly.setValue(false); } first = false; } diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp index f42dacd8b6f..630fcd44ac1 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp @@ -27,17 +27,17 @@ namespace sofa::core::behavior BaseMatrixLinearSystem::BaseMatrixLinearSystem() : Inherit1() , d_matrixSize(initData(&d_matrixSize, "matrixSize", "Size of the global matrix")) -, d_authorizeAssembly(initData(&d_authorizeAssembly, true, "authorizeAssembly", "Allows to assemble the system matrix")) +, d_enableAssembly(initData(&d_enableAssembly, true, "authorizeAssembly", "Allows to assemble the system matrix")) { d_matrixSize.setReadOnly(true); - d_authorizeAssembly.setReadOnly(true); - d_authorizeAssembly.setDisplayed(false); + d_enableAssembly.setReadOnly(true); + d_enableAssembly.setDisplayed(false); } void BaseMatrixLinearSystem::buildSystemMatrix(const core::MechanicalParams* mparams) { - if (d_authorizeAssembly.getValue()) + if (d_enableAssembly.getValue()) { preAssembleSystem(mparams); assembleSystem(mparams); diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h index 2b9ab50cd50..47462077153 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.h @@ -46,7 +46,7 @@ class SOFA_CORE_API BaseMatrixLinearSystem : public virtual core::objectmodel::B /// Size of the linear system Data< sofa::type::Vec2u > d_matrixSize; - Data< bool > d_authorizeAssembly; + Data< bool > d_enableAssembly; /// Returns the system matrix as a sofa::linearalgebra::BaseMatrix* virtual linearalgebra::BaseMatrix* getSystemBaseMatrix() const { return nullptr; } From 446edf939287d91dee7bfcb91d5bdfd95be75e7b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 16 Oct 2025 10:06:30 +0200 Subject: [PATCH 47/48] rename Data name --- .../sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl | 2 +- .../Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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 7865cd68b01..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 @@ -30,7 +30,7 @@ namespace sofa::component::linearsolver::direct template AsyncSparseLDLSolver::AsyncSparseLDLSolver() - : d_enableAssembly(initData(&d_enableAssembly, true, "authorizeAssembly", "Allow assembly of the linear system")) + : d_enableAssembly(initData(&d_enableAssembly, true, "enableAssembly", "Allow assembly of the linear system")) { } diff --git a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp index 630fcd44ac1..5a341c1277a 100644 --- a/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp +++ b/Sofa/framework/Core/src/sofa/core/behavior/BaseMatrixLinearSystem.cpp @@ -27,7 +27,7 @@ namespace sofa::core::behavior BaseMatrixLinearSystem::BaseMatrixLinearSystem() : Inherit1() , d_matrixSize(initData(&d_matrixSize, "matrixSize", "Size of the global matrix")) -, d_enableAssembly(initData(&d_enableAssembly, true, "authorizeAssembly", "Allows to assemble the system matrix")) +, d_enableAssembly(initData(&d_enableAssembly, true, "enableAssembly", "Allows to assemble the system matrix")) { d_matrixSize.setReadOnly(true); From cd7c3d36a6a51262baa604f9b387b3fc3f16ae31 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 16 Oct 2025 15:39:46 +0200 Subject: [PATCH 48/48] fix call to MechanicalMultiVectorPeqBaseVectorVisitor but actually it is not called --- .../component/linearsolver/iterative/MatrixLinearSolver.inl | 2 +- .../MechanicalMultiVectorPeqBaseVectorVisitor.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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 e1502d3413f..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 @@ -519,7 +519,7 @@ void MatrixLinearSolver::computeResidual(const core::ExecParams* sofa::core::behavior::MultiVecDeriv force(&vop, core::vec_id::write_access::force ); // force += rhs - // executeVisitor( MechanicalMultiVectorPeqBaseVectorVisitor(core::execparams::defaultInstance(), force, rhsVector, &(linearSystem.matrixAccessor)) ); + executeVisitor( MechanicalMultiVectorPeqBaseVectorVisitor(core::execparams::defaultInstance(), force, rhsVector, nullptr) ); } 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; }