[Diffusion] Explicit MechanicalObject link

Sofa/Component/Diffusion/src/sofa/component/diffusion/TetrahedronDiffusionFEMForceField.h

@@ -114,6 +114,7 @@ class TetrahedronDiffusionFEMForceField : public core::behavior::ForceField<Data
 
       /// Link to be set to the topology container in the component graph. 
       SingleLink<TetrahedronDiffusionFEMForceField<DataTypes>, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology;
+      SingleLink<TetrahedronDiffusionFEMForceField<DataTypes>, MechObject, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_mecaObj;
 
 
 protected:

Sofa/Component/Diffusion/src/sofa/component/diffusion/TetrahedronDiffusionFEMForceField.inl

@@ -128,6 +128,7 @@ TetrahedronDiffusionFEMForceField<DataTypes>::TetrahedronDiffusionFEMForceField(
       d_tagMeshMechanics(initData(&d_tagMeshMechanics, std::string("meca"),"tagMechanics","Tag of the Mechanical Object.")),
       d_drawConduc( initData(&d_drawConduc, (bool)false, "drawConduc","To display conductivity map."))
     , l_topology(initLink("topology", "link to the topology container"))
+    , l_mecaObj(initLink("mecaObj", "link to the 3D mechanical MechanicalObject"))
     , m_topology(nullptr)
 {
     this->f_listening.setValue(true);
@@ -198,12 +199,29 @@ void TetrahedronDiffusionFEMForceField<DataTypes>::init()
     }
 
     /// Get the mechanical object containing the mesh position in 3D
-    core::objectmodel::Tag mechanicalTag(d_tagMeshMechanics.getValue());
-    this->getContext()->get(mechanicalObject, mechanicalTag,sofa::core::objectmodel::BaseContext::SearchUp);
-    if (mechanicalObject==nullptr)
+    /// Priority is the explicit link, with a fallback to historical tag-based lookup.
+    if (!l_mecaObj.empty())
     {
-        msg_error() << "cannot find the mechanical object named '" << mechanicalObject << msgendl;
-        return;
+        mechanicalObject = l_mecaObj.get();
+        if (mechanicalObject == nullptr)
+        {
+            msg_error() << "No MechanicalState found at mecaObj link path '" << l_mecaObj.getLinkedPath() << "'.";
+            this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
+            return;
+        }
+    }
+    else
+    {
+        core::objectmodel::Tag mechanicalTag(d_tagMeshMechanics.getValue());
+        this->getContext()->get(mechanicalObject, mechanicalTag, sofa::core::objectmodel::BaseContext::SearchUp);
+        if (mechanicalObject == nullptr)
+        {
+            msg_error() << "No MechanicalState found. Set 'mecaObj' link, or provide a valid 'tagMechanics' (current value: '"
+                        << d_tagMeshMechanics.getValue() << "').";
+            this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
+            return;
+        }
+        msg_info() << "Using legacy tagMechanics lookup to find MechanicalState with tag '" << d_tagMeshMechanics.getValue() << "'.";
     }
 
     if(d_transverseAnisotropyRatio.getValue()!=1.0)

Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp

@@ -24,16 +24,241 @@
 #include <sofa/component/linearsolver/iterative/MatrixLinearSolver.inl>
 #include <sofa/core/ObjectFactory.h>
 #include <sofa/component/linearsystem/MatrixFreeSystem.h>
+#include <sofa/component/linearsystem/MatrixLinearSystem.h>
+#include <sofa/linearalgebra/CompressedRowSparseMatrix.h>
+#include <sofa/linearalgebra/FullVector.h>
 
 namespace sofa::component::linearsolver::iterative
 {
 
+using CRSd  = linearalgebra::CompressedRowSparseMatrix<SReal>;
+using CRS3d = linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>;
+using FVd   = linearalgebra::FullVector<SReal>;
+
+// Preconditioned CG on an assembled CRS matrix (both CRS<SReal> and CRS<Mat3x3>).
+//
+// Mirrors the generic PCGLinearSolver::solve() in the .inl, but drives the
+// preconditioner directly on the assembled matrix (precond->invert(M) once, then
+// precond->solve(M, z, r) per iteration) instead of the GraphScattered
+// setRHS/solveSystem/dispatch plumbing, which does not exist for a CRS
+// preconditioner such as SSORPreconditioner. The stopping test follows the
+// generic version: r_norm = <r, M^{-1} r> compared against tolerance * <b, b>.
+//
+// The body is inlined into each specialisation of solve() below: it relies on the
+// protected members (cgstep_*, d_graph, newton_iter, TempVectorContainer) of the
+// solver, so it cannot live in a free helper, and a member template would leak
+// into the GraphScattered specialisation which must keep the .inl implementation.
+
+// ---------------------------------------------------------------------------
+// CompressedRowSparseMatrix<SReal> specialisations
+// ---------------------------------------------------------------------------
+template<>
+void PCGLinearSolver<CRSd, FVd>::checkLinearSystem()
+{
+    this->template doCheckLinearSystem<
+        sofa::component::linearsystem::MatrixLinearSystem<CRSd, FVd>>();
+}
+
+template<>
+void PCGLinearSolver<CRSd, FVd>::ensureRequiredLinearSystemType()
+{
+    // MatrixLinearSystem<CRSd,FVd> is the expected assembled system — nothing to reject.
+}
+
+template<>
+void PCGLinearSolver<CRSd, FVd>::bwdInit()
+{
+    if (this->isComponentStateInvalid())
+        return;
+    // No PreconditionedMatrixFreeSystem linking is needed for the CRS pipeline.
+}
+
+template<>
+void PCGLinearSolver<CRSd, FVd>::solve(CRSd& M, FVd& x, FVd& b)
+{
+    using MLSolver = MatrixLinearSolver<CRSd, FVd>;
+
+    // Resolve the preconditioner (may be null -> plain CG).
+    MLSolver* precond = nullptr;
+    if (l_preconditioner.get() != nullptr && d_use_precond.getValue())
+    {
+        precond = dynamic_cast<MLSolver*>(l_preconditioner.get());
+        if (!precond)
+            msg_warning() << "Preconditioner '" << l_preconditioner->getName()
+                          << "' is not a MatrixLinearSolver<CRS,FullVector> "
+                          << "— running unpreconditioned.";
+        else
+            precond->invert(M); // factorise/refresh against the current matrix
+    }
+
+    std::map<std::string, sofa::type::vector<SReal>>& graph = *d_graph.beginEdit();
+    newton_iter++;
+    char name[256];
+    snprintf(name, sizeof(name), "Error %d", newton_iter);
+    sofa::type::vector<SReal>& graph_error = graph[std::string(name)];
+
+    const core::ExecParams* params = core::execparams::defaultInstance();
+    typename Inherit::TempVectorContainer vtmp(this, params, M, x, b);
+    FVd& r = *vtmp.createTempVector();
+    FVd& w = *vtmp.createTempVector();
+    FVd& s = *vtmp.createTempVector();
+    // createTempVector() default-constructs an empty FullVector; size it to the
+    // system before use, otherwise the preconditioner writes out of bounds.
+    r.resize(M.rowSize());
+    w.resize(M.rowSize());
+    s.resize(M.rowSize());
+
+    const SReal b_norm = b.dot(b);
+    const SReal tol    = d_tolerance.getValue() * b_norm;
+
+    r = M * x;
+    cgstep_beta(r, b, SReal(-1));   // r = b - M*x
+
+    if (precond) precond->solve(M, w, r);   // w = M^{-1} r
+    else         w = r;
+
+    SReal r_norm = r.dot(w);        // r . M^{-1} r
+    graph_error.push_back(r_norm / b_norm);
+
+    unsigned iter = 1;
+    while ((iter <= d_maxIter.getValue()) && (r_norm > tol))
+    {
+        s = M * w;
+        const SReal dtq   = w.dot(s);
+        const SReal alpha = r_norm / dtq;
+        cgstep_alpha(x, w,  alpha);   // x += alpha w
+        cgstep_alpha(r, s, -alpha);   // r -= alpha (M w)
+
+        if (precond) precond->solve(M, s, r);
+        else         s = r;
+
+        const SReal deltaOld = r_norm;
+        r_norm = r.dot(s);
+        graph_error.push_back(r_norm / b_norm);
+        const SReal beta = r_norm / deltaOld;
+        cgstep_beta(w, s, beta);      // w = s + beta w
+        iter++;
+    }
+
+    d_graph.endEdit();
+    vtmp.deleteTempVector(&r);
+    vtmp.deleteTempVector(&w);
+    vtmp.deleteTempVector(&s);
+    sofa::helper::AdvancedTimer::valSet("PCG iterations", iter);
+}
+
+// ---------------------------------------------------------------------------
+// CompressedRowSparseMatrix<Mat<3,3>> specialisations
+// ---------------------------------------------------------------------------
+template<>
+void PCGLinearSolver<CRS3d, FVd>::checkLinearSystem()
+{
+    this->template doCheckLinearSystem<
+        sofa::component::linearsystem::MatrixLinearSystem<CRS3d, FVd>>();
+}
+
+template<>
+void PCGLinearSolver<CRS3d, FVd>::ensureRequiredLinearSystemType()
+{
+    // MatrixLinearSystem<CRS3d,FVd> is the expected assembled system — nothing to reject.
+}
+
+template<>
+void PCGLinearSolver<CRS3d, FVd>::bwdInit()
+{
+    if (this->isComponentStateInvalid())
+        return;
+    // No PreconditionedMatrixFreeSystem linking is needed for the CRS pipeline.
+}
+
+template<>
+void PCGLinearSolver<CRS3d, FVd>::solve(CRS3d& M, FVd& x, FVd& b)
+{
+    using MLSolver = MatrixLinearSolver<CRS3d, FVd>;
+
+    // Resolve the preconditioner (may be null -> plain CG).
+    MLSolver* precond = nullptr;
+    if (l_preconditioner.get() != nullptr && d_use_precond.getValue())
+    {
+        precond = dynamic_cast<MLSolver*>(l_preconditioner.get());
+        if (!precond)
+            msg_warning() << "Preconditioner '" << l_preconditioner->getName()
+                          << "' is not a MatrixLinearSolver<CRS3,FullVector> "
+                          << "— running unpreconditioned.";
+        else
+            precond->invert(M); // factorise/refresh against the current matrix
+    }
+
+    std::map<std::string, sofa::type::vector<SReal>>& graph = *d_graph.beginEdit();
+    newton_iter++;
+    char name[256];
+    snprintf(name, sizeof(name), "Error %d", newton_iter);
+    sofa::type::vector<SReal>& graph_error = graph[std::string(name)];
+
+    const core::ExecParams* params = core::execparams::defaultInstance();
+    typename Inherit::TempVectorContainer vtmp(this, params, M, x, b);
+    FVd& r = *vtmp.createTempVector();
+    FVd& w = *vtmp.createTempVector();
+    FVd& s = *vtmp.createTempVector();
+    // createTempVector() default-constructs an empty FullVector; size it to the
+    // system before use, otherwise the preconditioner writes out of bounds.
+    r.resize(M.rowSize());
+    w.resize(M.rowSize());
+    s.resize(M.rowSize());
+
+    const SReal b_norm = b.dot(b);
+    const SReal tol    = d_tolerance.getValue() * b_norm;
+
+    r = M * x;
+    cgstep_beta(r, b, SReal(-1));   // r = b - M*x
+
+    if (precond) precond->solve(M, w, r);   // w = M^{-1} r
+    else         w = r;
+
+    SReal r_norm = r.dot(w);        // r . M^{-1} r
+    graph_error.push_back(r_norm / b_norm);
+
+    unsigned iter = 1;
+    while ((iter <= d_maxIter.getValue()) && (r_norm > tol))
+    {
+        s = M * w;
+        const SReal dtq   = w.dot(s);
+        const SReal alpha = r_norm / dtq;
+        cgstep_alpha(x, w,  alpha);   // x += alpha w
+        cgstep_alpha(r, s, -alpha);   // r -= alpha (M w)
+
+        if (precond) precond->solve(M, s, r);
+        else         s = r;
+
+        const SReal deltaOld = r_norm;
+        r_norm = r.dot(s);
+        graph_error.push_back(r_norm / b_norm);
+        const SReal beta = r_norm / deltaOld;
+        cgstep_beta(w, s, beta);      // w = s + beta w
+        iter++;
+    }
+
+    d_graph.endEdit();
+    vtmp.deleteTempVector(&r);
+    vtmp.deleteTempVector(&w);
+    vtmp.deleteTempVector(&s);
+    sofa::helper::AdvancedTimer::valSet("PCG iterations", iter);
+}
+
+// ---------------------------------------------------------------------------
+// Factory registration + explicit instantiations
+// ---------------------------------------------------------------------------
 void registerPCGLinearSolver(sofa::core::ObjectFactory* factory)
 {
-    factory->registerObjects(core::ObjectRegistrationData("Linear solver using the preconditioned conjugate gradient iterative algorithm.")
-        .add< PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >());
+    factory->registerObjects(core::ObjectRegistrationData(
+        "Linear solver using the preconditioned conjugate gradient iterative algorithm.")
+        .add< PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >()
+        .add< PCGLinearSolver<CRSd, FVd> >()
+        .add< PCGLinearSolver<CRS3d, FVd> >());
 }
 
 template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
+template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<CRSd, FVd>;
+template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<CRS3d, FVd>;
 
 } // namespace sofa::component::linearsolver::iterative

Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h

@@ -24,6 +24,9 @@
 
 #include <sofa/core/behavior/LinearSolver.h>
 #include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
+#include <sofa/linearalgebra/CompressedRowSparseMatrix.h>
+#include <sofa/linearalgebra/FullVector.h>
+#include <sofa/type/Mat.h>
 #include <sofa/helper/map.h>
 
 #include <cmath>
@@ -100,8 +103,80 @@ inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,compon
 template<>
 inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_alpha(Vector& x,Vector& p, Real alpha);
 
+// -----------------------------------------------------------------------------
+// CompressedRowSparseMatrix (CRS) specialisations.
+//
+// The generic (GraphScattered) implementation in the .inl targets a matrix-free
+// system and is incompatible with an assembled CRS system on three points, so
+// each is overridden in the .cpp:
+//   * checkLinearSystem()              -> the generic one auto-creates a
+//                                         PreconditionedMatrixFreeSystem<GraphScattered>;
+//                                         CRS needs a plain MatrixLinearSystem<CRS,FV>.
+//   * ensureRequiredLinearSystemType() -> the generic one rejects anything that is
+//                                         not a PreconditionedMatrixFreeSystem.
+//   * solve()                          -> the generic precond path drives the
+//                                         preconditioner through the GraphScattered
+//                                         setRHS/solveSystem/dispatch plumbing; for a
+//                                         CRS preconditioner (e.g. SSORPreconditioner)
+//                                         we call precond->invert(M)/precond->solve()
+//                                         directly on the assembled matrix.
+// This makes PCGLinearSolver usable with an assembled (CRS) linear system, hence
+// with LinearSolverConstraintCorrection, which requires one.
+// -----------------------------------------------------------------------------
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<SReal>,
+    linearalgebra::FullVector<SReal>>::solve(
+        linearalgebra::CompressedRowSparseMatrix<SReal>&,
+        linearalgebra::FullVector<SReal>&,
+        linearalgebra::FullVector<SReal>&);
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<SReal>,
+    linearalgebra::FullVector<SReal>>::checkLinearSystem();
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<SReal>,
+    linearalgebra::FullVector<SReal>>::ensureRequiredLinearSystemType();
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<SReal>,
+    linearalgebra::FullVector<SReal>>::bwdInit();
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>,
+    linearalgebra::FullVector<SReal>>::solve(
+        linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>&,
+        linearalgebra::FullVector<SReal>&,
+        linearalgebra::FullVector<SReal>&);
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>,
+    linearalgebra::FullVector<SReal>>::checkLinearSystem();
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>,
+    linearalgebra::FullVector<SReal>>::ensureRequiredLinearSystemType();
+
+template<>
+void PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>,
+    linearalgebra::FullVector<SReal>>::bwdInit();
+
 #if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)
-template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
+extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
+extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<SReal>,
+    linearalgebra::FullVector<SReal>>;
+extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<
+    linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal>>,
+    linearalgebra::FullVector<SReal>>;
 #endif // !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)
 
 } // namespace sofa::component::linearsolver::iterative