diff --git a/.github/workflows/nix.yml b/.github/workflows/nix.yml
new file mode 100644
index 00000000000..a1bff8fcec3
--- /dev/null
+++ b/.github/workflows/nix.yml
@@ -0,0 +1,21 @@
+name: "CI - Nix"
+
+on:
+  push:
+
+jobs:
+  nix:
+    runs-on: "${{ matrix.os }}-latest"
+    if: ${{ github.repository_owner == 'sofa-framework' }}
+    strategy:
+      matrix:
+        os: [ubuntu, macos]
+    steps:
+      - uses: actions/checkout@v4
+      - uses: cachix/install-nix-action@v27
+      # TODO: the "sofa" account on cachix does not exist yet
+      #- uses: cachix/cachix-action@v15
+        #with:
+          #name: sofa
+          #authToken: '${{ secrets.CACHIX_AUTH_TOKEN }}'
+      - run: nix build -L
diff --git a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.cpp b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.cpp
index e6932abfdae..5f694d58dd0 100644
--- a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.cpp
+++ b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.cpp
@@ -37,33 +37,37 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
 {
 public:
 
+    template<class T>
+    using SubsetIndices = typename BilateralLagrangianConstraint<T>::SubsetIndices;
+
+
     template<class T>
     static void bwdInit(BilateralLagrangianConstraint<T>& self) {
         if (!self.d_keepOrientDiff.getValue())
             return;
 
-        helper::WriteAccessor<Data<typename BilateralLagrangianConstraint<T>::VecDeriv > > wrest = self.d_restVector;
+        auto wrest = sofa::helper::getWriteAccessor(self.d_restVector);
 
-        if (wrest.size() > 0) {
-            msg_warning("BilateralLagrangianConstraintSpecialization") << "keepOrientationDifference is activated, rest_vector will be ignored! " ;
+        if (!wrest.empty()) {
+            msg_warning(&self) << "keepOrientationDifference is activated, rest_vector will be ignored! " ;
             wrest.resize(0);
         }
 
-        const typename BilateralLagrangianConstraint<T>::SubsetIndices& m1Indices = self.d_m1.getValue();
-        const typename BilateralLagrangianConstraint<T>::SubsetIndices& m2Indices = self.d_m2.getValue();
+        const SubsetIndices<T>& m1Indices = self.d_m1.getValue();
+        const SubsetIndices<T>& m2Indices = self.d_m2.getValue();
 
         const unsigned minp = std::min(m1Indices.size(),m2Indices.size());
 
-        const typename BilateralLagrangianConstraint<T>::DataVecCoord &d_x1 = *self.mstate1->read(core::ConstVecCoordId::position());
-        const typename BilateralLagrangianConstraint<T>::DataVecCoord &d_x2 = *self.mstate2->read(core::ConstVecCoordId::position());
+        const DataVecCoord_t<T> &d_x1 = *self.mstate1->read(core::ConstVecCoordId::position());
+        const DataVecCoord_t<T> &d_x2 = *self.mstate2->read(core::ConstVecCoordId::position());
 
-        const typename BilateralLagrangianConstraint<T>::VecCoord &x1 = d_x1.getValue();
-        const typename BilateralLagrangianConstraint<T>::VecCoord &x2 = d_x2.getValue();
+        const VecCoord_t<T> &x1 = d_x1.getValue();
+        const VecCoord_t<T> &x2 = d_x2.getValue();
 
         for (unsigned pid=0; pid<minp; pid++)
         {
-            const typename BilateralLagrangianConstraint<T>::Coord P = x1[m1Indices[pid]];
-            const typename BilateralLagrangianConstraint<T>::Coord Q = x2[m2Indices[pid]];
+            const Coord_t<T> P = x1[m1Indices[pid]];
+            const Coord_t<T> Q = x2[m2Indices[pid]];
 
             type::Quat<SReal> qP, qQ, dQP;
             qP = P.getOrientation();
@@ -73,7 +77,7 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
             dQP = qP.quatDiff(qQ, qP);
             dQP.normalize();
 
-            typename BilateralLagrangianConstraint<T>::Coord df;
+            Coord_t<T> df;
             df.getCenter() = Q.getCenter() - P.getCenter();
             df.getOrientation() = dQP;
             self.initialDifference.push_back(df);
@@ -86,12 +90,12 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
     static void getConstraintResolution(BilateralLagrangianConstraint<T>& self,
                                         const ConstraintParams* cParams,
                                         std::vector<ConstraintResolution*>& resTab,
-                                        unsigned int& offset, double tolerance)
+                                        unsigned int& offset, SReal tolerance)
     {
         SOFA_UNUSED(cParams);
-        const unsigned minp=std::min(self.d_m1.getValue().size(),
-                                     self.d_m2.getValue().size());
-        for (unsigned pid=0; pid<minp; pid++)
+        const unsigned minp = std::min(self.d_m1.getValue().size(),
+                                       self.d_m2.getValue().size());
+        for (unsigned pid = 0; pid < minp; pid++)
         {
             resTab[offset] = new BilateralConstraintResolution3Dof();
             offset += 3;
@@ -106,26 +110,26 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
     template <class T>
     static void buildConstraintMatrix(BilateralLagrangianConstraint<T>& self,
                                       const ConstraintParams* cParams,
-                                      typename BilateralLagrangianConstraint<T>::DataMatrixDeriv &c1_d,
-                                      typename BilateralLagrangianConstraint<T>::DataMatrixDeriv &c2_d,
+                                      DataMatrixDeriv_t<T> &c1_d,
+                                      DataMatrixDeriv_t<T> &c2_d,
                                       unsigned int &constraintId,
-                                      const typename BilateralLagrangianConstraint<T>::DataVecCoord &/*x1*/,
-                                      const typename BilateralLagrangianConstraint<T>::DataVecCoord &/*x2*/)
+                                      const DataVecCoord_t<T> &/*x1*/,
+                                      const DataVecCoord_t<T> &/*x2*/)
     {
         SOFA_UNUSED(cParams) ;
-        const typename BilateralLagrangianConstraint<T>::SubsetIndices& m1Indices = self.d_m1.getValue();
-        const typename BilateralLagrangianConstraint<T>::SubsetIndices& m2Indices = self.d_m2.getValue();
+        const SubsetIndices<T>& m1Indices = self.d_m1.getValue();
+        const SubsetIndices<T>& m2Indices = self.d_m2.getValue();
 
         unsigned minp = std::min(m1Indices.size(),m2Indices.size());
         self.cid.resize(minp);
 
-        typename BilateralLagrangianConstraint<T>::MatrixDeriv &c1 = *c1_d.beginEdit();
-        typename BilateralLagrangianConstraint<T>::MatrixDeriv &c2 = *c2_d.beginEdit();
+        auto c1 = sofa::helper::getWriteAccessor(c1_d);
+        auto c2 = sofa::helper::getWriteAccessor(c2_d);
 
-        const Vec<3, typename BilateralLagrangianConstraint<T>::Real> cx(1,0,0), cy(0,1,0), cz(0,0,1);
-        const Vec<3, typename BilateralLagrangianConstraint<T>::Real> vZero(0,0,0);
+        static constexpr Vec<3, Real_t<T>> cx(1,0,0), cy(0,1,0), cz(0,0,1);
+        static constexpr Vec<3, Real_t<T>> vZero(0,0,0);
 
-        for (unsigned pid=0; pid<minp; pid++)
+        for (unsigned pid = 0; pid < minp; ++pid)
         {
             int tm1 = m1Indices[pid];
             int tm2 = m2Indices[pid];
@@ -134,46 +138,43 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
             constraintId += 6;
 
             //Apply constraint for position
-            typename BilateralLagrangianConstraint<T>::MatrixDerivRowIterator c1_it = c1.writeLine(self.cid[pid]);
-            c1_it.addCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(-cx, vZero));
+            auto c1_it = c1->writeLine(self.cid[pid]);
+            c1_it.addCol(tm1, Deriv_t<T>(-cx, vZero));
 
-            typename BilateralLagrangianConstraint<T>::MatrixDerivRowIterator c2_it = c2.writeLine(self.cid[pid]);
-            c2_it.addCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(cx, vZero));
+            auto c2_it = c2->writeLine(self.cid[pid]);
+            c2_it.addCol(tm2, Deriv_t<T>(cx, vZero));
 
-            c1_it = c1.writeLine(self.cid[pid] + 1);
-            c1_it.setCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(-cy, vZero));
+            c1_it = c1->writeLine(self.cid[pid] + 1);
+            c1_it.setCol(tm1, Deriv_t<T>(-cy, vZero));
 
-            c2_it = c2.writeLine(self.cid[pid] + 1);
-            c2_it.setCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(cy, vZero));
+            c2_it = c2->writeLine(self.cid[pid] + 1);
+            c2_it.setCol(tm2, Deriv_t<T>(cy, vZero));
 
-            c1_it = c1.writeLine(self.cid[pid] + 2);
-            c1_it.setCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(-cz, vZero));
+            c1_it = c1->writeLine(self.cid[pid] + 2);
+            c1_it.setCol(tm1, Deriv_t<T>(-cz, vZero));
 
-            c2_it = c2.writeLine(self.cid[pid] + 2);
-            c2_it.setCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(cz, vZero));
+            c2_it = c2->writeLine(self.cid[pid] + 2);
+            c2_it.setCol(tm2, Deriv_t<T>(cz, vZero));
 
             //Apply constraint for orientation
-            c1_it = c1.writeLine(self.cid[pid] + 3);
-            c1_it.setCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(vZero, -cx));
+            c1_it = c1->writeLine(self.cid[pid] + 3);
+            c1_it.setCol(tm1, Deriv_t<T>(vZero, -cx));
 
-            c2_it = c2.writeLine(self.cid[pid] + 3);
-            c2_it.setCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(vZero, cx));
+            c2_it = c2->writeLine(self.cid[pid] + 3);
+            c2_it.setCol(tm2, Deriv_t<T>(vZero, cx));
 
-            c1_it = c1.writeLine(self.cid[pid] + 4);
-            c1_it.setCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(vZero, -cy));
+            c1_it = c1->writeLine(self.cid[pid] + 4);
+            c1_it.setCol(tm1, Deriv_t<T>(vZero, -cy));
 
-            c2_it = c2.writeLine(self.cid[pid] + 4);
-            c2_it.setCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(vZero, cy));
+            c2_it = c2->writeLine(self.cid[pid] + 4);
+            c2_it.setCol(tm2, Deriv_t<T>(vZero, cy));
 
-            c1_it = c1.writeLine(self.cid[pid] + 5);
-            c1_it.setCol(tm1, typename BilateralLagrangianConstraint<T>::Deriv(vZero, -cz));
+            c1_it = c1->writeLine(self.cid[pid] + 5);
+            c1_it.setCol(tm1, Deriv_t<T>(vZero, -cz));
 
-            c2_it = c2.writeLine(self.cid[pid] + 5);
-            c2_it.setCol(tm2, typename BilateralLagrangianConstraint<T>::Deriv(vZero, cz));
+            c2_it = c2->writeLine(self.cid[pid] + 5);
+            c2_it.setCol(tm2, Deriv_t<T>(vZero, cz));
         }
-
-        c1_d.endEdit();
-        c2_d.endEdit();
     }
 
 
@@ -190,41 +191,45 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
         const typename BilateralLagrangianConstraint<T>::SubsetIndices& m2Indices = self.d_m2.getValue();
 
         unsigned min = std::min(m1Indices.size(), m2Indices.size());
-        const  typename BilateralLagrangianConstraint<T>::VecDeriv& restVector = self.d_restVector.getValue();
-        self.dfree.resize(min);
+        const VecDeriv_t<T>& restVector = self.d_restVector.getValue();
+        self.m_violation.resize(min);
 
-        const  typename BilateralLagrangianConstraint<T>::VecCoord &x1 = d_x1.getValue();
-        const  typename BilateralLagrangianConstraint<T>::VecCoord &x2 = d_x2.getValue();
+        const VecCoord_t<T> &x1 = d_x1.getValue();
+        const VecCoord_t<T> &x2 = d_x2.getValue();
 
-        for (unsigned pid=0; pid<min; pid++)
+        for (unsigned pid = 0; pid < min; pid++)
         {
-            //typename BilateralLagrangianConstraint<T>::Coord dof1 = x1[m1Indices[pid]];
-            //typename BilateralLagrangianConstraint<T>::Coord dof2 = x2[m2Indices[pid]];
-            typename BilateralLagrangianConstraint<T>::Coord dof1;
+            Coord_t<T> dof1;
 
-             if (self.d_keepOrientDiff.getValue()) {
-                 const typename BilateralLagrangianConstraint<T>::Coord dof1c = x1[m1Indices[pid]];
+            if (self.d_keepOrientDiff.getValue())
+            {
+                const Coord_t<T> dof1c = x1[m1Indices[pid]];
 
-                 typename BilateralLagrangianConstraint<T>::Coord corr=self.initialDifference[pid];
-                 type::Quat<SReal> df = corr.getOrientation();
-                 type::Quat<SReal> o1 = dof1c.getOrientation();
-                 type::Quat<SReal> ro1 = o1 * df;
+                Coord_t<T> corr = self.initialDifference[pid];
+                type::Quat<Real_t<T>> df = corr.getOrientation();
+                type::Quat<Real_t<T>> o1 = dof1c.getOrientation();
+                type::Quat<Real_t<T>> ro1 = o1 * df;
 
-                 dof1.getCenter() = dof1c.getCenter() + corr.getCenter();
-                 dof1.getOrientation() = ro1;
-             } else
-                 dof1 = x1[m1Indices[pid]];
+                dof1.getCenter() = dof1c.getCenter() + corr.getCenter();
+                dof1.getOrientation() = ro1;
+            }
+            else
+            {
+                dof1 = x1[m1Indices[pid]];
+            }
 
-            const typename BilateralLagrangianConstraint<T>::Coord dof2 = x2[m2Indices[pid]];
+            const Coord_t<T> dof2 = x2[m2Indices[pid]];
 
-            getVCenter(self.dfree[pid]) = dof2.getCenter() - dof1.getCenter();
-            getVOrientation(self.dfree[pid]) =  dof1.rotate(self.q.angularDisplacement(dof2.getOrientation() ,
+            getVCenter(self.m_violation[pid]) = dof2.getCenter() - dof1.getCenter();
+            getVOrientation(self.m_violation[pid]) =  dof1.rotate(self.q.angularDisplacement(dof2.getOrientation() ,
                                                                                   dof1.getOrientation())); // angularDisplacement compute the rotation vector btw the two quaternions
             if (pid < restVector.size())
-                self.dfree[pid] -= restVector[pid];
+                self.m_violation[pid] -= restVector[pid];
 
-            for (unsigned int i=0 ; i<self.dfree[pid].size() ; i++)
-                v->set(self.cid[pid]+i, self.dfree[pid][i]);
+            for (unsigned int i = 0; i < self.m_violation[pid].size(); i++)
+            {
+                v->set(self.cid[pid]+i, self.m_violation[pid][i]);
+            }
         }
     }
 
@@ -235,9 +240,9 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
                            typename MyClass::Real /*contactDistance*/, int m1, int m2,
                            typename MyClass::Coord /*Pfree*/, typename MyClass::Coord /*Qfree*/, long /*id*/, typename MyClass::PersistentID /*localid*/)
     {
-        helper::WriteAccessor<Data<typename BilateralLagrangianConstraint<T>::SubsetIndices > > wm1 = self.d_m1;
-        helper::WriteAccessor<Data<typename BilateralLagrangianConstraint<T>::SubsetIndices > > wm2 = self.d_m2;
-        helper::WriteAccessor<Data<typename MyClass::VecDeriv > > wrest = self.d_restVector;
+        sofa::helper::WriteAccessor<Data<SubsetIndices<T>>> wm1 = self.d_m1;
+        sofa::helper::WriteAccessor<Data<SubsetIndices<T>>> wm2 = self.d_m2;
+        auto wrest = sofa::helper::WriteAccessor(self.d_restVector);
         wm1.push_back(m1);
         wm2.push_back(m2);
 
@@ -249,25 +254,30 @@ class BilateralLagrangianConstraintSpecialization<RigidImpl>
 
 };
 
+using RigidBilateralLagrangianConstraint = BilateralLagrangianConstraintSpecialization<RigidImpl>;
+
 
 template<> SOFA_COMPONENT_CONSTRAINT_LAGRANGIAN_MODEL_API
-void BilateralLagrangianConstraint<Rigid3Types>::init(){
+void BilateralLagrangianConstraint<Rigid3Types>::init()
+{
     unspecializedInit() ;
 }
 
 template<> SOFA_COMPONENT_CONSTRAINT_LAGRANGIAN_MODEL_API
-void BilateralLagrangianConstraint<Rigid3Types>::bwdInit() {
-    BilateralLagrangianConstraintSpecialization<RigidImpl>::bwdInit(*this);
+void BilateralLagrangianConstraint<Rigid3Types>::bwdInit()
+{
+    RigidBilateralLagrangianConstraint::bwdInit(*this);
 }
 
 template<> SOFA_COMPONENT_CONSTRAINT_LAGRANGIAN_MODEL_API
-void BilateralLagrangianConstraint<Rigid3Types>::getConstraintResolution(const ConstraintParams* cParams,
-                                                                           std::vector<ConstraintResolution*>& resTab,
-                                                                           unsigned int& offset)
+void BilateralLagrangianConstraint<Rigid3Types>::getConstraintResolution(
+    const ConstraintParams* cParams,
+    std::vector<ConstraintResolution*>& resTab,
+    unsigned int& offset)
 {
-    BilateralLagrangianConstraintSpecialization<RigidImpl>::getConstraintResolution(*this,
-                                                                                     cParams, resTab, offset,
-                                                                                     d_numericalTolerance.getValue()) ;
+    RigidBilateralLagrangianConstraint::getConstraintResolution(*this,
+        cParams, resTab, offset,
+        d_numericalTolerance.getValue());
 }
 
 template <> SOFA_COMPONENT_CONSTRAINT_LAGRANGIAN_MODEL_API
@@ -277,9 +287,9 @@ void BilateralLagrangianConstraint<Rigid3Types>::buildConstraintMatrix(const Con
                                                                          unsigned int &constraintId,
                                                                          const DataVecCoord &x1, const DataVecCoord &x2)
 {
-    BilateralLagrangianConstraintSpecialization<RigidImpl>::buildConstraintMatrix(*this,
-                                                                                   cParams, c1_d, c2_d, constraintId,
-                                                                                   x1, x2) ;
+    RigidBilateralLagrangianConstraint::buildConstraintMatrix(*this,
+        cParams, c1_d, c2_d, constraintId,
+        x1, x2);
 }
 
 
@@ -289,9 +299,9 @@ void BilateralLagrangianConstraint<Rigid3Types>::getConstraintViolation(const Co
                                                                           const DataVecCoord &d_x1, const DataVecCoord &d_x2,
                                                                           const DataVecDeriv &v1, const DataVecDeriv &v2)
 {
-    BilateralLagrangianConstraintSpecialization<RigidImpl>::getConstraintViolation(*this,
-                                                                                    cParams, v, d_x1, d_x2,
-                                                                                    v1, v2) ;
+    RigidBilateralLagrangianConstraint::getConstraintViolation(*this,
+        cParams, v, d_x1, d_x2,
+        v1, v2);
 }
 
 
@@ -312,9 +322,10 @@ void BilateralLagrangianConstraint<defaulttype::Rigid3Types>::addContact(Deriv n
                                                                            Coord Pfree, Coord Qfree,
                                                                            long id, PersistentID localid)
 {
-    BilateralLagrangianConstraintSpecialization<RigidImpl>::addContact(*this,
-                                                                        norm, P, Q, contactDistance, m1, m2, Pfree, Qfree,
-                                                                        id, localid) ;
+    RigidBilateralLagrangianConstraint::addContact(*this,
+                                                   norm, P, Q, contactDistance,
+                                                   m1, m2, Pfree, Qfree,
+                                                   id, localid);
 }
 
 void registerBilateralLagrangianConstraint(sofa::core::ObjectFactory* factory)
diff --git a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.h b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.h
index f986ad3321b..abc001ded95 100644
--- a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.h
+++ b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.h
@@ -45,7 +45,6 @@ namespace sofa::component::constraint::lagrangian::model
 using sofa::core::behavior::BaseConstraint ;
 using sofa::core::behavior::ConstraintResolution ;
 using sofa::core::behavior::PairInteractionConstraint ;
-using sofa::core::objectmodel::Data ;
 using sofa::core::ConstraintParams ;
 using sofa::core::ConstVecCoordId;
 
@@ -93,10 +92,11 @@ class BilateralLagrangianConstraint : public PairInteractionConstraint<DataTypes
     using DataSubsetIndices = sofa::core::topology::TopologySubsetIndices;
 
 protected:
-    std::vector<Deriv> dfree;
+    sofa::type::vector<Deriv> m_violation;
     Quat<SReal> q;
 
     std::vector<unsigned int> cid;
+
     SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_CONSTRAINT_LAGRANGIAN_MODEL()
     sofa::core::objectmodel::RenamedData<type::vector<Index> > m1;
 
@@ -114,7 +114,7 @@ class BilateralLagrangianConstraint : public PairInteractionConstraint<DataTypes
     Data<VecDeriv> d_restVector; ///< Relative position to maintain between attached points (optional)
     VecCoord initialDifference;
 
-    Data<double> d_numericalTolerance; ///< a real value specifying the tolerance during the constraint solving. (default=0.0001
+    Data<SReal> d_numericalTolerance; ///< a real value specifying the tolerance during the constraint solving. (default=0.0001
     Data<bool> d_activate; ///< control constraint activation (true by default)
     Data<bool> d_keepOrientDiff; ///< keep the initial difference in orientation (only for rigids)
 
@@ -128,7 +128,8 @@ class BilateralLagrangianConstraint : public PairInteractionConstraint<DataTypes
     BilateralLagrangianConstraint(MechanicalState* object) ;
     BilateralLagrangianConstraint();
 
-    virtual ~BilateralLagrangianConstraint(){}
+    ~BilateralLagrangianConstraint() override = default;
+
 public:
     void init() override;
 
@@ -137,22 +138,22 @@ class BilateralLagrangianConstraint : public PairInteractionConstraint<DataTypes
     void reinit() override;
 
     void buildConstraintMatrix(const ConstraintParams* cParams,
-                                       DataMatrixDeriv &c1, DataMatrixDeriv &c2,
-                                       unsigned int &cIndex,
-                                       const DataVecCoord &x1, const DataVecCoord &x2) override;
+                               DataMatrixDeriv& c1, DataMatrixDeriv& c2,
+                               unsigned int& cIndex,
+                               const DataVecCoord& x1, const DataVecCoord& x2) override;
 
     void getConstraintViolation(const ConstraintParams* cParams,
-                                        BaseVector *v,
-                                        const DataVecCoord &x1, const DataVecCoord &x2,
-                                        const DataVecDeriv &v1, const DataVecDeriv &v2) override;
+                                BaseVector* v,
+                                const DataVecCoord& x1, const DataVecCoord& x2,
+                                const DataVecDeriv& v1, const DataVecDeriv& v2) override;
 
     void getVelocityViolation(BaseVector *v,
                               const DataVecCoord &x1, const DataVecCoord &x2,
                               const DataVecDeriv &v1, const DataVecDeriv &v2);
 
     void getConstraintResolution(const ConstraintParams* cParams,
-                                         std::vector<ConstraintResolution*>& resTab,
-                                         unsigned int& offset) override;
+                                 std::vector<ConstraintResolution*>& resTab,
+                                 unsigned int& offset) override;
 
     void handleEvent(sofa::core::objectmodel::Event *event) override;
 
diff --git a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl
index ad18a6ed7f7..a486df45ecd 100644
--- a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl
+++ b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl
@@ -44,8 +44,8 @@ BilateralLagrangianConstraint<DataTypes>::BilateralLagrangianConstraint(Mechanic
     , d_m1(initData(&d_m1, "first_point","index of the constraint on the first model (object1)"))
     , d_m2(initData(&d_m2, "second_point","index of the constraint on the second model (object2)"))
     , d_restVector(initData(&d_restVector, "rest_vector","Relative position to maintain between attached points (optional)"))
-    , d_numericalTolerance(initData(&d_numericalTolerance, 0.0001, "numericalTolerance",
-                                    "a real value specifying the tolerance during the constraint solving. (optional, default=0.0001)") )
+    , d_numericalTolerance(initData(&d_numericalTolerance, 1e-4_sreal, "numericalTolerance",
+                                    "a real value specifying the tolerance during the constraint solving.") )
     , d_activate( initData(&d_activate, true, "activate", "control constraint activation (true by default)"))
     , d_keepOrientDiff(initData(&d_keepOrientDiff, false, "keepOrientationDifference", "keep the initial difference in orientation (only for rigids)"))
     , l_topology1(initLink("topology1", "link to the first topology container"))
@@ -81,6 +81,11 @@ void BilateralLagrangianConstraint<DataTypes>::unspecializedInit()
     assert(this->mstate1);
     assert(this->mstate2);
 
+    if (!this->mstate1 || !this->mstate2)
+    {
+        this->d_componentState.setValue(core::objectmodel::ComponentState::Invalid);
+    }
+
     prevForces.clear();
 }
 
@@ -120,8 +125,10 @@ void BilateralLagrangianConstraint<DataTypes>::reinit()
 
 
 template<class DataTypes>
-void BilateralLagrangianConstraint<DataTypes>::buildConstraintMatrix(const ConstraintParams*, DataMatrixDeriv &c1_d, DataMatrixDeriv &c2_d, unsigned int &constraintId
-                                                                      , const DataVecCoord &/*x1*/, const DataVecCoord &/*x2*/)
+void BilateralLagrangianConstraint<DataTypes>::buildConstraintMatrix(
+    const ConstraintParams*, DataMatrixDeriv& c1_d, DataMatrixDeriv& c2_d,
+    unsigned int& constraintId,
+    const DataVecCoord&/*x1*/, const DataVecCoord&/*x2*/)
 {
     if (!d_activate.getValue())
         return;
@@ -133,12 +140,12 @@ void BilateralLagrangianConstraint<DataTypes>::buildConstraintMatrix(const Const
     const SubsetIndices& m1Indices = d_m1.getValue();
     const SubsetIndices& m2Indices = d_m2.getValue();
 
-    MatrixDeriv &c1 = *c1_d.beginEdit();
-    MatrixDeriv &c2 = *c2_d.beginEdit();
+    auto c1 = sofa::helper::getWriteAccessor(c1_d);
+    auto c2 = sofa::helper::getWriteAccessor(c2_d);
 
     cid.resize(minp);
 
-    for (unsigned pid=0; pid<minp; pid++)
+    for (unsigned pid = 0; pid < minp; ++pid)
     {
         int tm1 = m1Indices[pid];
         int tm2 = m2Indices[pid];
@@ -148,35 +155,32 @@ void BilateralLagrangianConstraint<DataTypes>::buildConstraintMatrix(const Const
         cid[pid] = constraintId;
         constraintId += 3;
 
-        MatrixDerivRowIterator c1_it = c1.writeLine(cid[pid]);
+        MatrixDerivRowIterator c1_it = c1->writeLine(cid[pid]);
         c1_it.addCol(tm1, -cx);
 
-        MatrixDerivRowIterator c2_it = c2.writeLine(cid[pid]);
+        MatrixDerivRowIterator c2_it = c2->writeLine(cid[pid]);
         c2_it.addCol(tm2, cx);
 
-        c1_it = c1.writeLine(cid[pid] + 1);
+        c1_it = c1->writeLine(cid[pid] + 1);
         c1_it.setCol(tm1, -cy);
 
-        c2_it = c2.writeLine(cid[pid] + 1);
+        c2_it = c2->writeLine(cid[pid] + 1);
         c2_it.setCol(tm2, cy);
 
-        c1_it = c1.writeLine(cid[pid] + 2);
+        c1_it = c1->writeLine(cid[pid] + 2);
         c1_it.setCol(tm1, -cz);
 
-        c2_it = c2.writeLine(cid[pid] + 2);
+        c2_it = c2->writeLine(cid[pid] + 2);
         c2_it.setCol(tm2, cz);
     }
-
-    c1_d.endEdit();
-    c2_d.endEdit();
 }
 
 
 template<class DataTypes>
 void BilateralLagrangianConstraint<DataTypes>::getConstraintViolation(const ConstraintParams* cParams,
-                                                                       BaseVector *v,
-                                                                       const DataVecCoord &d_x1, const DataVecCoord &d_x2
-                                                                       , const DataVecDeriv & d_v1, const DataVecDeriv & d_v2)
+    BaseVector* v,
+    const DataVecCoord& d_x1, const DataVecCoord& d_x2,
+    const DataVecDeriv& d_v1, const DataVecDeriv& d_v2)
 {
     if (!d_activate.getValue()) return;
 
@@ -196,18 +200,20 @@ void BilateralLagrangianConstraint<DataTypes>::getConstraintViolation(const Cons
     const VecCoord &x1 = d_x1.getValue();
     const VecCoord &x2 = d_x2.getValue();
 
-    dfree.resize(minp);
+    m_violation.resize(minp);
 
-    for (unsigned pid=0; pid<minp; pid++)
+    for (unsigned pid = 0; pid < minp; pid++)
     {
-        dfree[pid] = x2[m2Indices[pid]] - x1[m1Indices[pid]];
+        m_violation[pid] = x2[m2Indices[pid]] - x1[m1Indices[pid]];
 
         if (pid < restVector.size())
-            dfree[pid] -= restVector[pid];
+        {
+            m_violation[pid] -= restVector[pid];
+        }
 
-        v->set(cid[pid]  , dfree[pid][0]);
-        v->set(cid[pid]+1, dfree[pid][1]);
-        v->set(cid[pid]+2, dfree[pid][2]);
+        v->set(cid[pid]  , m_violation[pid][0]);
+        v->set(cid[pid]+1, m_violation[pid][1]);
+        v->set(cid[pid]+2, m_violation[pid][2]);
     }
 }
 
@@ -235,7 +241,7 @@ void BilateralLagrangianConstraint<DataTypes>::getVelocityViolation(BaseVector *
     auto pos2 = this->getMState2()->readPositions();
 
     const SReal dt = this->getContext()->getDt();
-    const SReal invDt = SReal(1.0) / dt;
+    const SReal invDt = 1_sreal / dt;
 
     for (unsigned pid=0; pid<minp; ++pid)
     {
diff --git a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/SlidingLagrangianConstraint.h b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/SlidingLagrangianConstraint.h
index 2aec761548d..4c9e9aa5c9c 100644
--- a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/SlidingLagrangianConstraint.h
+++ b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/SlidingLagrangianConstraint.h
@@ -66,7 +66,7 @@ class SlidingLagrangianConstraint : public core::behavior::PairInteractionConstr
     SlidingLagrangianConstraint(MechanicalState* object);
     SlidingLagrangianConstraint(MechanicalState* object1, MechanicalState* object2);
 
-    virtual ~SlidingLagrangianConstraint(){}
+    ~SlidingLagrangianConstraint() override {}
 
 
 
diff --git a/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.cpp b/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.cpp
index cf3d4e1a770..e33cbaa1b0b 100644
--- a/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.cpp
+++ b/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.cpp
@@ -60,8 +60,10 @@ void clearMultiVecId(sofa::core::objectmodel::BaseContext* ctx, const sofa::core
 
 }
 
+static constexpr GenericConstraintSolver::ResolutionMethod defaultResolutionMethod("ProjectedGaussSeidel");
+
 GenericConstraintSolver::GenericConstraintSolver()
-    : d_resolutionMethod( initData(&d_resolutionMethod, "resolutionMethod", "Method used to solve the constraint problem, among: \"ProjectedGaussSeidel\", \"UnbuiltGaussSeidel\" or \"for NonsmoothNonlinearConjugateGradient\""))
+    : d_resolutionMethod( initData(&d_resolutionMethod, defaultResolutionMethod, "resolutionMethod", ("Method used to solve the constraint problem\n" + ResolutionMethod::dataDescription()).c_str()))
     , d_maxIt(initData(&d_maxIt, 1000, "maxIterations", "maximal number of iterations of the Gauss-Seidel algorithm"))
     , d_tolerance(initData(&d_tolerance, 0.001_sreal, "tolerance", "residual error threshold for termination of the Gauss-Seidel algorithm"))
     , d_sor(initData(&d_sor, 1.0_sreal, "sor", "Successive Over Relaxation parameter (0-2)"))
@@ -86,10 +88,6 @@ GenericConstraintSolver::GenericConstraintSolver()
     , current_cp(&m_cpBuffer[0])
     , last_cp(nullptr)
 {
-    sofa::helper::OptionsGroup m_newoptiongroup{"ProjectedGaussSeidel","UnbuiltGaussSeidel", "NonsmoothNonlinearConjugateGradient"};
-    m_newoptiongroup.setSelectedItem("ProjectedGaussSeidel");
-    d_resolutionMethod.setValue(m_newoptiongroup);
-
     addAlias(&d_maxIt, "maxIt");
 
     d_graphErrors.setWidget("graph");
@@ -159,7 +157,8 @@ void GenericConstraintSolver::init()
 
     if(d_newtonIterations.isSet())
     {
-        if (d_resolutionMethod.getValue().getSelectedId() != 2)
+        static constexpr ResolutionMethod NonsmoothNonlinearConjugateGradient("NonsmoothNonlinearConjugateGradient");
+        if (d_resolutionMethod.getValue() != NonsmoothNonlinearConjugateGradient)
         {
             msg_warning() << "data \"newtonIterations\" is not only taken into account when using the NonsmoothNonlinearConjugateGradient solver";
         }
@@ -225,15 +224,15 @@ bool GenericConstraintSolver::buildSystem(const core::ConstraintParams *cParams,
     }
 
     // Resolution depending on the method selected
-    switch ( d_resolutionMethod.getValue().getSelectedId() )
+    switch ( d_resolutionMethod.getValue() )
     {
-        case 0: // ProjectedGaussSeidel
-        case 2: // NonsmoothNonlinearConjugateGradient
+        case ResolutionMethod("ProjectedGaussSeidel"):
+        case ResolutionMethod("NonsmoothNonlinearConjugateGradient"):
         {
             buildSystem_matrixAssembly(cParams);
             break;
         }
-        case 1: // UnbuiltGaussSeidel
+        case ResolutionMethod("UnbuiltGaussSeidel"):
         {
             buildSystem_matrixFree(numConstraints);
             break;
@@ -429,10 +428,9 @@ bool GenericConstraintSolver::solveSystem(const core::ConstraintParams * /*cPara
 
 
     // Resolution depending on the method selected
-    switch ( d_resolutionMethod.getValue().getSelectedId() )
+    switch ( d_resolutionMethod.getValue())
     {
-        // ProjectedGaussSeidel
-        case 0: {
+        case ResolutionMethod("ProjectedGaussSeidel"): {
             if (notMuted())
             {
                 std::stringstream tmp;
@@ -445,14 +443,12 @@ bool GenericConstraintSolver::solveSystem(const core::ConstraintParams * /*cPara
             current_cp->gaussSeidel(0, this);
             break;
         }
-        // UnbuiltGaussSeidel
-        case 1: {
+        case ResolutionMethod("UnbuiltGaussSeidel"): {
             SCOPED_TIMER_VARNAME(unbuiltGaussSeidelTimer, "ConstraintsUnbuiltGaussSeidel");
             current_cp->unbuiltGaussSeidel(0, this);
             break;
         }
-        // NonsmoothNonlinearConjugateGradient
-        case 2: {
+        case ResolutionMethod("NonsmoothNonlinearConjugateGradient"): {
             current_cp->NNCG(this, d_newtonIterations.getValue());
             break;
         }
diff --git a/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.h b/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.h
index 7d752048930..cd3c770fd97 100644
--- a/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.h
+++ b/Sofa/Component/Constraint/Lagrangian/Solver/src/sofa/component/constraint/lagrangian/solver/GenericConstraintSolver.h
@@ -33,6 +33,8 @@
 #include <sofa/component/constraint/lagrangian/solver/visitors/MechanicalGetConstraintResolutionVisitor.h>
 
 #include <sofa/core/objectmodel/RenamedData.h>
+#include <sofa/helper/SelectableItem.h>
+
 
 namespace sofa::component::constraint::lagrangian::solver
 {
@@ -60,7 +62,13 @@ class SOFA_COMPONENT_CONSTRAINT_LAGRANGIAN_SOLVER_API GenericConstraintSolver :
     ConstraintProblem* getConstraintProblem() override;
     void lockConstraintProblem(sofa::core::objectmodel::BaseObject* from, ConstraintProblem* p1, ConstraintProblem* p2 = nullptr) override;
 
-    Data< sofa::helper::OptionsGroup > d_resolutionMethod; ///< Method used to solve the constraint problem, among: "ProjectedGaussSeidel", "UnbuiltGaussSeidel" or "for NonsmoothNonlinearConjugateGradient"
+    MAKE_SELECTABLE_ITEMS(ResolutionMethod,
+        sofa::helper::Item{"ProjectedGaussSeidel", "Projected Gauss-Seidel"},
+        sofa::helper::Item{"UnbuiltGaussSeidel", "Gauss-Seidel where the matrix is not assembled"},
+        sofa::helper::Item{"NonsmoothNonlinearConjugateGradient", "Non-smooth non-linear conjugate gradient"}
+    );
+
+    Data< ResolutionMethod > d_resolutionMethod; ///< Method used to solve the constraint problem, among: "ProjectedGaussSeidel", "UnbuiltGaussSeidel" or "for NonsmoothNonlinearConjugateGradient"
 
     SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_CONSTRAINT_LAGRANGIAN_SOLVER()
     sofa::core::objectmodel::RenamedData<int> maxIt;
diff --git a/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt b/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt
index e94666809ec..50514787466 100644
--- a/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt
+++ b/Sofa/Component/LinearSolver/Iterative/CMakeLists.txt
@@ -18,6 +18,8 @@ set(HEADER_FILES
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MinResLinearSolver.inl
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/ShewchukPCGLinearSolver.h
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/ShewchukPCGLinearSolver.inl
+    ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.h
+    ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.inl
 )
 
 set(SOURCE_FILES
@@ -29,7 +31,7 @@ set(SOURCE_FILES
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MatrixLinearSolver.cpp
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MatrixLinearSystem[GraphScattered].cpp
     ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/MinResLinearSolver.cpp
-    ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/ShewchukPCGLinearSolver.cpp
+    ${SOFACOMPONENTLINEARSOLVERITERATIVE_SOURCE_DIR}/PCGLinearSolver.cpp
 )
 
 sofa_find_package(Sofa.Simulation.Core REQUIRED)
diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/GraphScatteredTypes.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/GraphScatteredTypes.cpp
index 51dacd5f02d..f9aa5f3a656 100644
--- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/GraphScatteredTypes.cpp
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/GraphScatteredTypes.cpp
@@ -34,7 +34,10 @@ void GraphScatteredMatrix::apply(GraphScatteredVector& res, GraphScatteredVector
 {
     // matrix-vector product through visitors
     parent->propagateDxAndResetDf(x,res);
-    parent->addMBKdx(res,parent->mparams.mFactor(),parent->mparams.bFactor(),parent->mparams.kFactor(), false); // df = (m M + b B + k K) dx
+    parent->addMBKdx(res,
+        core::MatricesFactors::M(parent->mparams.mFactor()),
+        core::MatricesFactors::B(parent->mparams.bFactor()),
+        core::MatricesFactors::K(parent->mparams.kFactor()), false); // df = (m M + b B + k K) dx
 
     // filter the product to take the constraints into account
     //
diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.cpp b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp
similarity index 82%
rename from Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.cpp
rename to Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp
index 5cd240c90db..72f4df45bed 100644
--- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.cpp
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp
@@ -19,8 +19,8 @@
 *                                                                             *
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
-#define SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_SHEWCHUKPCGLINEARSOLVER_CPP
-#include <sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.inl>
+#define SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP
+#include <sofa/component/linearsolver/iterative/PCGLinearSolver.inl>
 #include <sofa/component/linearsolver/iterative/MatrixLinearSolver.inl>
 #include <sofa/core/ObjectFactory.h>
 #include <sofa/component/linearsystem/MatrixFreeSystem.h>
@@ -28,12 +28,12 @@
 namespace sofa::component::linearsolver::iterative
 {
 
-void registerShewchukPCGLinearSolver(sofa::core::ObjectFactory* factory)
+void registerPCGLinearSolver(sofa::core::ObjectFactory* factory)
 {
     factory->registerObjects(core::ObjectRegistrationData("Linear system solver using the Shewchuk conjugate gradient iterative algorithm.")
-        .add< ShewchukPCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >());
+        .add< PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >());
 }
 
-template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API ShewchukPCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
+template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
 
 } // namespace sofa::component::linearsolver::iterative
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
new file mode 100644
index 00000000000..1506256becd
--- /dev/null
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h
@@ -0,0 +1,123 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/component/linearsolver/iterative/config.h>
+
+#include <sofa/core/behavior/LinearSolver.h>
+#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
+#include <sofa/helper/map.h>
+
+#include <cmath>
+
+#include <sofa/core/objectmodel/RenamedData.h>
+
+namespace sofa::component::linearsolver::iterative
+{
+
+/// Linear system solver using the conjugate gradient iterative algorithm
+template<class TMatrix, class TVector>
+class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector>
+{
+public:
+    SOFA_CLASS(SOFA_TEMPLATE2(PCGLinearSolver,TMatrix,TVector),SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector));
+
+    typedef TMatrix Matrix;
+    typedef TVector Vector;
+    typedef sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector> Inherit;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<unsigned> f_maxIter;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<double> f_tolerance;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<bool> f_use_precond;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<unsigned> f_update_step;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<bool> f_build_precond;
+
+    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
+    sofa::core::objectmodel::RenamedData<std::map < std::string, sofa::type::vector<double> > >  f_graph;
+
+
+    Data<unsigned> d_maxIter; ///< Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop
+    Data<double> d_tolerance; ///< Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm)
+    Data<bool> d_use_precond; ///< Use a preconditioner
+    SingleLink<PCGLinearSolver, sofa::core::behavior::LinearSolver, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_preconditioner; ///< Link towards the linear solver used to precondition the conjugate gradient
+    Data<unsigned> d_update_step; ///< Number of steps before the next refresh of preconditioners
+    Data<bool> d_build_precond; ///< Build the preconditioners, if false build the preconditioner only at the initial step
+    Data<std::map < std::string, sofa::type::vector<double> > > d_graph; ///< Graph of residuals at each iteration
+
+protected:
+    PCGLinearSolver();
+
+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;
+    bool first;
+    int newton_iter;
+
+protected:
+    /// This method is separated from the rest to be able to use custom/optimized versions depending on the types of vectors.
+    /// It computes: p = p*beta + r
+    inline void cgstep_beta(Vector& p, Vector& r, double beta);
+    /// 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);
+
+    void handleEvent(sofa::core::objectmodel::Event* event) override;
+
+
+};
+
+template<class TMatrix, class TVector>
+inline void PCGLinearSolver<TMatrix,TVector>::cgstep_beta(Vector& p, Vector& r, double beta)
+{
+    p *= beta;
+    p += r; //z;
+}
+
+template<class TMatrix, class TVector>
+inline void PCGLinearSolver<TMatrix,TVector>::cgstep_alpha(Vector& x,Vector& p,double alpha)
+{
+    x.peq(p,alpha);                 // x = x + alpha p
+}
+
+template<>
+inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_beta(Vector& p, Vector& r, double beta);
+
+template<>
+inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_alpha(Vector& x,Vector& p,double alpha);
+
+#if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)
+template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
+#endif // !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)
+
+} // namespace sofa::component::linearsolver::iterative
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
new file mode 100644
index 00000000000..ff9da2bb46e
--- /dev/null
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl
@@ -0,0 +1,247 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/component/linearsolver/iterative/PCGLinearSolver.h>
+
+#include <sofa/core/behavior/LinearSolver.h>
+#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
+#include <sofa/simulation/AnimateBeginEvent.h>
+#include <sofa/helper/map.h>
+#include <sofa/helper/AdvancedTimer.h>
+#include <sofa/helper/ScopedAdvancedTimer.h>
+ 
+#include <cmath>
+
+namespace sofa::component::linearsolver::iterative
+{
+
+template<class TMatrix, class TVector>
+PCGLinearSolver<TMatrix,TVector>::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_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)
+{
+    d_graph.setWidget("graph");
+    first = true;
+    this->f_listening.setValue(true);
+
+    f_maxIter.setOriginalData(&d_maxIter);
+    f_tolerance.setOriginalData(&d_tolerance);
+    f_use_precond.setOriginalData(&d_use_precond);
+    f_update_step.setOriginalData(&d_update_step);
+    f_build_precond.setOriginalData(&d_build_precond);
+    f_graph.setOriginalData(&d_graph);
+
+}
+
+template<class TMatrix, class TVector>
+void PCGLinearSolver<TMatrix,TVector>::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.";
+    }
+    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);
+            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);
+                return;
+            }
+            else
+            {
+                msg_info() << "Preconditioner path used: '" << l_preconditioner.getLinkedPath() << "'";
+            }
+        }
+    }
+
+    first = true;
+    sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid);
+}
+
+template<class TMatrix, class TVector>
+void PCGLinearSolver<TMatrix,TVector>::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<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_beta(Vector& p, Vector& r, double beta)
+{
+    p.eq(r,p,beta); // p = p*beta + r
+}
+
+template<>
+inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_alpha(Vector& x, Vector& p, double alpha)
+{
+    x.peq(p,alpha);                 // x = x + alpha p
+}
+
+template<class Matrix, class Vector>
+void PCGLinearSolver<Matrix,Vector>::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<double> >& graph = * d_graph.beginEdit();
+        graph.clear();
+    }
+}
+
+
+template<class TMatrix, class TVector>
+void PCGLinearSolver<TMatrix,TVector>::solve (Matrix& M, Vector& x, Vector& b)
+{
+    SCOPED_TIMER_VARNAME(solveTimer, "PCGLinearSolver::solve");
+
+    std::map < std::string, sofa::type::vector<double> >& graph = * d_graph.beginEdit();
+//    sofa::type::vector<double>& graph_error = graph["Error"];
+
+    newton_iter++;
+    char name[256];
+    sprintf(name,"Error %d",newton_iter);
+    sofa::type::vector<double>& graph_error = graph[std::string(name)];
+
+    const core::ExecParams* params = core::execparams::defaultInstance();
+    typename Inherit::TempVectorContainer vtmp(this, params, M, x, b);
+    Vector& r = *vtmp.createTempVector();
+    Vector& w = *vtmp.createTempVector();
+    Vector& s = *vtmp.createTempVector();
+
+    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;
+
+    r = M * x;
+    cgstep_beta(r,b,-1);// r = -1 * r + b  =   b - (M * x)
+
+    if (apply_precond)
+    {
+        SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
+        l_preconditioner.get()->setSystemLHVector(w);
+        l_preconditioner.get()->setSystemRHVector(r);
+        l_preconditioner.get()->solveSystem();
+    }
+    else
+    {
+        w = r;
+    }
+
+    double 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;
+
+        cgstep_alpha(x,w,alpha);//for(int i=0; i<n; i++) x[i] += alpha * d[i];
+        cgstep_alpha(r,s,-alpha);//for (int i=0; i<n; i++) r[i] = r[i] - alpha * q[i];
+
+        if (apply_precond)
+        {
+            SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
+            l_preconditioner.get()->setSystemLHVector(s);
+            l_preconditioner.get()->setSystemRHVector(r);
+            l_preconditioner.get()->solveSystem();
+        }
+        else
+        {
+            s = r;
+        }
+
+        const double deltaOld = r_norm;
+        r_norm = r.dot(s);
+        graph_error.push_back(r_norm/b_norm);
+
+        double beta = r_norm / deltaOld;
+
+        cgstep_beta(w,s,beta);//for (int i=0; i<n; i++) d[i] = r[i] + beta * d[i];
+
+        iter++;
+    }
+
+    d_graph.endEdit();
+
+    vtmp.deleteTempVector(&r);
+    vtmp.deleteTempVector(&w);
+    vtmp.deleteTempVector(&s);
+
+    sofa::helper::AdvancedTimer::valSet("PCG iterations", iter);
+}
+
+} // namespace sofa::component::linearsolver::iterative
diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.h b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.h
index a77daf51db7..3b0edb3e7f7 100644
--- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.h
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.h
@@ -20,104 +20,13 @@
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
 #pragma once
-#include <sofa/component/linearsolver/iterative/config.h>
 
-#include <sofa/core/behavior/LinearSolver.h>
-#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
-#include <sofa/helper/map.h>
+#include <sofa/component/linearsolver/iterative/PCGLinearSolver.h>
 
-#include <cmath>
-
-#include <sofa/core/objectmodel/RenamedData.h>
+SOFA_HEADER_DEPRECATED("v24.12", "v25.12", "sofa/component/linearsolver/iterative/PCGLinearSolver.h")
 
 namespace sofa::component::linearsolver::iterative
 {
-
-/// Linear system solver using the conjugate gradient iterative algorithm
-template<class TMatrix, class TVector>
-class ShewchukPCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector>
-{
-public:
-    SOFA_CLASS(SOFA_TEMPLATE2(ShewchukPCGLinearSolver,TMatrix,TVector),SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector));
-
-    typedef TMatrix Matrix;
-    typedef TVector Vector;
-    typedef sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector> Inherit;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<unsigned> f_maxIter;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<double> f_tolerance;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<bool> f_use_precond;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<unsigned> f_update_step;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<bool> f_build_precond;
-
-    SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_LINEARSOLVER_ITERATIVE()
-    sofa::core::objectmodel::RenamedData<std::map < std::string, sofa::type::vector<double> > >  f_graph;
-
-
-    Data<unsigned> d_maxIter; ///< Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop
-    Data<double> d_tolerance; ///< Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm)
-    Data<bool> d_use_precond; ///< Use a preconditioner
-    SingleLink<ShewchukPCGLinearSolver, sofa::core::behavior::LinearSolver, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_preconditioner; ///< Link towards the linear solver used to precondition the conjugate gradient
-    Data<unsigned> d_update_step; ///< Number of steps before the next refresh of preconditioners
-    Data<bool> d_build_precond; ///< Build the preconditioners, if false build the preconditioner only at the initial step
-    Data<std::map < std::string, sofa::type::vector<double> > > d_graph; ///< Graph of residuals at each iteration
-
-protected:
-    ShewchukPCGLinearSolver();
-
-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;
-    bool first;
-    int newton_iter;
-
-protected:
-    /// This method is separated from the rest to be able to use custom/optimized versions depending on the types of vectors.
-    /// It computes: p = p*beta + r
-    inline void cgstep_beta(Vector& p, Vector& r, double beta);
-    /// 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);
-
-    void handleEvent(sofa::core::objectmodel::Event* event) override;
-
-
-};
-
-template<class TMatrix, class TVector>
-inline void ShewchukPCGLinearSolver<TMatrix,TVector>::cgstep_beta(Vector& p, Vector& r, double beta)
-{
-    p *= beta;
-    p += r; //z;
-}
-
 template<class TMatrix, class TVector>
-inline void ShewchukPCGLinearSolver<TMatrix,TVector>::cgstep_alpha(Vector& x,Vector& p,double alpha)
-{
-    x.peq(p,alpha);                 // x = x + alpha p
+using ShewchukPCGLinearSolver SOFA_ATTRIBUTE_DEPRECATED("v24.12", "v25.12", "ShewchukPCGLinearSolver has been renamed to PCGLinearSolver") = PCGLinearSolver<TMatrix, TVector>;
 }
-
-template<>
-inline void ShewchukPCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_beta(Vector& p, Vector& r, double beta);
-
-template<>
-inline void ShewchukPCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_alpha(Vector& x,Vector& p,double alpha);
-
-#if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_SHEWCHUKPCGLINEARSOLVER_CPP)
-template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API ShewchukPCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
-#endif // !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_SHEWCHUKPCGLINEARSOLVER_CPP)
-
-} // namespace sofa::component::linearsolver::iterative
diff --git a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.inl b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.inl
index 68eb09ae9d0..8345764ae6c 100644
--- a/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.inl
+++ b/Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.inl
@@ -20,228 +20,7 @@
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
 #pragma once
-#include <sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.h>
 
-#include <sofa/core/behavior/LinearSolver.h>
-#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
-#include <sofa/simulation/AnimateBeginEvent.h>
-#include <sofa/helper/map.h>
-#include <sofa/helper/AdvancedTimer.h>
-#include <sofa/helper/ScopedAdvancedTimer.h>
+#include <sofa/component/linearsolver/iterative/PCGLinearSolver.inl>
 
-#include <cmath>
-
-namespace sofa::component::linearsolver::iterative
-{
-
-template<class TMatrix, class TVector>
-ShewchukPCGLinearSolver<TMatrix,TVector>::ShewchukPCGLinearSolver()
-    : d_maxIter(initData(&d_maxIter, (unsigned)25, "iterations", "Maximum number of iterations after which the iterative descent of the Conjugate Gradient must stop") )
-    , d_tolerance(initData(&d_tolerance, 1e-5, "tolerance", "Desired accuracy of the Conjugate Gradient solution evaluating: |r|²/|b|² (ratio of current residual norm over initial residual norm)") )
-    , d_use_precond(initData(&d_use_precond, true, "use_precond", "Use a preconditioner") )
-    , l_preconditioner(initLink("preconditioner", "Link towards the linear solver used to precondition the conjugate gradient"))
-    , d_update_step(initData(&d_update_step, (unsigned)1, "update_step", "Number of steps before the next refresh of preconditioners") )
-    , d_build_precond(initData(&d_build_precond, true, "build_precond", "Build the preconditioners, if false build the preconditioner only at the initial step") )
-    , d_graph(initData(&d_graph, "graph", "Graph of residuals at each iteration") )
-    , next_refresh_step(0)
-    , newton_iter(0)
-{
-    d_graph.setWidget("graph");
-    first = true;
-    this->f_listening.setValue(true);
-
-    f_maxIter.setOriginalData(&d_maxIter);
-    f_tolerance.setOriginalData(&d_tolerance);
-    f_use_precond.setOriginalData(&d_use_precond);
-    f_update_step.setOriginalData(&d_update_step);
-    f_build_precond.setOriginalData(&d_build_precond);
-    f_graph.setOriginalData(&d_graph);
-
-}
-
-template<class TMatrix, class TVector>
-void ShewchukPCGLinearSolver<TMatrix,TVector>::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.";
-    }
-    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);
-            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);
-                return;
-            }
-            else
-            {
-                msg_info() << "Preconditioner path used: '" << l_preconditioner.getLinkedPath() << "'";
-            }
-        }
-    }
-
-    first = true;
-    sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid);
-}
-
-template<class TMatrix, class TVector>
-void ShewchukPCGLinearSolver<TMatrix,TVector>::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 ShewchukPCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_beta(Vector& p, Vector& r, double beta)
-{
-    p.eq(r,p,beta); // p = p*beta + r
-}
-
-template<>
-inline void ShewchukPCGLinearSolver<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector>::cgstep_alpha(Vector& x, Vector& p, double alpha)
-{
-    x.peq(p,alpha);                 // x = x + alpha p
-}
-
-template<class Matrix, class Vector>
-void ShewchukPCGLinearSolver<Matrix,Vector>::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<double> >& graph = * d_graph.beginEdit();
-        graph.clear();
-    }
-}
-
-
-template<class TMatrix, class TVector>
-void ShewchukPCGLinearSolver<TMatrix,TVector>::solve (Matrix& M, Vector& x, Vector& b)
-{
-    SCOPED_TIMER_VARNAME(solveTimer, "PCGLinearSolver::solve");
-
-    std::map < std::string, sofa::type::vector<double> >& graph = * d_graph.beginEdit();
-//    sofa::type::vector<double>& graph_error = graph["Error"];
-
-    newton_iter++;
-    char name[256];
-    sprintf(name,"Error %d",newton_iter);
-    sofa::type::vector<double>& graph_error = graph[std::string(name)];
-
-    const core::ExecParams* params = core::execparams::defaultInstance();
-    typename Inherit::TempVectorContainer vtmp(this, params, M, x, b);
-    Vector& r = *vtmp.createTempVector();
-    Vector& w = *vtmp.createTempVector();
-    Vector& s = *vtmp.createTempVector();
-
-    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;
-
-    r = M * x;
-    cgstep_beta(r,b,-1);// r = -1 * r + b  =   b - (M * x)
-
-    if (apply_precond)
-    {
-        SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
-        l_preconditioner.get()->setSystemLHVector(w);
-        l_preconditioner.get()->setSystemRHVector(r);
-        l_preconditioner.get()->solveSystem();
-    }
-    else
-    {
-        w = r;
-    }
-
-    double 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;
-
-        cgstep_alpha(x,w,alpha);//for(int i=0; i<n; i++) x[i] += alpha * d[i];
-        cgstep_alpha(r,s,-alpha);//for (int i=0; i<n; i++) r[i] = r[i] - alpha * q[i];
-
-        if (apply_precond)
-        {
-            SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
-            l_preconditioner.get()->setSystemLHVector(s);
-            l_preconditioner.get()->setSystemRHVector(r);
-            l_preconditioner.get()->solveSystem();
-        }
-        else
-        {
-            s = r;
-        }
-
-        const double deltaOld = r_norm;
-        r_norm = r.dot(s);
-        graph_error.push_back(r_norm/b_norm);
-
-        double beta = r_norm / deltaOld;
-
-        cgstep_beta(w,s,beta);//for (int i=0; i<n; i++) d[i] = r[i] + beta * d[i];
-
-        iter++;
-    }
-
-    d_graph.endEdit();
-
-    vtmp.deleteTempVector(&r);
-    vtmp.deleteTempVector(&w);
-    vtmp.deleteTempVector(&s);
-
-    sofa::helper::AdvancedTimer::valSet("PCG iterations", iter);
-}
-
-} // namespace sofa::component::linearsolver::iterative
+SOFA_HEADER_DEPRECATED("v24.12", "v25.12", "sofa/component/linearsolver/iterative/PCGLinearSolver.inl")
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 c29da966335..f6fa6c78938 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
@@ -33,7 +33,7 @@ namespace sofa::component::linearsolver::iterative
 
 extern void registerCGLinearSolver(sofa::core::ObjectFactory* factory);
 extern void registerMinResLinearSolver(sofa::core::ObjectFactory* factory);
-extern void registerShewchukPCGLinearSolver(sofa::core::ObjectFactory* factory);
+extern void registerPCGLinearSolver(sofa::core::ObjectFactory* factory);
 
 extern "C" {
     SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule();
@@ -62,7 +62,7 @@ void registerObjects(sofa::core::ObjectFactory* factory)
     registerCGLinearSolver(factory);
     linearsystem::registerMatrixFreeSystemGraphScattered(factory);
     registerMinResLinearSolver(factory);
-    registerShewchukPCGLinearSolver(factory);
+    registerPCGLinearSolver(factory);
 }
 
 void init()
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.h
index 8f5b998cd95..2ebb82cf83b 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.h
@@ -30,7 +30,7 @@
 namespace sofa::component::mapping::linear
 {
 
-/// concatanate several entire mechanical states together
+/// Concatenate several mechanical states together
 template <class TIn, class TOut>
 class IdentityMultiMapping : public LinearMultiMapping<TIn, TOut>
 {
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.h
index 066888c34e5..75fe2fe440c 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.h
@@ -34,10 +34,6 @@
 namespace sofa::component::mapping::linear
 {
 
-/**
- * @class SubsetMapping
- * @brief Compute a subset of input points
- */
 template <class TIn, class TOut>
 class SubsetMultiMapping : public LinearMultiMapping<TIn, TOut>
 {
@@ -48,28 +44,6 @@ class SubsetMultiMapping : public LinearMultiMapping<TIn, TOut>
     typedef TIn In;
     typedef TOut Out;
 
-    typedef typename Out::VecCoord OutVecCoord;
-    typedef typename Out::VecDeriv OutVecDeriv;
-    typedef typename Out::Coord OutCoord;
-    typedef typename Out::Deriv OutDeriv;
-    typedef typename Out::MatrixDeriv OutMatrixDeriv;
-    typedef typename In::VecCoord InVecCoord;
-    typedef typename In::VecDeriv InVecDeriv;
-    typedef typename In::Coord InCoord;
-    typedef typename In::Deriv InDeriv;
-    typedef typename In::MatrixDeriv InMatrixDeriv;
-    typedef typename InCoord::value_type Real;
-    typedef typename type::vector<const InVecCoord*> vecConstInVecCoord;
-    typedef typename type::vector<OutVecCoord*> vecOutVecCoord;
-
-    typedef Data<InVecCoord> InDataVecCoord;
-    typedef Data<InVecDeriv> InDataVecDeriv;
-    typedef Data<InMatrixDeriv> InDataMatrixDeriv;
-
-    typedef Data<OutVecCoord> OutDataVecCoord;
-    typedef Data<OutVecDeriv> OutDataVecDeriv;
-    typedef Data<OutMatrixDeriv> OutDataMatrixDeriv;
-
     /// Correspondence array
     typedef core::topology::BaseMeshTopology::SetIndex IndexArray;
 
@@ -81,12 +55,12 @@ class SubsetMultiMapping : public LinearMultiMapping<TIn, TOut>
     void addPoint(int fromModel, int index);
 
     // usual Mapping API
-    void apply(const core::MechanicalParams* mparams, const type::vector<OutDataVecCoord*>& dataVecOutPos, const type::vector<const InDataVecCoord*>& dataVecInPos) override;
-    void applyJ(const core::MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel) override;
-    void applyJT(const core::MechanicalParams* mparams, const type::vector<InDataVecDeriv*>& dataVecOutForce, const type::vector<const OutDataVecDeriv*>& dataVecInForce) override;
+    void apply(const core::MechanicalParams* mparams, const type::vector<DataVecCoord_t<Out>*>& dataVecOutPos, const type::vector<const DataVecCoord_t<In>*>& dataVecInPos) override;
+    void applyJ(const core::MechanicalParams* mparams, const type::vector<DataVecDeriv_t<Out>*>& dataVecOutVel, const type::vector<const DataVecDeriv_t<In>*>& dataVecInVel) override;
+    void applyJT(const core::MechanicalParams* mparams, const type::vector<DataVecDeriv_t<In>*>& dataVecOutForce, const type::vector<const DataVecDeriv_t<Out>*>& dataVecInForce) override;
     void applyDJT(const core::MechanicalParams* /*mparams*/, core::MultiVecDerivId /*inForce*/, core::ConstMultiVecDerivId /*outForce*/) override {}
 
-    void applyJT( const core::ConstraintParams* cparams, const type::vector< InDataMatrixDeriv* >& dataMatOutConst, const type::vector< const OutDataMatrixDeriv* >& dataMatInConst ) override;
+    void applyJT( const core::ConstraintParams* cparams, const type::vector< DataMatrixDeriv_t<In>* >& dataMatOutConst, const type::vector< const DataMatrixDeriv_t<Out>* >& dataMatInConst ) override;
 
     /// Experimental API used to handle multimappings in matrix assembly. Returns pointers to matrices associated with parent states, consistently with  getFrom().
     virtual const type::vector<sofa::linearalgebra::BaseMatrix*>* getJs() override;
@@ -99,7 +73,7 @@ class SubsetMultiMapping : public LinearMultiMapping<TIn, TOut>
 protected :
 
     SubsetMultiMapping();
-    virtual ~SubsetMultiMapping();
+    virtual ~SubsetMultiMapping() override;
 
     type::vector<linearalgebra::BaseMatrix*> baseMatrices;      ///< Jacobian of the mapping, in a vector
 
@@ -112,7 +86,6 @@ extern template class SOFA_COMPONENT_MAPPING_LINEAR_API SubsetMultiMapping< defa
 extern template class SOFA_COMPONENT_MAPPING_LINEAR_API SubsetMultiMapping< defaulttype::Vec1Types, defaulttype::Vec1Types >;
 extern template class SOFA_COMPONENT_MAPPING_LINEAR_API SubsetMultiMapping< defaulttype::Rigid3Types, defaulttype::Rigid3Types >;
 extern template class SOFA_COMPONENT_MAPPING_LINEAR_API SubsetMultiMapping< defaulttype::Rigid3Types, defaulttype::Vec3Types >;
-
 #endif
 
 } // namespace sofa::component::mapping::linear
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.inl
index 261b0f853fd..62e1d9dd422 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMultiMapping.inl
@@ -43,25 +43,29 @@ void SubsetMultiMapping<TIn, TOut>::init()
     static constexpr auto Nout = TOut::deriv_total_size;
 
 
-    for( unsigned i=0; i<baseMatrices.size(); i++ )
+    for (unsigned i = 0; i < baseMatrices.size(); i++)
+    {
         delete baseMatrices[i];
+    }
 
     typedef linearalgebra::EigenSparseMatrix<TIn,TOut> Jacobian;
     baseMatrices.resize( this->getFrom().size() );
     type::vector<Jacobian*> jacobians( this->getFrom().size() );
-    for(unsigned i=0; i<baseMatrices.size(); i++ )
+    for (unsigned i = 0; i < baseMatrices.size(); i++)
     {
         baseMatrices[i] = jacobians[i] = new linearalgebra::EigenSparseMatrix<TIn,TOut>;
         jacobians[i]->resize(Nout*indexPairSize,Nin*this->fromModels[i]->getSize() ); // each jacobian has the same number of rows
     }
 
+    const auto& indexPairsValue = d_indexPairs.getValue();
+
     // fill the Jacobians
-    for(unsigned i=0; i<indexPairSize; i++)
+    for (unsigned i = 0; i < indexPairSize; i++)
     {
-        unsigned parent = d_indexPairs.getValue()[i * 2];
+        unsigned parent = indexPairsValue[i * 2];
         Jacobian* jacobian = jacobians[parent];
-        unsigned bcol = d_indexPairs.getValue()[i * 2 + 1];  // parent particle
-        for(unsigned k=0; k<Nout; k++ )
+        unsigned bcol = indexPairsValue[i * 2 + 1];  // parent particle
+        for (unsigned k = 0; k < Nout; k++)
         {
             unsigned row = i*Nout + k;
 
@@ -70,22 +74,24 @@ void SubsetMultiMapping<TIn, TOut>::init()
         }
     }
     // finalize the Jacobians
-    for(unsigned i=0; i<baseMatrices.size(); i++ )
-        baseMatrices[i]->compress();
+    for (linearalgebra::BaseMatrix* baseMat : baseMatrices)
+    {
+        baseMat->compress();
+    }
 }
 
 template <class TIn, class TOut>
 SubsetMultiMapping<TIn, TOut>::SubsetMultiMapping()
-    : Inherit()
-    , d_indexPairs(initData(&d_indexPairs, type::vector<unsigned>(), "indexPairs", "list of couples (parent index + index in the parent)"))
+    : d_indexPairs(initData(&d_indexPairs, type::vector<unsigned>(), "indexPairs",
+        "list of couples (parent index + index in the parent)"))
 {
-        indexPairs.setOriginalData(&d_indexPairs);
+    indexPairs.setOriginalData(&d_indexPairs);
 }
 
 template <class TIn, class TOut>
 SubsetMultiMapping<TIn, TOut>::~SubsetMultiMapping()
 {
-    for(unsigned i=0; i<baseMatrices.size(); i++ )
+    for (unsigned i = 0; i < baseMatrices.size(); i++)
     {
         delete baseMatrices[i];
     }
@@ -100,15 +106,18 @@ const type::vector<sofa::linearalgebra::BaseMatrix*>* SubsetMultiMapping<TIn, TO
 template <class TIn, class TOut>
 void SubsetMultiMapping<TIn, TOut>::addPoint( const core::BaseState* from, int index)
 {
-
     // find the index of the parent state
     unsigned i;
-    for(i=0; i<this->fromModels.size(); i++)
-        if(this->fromModels.get(i)==from )
+    for (i = 0; i < this->fromModels.size(); i++)
+    {
+        if (this->fromModels.get(i) == from)
+        {
             break;
-    if(i==this->fromModels.size())
+        }
+    }
+    if (i == this->fromModels.size())
     {
-        msg_error() << "SubsetMultiMapping<TIn, TOut>::addPoint, parent " << from->getName() << " not found !";
+        msg_error() << "addPoint, parent " << from->getName() << " not found !";
         assert(0);
     }
 
@@ -119,56 +128,52 @@ template <class TIn, class TOut>
 void SubsetMultiMapping<TIn, TOut>::addPoint( int from, int index)
 {
     assert((size_t)from < this->fromModels.size());
-    type::vector<unsigned>& indexPairsVector = *d_indexPairs.beginEdit();
+    auto indexPairsVector = sofa::helper::getWriteAccessor(d_indexPairs);
     indexPairsVector.push_back(from);
     indexPairsVector.push_back(index);
-    d_indexPairs.endEdit();
 }
 
-
-template <class TIn, class TOut>
-void SubsetMultiMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams, const type::vector<OutDataVecCoord*>& dataVecOutPos, const type::vector<const InDataVecCoord*>& dataVecInPos)
+template <class DataVecOut, class DataVecIn>
+void apply_impl(
+    const type::vector<DataVecOut*>& dataVecOut,
+    const type::vector<const DataVecIn*>& dataVecIn,
+    const type::vector<unsigned>& indexPairs)
 {
-    SOFA_UNUSED(mparams);
-
-    OutVecCoord& out = *(dataVecOutPos[0]->beginEdit());
+    auto out = sofa::helper::getWriteOnlyAccessor(*dataVecOut[0]);
 
-    for(unsigned i=0; i<out.size(); i++)
+    for (unsigned i = 0; i < out.size(); i++)
     {
-        const InDataVecCoord* inPosPtr = dataVecInPos[d_indexPairs.getValue()[i * 2]];
-        const InVecCoord& inPos = (*inPosPtr).getValue();
-
-        core::eq( out[i], inPos[d_indexPairs.getValue()[i * 2 + 1]] );
-    }
+        const unsigned stateId = indexPairs[i * 2];
+        const unsigned coordId = indexPairs[i * 2 + 1];
 
-    dataVecOutPos[0]->endEdit();
+        const DataVecIn* inPosPtr = dataVecIn[stateId];
+        const typename DataVecIn::value_type& inPos = inPosPtr->getValue();
 
+        core::eq( out[i], inPos[coordId] );
+    }
 }
 
 template <class TIn, class TOut>
-void SubsetMultiMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel)
+void SubsetMultiMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams, const type::vector<DataVecCoord_t<Out>*>& dataVecOutPos, const type::vector<const DataVecCoord_t<In>*>& dataVecInPos)
 {
     SOFA_UNUSED(mparams);
+    apply_impl( dataVecOutPos, dataVecInPos, d_indexPairs.getValue());
+}
 
-    OutVecDeriv& out = *(dataVecOutVel[0]->beginEdit());
-
-    for(unsigned i=0; i<out.size(); i++)
-    {
-        const InDataVecDeriv* inDerivPtr = dataVecInVel[d_indexPairs.getValue()[i * 2]];
-        const InVecDeriv& inDeriv = (*inDerivPtr).getValue();
-        core::eq( out[i], inDeriv[d_indexPairs.getValue()[i * 2 + 1]] );
-    }
-
-    dataVecOutVel[0]->endEdit();
+template <class TIn, class TOut>
+void SubsetMultiMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mparams, const type::vector<DataVecDeriv_t<Out>*>& dataVecOutVel, const type::vector<const DataVecDeriv_t<In>*>& dataVecInVel)
+{
+    SOFA_UNUSED(mparams);
+    apply_impl( dataVecOutVel, dataVecInVel, d_indexPairs.getValue());
 }
 
 template <class TIn, class TOut>
-void SubsetMultiMapping<TIn, TOut>::applyJT( const core::ConstraintParams* /*cparams*/, const type::vector< InDataMatrixDeriv* >& dOut, const type::vector< const OutDataMatrixDeriv* >& dIn)
+void SubsetMultiMapping<TIn, TOut>::applyJT( const core::ConstraintParams* /*cparams*/, const type::vector< DataMatrixDeriv_t<In>* >& dOut, const type::vector< const DataMatrixDeriv_t<Out>* >& dIn)
 {
     type::vector<unsigned>  indexP = d_indexPairs.getValue();
 
     // hypothesis: one child only:
-    const OutMatrixDeriv& in = dIn[0]->getValue();
+    const MatrixDeriv_t<Out>& in = dIn[0]->getValue();
 
     if (dOut.size() != this->fromModels.size())
     {
@@ -176,13 +181,13 @@ void SubsetMultiMapping<TIn, TOut>::applyJT( const core::ConstraintParams* /*cpa
         return;
     }
 
-    typename OutMatrixDeriv::RowConstIterator rowItEnd = in.end();
+    typename MatrixDeriv_t<Out>::RowConstIterator rowItEnd = in.end();
     // loop on the constraints defined on the child of the mapping
-    for (typename OutMatrixDeriv::RowConstIterator rowIt = in.begin(); rowIt != rowItEnd; ++rowIt)
+    for (typename MatrixDeriv_t<Out>::RowConstIterator rowIt = in.begin(); rowIt != rowItEnd; ++rowIt)
     {
 
-        typename OutMatrixDeriv::ColConstIterator colIt = rowIt.begin();
-        typename OutMatrixDeriv::ColConstIterator colItEnd = rowIt.end();
+        typename MatrixDeriv_t<Out>::ColConstIterator colIt = rowIt.begin();
+        typename MatrixDeriv_t<Out>::ColConstIterator colItEnd = rowIt.end();
 
 
         // A constraint can be shared by several nodes,
@@ -194,38 +199,40 @@ void SubsetMultiMapping<TIn, TOut>::applyJT( const core::ConstraintParams* /*cpa
         {
             unsigned int index_parent = indexP[colIt.index() * 2]; // 0 or 1 (for now...)
             // writeLine provide an iterator on the line... if this line does not exist, the line is created:
-            typename InMatrixDeriv::RowIterator o = dOut[index_parent]->beginEdit()->writeLine(rowIt.index());
+            typename MatrixDeriv_t<In>::RowIterator o = dOut[index_parent]->beginEdit()->writeLine(rowIt.index());
             dOut[index_parent]->endEdit();
 
             // for each col of the constraint direction, it adds a col in the corresponding parent's constraint direction
             if(d_indexPairs.getValue()[colIt.index() * 2 + 1] < (unsigned int)this->fromModels[index_parent]->getSize())
             {
-                InDeriv tmp; core::eq( tmp, colIt.val() );
+                Deriv_t<In> tmp;
+                core::eq( tmp, colIt.val() );
                 o.addCol(indexP[colIt.index() * 2 + 1], tmp);
             }
             ++colIt;
         }
-
-
-
     }
 }
 
 
 template <class TIn, class TOut>
-void SubsetMultiMapping<TIn, TOut>::applyJT(const core::MechanicalParams* mparams, const type::vector<InDataVecDeriv*>& dataVecOutForce, const type::vector<const OutDataVecDeriv*>& dataVecInForce)
+void SubsetMultiMapping<TIn, TOut>::applyJT(const core::MechanicalParams* mparams, const type::vector<DataVecDeriv_t<In>*>& dataVecOutForce, const type::vector<const DataVecDeriv_t<Out>*>& dataVecInForce)
 {
     SOFA_UNUSED(mparams);
 
-    const OutDataVecDeriv* cderData = dataVecInForce[0];
-    const OutVecDeriv& cder = cderData->getValue();
+    const DataVecDeriv_t<Out>* cderData = dataVecInForce[0];
+    const VecDeriv_t<Out>& cder = cderData->getValue();
+
+    const auto& indexPairsValue = d_indexPairs.getValue();
 
-    for(unsigned i=0; i<cder.size(); i++)
+    for (unsigned i = 0; i < cder.size(); i++)
     {
-        InDataVecDeriv* inDerivPtr = dataVecOutForce[d_indexPairs.getValue()[i * 2]];
-        InVecDeriv& inDeriv = *(*inDerivPtr).beginEdit();
-        core::peq(inDeriv[d_indexPairs.getValue()[i * 2 + 1]], cder[i] );
-        (*inDerivPtr).endEdit();
+        const unsigned stateId = indexPairsValue[i * 2];
+        const unsigned coordId = indexPairsValue[i * 2 + 1];
+
+        DataVecDeriv_t<In>* inDerivPtr = dataVecOutForce[stateId];
+        auto inDeriv = sofa::helper::getWriteAccessor(*inDerivPtr);
+        core::peq(inDeriv[coordId], cder[i] );
     }
 }
 
diff --git a/Sofa/Component/Mapping/Linear/tests/SubsetMultiMapping_test.cpp b/Sofa/Component/Mapping/Linear/tests/SubsetMultiMapping_test.cpp
index b70504770a6..b9fec71addd 100644
--- a/Sofa/Component/Mapping/Linear/tests/SubsetMultiMapping_test.cpp
+++ b/Sofa/Component/Mapping/Linear/tests/SubsetMultiMapping_test.cpp
@@ -21,75 +21,39 @@
 ******************************************************************************/
 #include <sofa/simulation/graph/DAGSimulation.h>
 #include <sofa/defaulttype/VecTypes.h>
-//#include <sofa/defaulttype/RigidTypes.h>
+#include <sofa/defaulttype/RigidTypes.h>
 #include <sofa/component/mapping/linear/SubsetMultiMapping.h>
 
 #include <sofa/component/mapping/testing/MultiMappingTestCreation.h>
 
 namespace sofa {
-namespace {
 
-using namespace core;
-using namespace component;
-using type::Vec;
-using type::Mat;
-using type::vector;
+using component::mapping::linear::SubsetMultiMapping;
 
-
-/**  Test suite for SubsetMultiMapping.
-  */
+/**
+ * Test suite for SubsetMultiMapping.
+ */
 template <typename _SubsetMultiMapping>
 struct SubsetMultiMappingTest : public MultiMapping_test<_SubsetMultiMapping>
 {
+    using In = typename _SubsetMultiMapping::In;
+    using Out = typename _SubsetMultiMapping::Out;
 
-    typedef _SubsetMultiMapping SubsetMultiMapping;
-
-    typedef typename SubsetMultiMapping::In InDataTypes;
-    typedef typename InDataTypes::VecCoord InVecCoord;
-    typedef typename InDataTypes::VecDeriv InVecDeriv;
-    typedef typename InDataTypes::Coord InCoord;
-    typedef typename InDataTypes::Deriv InDeriv;
-    typedef statecontainer::MechanicalObject<InDataTypes> InMechanicalObject;
-    typedef typename InMechanicalObject::ReadVecCoord  ReadInVecCoord;
-    typedef typename InMechanicalObject::WriteVecCoord WriteInVecCoord;
-    typedef typename InMechanicalObject::WriteVecDeriv WriteInVecDeriv;
-    typedef typename InDataTypes::Real InReal;
-    typedef Mat<InDataTypes::spatial_dimensions,InDataTypes::spatial_dimensions,InReal> RotationMatrix;
-
-
-    typedef typename SubsetMultiMapping::Out OutDataTypes;
-    typedef typename OutDataTypes::VecCoord OutVecCoord;
-    typedef typename OutDataTypes::VecDeriv OutVecDeriv;
-    typedef typename OutDataTypes::Coord OutCoord;
-    typedef typename OutDataTypes::Deriv OutDeriv;
-    typedef statecontainer::MechanicalObject<OutDataTypes> OutMechanicalObject;
-    typedef typename OutMechanicalObject::WriteVecCoord WriteOutVecCoord;
-    typedef typename OutMechanicalObject::WriteVecDeriv WriteOutVecDeriv;
-    typedef typename OutMechanicalObject::ReadVecCoord ReadOutVecCoord;
-    typedef typename OutMechanicalObject::ReadVecDeriv ReadOutVecDeriv;
-
-    void SetUp() override
-    {
-    }
-
-    /** @name Test_Cases
-      For each of these cases, we can test if the mapping work
-      */
-    ///@{
-    /** Two parent particles, two children
-    */
+    /**
+     * Two parent particles, two children
+     */
     bool test_two_parents_one_child()
     {
-        const int NP = 2;
+        static constexpr int NP = 2;
         this->setupScene(NP); // NP parents, 1 child
-        SubsetMultiMapping* smm = static_cast<SubsetMultiMapping*>( this->mapping );
+        _SubsetMultiMapping* smm = static_cast<_SubsetMultiMapping*>( this->mapping );
 
         // parent positions
-        vector< InVecCoord > incoords(NP);
-        for( int i=0; i<NP; i++ )
+        type::vector< VecCoord_t<In> > incoords(NP);
+        for (int i = 0; i < NP; i++)
         {
             incoords[i].resize(1);
-            InDataTypes::set( incoords[i][0], i+1.,-2., 3. );
+            In::set( incoords[i][0], i+1.,-2., 3. );
         }
 
         // subset
@@ -97,34 +61,28 @@ struct SubsetMultiMappingTest : public MultiMapping_test<_SubsetMultiMapping>
         smm->addPoint(smm->getMechFrom()[1],0);  // parent, index in parent
 
         // expected child positions
-        OutVecCoord outcoords(2);
-        OutDataTypes::set( outcoords[0], 1.  , -2., 3. );
-        OutDataTypes::set( outcoords[1], 1+1., -2., 3. );
+        VecCoord_t<Out> outcoords(2);
+        Out::set( outcoords[0], 1.  , -2., 3. );
+        Out::set( outcoords[1], 1+1., -2., 3. );
 
         return this->runTest(incoords,outcoords);
     }
-
-
-    ///@}
-
-
 };
 
-
-// Define the list of types to instantiate. We do not necessarily need to test all combinations.
 using ::testing::Types;
 typedef Types<
-mapping::linear::SubsetMultiMapping<defaulttype::Rigid3Types,defaulttype::Rigid3Types>,
-mapping::linear::SubsetMultiMapping<defaulttype::Vec3Types,defaulttype::Vec3Types>
-> DataTypes; // the types to instantiate.
+    SubsetMultiMapping<defaulttype::Vec3Types,defaulttype::Vec3Types>,
+    SubsetMultiMapping<defaulttype::Vec2Types,defaulttype::Vec2Types>,
+    SubsetMultiMapping<defaulttype::Vec1Types,defaulttype::Vec1Types>,
+    SubsetMultiMapping<defaulttype::Rigid3Types,defaulttype::Rigid3Types>,
+    SubsetMultiMapping<defaulttype::Rigid3Types,defaulttype::Vec3Types>
+> DataTypes;
 
-// Test suite for all the instantiations
 TYPED_TEST_SUITE(SubsetMultiMappingTest, DataTypes);
-// first test case
+
 TYPED_TEST( SubsetMultiMappingTest , two_parents_one_child )
 {
     ASSERT_TRUE(this->test_two_parents_one_child());
 }
 
-} // namespace
 } // namespace sofa
diff --git a/Sofa/Component/Mapping/NonLinear/tests/RigidMapping_test.cpp b/Sofa/Component/Mapping/NonLinear/tests/RigidMapping_test.cpp
index 46a53cb53ef..630add17ad0 100644
--- a/Sofa/Component/Mapping/NonLinear/tests/RigidMapping_test.cpp
+++ b/Sofa/Component/Mapping/NonLinear/tests/RigidMapping_test.cpp
@@ -90,7 +90,7 @@ struct RigidMappingTest : public sofa::mapping_test::Mapping_test<_RigidMapping>
         if constexpr ( InDataTypes::spatial_dimensions != 3 )
         {
             // RigidMapping::getK is not yet implemented for 2D rigids
-            this->flags &= ~Inherit::TEST_getK;
+            this->setTestExecution(Inherit::TEST_getK, false);
         }
     }
 
diff --git a/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MappingTestCreation.h b/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MappingTestCreation.h
index 3a430a8deed..075524fe651 100644
--- a/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MappingTestCreation.h
+++ b/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MappingTestCreation.h
@@ -22,6 +22,7 @@
 #pragma once
 
 #include <sofa/core/BaseLocalMappingMatrix.h>
+#include <sofa/simpleapi/SimpleApi.h>
 #include <sofa/testing/BaseSimulationTest.h>
 using sofa::testing::BaseSimulationTest;
 
@@ -30,7 +31,6 @@ using sofa::testing::NumericTest;
 
 #include <sofa/simulation/VectorOperations.h>
 
-#include <sofa/linearalgebra/FullVector.h>
 #include <sofa/linearalgebra/EigenSparseMatrix.h>
 #include <sofa/component/statecontainer/MechanicalObject.h>
 #include <sofa/simulation/graph/DAGSimulation.h>
@@ -38,16 +38,17 @@ using sofa::testing::NumericTest;
 
 #include <sofa/helper/logging/Messaging.h>
 
-namespace sofa::mapping_test {
+namespace sofa::mapping_test
+{
 
 
 /** @brief Base class for the Mapping tests, with helpers to automatically test applyJ, applyJT, applyDJT and getJs using finite differences.
 
   Specific test cases can be created using a derived class instantiated on the mapping class to test,
-  and calling function runTest( const InVecCoord& parentInit,
-                  const OutVecCoord& childInit,
-                  const InVecCoord parentNew,
-                  const OutVecCoord expectedChildNew);
+  and calling function runTest( const VecCoord_t<In>& parentInit,
+                  const VecCoord_t<Out>& childInit,
+                  const VecCoord_t<In> parentNew,
+                  const VecCoord_t<Out> expectedChildNew);
 
 
   This function compares the actual output positions with the expected ones, then automatically tests the methods related to
@@ -65,41 +66,16 @@ namespace sofa::mapping_test {
 
   @author François Faure @date 2013
   */
-
 template< class _Mapping>
 struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::In::Real>
 {
-    typedef _Mapping Mapping;
-    typedef typename Mapping::In In;
-    typedef component::statecontainer::MechanicalObject<In> InDOFs;
-    typedef typename InDOFs::Real  Real;
-    typedef typename InDOFs::Coord  InCoord;
-    typedef typename InDOFs::Deriv  InDeriv;
-    typedef typename InDOFs::VecCoord  InVecCoord;
-    typedef typename InDOFs::VecDeriv  InVecDeriv;
-    typedef typename InDOFs::ReadVecCoord  ReadInVecCoord;
-    typedef typename InDOFs::WriteVecCoord WriteInVecCoord;
-    typedef typename InDOFs::ReadVecDeriv  ReadInVecDeriv;
-    typedef typename InDOFs::WriteVecDeriv WriteInVecDeriv;
-    typedef typename InDOFs::MatrixDeriv  InMatrixDeriv;
-    typedef Data<InVecCoord> InDataVecCoord;
-    typedef Data<InVecDeriv> InDataVecDeriv;
-    typedef Data<InMatrixDeriv> InDataMatrixDeriv;
-
-    typedef typename Mapping::Out Out;
-    typedef component::statecontainer::MechanicalObject<Out> OutDOFs;
-    typedef typename OutDOFs::Coord     OutCoord;
-    typedef typename OutDOFs::Deriv     OutDeriv;
-    typedef typename OutDOFs::VecCoord  OutVecCoord;
-    typedef typename OutDOFs::VecDeriv  OutVecDeriv;
-    typedef typename OutDOFs::ReadVecCoord  ReadOutVecCoord;
-    typedef typename OutDOFs::WriteVecCoord WriteOutVecCoord;
-    typedef typename OutDOFs::ReadVecDeriv  ReadOutVecDeriv;
-    typedef typename OutDOFs::WriteVecDeriv WriteOutVecDeriv;
-    typedef typename OutDOFs::MatrixDeriv  OutMatrixDeriv;
-    typedef Data<OutVecCoord> OutDataVecCoord;
-    typedef Data<OutVecDeriv> OutDataVecDeriv;
-    typedef Data<OutMatrixDeriv> OutDataMatrixDeriv;
+    using Mapping = _Mapping;
+
+    using In = typename Mapping::In;
+    using InDOFs = component::statecontainer::MechanicalObject<In>;
+
+    using Out = typename Mapping::Out;
+    using OutDOFs = component::statecontainer::MechanicalObject<Out>;
 
     typedef linearalgebra::EigenSparseMatrix<In,Out> EigenSparseMatrix;
 
@@ -107,10 +83,9 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
     typename InDOFs::SPtr  inDofs;  ///< mapping input
     typename OutDOFs::SPtr outDofs; ///< mapping output
     simulation::Node::SPtr root;         ///< Root of the scene graph, created by the constructor an re-used in the tests
-    simulation::Simulation* simulation;  ///< created by the constructor an re-used in the tests
-    std::pair<Real,Real> deltaRange; ///< The minimum and maximum magnitudes of the change of each scalar value of the small displacement is perturbation * numeric_limits<Real>::epsilon. This epsilon is 1.19209e-07 for float and 2.22045e-16 for double.
-    Real errorMax;     ///< The test is successful if the (infinite norm of the) difference is less than  errorMax * numeric_limits<Real>::epsilon
-    Real errorFactorDJ;     ///< The test for geometric stiffness is successful if the (infinite norm of the) difference is less than  errorFactorDJ * errorMax * numeric_limits<Real>::epsilon
+    std::pair<Real_t<In>, Real_t<In>> deltaRange; ///< The minimum and maximum magnitudes of the change of each scalar value of the small displacement is perturbation * numeric_limits<Real>::epsilon. This epsilon is 1.19209e-07 for float and 2.22045e-16 for double.
+    Real_t<In> errorMax;     ///< The test is successful if the (infinite norm of the) difference is less than  errorMax * numeric_limits<Real>::epsilon
+    Real_t<In> errorFactorDJ;     ///< The test for geometric stiffness is successful if the (infinite norm of the) difference is less than  errorFactorDJ * errorMax * numeric_limits<Real>::epsilon
 
 
     static constexpr unsigned char TEST_getJs = 1; ///< testing getJs used in assembly API
@@ -122,6 +97,25 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
     static constexpr unsigned char TEST_GEOMETRIC_STIFFNESS = TEST_applyDJT | TEST_getK | TEST_buildGeometricStiffnessMatrix; ///< testing functions used in assembly API getJS getKS
     unsigned char flags; ///< testing options. (all by default). To be used with precaution. Please implement the missing API in the mapping rather than not testing it.
 
+    /// Returns true if the test of the method corresponding to the flag parameter must be executed, false otherwise
+    bool isTestExecuted(const unsigned char testFlag) const
+    {
+        return flags & testFlag;
+    }
+
+    /// Defines if a test of the method corresponding to the flag parameter must be executed or not.
+    void setTestExecution(const unsigned char testFlag, const bool doTheTest)
+    {
+        if (doTheTest)
+        {
+            flags |= testFlag;  // Set the flag
+        }
+        else
+        {
+            flags &= ~testFlag; // Unset the flag
+        }
+    }
+
 
     Mapping_test()
         : deltaRange(1,1000)
@@ -129,10 +123,8 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
         , errorFactorDJ(1)
         , flags(TEST_ASSEMBLY_API | TEST_GEOMETRIC_STIFFNESS)
     {
-        simulation = sofa::simulation::getSimulation();
-
         /// Parent node
-        root = simulation->createNewGraph("root");
+        root = simpleapi::createRootNode(simulation::getSimulation(), "root");
 
         inDofs = core::objectmodel::New<InDOFs>();
         root->addObject(inDofs);
@@ -149,16 +141,15 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
     }
 
     Mapping_test(std::string fileName)
-        : simulation(sofa::simulation::getSimulation()),
-          deltaRange(1, 1000),
+        : deltaRange(1, 1000),
           errorMax(100),
           errorFactorDJ(1),
           flags(TEST_ASSEMBLY_API|TEST_GEOMETRIC_STIFFNESS)
     {
-        assert(simulation);
+        assert(simulation::getSimulation());
 
         /// Load the scene
-        root = simulation->createNewGraph("root");
+        root = simpleapi::createRootNode(simulation::getSimulation(), "root");
         root = sofa::simulation::node::load(fileName.c_str(), false);
 
         // InDofs
@@ -180,34 +171,8 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
         
     }
 
-
-    /** Returns OutCoord subtraction a-b */
-    virtual OutDeriv difference( const OutCoord& a, const OutCoord& b )
-    {
-        return Out::coordDifference(a,b);
-    }
-
-    virtual OutVecDeriv difference( const OutVecDeriv& a, const OutVecDeriv& b )
-    {
-        if( a.size()!=b.size() ){
-            ADD_FAILURE() << "OutVecDeriv have different sizes";
-            return OutVecDeriv();
-        }
-
-        OutVecDeriv c(a.size());
-        for (size_t i=0; i<a.size() ; ++i)
-        {
-            c[i] = a[i]-b[i];
-        }
-        return c;
-    }
-
-    /** Possible child force pre-treatment, does nothing by default
-      */
-    virtual OutVecDeriv preTreatment( const OutVecDeriv& f ) { return f; }
-
-
-    /** Test the mapping using the given values and small changes.
+    /**
+     * Test the mapping using the given values and small changes.
      * Return true in case of success, if all errors are below maxError*epsilon.
      * The mapping is initialized using the two first parameters,
      * then a new parent position is applied,
@@ -222,149 +187,73 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
      *\param parentNew new parent position
      *\param expectedChildNew expected position of the child corresponding to the new parent position
      */
-    virtual bool runTest( const InVecCoord& parentInit,
-                          const OutVecCoord& childInit,
-                          const InVecCoord& parentNew,
-                          const OutVecCoord& expectedChildNew)
+    virtual bool runTest( const VecCoord_t<In>& parentInit,
+                          const VecCoord_t<Out>& childInit,
+                          const VecCoord_t<In>& parentNew,
+                          const VecCoord_t<Out>& expectedChildNew)
     {
-        if( deltaRange.second / errorMax <= sofa::testing::g_minDeltaErrorRatio )
-            ADD_FAILURE() << "The comparison threshold is too large for the finite difference delta";
-
-        if( !(flags & TEST_getJs) )          msg_warning("MappingTest") << "getJs is not tested";
-        if( !(flags & TEST_getK) )           msg_warning("MappingTest") << "getK is not tested";
-        if( !(flags & TEST_applyJT_matrix) ) msg_warning("MappingTest") << "applyJT on matrices is not tested";
-        if( !(flags & TEST_applyDJT) )       msg_warning("MappingTest") << "applyDJT is not tested";
+        checkComparisonThreshold();
+        warnMissingTests();
 
+        const auto errorThreshold = this->epsilon() * errorMax;
 
-        const Real errorThreshold = this->epsilon()*errorMax;
+        using EigenSparseMatrix = linearalgebra::EigenSparseMatrix<In, Out>;
 
-        typedef linearalgebra::EigenSparseMatrix<In,Out> EigenSparseMatrix;
         core::MechanicalParams mparams;
         mparams.setKFactor(1.0);
         mparams.setSupportOnlySymmetricMatrix(false);
-        inDofs->resize(parentInit.size());
-        WriteInVecCoord xin = inDofs->writePositions();
-        sofa::testing::copyToData(xin,parentInit); // xin = parentInit
-
-        outDofs->resize(childInit.size());
-        WriteOutVecCoord xout = outDofs->writePositions();
-        sofa::testing::copyToData(xout,childInit);
-
-        /// Init based on parentInit
-        sofa::simulation::node::initRoot(root.get());
-
-        /// Updated to parentNew
-        sofa::testing::copyToData(xin,parentNew);
-        mapping->apply(&mparams, core::VecCoordId::position(), core::VecCoordId::position());
-        mapping->applyJ(&mparams, core::VecDerivId::velocity(), core::VecDerivId::velocity());
 
         /// test apply: check if the child positions are the expected ones
-        bool succeed=true;
+        bool succeed = testMappingPositionVelocity(
+            parentInit, childInit, parentNew, expectedChildNew, errorThreshold,
+            mparams);
 
-        if (expectedChildNew.size() != xout.size()) {
-            ADD_FAILURE() << "Size of output dofs is wrong: " << xout.size() << " expected: " << expectedChildNew.size();
-            succeed = false;
-        }
-        for( unsigned i=0; i<xout.size(); i++ )
-        {
-            if( !this->isSmall( difference(xout[i],expectedChildNew[i]).norm(), errorMax ) ) {
-                ADD_FAILURE() << "Position of mapped particle " << i << " is wrong: \n" << xout[i] <<"\nexpected: \n" << expectedChildNew[i]
-                              <<  "\ndifference should be less than " << errorThreshold << " (" << difference(xout[i],expectedChildNew[i]).norm() << ")" << std::endl;
-                succeed = false;
-            }
-        }
-
-        /// test applyJ and everything related to Jacobians
-        size_t Np = inDofs->getSize(), Nc=outDofs->getSize();
-
-        InVecCoord xp(Np),xp1(Np);
-        InVecDeriv vp(Np),fp(Np);
-        InVecDeriv dfp_a(Np);
-        InVecDeriv dfp_b(Np);
-        InVecDeriv fp2(Np);
-        OutVecCoord xc(Nc),xc1(Nc);
-        OutVecDeriv vc(Nc),fc(Nc);
+        const std::size_t sizeIn = inDofs->getSize();
+        const std::size_t sizeOut = outDofs->getSize();
 
         // get position data
-        sofa::testing::copyFromData( xp, inDofs->readPositions() );
-        sofa::testing::copyFromData( xc, outDofs->readPositions() ); // positions and have already been propagated
+        VecCoord_t<Out> positionOut = outDofs->readPositions().ref();
 
         // set random child forces and propagate them to the parent
-        for( unsigned i=0; i<Nc; i++ ){
-            fc[i] = Out::randomDeriv( 0.1, 1.0 );
-        }
-        fp2.fill( InDeriv() );
-        WriteInVecDeriv fin = inDofs->writeForces();
-        sofa::testing::copyToData( fin, fp2 );  // reset parent forces before accumulating child forces
+        VecDeriv_t<Out> forceOut = generateRandomVecDeriv<Out>(sizeOut, 0.1, 1.);
 
-        WriteOutVecDeriv fout = outDofs->writeForces();
-        sofa::testing::copyToData( fout, fc );
-        mapping->applyJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
-        sofa::testing::copyFromData( fp, inDofs->readForces() );
+        VecDeriv_t<In> forceIn;
+        computeForceInFromForceOut(mparams, forceIn, forceOut);
 
         // set small parent velocities and use them to update the child
-        for( unsigned i=0; i<Np; i++ ){
-            vp[i] = In::randomDeriv(this->epsilon() * deltaRange.first, this->epsilon() * deltaRange.second );
-        }
+        const VecDeriv_t<In> velocityIn = generateRandomVecDeriv<In>(sizeIn,
+            this->epsilon() * deltaRange.first,
+            this->epsilon() * deltaRange.second);
 
-        for( unsigned i=0; i<Np; i++ ){             // and small displacements
-            xp1[i] = xp[i] + vp[i];
-        }
+        const VecCoord_t<In> perturbedPositionIn = computePerturbedPositions(sizeIn, velocityIn);
 
-        // propagate small velocity
-        WriteInVecDeriv vin = inDofs->writeVelocities();
-        sofa::testing::copyToData( vin, vp );
-        mapping->applyJ( &mparams, core::VecDerivId::velocity(), core::VecDerivId::velocity() );
-        ReadOutVecDeriv vout = outDofs->readVelocities();
-        sofa::testing::copyFromData( vc, vout);
+        VecDeriv_t<Out> velocityOut;
+        computeVelocityOutFromVelocityIn(mparams, velocityOut, velocityIn);
 
         // apply geometric stiffness
         inDofs->vRealloc( &mparams, core::VecDerivId::dx() ); // dx is not allocated by default
-        WriteInVecDeriv dxin = inDofs->writeDx();
-        sofa::testing::copyToData( dxin, vp );
-        dfp_a.fill( InDeriv() );
-        dfp_b.fill( InDeriv() );
-        sofa::testing::copyToData( fin, dfp_a );
-        sofa::testing::copyToData( fin, dfp_b );
-
-        //calling applyDJT without updateK before
-        mapping->applyDJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
-        sofa::testing::copyFromData( dfp_a, inDofs->readForces() ); // fp + df due to geometric stiffness
-        sofa::testing::copyToData( fin, fp2 ); //reset forces
+        inDofs->writeDx().wref() = velocityIn;
 
-        //calling applyDJT after updateK
-        mapping->updateK( &mparams, core::ConstVecDerivId::force() ); // updating stiffness matrix for the current state and force
-        mapping->applyDJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
-        sofa::testing::copyFromData( dfp_b, inDofs->readForces() ); // fp + df due to geometric stiffness
+        const VecDeriv_t<In> dfp_withoutUpdateK = applyDJT(mparams, false);
+        const VecDeriv_t<In> dfp_withUpdateK = applyDJT(mparams, true);
 
         // Jacobian will be obsolete after applying new positions
-        if( flags & TEST_getJs )
+        if( isTestExecuted(TEST_getJs) )
         {
             EigenSparseMatrix* J = this->getMatrix<EigenSparseMatrix>(mapping->getJs());
-            //        cout<<"J = "<< endl << *J << endl;
-            OutVecDeriv Jv(Nc);
-            J->mult(Jv,vp);
-
-            // ================ test applyJT()
-            InVecDeriv jfc( (long)Np,InDeriv());
-            J->addMultTranspose(jfc,fc);
-            if( this->vectorMaxDiff(jfc,fp)>errorThreshold ){
-                succeed = false;
-                ADD_FAILURE() << "applyJT test failed, difference should be less than " << errorThreshold << " (" << this->vectorMaxDiff(jfc,fp) << ")" << std::endl
-                              << "jfc = " << jfc << std::endl<<" fp = " << fp << std::endl;
-            }
-            // ================ test getJs()
-            // check that J.vp = vc
-            if(this->vectorMaxDiff(Jv,vc)>errorThreshold ){
-                succeed = false;
-                std::cout<<"vp = " << vp << std::endl;
-                std::cout<<"Jvp = " << Jv << std::endl;
-                std::cout<<"vc  = " << vc << std::endl;
-                ADD_FAILURE() << "getJs() test failed"<<std::endl<<"vp = " << vp << std::endl<<"Jvp = " << Jv << std::endl <<"vc  = " << vc << std::endl;
-            }
+
+            // forceIn has been computed using applyJT
+            // The following tests that forceIn is also the result of applying
+            // the transposed jacobian matrix
+            succeed &= checkJacobianMatrixTranspose(J, forceOut, forceIn, errorThreshold);
+
+            // velocityOut has been computed using applyJ
+            // The following tests that velocityOut is also the result of applying
+            // the jacobian matrix
+            succeed &= checkJacobianMatrix(J, velocityIn, velocityOut, errorThreshold);
         }
 
-        if( flags & TEST_applyJT_matrix )
+        if( isTestExecuted(TEST_applyJT_matrix) )
         {
             // TODO test applyJT on matrices
             // basic idea build a out random matrix  e.g. with helper::drand(100.0)
@@ -382,201 +271,439 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
         }
 
         // compute parent forces from pre-treated child forces (in most cases, the pre-treatment does nothing)
-        // the pre-treatement can be useful to be able to compute 2 comparable results of applyJT with a small displacement to test applyDJT
-        fp.fill( InDeriv() );
-        sofa::testing::copyToData( fin, fp );  // reset parent forces before accumulating child forces
-        sofa::testing::copyToData( fout, preTreatment(fc) );
-        mapping->applyJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
-        sofa::testing::copyFromData( fp, inDofs->readForces() );
-
-
+        // the pre-treatment can be useful to be able to compute 2 comparable results of applyJT with a small displacement to test applyDJT
+        computeForceInFromForceOut(mparams, forceIn, preTreatment(forceOut));
 
         ///////////////////////////////
 
-
         // propagate small displacement
-        WriteInVecCoord pin (inDofs->writePositions());
-        sofa::testing::copyToData( pin, xp1 );
+        inDofs->writePositions().wref() = perturbedPositionIn;
+        succeed &= testApplyJonPosition(mparams, positionOut, velocityOut, errorThreshold);
 
-        mapping->apply ( &mparams, core::VecCoordId::position(), core::VecCoordId::position() );
-        ReadOutVecCoord pout = outDofs->readPositions();
-        sofa::testing::copyFromData( xc1, pout );
+        const VecDeriv_t<In> forceChange = computeForceChange(mparams, sizeIn, forceOut, forceIn);
 
-        // ================ test applyJ: compute the difference between propagated displacements and velocities
-        OutVecDeriv dxc(Nc);
-        for(unsigned i=0; i<Nc; i++ ){
-            dxc[i] = difference( xc1[i], xc[i] );
+        if( isTestExecuted(TEST_applyDJT) )
+        {
+            succeed &= checkApplyDJT(dfp_withoutUpdateK, forceChange, errorThreshold, false);
+            succeed &= checkApplyDJT(dfp_withUpdateK, forceChange, errorThreshold, true);
         }
 
-        if(this->vectorMaxAbs(difference(dxc,vc))>errorThreshold ){
-            succeed = false;
-            ADD_FAILURE() << "applyJ test failed: the difference between child position change and child velocity (dt=1) "<< this->vectorMaxAbs(difference(dxc,vc))<<" should be less than  " << errorThreshold << std::endl
-                          << "position change = " << dxc << std::endl
-                          << "velocity        = " << vc << std::endl;
+        if( isTestExecuted(TEST_getK))
+        {
+            succeed &= testGetK(sizeIn, velocityIn, forceChange, errorThreshold);
         }
 
+        if( isTestExecuted(TEST_buildGeometricStiffnessMatrix) )
+        {
+            succeed &= testBuildGeometricStiffnessMatrix(sizeIn, velocityIn, forceChange, errorThreshold);
+        }
 
+        if(!succeed)
+        {
+            ADD_FAILURE() << "Failed Seed number = " << this->seed << std::endl;
+        }
+        return succeed;
+    }
 
-        // update parent force based on the same child forces
-        fp2.fill( InDeriv() );
-        sofa::testing::copyToData( fin, fp2 );  // reset parent forces before accumulating child forces
-        sofa::testing::copyToData( fout, preTreatment(fc) );
-        mapping->applyJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
-        sofa::testing::copyFromData( fp2, inDofs->readForces() );
 
-        InVecDeriv fp12(Np);
-        for(unsigned i=0; i<Np; i++){
-            fp12[i] = fp2[i]-fp[i];       // fp2 - fp
+    /** Test the mapping using the given values and small changes.
+     * Return true in case of success, if all errors are below maxError*epsilon.
+     * The mapping is initialized using the first parameter,
+     * @warning this version supposes the mapping initialization does not depend on child positions
+     * otherwise, use the runTest functions with 4 parameters
+     * the child position is computed from parent position and compared with the expected one.
+     * Additionally, the Jacobian-related methods are tested using finite differences.
+     *
+     *\param parent parent position
+     *\param expectedChild expected position of the child corresponding to the parent position
+     */
+    virtual bool runTest( const VecCoord_t<In>& parent,
+                          const VecCoord_t<Out> expectedChild)
+    {
+        VecCoord_t<Out> childInit( expectedChild.size() ); // start with null child
+        return runTest( parent, childInit, parent, expectedChild );
+    }
+
+    ~Mapping_test() override
+    {
+        if (root != nullptr)
+        {
+            sofa::simulation::node::unload(root);
         }
+    }
+
+protected:
+
 
 
-        // ================ test applyDJT()
-        if( flags & TEST_applyDJT )
+    /** Returns OutCoord substraction a-b */
+    virtual Deriv_t<Out> difference( const Coord_t<Out>& a, const Coord_t<Out>& b )
+    {
+        return Out::coordDifference(a,b);
+    }
+
+    virtual VecDeriv_t<Out> difference( const VecDeriv_t<Out>& a, const VecDeriv_t<Out>& b )
+    {
+        if (a.size() != b.size())
         {
-            //we test the function applyDJT in two different contexts because
-            //the implementation usually depends on the context
-            if (this->vectorMaxDiff(dfp_a,fp12) > errorThreshold * errorFactorDJ)
-            {
-                succeed = false;
-                ADD_FAILURE() << "applyDJT (no call to updateK) test failed" << std::endl
-                    << "dfp    = " << dfp_a << std::endl
-                    << "fp2-fp = " << fp12 << std::endl
-                    << "error threshold = " << errorThreshold * errorFactorDJ << std::endl;
-            }
-            if (this->vectorMaxDiff(dfp_b, fp12) > errorThreshold * errorFactorDJ)
-            {
-                succeed = false;
-                ADD_FAILURE() << "applyDJT (after call to updateK) test failed" << std::endl
-                    << "dfp    = " << dfp_b << std::endl
-                    << "fp2-fp = " << fp12 << std::endl
-                    << "error threshold = " << errorThreshold * errorFactorDJ << std::endl;
-            }
+            ADD_FAILURE() << "VecDeriv_t<Out> have different sizes";
+            return {};
         }
 
+        VecDeriv_t<Out> c;
+        c.reserve(a.size());
+        std::transform(a.begin(), a.end(), b.begin(), std::back_inserter(c), std::minus());
 
-        // ================ test getK()
-        if( flags & TEST_getK )
+        return c;
+    }
+
+    /** Possible child force pre-treatment, does nothing by default
+      */
+    virtual VecDeriv_t<Out> preTreatment( const VecDeriv_t<Out>& f )
+    {
+        return f;
+    }
+
+
+    void checkComparisonThreshold()
+    {
+        if( deltaRange.second / errorMax <= sofa::testing::g_minDeltaErrorRatio )
         {
-            InVecDeriv Kv(Np);
+            ADD_FAILURE() << "The comparison threshold is too large for the finite difference delta";
+        }
+    }
 
-            const linearalgebra::BaseMatrix* bk = mapping->getK();
+    void warnMissingTests() const
+    {
+        msg_warning_when(!isTestExecuted(TEST_getJs), "MappingTest") << "getJs is not tested";
+        msg_warning_when(!isTestExecuted(TEST_getK)           , "MappingTest") << "getK is not tested";
+        msg_warning_when(!isTestExecuted(TEST_applyJT_matrix) , "MappingTest") << "applyJT on matrices is not tested";
+        msg_warning_when(!isTestExecuted(TEST_applyDJT)       , "MappingTest") << "applyDJT is not tested";
+        msg_warning_when(!isTestExecuted(TEST_buildGeometricStiffnessMatrix) , "MappingTest") << "buildGeometricStiffnessMatrix is not tested";
+    }
 
-            // K can be null or empty for linear mappings
-            // still performing the test with a null Kv vector to check if the mapping is really linear
+    bool testMappingPositionVelocity(const VecCoord_t<In>& parentInit,
+                                     const VecCoord_t<Out>& childInit,
+                                     const VecCoord_t<In>& parentNew,
+                                     const VecCoord_t<Out>& expectedChildNew,
+                                     const Real_t<In> errorThreshold,
+                                     core::MechanicalParams mparams)
+    {
+        helper::WriteAccessor positionAccessorIn = inDofs->writePositions();
+        positionAccessorIn.wref() = parentInit;
 
-            if( bk != nullptr ){
+        helper::WriteAccessor positionAccessorOut = outDofs->writePositions();
+        positionAccessorOut.wref() = childInit;
 
-                typedef linearalgebra::EigenSparseMatrix<In,In> EigenSparseKMatrix;
-                const EigenSparseKMatrix* K = dynamic_cast<const EigenSparseKMatrix*>(bk);
-                if( K == nullptr ){
-                    succeed = false;
-                    ADD_FAILURE() << "getK returns a matrix of non-EigenSparseMatrix type";
-                    // TODO perform a slow conversion with a big warning rather than a failure?
-                }
+        /// Init based on parentInit
+        sofa::simulation::node::initRoot(root.get());
 
-                if( K->compressedMatrix.nonZeros() ) K->mult(Kv,vp);
-            }
+        /// Updated to parentNew
+        positionAccessorIn.wref() = parentNew;
+        mapping->apply(&mparams, core::VecCoordId::position(), core::VecCoordId::position());
+        mapping->applyJ(&mparams, core::VecDerivId::velocity(), core::VecDerivId::velocity());
 
-            // check that K.vp = dfp
-            if(this->vectorMaxDiff(Kv,fp12)>errorThreshold*errorFactorDJ ){
-                succeed = false;
-                ADD_FAILURE() << "K test failed, difference should be less than " << errorThreshold*errorFactorDJ  << std::endl
-                              << "Kv    = " << Kv << std::endl
-                              << "dfp = " << fp12 << std::endl;
-            }
-        }
+        bool succeed = true;
 
-        if( flags & TEST_buildGeometricStiffnessMatrix )
+        if (expectedChildNew.size() != positionAccessorOut.size())
         {
-            core::GeometricStiffnessMatrix testGeometricStiffness;
-
-            struct GeometricStiffnessAccumulator : core::MappingMatrixAccumulator
-            {
-                void add(sofa::SignedIndex row, sofa::SignedIndex col, float value) override
-                {
-                    assembledMatrix.add(row, col, value);
-                }
-                void add(sofa::SignedIndex row, sofa::SignedIndex col, double value) override
-                {
-                    assembledMatrix.add(row, col, value);
-                }
-
-                linearalgebra::EigenSparseMatrix<In,In> assembledMatrix;
-            } accumulator;
-
-            testGeometricStiffness.setMatrixAccumulator(&accumulator, mapping->getFromModel(), mapping->getFromModel());
-
-            accumulator.assembledMatrix.resize(mapping->getFromModel()->getSize() * In::deriv_total_size, mapping->getFromModel()->getSize() * In::deriv_total_size);
-            mapping->buildGeometricStiffnessMatrix(&testGeometricStiffness);
-            accumulator.assembledMatrix.compress();
-
-            InVecDeriv Kv(Np);
-            if( accumulator.assembledMatrix.compressedMatrix.nonZeros() )
-            {
-                accumulator.assembledMatrix.mult(Kv,vp);
-            }
+            ADD_FAILURE() << "Size of output dofs is wrong: " << positionAccessorOut.size() << " expected: " << expectedChildNew.size();
+            succeed = false;
+        }
 
-            // check that K.vp = dfp
-            if (this->vectorMaxDiff(Kv,fp12) > errorThreshold*errorFactorDJ )
+        for (unsigned i = 0; i < positionAccessorOut.size(); ++i)
+        {
+            if (!this->isSmall(difference(positionAccessorOut[i], expectedChildNew[i]).norm(), errorMax))
             {
+                ADD_FAILURE() << "Position of mapped particle " << i << " is wrong: \n" << positionAccessorOut[i] <<"\nexpected: \n" << expectedChildNew[i]
+                        <<  "\ndifference should be less than " << errorThreshold << " (" << difference(positionAccessorOut[i],expectedChildNew[i]).norm() << ")" << std::endl;
                 succeed = false;
-                ADD_FAILURE() << "buildGeometricStiffnessMatrix test failed, difference should be less than " << errorThreshold*errorFactorDJ  << std::endl
-                              << "Kv    = " << Kv << std::endl
-                              << "dfp = " << fp12 << std::endl;
             }
         }
 
-        if(!succeed)
-        { ADD_FAILURE() << "Failed Seed number = " << this->seed << std::endl;}
         return succeed;
     }
 
+    template<class DataTypes>
+    VecDeriv_t<DataTypes> generateRandomVecDeriv(const std::size_t size, const Real_t<DataTypes> minMagnitude, const Real_t<DataTypes> maxMagnitude)
+    {
+        VecDeriv_t<DataTypes> randomForce;
+        randomForce.reserve(size);
+        for (std::size_t i = 0; i < size; i++)
+        {
+            randomForce.push_back(DataTypes::randomDeriv(minMagnitude, maxMagnitude));
+        }
+        return randomForce;
+    }
 
-    /** Test the mapping using the given values and small changes.
-     * Return true in case of success, if all errors are below maxError*epsilon.
-     * The mapping is initialized using the first parameter,
-     * @warning this version supposes the mapping initialization does not depend on child positions
-     * otherwise, use the runTest functions with 4 parameters
-     * the child position is computed from parent position and compared with the expected one.
-     * Additionally, the Jacobian-related methods are tested using finite differences.
-     *
-     *\param parent parent position
-     *\param expectedChild expected position of the child corresponding to the parent position
-     */
-    virtual bool runTest( const InVecCoord& parent,
-                          const OutVecCoord expectedChild)
+    void computeForceInFromForceOut(core::MechanicalParams mparams, VecDeriv_t<In>& forceIn, const VecDeriv_t<Out>& forceOut)
     {
-        OutVecCoord childInit( expectedChild.size() ); // start with null child
-        return runTest( parent, childInit, parent, expectedChild );
+        inDofs->writeForces()->fill(Deriv_t<In>());  // reset parent forces before accumulating child forces
+
+        outDofs->writeForces().wref() = forceOut;
+        mapping->applyJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
+        forceIn = inDofs->readForces().ref();
     }
 
-    ~Mapping_test() override
+    void computeVelocityOutFromVelocityIn(core::MechanicalParams mparams, VecDeriv_t<Out>& velocityOut, const VecDeriv_t<In>& velocityIn)
     {
-        if (root!=nullptr)
-            sofa::simulation::node::unload(root);
+        inDofs->writeVelocities().wref() = velocityIn;
+        mapping->applyJ( &mparams, core::VecDerivId::velocity(), core::VecDerivId::velocity() );
+        velocityOut = outDofs->readVelocities().ref();
+    }
+
+    const VecDeriv_t<In>& applyDJT(core::MechanicalParams mparams, bool updateK)
+    {
+        inDofs->writeForces()->fill(Deriv_t<In>()); //reset force
+
+        if (updateK)
+        {
+            mapping->updateK( &mparams, core::ConstVecDerivId::force() ); // updating stiffness matrix for the current state and force
+        }
+        mapping->applyDJT( &mparams, core::VecDerivId::force(), core::VecDerivId::force() );
+        return inDofs->readForces().ref();
+    }
+
+    [[nodiscard]] bool checkApplyDJT(const VecDeriv_t<In>& dfp, const VecDeriv_t<In>& fp12, Real_t<In> errorThreshold, bool updateK)
+    {
+        if (this->vectorMaxDiff(dfp, fp12) > errorThreshold * errorFactorDJ)
+        {
+            const std::string updateKString = updateK ? "after call to updateK" : "no call to updateK";
+            ADD_FAILURE() << "applyDJT (" << updateKString << ") test failed" << std::endl
+                << "dfp    = " << dfp << std::endl
+                << "fp2-fp = " << fp12 << std::endl
+                << "error threshold = " << errorThreshold * errorFactorDJ << std::endl;
+            return false;
+        }
+        return true;
+    }
+
+    VecCoord_t<In> computePerturbedPositions(const std::size_t sizeIn, const VecDeriv_t<In> velocityIn)
+    {
+        VecCoord_t<In> perturbedPositionIn;
+        perturbedPositionIn.reserve(sizeIn);
+        const helper::ReadAccessor positionIn = inDofs->readPositions();
+        for (std::size_t i = 0; i < sizeIn; ++i)
+        {
+            perturbedPositionIn.push_back(positionIn[i] + velocityIn[i]);
+        }
+        return perturbedPositionIn;
     }
 
-protected:
 
     /// Get one EigenSparseMatrix out of a list. Error if not one single matrix in the list.
     template<class EigenSparseMatrixType>
     static EigenSparseMatrixType* getMatrix(const type::vector<sofa::linearalgebra::BaseMatrix*>* matrices)
     {
-        if( !matrices ){
+        if( !matrices )
+        {
             ADD_FAILURE()<< "Matrix list is nullptr (API for assembly is not implemented)";
         }
-        if( matrices->size() != 1 ){
-            ADD_FAILURE()<< "Matrix list should have size == 1 in simple mappings";
+
+        if( matrices->empty() )
+        {
+            ADD_FAILURE()<< "Matrix list is empty";
+            return nullptr;
+        }
+
+        if( matrices->size() != 1 )
+        {
+            ADD_FAILURE()<< "Matrix list should have size == 1 in simple mappings (current size = " << matrices->size() << ")";
         }
         EigenSparseMatrixType* ei = dynamic_cast<EigenSparseMatrixType*>((*matrices)[0] );
-        if( ei == nullptr ){
+        if( ei == nullptr )
+        {
             ADD_FAILURE() << "getJs returns a matrix of non-EigenSparseMatrix type";
             // TODO perform a slow conversion with a big warning rather than a failure?
         }
         return ei;
     }
 
+    bool checkJacobianMatrixTranspose(
+        EigenSparseMatrix* jacobianMatrix,
+        const VecDeriv_t<Out>& forceOut,
+        const VecDeriv_t<In>& expectedForceIn,
+        Real_t<In> errorThreshold)
+    {
+        VecDeriv_t<In> computedForceIn(expectedForceIn.size(), Deriv_t<In>());
+
+        //computedForceIn += J^T * forceOut
+        jacobianMatrix->addMultTranspose(computedForceIn, forceOut);
 
+        const auto diff = this->vectorMaxDiff(computedForceIn, expectedForceIn);
+        if (diff > errorThreshold)
+        {
+            ADD_FAILURE() <<
+                "getJs is not consistent with applyJT, difference should be "
+                "less than " << errorThreshold << " (" << diff << ")" << std::endl
+                << "computedForceIn = " << computedForceIn << std::endl
+                << "expectedForceIn = " << expectedForceIn << std::endl;
+            return false;
+        }
+        return true;
+    }
+
+    bool checkJacobianMatrix(
+        EigenSparseMatrix* jacobianMatrix,
+        const VecDeriv_t<In>& velocityIn,
+        const VecDeriv_t<Out>& expectedVelocityOut,
+        Real_t<In> errorThreshold)
+    {
+        VecDeriv_t<Out> computedVelocityOut(expectedVelocityOut.size());
+
+        jacobianMatrix->mult(computedVelocityOut, velocityIn);
+
+        const auto diff = this->vectorMaxDiff(computedVelocityOut, expectedVelocityOut);
+        if (diff > errorThreshold)
+        {
+            ADD_FAILURE() <<
+                "getJs is not consistent with applyJ, difference should be "
+                "less than " << errorThreshold << " (" << diff << ")" << std::endl
+                << "velocityIn = " << velocityIn << std::endl
+                << "computedVelocityOut = " << computedVelocityOut << std::endl
+                << "expectedVelocityOut = " << expectedVelocityOut << std::endl;
+            return false;
+        }
+
+        return true;
+    }
+
+    bool testApplyJonPosition(
+        core::MechanicalParams mparams,
+        const VecCoord_t<Out>& positionOut,
+        const VecDeriv_t<Out>& expectedVelocityOut,
+        Real_t<In> errorThreshold
+        )
+    {
+        mapping->apply ( &mparams, core::VecCoordId::position(), core::VecCoordId::position() );
+        const VecCoord_t<Out>& positionOut1 = outDofs->readPositions();
+
+        const auto sizeOut = positionOut.size();
+        VecDeriv_t<Out> dxOut(sizeOut);
+        for (unsigned i = 0; i < sizeOut; i++)
+        {
+            dxOut[i] = difference(positionOut1[i], positionOut[i]);
+        }
+
+        const auto diff = this->vectorMaxAbs(difference(dxOut, expectedVelocityOut));
+        if (diff > errorThreshold)
+        {
+            ADD_FAILURE() << "applyJ test failed: the difference between child "
+                "position change and child velocity (dt=1) " << diff <<
+                " should be less than  " << errorThreshold << std::endl
+                << "position change = " << dxOut << std::endl
+                << "expectedVelocityOut = " << expectedVelocityOut << std::endl;
+            return false;
+        }
+
+        return true;
+    }
+
+    VecDeriv_t<In> computeForceChange(core::MechanicalParams mparams, const std::size_t sizeIn, VecDeriv_t<Out> forceOut, VecDeriv_t<In> forceIn)
+    {
+        VecDeriv_t<In> forceChange;
+        forceChange.reserve(sizeIn);
+
+        // apply has been called, therefore parent force must be updated
+        // based on the same child forces
+        VecDeriv_t<In> forceIn2;
+        computeForceInFromForceOut(mparams, forceIn2, preTreatment(forceOut));
+
+        for (unsigned i = 0; i < sizeIn; ++i)
+        {
+            forceChange.push_back(forceIn2[i] - forceIn[i]);
+        }
+
+        return forceChange;
+    }
+
+    bool testGetK(const std::size_t& sizeIn,
+        const VecDeriv_t<In>& velocityIn,
+        const VecDeriv_t<In>& forceChange,
+        Real_t<In> errorThreshold)
+    {
+        VecDeriv_t<In> Kv(sizeIn);
+
+        const linearalgebra::BaseMatrix* bk = mapping->getK();
+
+        // K can be null or empty for linear mappings
+        // still performing the test with a null Kv vector to check if the mapping is really linear
+
+        if( bk != nullptr )
+        {
+            typedef linearalgebra::EigenSparseMatrix<In,In> EigenSparseKMatrix;
+            const EigenSparseKMatrix* K = dynamic_cast<const EigenSparseKMatrix*>(bk);
+            if( K == nullptr )
+            {
+                ADD_FAILURE() << "getK returns a matrix of non-EigenSparseMatrix type";
+                // TODO perform a slow conversion with a big warning rather than a failure?
+                return false;
+            }
+
+            if( K->compressedMatrix.nonZeros() )
+            {
+                K->mult(Kv,velocityIn);
+            }
+        }
+
+        // check that K.vp = dfp
+        if (this->vectorMaxDiff(Kv, forceChange) > errorThreshold * errorFactorDJ)
+        {
+            ADD_FAILURE() << "K test failed, difference should be less than " << errorThreshold*errorFactorDJ  << std::endl
+                          << "Kv    = " << Kv << std::endl
+                          << "dfp = " << forceChange << std::endl;
+            return false;
+        }
+
+        return true;
+    }
+
+    bool testBuildGeometricStiffnessMatrix(
+        std::size_t sizeIn,
+        const VecDeriv_t<In>& velocityIn,
+        const VecDeriv_t<In>& forceChange,
+        Real_t<In> errorThreshold
+        )
+    {
+        core::GeometricStiffnessMatrix testGeometricStiffness;
+
+        struct GeometricStiffnessAccumulator : core::MappingMatrixAccumulator
+        {
+            void add(sofa::SignedIndex row, sofa::SignedIndex col, float value) override
+            {
+                assembledMatrix.add(row, col, value);
+            }
+            void add(sofa::SignedIndex row, sofa::SignedIndex col, double value) override
+            {
+                assembledMatrix.add(row, col, value);
+            }
+
+            linearalgebra::EigenSparseMatrix<In,In> assembledMatrix;
+        } accumulator;
+
+        testGeometricStiffness.setMatrixAccumulator(&accumulator, mapping->getFromModel(), mapping->getFromModel());
+
+        accumulator.assembledMatrix.resize(mapping->getFromModel()->getSize() * In::deriv_total_size, mapping->getFromModel()->getSize() * In::deriv_total_size);
+        mapping->buildGeometricStiffnessMatrix(&testGeometricStiffness);
+        accumulator.assembledMatrix.compress();
+
+        VecDeriv_t<In> Kv(sizeIn);
+        if( accumulator.assembledMatrix.compressedMatrix.nonZeros() )
+        {
+            accumulator.assembledMatrix.mult(Kv,velocityIn);
+        }
+
+        // check that K.vp = dfp
+        if (this->vectorMaxDiff(Kv,forceChange) > errorThreshold*errorFactorDJ )
+        {
+            ADD_FAILURE() << "buildGeometricStiffnessMatrix test failed, difference should be less than " << errorThreshold*errorFactorDJ  << std::endl
+                          << "Kv    = " << Kv << std::endl
+                          << "dfp = " << forceChange << std::endl;
+            return false;
+        }
+
+        return true;
+    }
 
 };
 
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 e89ed41245e..fd44c92cb7c 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
@@ -145,8 +145,9 @@ void EulerImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa::
             msg_info() << "EulerImplicitSolver, f = " << f;
 
             // add the change of force due to stiffness + Rayleigh damping
-            mop.addMBKv(b, -d_rayleighMass.getValue(), 0,
-                        h + d_rayleighStiffness.getValue()); // b =  f + ( rm M + (h+rs) K ) v
+            mop.addMBKv(b, core::MatricesFactors::M(-d_rayleighMass.getValue()),
+                        core::MatricesFactors::B(0),
+                        core::MatricesFactors::K(h + d_rayleighStiffness.getValue())); // b =  f + ( rm M + (h+rs) K ) v
 
             // integration over a time step
             b.teq(h *
@@ -162,19 +163,10 @@ void EulerImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa::
 
     {
         SCOPED_TIMER("setSystemMBKMatrix");
-        SReal mFact, kFact, bFact;
-        if (firstOrder)
-        {
-            mFact = 1;
-            bFact = 0;
-            kFact = -h * tr;
-        }
-        else
-        {
-            mFact = 1 + tr * h * d_rayleighMass.getValue();
-            bFact = -tr * h;
-            kFact = -tr * h * (h + d_rayleighStiffness.getValue());
-        }
+        const core::MatricesFactors::M mFact (firstOrder ? 1 : 1 + tr * h * d_rayleighMass.getValue());
+        const core::MatricesFactors::B bFact (firstOrder ? 0 : -tr * h);
+        const core::MatricesFactors::K kFact (firstOrder ? -h * tr : -tr * h * (h + d_rayleighStiffness.getValue()));
+
         mop.setSystemMBKMatrix(mFact, bFact, kFact, l_linearSolver.get());
 
 #ifdef SOFA_DUMP_VISITOR_INFO
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 34e6f35298c..cf35ceb506e 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
@@ -113,11 +113,15 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
         mop.propagateDx(a);
 
         // b += (-h (1-\gamma)(r_M M + r_K K) - h^2/2 (1-2\beta) K ) a
-        mop.addMBKdx(b, -h*(1-gamma)*rM, h*(1-gamma), h*(1-gamma)*rK + h*h*(1-2*beta)/2.0,true,true);
+        mop.addMBKdx(b,
+            core::MatricesFactors::M(-h*(1-gamma)*rM),
+            core::MatricesFactors::B(h*(1-gamma)),
+            core::MatricesFactors::K(h*(1-gamma)*rK + h*h*(1-2*beta)/2.0),
+            true,true);
     }
 
     // b += -h K v
-    mop.addMBKv(b, -rM, 1,rK+h);
+    mop.addMBKv(b, core::MatricesFactors::M(-rM), core::MatricesFactors::B(1), core::MatricesFactors::K(rK+h));
 
     msg_info() << "b = " << b;
     mop.projectResponse(b);          // b is projected to the constrained space
@@ -126,9 +130,9 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
 
     // 3. Solve system of equations on a_{t+h}
 
-    const SReal mFact = 1 + h * gamma * rM;
-    const SReal bFact = (-h) * gamma;
-    const SReal kFact = -h * h * beta - h * rK * gamma;
+    const core::MatricesFactors::M mFact ( 1 + h * gamma * rM);
+    const core::MatricesFactors::B bFact ( (-h) * gamma);
+    const core::MatricesFactors::K kFact ( -h * h * beta - h * rK * gamma);
 
     mop.setSystemMBKMatrix(mFact, bFact, kFact, l_linearSolver.get());
     l_linearSolver->setSystemLHVector(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 48d6750452f..61f3eba93a2 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
@@ -45,29 +45,29 @@ StaticSolver::StaticSolver()
     : d_newton_iterations(initData(&d_newton_iterations,
             (unsigned) 1,
             "newton_iterations",
-            "Number of Netwon iterations between each load increments (normally, one load increment per simulation time-step."))
+            "Number of Newton iterations between each load increments (normally, one load increment per simulation time-step."))
     , d_absolute_correction_tolerance_threshold(initData(&d_absolute_correction_tolerance_threshold,
             1e-5_sreal,
             "absolute_correction_tolerance_threshold",
-            "Convergence criterion of the norm |du| under which the Netwon iterations stop"))
+            "Convergence criterion of the norm |du| under which the Newton iterations stop"))
     , d_relative_correction_tolerance_threshold(initData(&d_relative_correction_tolerance_threshold,
             1e-5_sreal,
             "relative_correction_tolerance_threshold",
-            "Convergence criterion regarding the ratio |du| / |U| under which the Netwon iterations stop"))
+            "Convergence criterion regarding the ratio |du| / |U| under which the Newton iterations stop"))
     , d_absolute_residual_tolerance_threshold( initData(&d_absolute_residual_tolerance_threshold,
             1e-5_sreal,
             "absolute_residual_tolerance_threshold",
-            "Convergence criterion of the norm |R| under which the Netwon iterations stop."
+            "Convergence criterion of the norm |R| under which the Newton iterations stop."
             "Use a negative value to disable this criterion"))
     , d_relative_residual_tolerance_threshold( initData(&d_relative_residual_tolerance_threshold,
             1e-5_sreal,
             "relative_residual_tolerance_threshold",
-            "Convergence criterion regarding the ratio |R|/|R0| under which the Netwon iterations stop."
+            "Convergence criterion regarding the ratio |R|/|R0| under which the Newton iterations stop."
             "Use a negative value to disable this criterion"))
     , d_should_diverge_when_residual_is_growing( initData(&d_should_diverge_when_residual_is_growing,
             false,
             "should_diverge_when_residual_is_growing",
-            "Boolean stopping Netwon iterations when the residual is greater than the one from the previous iteration"))
+            "Boolean stopping Newton iterations when the residual is greater than the one from the previous iteration"))
 {}
 
 void StaticSolver::solve(const sofa::core::ExecParams* params, SReal dt, sofa::core::MultiVecCoordId xResult, sofa::core::MultiVecDerivId vResult)
@@ -212,7 +212,10 @@ void StaticSolver::solve(const sofa::core::ExecParams* params, SReal dt, sofa::c
             //       is called on every BaseProjectiveConstraintSet objects. An example of such constraint set is the
             //       FixedProjectiveConstraint. In this case, it will set to 0 every column (_, i) and row (i, _) of the assembled
             //       matrix for the ith degree of freedom.
-            mop.setSystemMBKMatrix(0, 0, -1, l_linearSolver.get());
+            mop.setSystemMBKMatrix(
+                core::MatricesFactors::M(0),
+                core::MatricesFactors::B(0),
+                core::MatricesFactors::K(-1), l_linearSolver.get());
         }
 
         // Part II. Solve the unknown increment.
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 173b2523a89..782466879b9 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
@@ -215,7 +215,9 @@ void VariationalSymplecticSolver::solve(const core::ExecParams* params, SReal dt
 
 			// corresponds to do b+=-K*res, where res=res(i-1)=q(k,i-1)-q(k,0)
 			mop.propagateDx(res);
-			mop.addMBKdx(b,0,0,-1.0);
+			mop.addMBKdx(b,core::MatricesFactors::M(0),
+			    core::MatricesFactors::B(0),
+			    core::MatricesFactors::K(-1.0));
 
 
 			mop.projectResponse(b);
@@ -224,9 +226,9 @@ void VariationalSymplecticSolver::solve(const core::ExecParams* params, SReal dt
 			// add left term : matrix=-K+4/h^(2)M, but with dampings rK and rM
 		    {
 			    SCOPED_TIMER("MBKBuild");
-                const SReal mFact = 4.0 / (h * h) + 4 * rM / h;
-			    const SReal kFact = -1.0 - 4 * rK / h;
-                mop.setSystemMBKMatrix(mFact, 0, kFact, l_linearSolver.get());
+                const core::MatricesFactors::M mFact ( 4.0 / (h * h) + 4 * rM / h );
+			    const core::MatricesFactors::K kFact ( -1.0 - 4 * rK / h );
+                mop.setSystemMBKMatrix(mFact, core::MatricesFactors::B(0), kFact, l_linearSolver.get());
 		    }
 
             {
@@ -312,7 +314,7 @@ void VariationalSymplecticSolver::solve(const core::ExecParams* params, SReal dt
             MultiVecDeriv b(&vop);
             b.clear();
             // Mass matrix
-            mop.setSystemMBKMatrix(1, 0, 0, l_linearSolver.get());
+            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);
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 a1a0ecc83b2..b9ccad03a69 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
@@ -328,7 +328,10 @@ void EulerExplicitSolver::assembleSystemMatrix(sofa::simulation::common::Mechani
     //       is called on every BaseProjectiveConstraintSet objects. An example of such constraint set is the
     //       FixedProjectiveConstraint. In this case, it will set to 0 every column (_, i) and row (i, _) of the assembled
     //       matrix for the ith degree of freedom.
-    mop->setSystemMBKMatrix(1, 0, 0, l_linearSolver.get());
+    mop->setSystemMBKMatrix(
+        core::MatricesFactors::M(1),
+        core::MatricesFactors::B(0),
+        core::MatricesFactors::K(0), l_linearSolver.get());
 }
 
 void EulerExplicitSolver::solveSystem(core::MultiVecDerivId solution, core::MultiVecDerivId rhs) const
diff --git a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.cpp b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.cpp
index 2249b18d5e1..d126f8ac906 100644
--- a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.cpp
+++ b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.cpp
@@ -56,8 +56,8 @@ DataDisplay::DataDisplay()
     , state(nullptr)
     , m_topology(nullptr)
     , l_topology(initLink("topology", "link to the topology container"))
-    , oldMin(0)
-    , oldMax(0)
+    , m_oldMin(std::numeric_limits<Real>::max())
+    , m_oldMax(std::numeric_limits<Real>::lowest())
 {
     this->addAlias(&f_triangleData,"cellData"); // backward compatibility
     d_currentMin.setReadOnly(true);
@@ -141,9 +141,8 @@ void DataDisplay::doDrawVisual(const core::visual::VisualParams* vparams)
     }
 
     // Range for points
-    Real min ;
-    Real max ;
-    min = max = 0;
+    Real min = std::numeric_limits<Real>::max();
+    Real max = std::numeric_limits<Real>::lowest();
     if (bDrawPointData) {
         VecPointData::const_iterator i = ptData.begin();
         min = *i;
@@ -212,12 +211,12 @@ void DataDisplay::doDrawVisual(const core::visual::VisualParams* vparams)
     }
 
 
-    if (max > oldMax) oldMax = max;
-    if (min < oldMin) oldMin = min;
+    if (max > m_oldMax) m_oldMax = max;
+    if (min < m_oldMin) m_oldMin = min;
 
     if (f_maximalRange.getValue()) {
-        max = oldMax;
-        min = oldMin;
+        max = m_oldMax;
+        min = m_oldMin;
     }
     d_currentMin.setValue(min);
     d_currentMax.setValue(max);
diff --git a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.h b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.h
index 91491ef5d43..91bb25b1acb 100644
--- a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.h
+++ b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/DataDisplay.h
@@ -66,8 +66,6 @@ class SOFA_GL_COMPONENT_RENDERING3D_API DataDisplay : public core::visual::Visua
     /// Link to be set to the topology container in the component graph.
     SingleLink <DataDisplay, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology;
 
-    Real oldMin, oldMax;
-
     void init() override;
     void doDrawVisual(const core::visual::VisualParams* vparams) override;
     void doUpdateVisual(const core::visual::VisualParams* vparams) override;
@@ -80,6 +78,9 @@ class SOFA_GL_COMPONENT_RENDERING3D_API DataDisplay : public core::visual::Visua
     type::vector<type::Vec3f> m_normals;
 
     DataDisplay();
+
+    Real m_oldMin, m_oldMax;
+
 };
 
 } // namespace sofa::gl::component::rendering3d
diff --git a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.cpp b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.cpp
index 41d46622f1a..f2bba6ceb04 100644
--- a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.cpp
+++ b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.cpp
@@ -32,20 +32,15 @@ int OglSceneFrameClass = core::RegisterObject("Display a frame at the corner of
 
 using namespace sofa::defaulttype;
 
+static constexpr OglSceneFrame::Alignment defaultAlignment("BottomRight");
+static constexpr OglSceneFrame::Style defaultStyle("Cylinders");
+
 OglSceneFrame::OglSceneFrame()
     : d_drawFrame(initData(&d_drawFrame, true,  "draw", "Display the frame or not"))
-    , d_style(initData(&d_style, "style", "Style of the frame"))
-    , d_alignment(initData(&d_alignment, "alignment", "Alignment of the frame in the view"))
+    , d_style(initData(&d_style, defaultStyle, "style", ("Style of the frame\n" + Style::dataDescription()).c_str()))
+    , d_alignment(initData(&d_alignment, defaultAlignment, "alignment", ("Alignment of the frame in the view\n" + Alignment::dataDescription()).c_str()))
     , d_viewportSize(initData(&d_viewportSize, 150, "viewportSize", "Size of the viewport where the frame is rendered"))
-{
-    sofa::helper::OptionsGroup styleOptions{"Arrows", "Cylinders", "CubeCones"};
-    styleOptions.setSelectedItem(1);
-    d_style.setValue(styleOptions);
-
-    sofa::helper::OptionsGroup alignmentOptions{"BottomLeft", "BottomRight", "TopRight", "TopLeft"};
-    alignmentOptions.setSelectedItem(1);
-    d_alignment.setValue(alignmentOptions);
-}
+{}
 
 void OglSceneFrame::drawArrows(const core::visual::VisualParams* vparams)
 {
@@ -112,22 +107,22 @@ void OglSceneFrame::doDrawVisual(const core::visual::VisualParams* vparams)
 
     const auto viewportSize = d_viewportSize.getValue();
 
-    switch(d_alignment.getValue().getSelectedId())
+    switch(d_alignment.getValue())
     {
-        case 0: //BottomLeft
+        case Alignment("BottomLeft"):
         default:
             glViewport(0,0,viewportSize,viewportSize);
             glScissor(0,0,viewportSize,viewportSize);
             break;
-        case 1: //BottomRight
+        case Alignment("BottomRight"):
             glViewport(viewport[2]-viewportSize,0,viewportSize,viewportSize);
             glScissor(viewport[2]-viewportSize,0,viewportSize,viewportSize);
             break;
-        case 2: //TopRight
+        case Alignment("TopRight"):
             glViewport(viewport[2]-viewportSize,viewport[3]-viewportSize,viewportSize,viewportSize);
             glScissor(viewport[2]-viewportSize,viewport[3]-viewportSize,viewportSize,viewportSize);
             break;
-        case 3: //TopLeft
+        case Alignment("TopLeft"):
             glViewport(0,viewport[3]-viewportSize,viewportSize,viewportSize);
             glScissor(0,viewport[3]-viewportSize,viewportSize,viewportSize);
             break;
@@ -157,18 +152,18 @@ void OglSceneFrame::doDrawVisual(const core::visual::VisualParams* vparams)
 
     vparams->drawTool()->disableLighting();
 
-    switch (d_style.getValue().getSelectedId())
+    switch (d_style.getValue())
     {
-    case 0:
+    case Style("Arrows"):
     default:
         drawArrows(vparams);
         break;
 
-    case 1:
+    case Style("Cylinders"):
         drawCylinders(vparams);
         break;
 
-    case 2:
+    case Style("CubeCones"):
         drawCubeCones(vparams);
         break;
     }
diff --git a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.h b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.h
index c6cdd2288bb..d834c9584d1 100644
--- a/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.h
+++ b/Sofa/GL/Component/Rendering3D/src/sofa/gl/component/rendering3d/OglSceneFrame.h
@@ -25,6 +25,8 @@
 #include <sofa/core/visual/VisualModel.h>
 #include <sofa/core/visual/VisualParams.h>
 #include <sofa/helper/OptionsGroup.h>
+#include <sofa/helper/SelectableItem.h>
+
 
 namespace sofa::gl::component::rendering3d
 {
@@ -38,8 +40,23 @@ class SOFA_GL_COMPONENT_RENDERING3D_API OglSceneFrame : public core::visual::Vis
     typedef core::visual::VisualParams::Viewport Viewport;
 
     Data<bool> d_drawFrame; ///< Display the frame or not
-    Data<sofa::helper::OptionsGroup> d_style; ///< Style of the frame
-    Data<sofa::helper::OptionsGroup> d_alignment; ///< Alignment of the frame in the view
+
+    MAKE_SELECTABLE_ITEMS(Style,
+        sofa::helper::Item{"Arrows", "The frame is composed of arrows"},
+        sofa::helper::Item{"Cylinders", "The frame is composed of cylinders"},
+        sofa::helper::Item{"CubeCones", "The frame is composed of cubes and cones"},
+    );
+
+    Data<Style> d_style; ///< Style of the frame
+
+    MAKE_SELECTABLE_ITEMS(Alignment,
+        sofa::helper::Item{"BottomLeft", "The scene frame is displayed in the bottom-left corner"},
+        sofa::helper::Item{"BottomRight", "The scene frame is displayed in the bottom-right corner"},
+        sofa::helper::Item{"TopRight", "The scene frame is displayed in the top-right corner"},
+        sofa::helper::Item{"TopLeft", "The scene frame is displayed in the top-left corner"}
+    );
+
+    Data<Alignment> d_alignment; ///< Alignment of the frame in the view
     Data<int> d_viewportSize; ///< Size of the viewport where the frame is rendered
 
     OglSceneFrame();
diff --git a/Sofa/GUI/Qt/CMakeLists.txt b/Sofa/GUI/Qt/CMakeLists.txt
index 2989ce4279f..85ebe4e49ee 100644
--- a/Sofa/GUI/Qt/CMakeLists.txt
+++ b/Sofa/GUI/Qt/CMakeLists.txt
@@ -24,8 +24,11 @@ endif()
 
 if (Qt6Core_FOUND)
     message("${PROJECT_NAME}: will use Qt6")
-    sofa_find_package(Qt6 COMPONENTS Gui GuiTools Widgets WidgetsTools OpenGLWidgets REQUIRED)
-    set(SOFA_GUI_QT_TARETS ${SOFA_GUI_QT_TARGETS} Qt::Core Qt::Gui Qt::Widgets Qt::OpenGLWidgets )
+    find_package(Qt6 COMPONENTS Gui GuiTools Widgets WidgetsTools OpenGLWidgets REQUIRED)
+    # GuiTools & WidgetsTools does not define a Qt6::component CMake target
+    # So we can't look for those target to know how we should define SOFA_GUI_QT_HAVE_QT6
+    sofa_find_package(Qt6 COMPONENTS Gui Widgets OpenGLWidgets REQUIRED)
+    set(SOFA_GUI_QT_TARGETS ${SOFA_GUI_QT_TARGETS} Qt::Core Qt::Gui Qt::Widgets Qt::OpenGLWidgets )
 elseif (Qt5Core_FOUND)
     message("${PROJECT_NAME}: will use Qt5 (deprecated)")
     sofa_find_package(Qt5 COMPONENTS Core Gui OpenGL REQUIRED)
diff --git a/Sofa/GUI/Qt/src/sofa/gui/qt/DataWidget.h b/Sofa/GUI/Qt/src/sofa/gui/qt/DataWidget.h
index 52b9c8f5313..7b068fedf58 100644
--- a/Sofa/GUI/Qt/src/sofa/gui/qt/DataWidget.h
+++ b/Sofa/GUI/Qt/src/sofa/gui/qt/DataWidget.h
@@ -177,11 +177,37 @@ class TDataWidget : public DataWidget
     template <class RealObject>
     static RealObject* create( RealObject*, CreatorArgument& arg)
     {
+        RealObject* obj = nullptr;
+
         typename RealObject::MyTData* realData = dynamic_cast< typename RealObject::MyTData* >(arg.data);
-        if (!realData) return nullptr;
+        if (!realData)
+        {
+            if constexpr (std::is_constructible_v<RealObject, QWidget*, const char*, core::objectmodel::BaseData*, const T*>)
+            {
+                if (arg.data)
+                {
+                    const void* rawPtr = arg.data->getValueVoidPtr();
+
+                    //note that this cast is not type-safe. You need to check
+                    //later in createWidget that the baseData is the expected
+                    //type. You can use getValueTypeString for example
+                    if (const T* castedPtr = static_cast<const T*>(rawPtr))
+                    {
+                        obj = new RealObject(arg.parent, arg.name.c_str(), arg.data, castedPtr);
+                    }
+                }
+            }
+        }
         else
         {
-            RealObject* obj = new RealObject(arg.parent,arg.name.c_str(), realData);
+            if constexpr (std::is_constructible_v<RealObject, QWidget*, const char*, MyTData*>)
+            {
+                obj = new RealObject(arg.parent, arg.name.c_str(), realData);
+            }
+        }
+
+        if (obj)
+        {
             if( !obj->createWidgets() )
             {
                 delete obj;
@@ -191,13 +217,20 @@ class TDataWidget : public DataWidget
             {
                 obj->setDataReadOnly(arg.readOnly);
             }
-            return obj;
         }
+        return obj;
+    }
 
+    TDataWidget(QWidget* parent,const char* name, core::objectmodel::BaseData* d, const T* object):
+        DataWidget(parent, name, d),
+        Tdata(nullptr)
+    {
+        SOFA_UNUSED(object);
     }
 
     TDataWidget(QWidget* parent,const char* name, MyTData* d):
         DataWidget(parent,name,d),Tdata(d) {}
+
     /// Accessor function. Gives you the actual data instead
     /// of a BaseData pointer of it like in getBaseData().
     sofa::core::objectmodel::Data<T>* getData() {return Tdata;}
@@ -209,7 +242,7 @@ class TDataWidget : public DataWidget
         Tdata = d;
     }
 protected:
-    MyTData* Tdata;
+    MyTData* Tdata { nullptr };
 };
 
 
diff --git a/Sofa/GUI/Qt/src/sofa/gui/qt/DisplayFlagsDataWidget.cpp b/Sofa/GUI/Qt/src/sofa/gui/qt/DisplayFlagsDataWidget.cpp
index b2b609088c2..e09f283ef12 100644
--- a/Sofa/GUI/Qt/src/sofa/gui/qt/DisplayFlagsDataWidget.cpp
+++ b/Sofa/GUI/Qt/src/sofa/gui/qt/DisplayFlagsDataWidget.cpp
@@ -98,7 +98,8 @@ DisplayFlagWidget::DisplayFlagWidget(QWidget* parent, const char* name,  Qt::Win
     itemShowFlag[NORMALS]   = new QTreeWidgetItem(itemShowOptions, itemShowFlag[WIREFRAME]);
     this->setTreeWidgetCheckable(itemShowFlag[NORMALS], "Normals");
     this->addTopLevelItem(itemShowAll);
-    itemShowAll->addChild(itemShowVisual); itemShowAll->setExpanded(true);
+    itemShowAll->setExpanded(true);
+    itemShowAll->addChild(itemShowVisual);
     itemShowVisual->addChild(itemShowFlag[VISUALMODELS]);
     itemShowAll->addChild(itemShowBehavior);
     itemShowBehavior->addChild(itemShowFlag[BEHAVIORMODELS]);
@@ -112,7 +113,8 @@ DisplayFlagWidget::DisplayFlagWidget(QWidget* parent, const char* name,  Qt::Win
     itemShowMapping->addChild(itemShowFlag[MAPPINGS]);
     itemShowMapping->addChild(itemShowFlag[MECHANICALMAPPINGS]);
 
-    this->addTopLevelItem(itemShowOptions); itemShowOptions->setExpanded(true);
+    this->addTopLevelItem(itemShowOptions);
+    itemShowOptions->setExpanded(true);
     itemShowOptions->addChild(itemShowFlag[RENDERING]);
     itemShowOptions->addChild(itemShowFlag[WIREFRAME]);
     itemShowOptions->addChild(itemShowFlag[NORMALS]);
diff --git a/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.cpp b/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.cpp
index 0f051db8bf3..7942456a333 100644
--- a/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.cpp
+++ b/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.cpp
@@ -130,6 +130,15 @@ Creator<DataWidgetFactory, SimpleDataWidget< Mat<6,6,double> > > DWClass_Mat66d(
 Creator<DataWidgetFactory, SimpleDataWidget< sofa::core::objectmodel::TagSet > > DWClass_TagSet("default",true);
 
 
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<1,float>> > > DWClass_CRSCVec1f("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<2,float>> > > DWClass_CRSCVec2f("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<3,float>> > > DWClass_CRSCVec3f("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<6,float>> > > DWClass_CRSCVec6f("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<1,double>> > > DWClass_CRSCVec1d("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<2,double>> > > DWClass_CRSCVec2d("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<3,double>> > > DWClass_CRSCVec3d("default",true);
+Creator<DataWidgetFactory, SimpleDataWidget< sofa::linearalgebra::CompressedRowSparseMatrixConstraint<Vec<6,double>> > > DWClass_CRSCVec6d("default",true);
+
 ////////////////////////////////////////////////////////////////
 /// OptionsGroup support
 ////////////////////////////////////////////////////////////////
@@ -222,6 +231,122 @@ void RadioDataWidget::writeToData()
     this->getData()->setValue(m_radiotrick);
 }
 
+Creator<DataWidgetFactory, SelectableItemWidget> DWClass_SelectableItem("default",true);
+
+SelectableItemWidget::SelectableItemWidget(QWidget* parent, const char* name,
+    core::BaseData* m_data, const helper::BaseSelectableItem* item)
+: TDataWidget(parent, name, m_data, item)
+, m_selectableItem(item)
+{}
+
+bool SelectableItemWidget::createWidgets()
+{
+    if (Tdata && Tdata->getValueTypeString() != "SelectableItem" ||
+        baseData && baseData->getValueTypeString() != "SelectableItem")
+    {
+        return false;
+    }
+
+    QVBoxLayout* layout = new QVBoxLayout(this);
+    static constexpr unsigned int LIMIT_NUM_BUTTON = 4;
+
+    assert(m_selectableItem);
+    const std::size_t nbItems = m_selectableItem->getNumberOfItems();
+    m_buttonMode = nbItems < LIMIT_NUM_BUTTON;
+    const auto* items = m_selectableItem->getItemsData();
+
+    const auto getItem = [](const sofa::helper::Item* item)
+    {
+        std::stringstream ss;
+        ss << item->key;
+        if (!item->description.empty())
+        {
+            ss << " (" << item->description << ")";
+        }
+        return ss.str();
+    };
+
+    if (m_buttonMode)
+    {
+        m_buttonList = new QButtonGroup(this);
+
+        for (std::size_t i = 0; i < nbItems; i++)
+        {
+            const helper::Item* item_i = items + i;
+
+            QRadioButton * m_radiobutton = new QRadioButton(QString(getItem(item_i).c_str()), this);
+            if (i == m_selectableItem->getSelectedId())
+            {
+                m_radiobutton->setChecked(true);
+            }
+            layout->addWidget(m_radiobutton);
+            m_buttonList->addButton(m_radiobutton, i);
+        }
+        connect(m_buttonList, SIGNAL(buttonClicked(int)), this, SLOT(setWidgetDirty())) ;
+    }
+    else
+    {
+        m_comboList=new QComboBox(this);
+        m_comboList->setSizePolicy(QSizePolicy::Expanding, QSizePolicy::Minimum);
+
+        QStringList list;
+        for (std::size_t i = 0; i < nbItems; i++)
+        {
+            const helper::Item* item_i = items + i;
+            list << getItem(item_i).c_str();
+        }
+
+        m_comboList->insertItems(0, list);
+
+        m_comboList->setCurrentIndex(m_selectableItem->getSelectedId());
+
+        connect(m_comboList, SIGNAL(activated(int)), this, SLOT(setWidgetDirty()));
+        layout->addWidget(m_comboList);
+
+    }
+
+    return true;
+}
+
+void SelectableItemWidget::setDataReadOnly(const bool readOnly)
+{
+    if (m_buttonMode)
+    {
+        const QList<QAbstractButton *> buttons = m_buttonList->buttons();
+        for (auto& button : buttons)
+        {
+            button->setEnabled(!readOnly);
+        }
+    }
+    else
+    {
+        m_comboList->setEnabled(!readOnly);
+    }
+}
+
+void SelectableItemWidget::readFromData()
+{
+    if (m_buttonMode)
+    {
+        m_buttonList->button(m_selectableItem->getSelectedId())->setChecked(true);
+    }
+    else
+    {
+        m_comboList->setCurrentIndex(m_selectableItem->getSelectedId());
+    }
+}
+
+void SelectableItemWidget::writeToData()
+{
+    if (m_buttonMode)
+    {
+        const_cast<helper::BaseSelectableItem*>(m_selectableItem)->setSelectedId(static_cast<std::size_t>(m_buttonList->checkedId()));
+    }
+    else
+    {
+        const_cast<helper::BaseSelectableItem*>(m_selectableItem)->setSelectedId(static_cast<std::size_t>(m_comboList->currentIndex()));
+    }
+}
 
 
 } //namespace sofa::gui::qt
diff --git a/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.h b/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.h
index 8da330398a4..3583f2c29a8 100644
--- a/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.h
+++ b/Sofa/GUI/Qt/src/sofa/gui/qt/SimpleDataWidget.h
@@ -27,6 +27,7 @@
 //#include <sofa/defaulttype/RigidTypes.h>
 #include <sofa/type/fixed_array.h>
 #include <sofa/core/topology/Topology.h>
+#include <sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h>
 #include "WDoubleLineEdit.h"
 #include <climits>
 
@@ -932,6 +933,45 @@ template<sofa::Size L, sofa::Size C, class T>
 class data_widget_container < sofa::type::Mat<L, C, T> > : public fixed_grid_data_widget_container < sofa::type::Mat<L, C, T> >
 {};
 
+////////////////////////////////////////////////////////////////
+/// sofa::linearalgebra::CompressedRowSparseMatrixConstraint support
+////////////////////////////////////////////////////////////////
+
+template<typename TBlock>
+class data_widget_trait <sofa::linearalgebra::CompressedRowSparseMatrixConstraint<TBlock>>
+{
+public:
+    typedef sofa::linearalgebra::CompressedRowSparseMatrixConstraint<TBlock> data_type;
+    typedef QLineEdit Widget;
+    static Widget* create(QWidget* parent, const data_type& /*d*/)
+    {
+        Widget* w = new Widget(parent);
+        w->setSizePolicy(QSizePolicy::Expanding, QSizePolicy::Minimum);
+        return w;
+    }
+    static void readFromData(Widget* w, const data_type& d)
+    {
+        std::ostringstream oss;
+        d.prettyPrint(oss);
+        w->setText(QString(oss.str().c_str()));
+    }
+    static void writeToData(Widget* /* w */, const data_type& /* d */)
+    {
+        // not supported by this type
+        // TODO: CompressedRowSparseMatrixConstraint needs a parser for its pretty output
+    }
+    static void setReadOnly(Widget* w, bool readOnly)
+    {
+        w->setEnabled(!readOnly);
+        w->setReadOnly(readOnly);
+    }
+    static void connectChanged(Widget* w, DataWidget* datawidget)
+    {
+        datawidget->connect(w, SIGNAL( textChanged(const QString&) ), datawidget, SLOT(setWidgetDirty()) );
+    }
+    
+};
+
 ////////////////////////////////////////////////////////////////
 /// OptionsGroup support
 ////////////////////////////////////////////////////////////////
@@ -946,7 +986,7 @@ public :
     ///a widget for a that particular data type.
     RadioDataWidget(QWidget* parent, const char* name,
             core::objectmodel::Data<sofa::helper::OptionsGroup >* m_data)
-        : TDataWidget<sofa::helper::OptionsGroup >(parent,name,m_data) {};
+        : TDataWidget<sofa::helper::OptionsGroup >(parent,name,m_data) {}
 
     ///In this method we  create the widgets and perform the signal / slots connections.
     virtual bool createWidgets();
@@ -964,5 +1004,28 @@ public :
     bool buttonMode;
 };
 
+class SelectableItemWidget final : public TDataWidget<helper::BaseSelectableItem>
+{
+    Q_OBJECT
+public :
+
+    SelectableItemWidget(QWidget* parent, const char* name,
+            core::BaseData* m_data, const helper::BaseSelectableItem* item);
+
+    bool createWidgets() override;
+    void setDataReadOnly(bool readOnly) override;
+
+protected:
+    void readFromData() override;
+
+    void writeToData() override;
+
+    QButtonGroup *m_buttonList { nullptr };
+    QComboBox    *m_comboList { nullptr };
+    bool m_buttonMode { false };
+
+    const helper::BaseSelectableItem* m_selectableItem { nullptr };
+};
+
 
 } //namespace sofa::gui::qt
diff --git a/Sofa/framework/Config/src/sofa/config.h.in b/Sofa/framework/Config/src/sofa/config.h.in
index df763c20356..ccd35844237 100644
--- a/Sofa/framework/Config/src/sofa/config.h.in
+++ b/Sofa/framework/Config/src/sofa/config.h.in
@@ -69,12 +69,6 @@ constexpr SReal operator"" _sreal(unsigned long long int integer)
 #  endif
 #endif
 
-// snprintf() has been provided since MSVC++ 14 (Visual Studio 2015).  For other
-// versions, it is simply #defined to _snprintf().
-#if (defined(_MSC_VER) && _MSC_VER < 1900)
-#  define snprintf _snprintf
-#endif
-
 #ifdef _WIN32
 #define WIN32_LEAN_AND_MEAN
 #  include <windows.h>
@@ -223,8 +217,8 @@ class DeprecatedAndRemoved {};
 
 /**********************************************/
 
-#define SOFA_DECL_CLASS(name) // extern "C" { int sofa_concat(class_,name) = 0; }
-#define SOFA_LINK_CLASS(name) // extern "C" { extern int sofa_concat(class_,name); int sofa_concat(link_,name) = sofa_concat(class_,name); }
+#define SOFA_DECL_CLASS(name) SOFA_PRAGMA_WARNING("The macro SOFA_DECL_CLASS is no longer used")
+#define SOFA_LINK_CLASS(name) SOFA_PRAGMA_WARNING("The macro SOFA_LINK_CLASS is no longer used")
 
 // Prevent compiler warnings about 'unused variables'.
 // This should be used when a parameter name is needed (e.g. for
@@ -238,12 +232,6 @@ class DeprecatedAndRemoved {};
     #define SOFA_CLASS_METHOD ( std::string(this->getClassName()) + "::" + __func__ + " " )
 #endif
 
-
-// The SOFA_EXTERN_TEMPLATE macro was used to control the use of extern templates in Sofa.
-// It has been cleaned out in 41e0ab98bbb6e22b053b24e7bbdd31c8636336c9 "[ALL] Remove SOFA_EXTERN_TEMPLATE".
-// Macro definition is kept to avoid breaking all external plugins.
-#define SOFA_EXTERN_TEMPLATE SOFA_PRAGMA_ERROR("The usage of SOFA_EXTERN_TEMPLATE is now disabled. See commit 41e0ab9 for more information")
-
 namespace sofa
 {
 
diff --git a/Sofa/framework/Core/CMakeLists.txt b/Sofa/framework/Core/CMakeLists.txt
index e490b7af3bc..4910d364c28 100644
--- a/Sofa/framework/Core/CMakeLists.txt
+++ b/Sofa/framework/Core/CMakeLists.txt
@@ -31,6 +31,7 @@ set(HEADER_FILES
     ${SRC_ROOT}/Mapping.h
     ${SRC_ROOT}/Mapping.inl
     ${SRC_ROOT}/MappingHelper.h
+    ${SRC_ROOT}/MatricesFactors.h
     ${SRC_ROOT}/MatrixAccumulator.h
     ${SRC_ROOT}/MechanicalParams.h
     ${SRC_ROOT}/MechanicalStatesMatrixAccumulators.h
diff --git a/Sofa/framework/Core/src/sofa/core/MatricesFactors.h b/Sofa/framework/Core/src/sofa/core/MatricesFactors.h
new file mode 100644
index 00000000000..50a6a6561d7
--- /dev/null
+++ b/Sofa/framework/Core/src/sofa/core/MatricesFactors.h
@@ -0,0 +1,44 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/type/StrongType.h>
+#include <sofa/config.h>
+
+namespace sofa::core
+{
+
+/**
+ * Contains a strong type for each of the 3 main matrices
+ */
+struct MatricesFactors
+{
+    // A strong type for the mass matrix factor
+    using M = sofa::type::StrongType<SReal, struct MFactorTag, sofa::type::functionality::Arithmetic>;
+
+    // A strong type for the damping matrix factor
+    using B = sofa::type::StrongType<SReal, struct BFactorTag, sofa::type::functionality::Arithmetic>;
+
+    // A strong type for the stiffness matrix factor
+    using K = sofa::type::StrongType<SReal, struct KFactorTag, sofa::type::functionality::Arithmetic>;
+};
+
+}
diff --git a/Sofa/framework/Core/src/sofa/core/objectmodel/BaseClassNameHelper.h b/Sofa/framework/Core/src/sofa/core/objectmodel/BaseClassNameHelper.h
index 2ae1e5d194b..495389d87ae 100644
--- a/Sofa/framework/Core/src/sofa/core/objectmodel/BaseClassNameHelper.h
+++ b/Sofa/framework/Core/src/sofa/core/objectmodel/BaseClassNameHelper.h
@@ -27,7 +27,7 @@
 
 namespace sofa::defaulttype
 {
-    template<class T> struct DataTypeInfo;
+    template<class T, typename Enable> struct DataTypeInfo;
 }
 
 namespace sofa::core::objectmodel
@@ -78,7 +78,7 @@ class HasDataTypeInfo
     typedef char YesType[1];
     typedef char NoType[2];
 
-    template<typename C> static YesType& test( decltype (&sofa::defaulttype::DataTypeInfo<C>::name) );
+    template<typename C> static YesType& test( decltype (&sofa::defaulttype::DataTypeInfo<C, void>::name) );
     template<typename C> static NoType& test(...);
 
 public:
@@ -264,7 +264,7 @@ std::string GetSofaTypeTemplateName(const std::string prefix)
     if constexpr (HasName<T>::value )
             return prefix + T::Name();
     else if constexpr (HasDataTypeInfo<T>::value )
-            return prefix + sofa::defaulttype::DataTypeInfo<T>::name();
+            return prefix + sofa::defaulttype::DataTypeInfo<T, void>::name();
     else
         return prefix + sofa::helper::NameDecoder::decodeTypeName(typeid(T));
 }
diff --git a/Sofa/framework/Core/src/sofa/core/objectmodel/BaseData.cpp b/Sofa/framework/Core/src/sofa/core/objectmodel/BaseData.cpp
index 3d16511118f..d449b0b7a7e 100644
--- a/Sofa/framework/Core/src/sofa/core/objectmodel/BaseData.cpp
+++ b/Sofa/framework/Core/src/sofa/core/objectmodel/BaseData.cpp
@@ -128,9 +128,14 @@ bool BaseData::setParent(BaseData* parent, const std::string& path)
         addInput(parent);
         BaseData::setDirtyValue();
         m_counter++;
-        m_isSet = true;
-    }else if (!path.empty())
+    }
+    // the referenced parent data has not been created yet but a path has been given
+    else if (!path.empty())
+    {
         parentData.setPath(path);
+    }
+
+    m_isSet = true;
 
     return true;
 }
diff --git a/Sofa/framework/Core/test/objectmodel/DataLink_test.cpp b/Sofa/framework/Core/test/objectmodel/DataLink_test.cpp
index 5224554b18f..34cd074c650 100644
--- a/Sofa/framework/Core/test/objectmodel/DataLink_test.cpp
+++ b/Sofa/framework/Core/test/objectmodel/DataLink_test.cpp
@@ -40,7 +40,6 @@ struct DataLink_test: public BaseTest
     Data<Vec3f> dataVec3f;
     Data<Vec3d> dataVec3d;
 
-    /// This method is defined in gtest framework to setting the test up.
     void SetUp() override
     {
         /// Setup the data and create a link between the two data
@@ -58,7 +57,7 @@ struct DataLink_test: public BaseTest
     }
 };
 
-/// This test check that the setting/unsetting mechanisme when the value is changed is working
+/// This test check that the setting/unsetting mechanism when the value is changed is working
 /// Currently in Sofa the parenting link is not broken if the value is written.
 TEST_F(DataLink_test, UnsetByValue)
 {
@@ -79,7 +78,19 @@ TEST_F(DataLink_test, SetParentOfDifferentType)
     ASSERT_FLOAT_EQ(dataVec3f.getValue().z(), 3.0f);
 }
 
-/// This test check that the setting/unsetting mechanisme is working
+TEST_F(DataLink_test, SetParentOfDifferentType2)
+{
+    Data<sofa::type::vector<sofa::type::Vec1d>> dataVecVec1d;
+    Data<sofa::type::vector<SReal>> dataVecSReal;
+
+    ASSERT_TRUE(dataVecSReal.setParent(&dataVecVec1d));
+    ASSERT_TRUE(dataVecVec1d.getParent() == nullptr);
+    dataVecVec1d.setValue({sofa::type::Vec1d(1.0), sofa::type::Vec1d(2.0)});
+    ASSERT_DOUBLE_EQ(dataVecSReal.getValue()[0], 1.);
+    ASSERT_DOUBLE_EQ(dataVecSReal.getValue()[1], 2.);
+}
+
+/// This test check that the setting/unsetting mechanism is working
 TEST_F(DataLink_test, Set)
 {
     ASSERT_EQ(data1.getParent(),nullptr);
@@ -89,7 +100,7 @@ TEST_F(DataLink_test, Set)
     ASSERT_EQ(data2.getParent(), &data1);
 }
 
-/// This test check that the setting/unsetting mechanisme is working
+/// This test check that the setting/unsetting mechanism is working
 TEST_F(DataLink_test, Unset)
 {
     ASSERT_EQ(data1.getParent(),nullptr);
diff --git a/Sofa/framework/DefaultType/CMakeLists.txt b/Sofa/framework/DefaultType/CMakeLists.txt
index b35dc3e60ee..c0142271c3e 100644
--- a/Sofa/framework/DefaultType/CMakeLists.txt
+++ b/Sofa/framework/DefaultType/CMakeLists.txt
@@ -39,6 +39,7 @@ set(HEADER_FILES
     ${SRC_ROOT}/typeinfo/TypeInfo_RGBAColor.h
     ${SRC_ROOT}/typeinfo/TypeInfo_Set.h
     ${SRC_ROOT}/typeinfo/TypeInfo_Scalar.h
+    ${SRC_ROOT}/typeinfo/TypeInfo_SelectableItem.h
     ${SRC_ROOT}/typeinfo/TypeInfo_Text.h
     ${SRC_ROOT}/typeinfo/TypeInfo_Vector.h
     ${SRC_ROOT}/typeinfo/TypeInfo_Vec.h
diff --git a/Sofa/framework/DefaultType/src/sofa/defaulttype/DataTypeInfo.h b/Sofa/framework/DefaultType/src/sofa/defaulttype/DataTypeInfo.h
index 79562184484..302989871a3 100644
--- a/Sofa/framework/DefaultType/src/sofa/defaulttype/DataTypeInfo.h
+++ b/Sofa/framework/DefaultType/src/sofa/defaulttype/DataTypeInfo.h
@@ -35,6 +35,7 @@
 #include <sofa/defaulttype/typeinfo/TypeInfo_Mat.h>
 #include <sofa/defaulttype/typeinfo/TypeInfo_Quat.h>
 #include <sofa/defaulttype/typeinfo/TypeInfo_Scalar.h>
+#include <sofa/defaulttype/typeinfo/TypeInfo_SelectableItem.h>
 #include <sofa/defaulttype/typeinfo/TypeInfo_Set.h>
 #include <sofa/defaulttype/typeinfo/TypeInfo_Text.h>
 #include <sofa/defaulttype/typeinfo/TypeInfo_Vec.h>
diff --git a/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/DataTypeInfo.h b/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/DataTypeInfo.h
index fb34f4d3cf3..1bfb0913c93 100644
--- a/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/DataTypeInfo.h
+++ b/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/DataTypeInfo.h
@@ -67,7 +67,7 @@ namespace sofa::defaulttype
     \see AbstractTypeInfo provides similar mechanisms to manipulate Data objects
     generically in non-template code.
 */
-template<class TDataType>
+template<class TDataType, typename Enable = void>
 struct DataTypeInfo;
 
 template<class TDataType>
@@ -140,7 +140,7 @@ struct DefaultDataTypeInfo
     static const std::string GetTypeName() { return sofa::helper::NameDecoder::decodeTypeName(typeid(DataType)); }
 };
 
-template<class TDataType>
+template<class TDataType, typename Enable>
 struct DataTypeInfo : public DefaultDataTypeInfo<TDataType>
 {
 };
diff --git a/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/TypeInfo_SelectableItem.h b/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/TypeInfo_SelectableItem.h
new file mode 100644
index 00000000000..cd7e2f5449b
--- /dev/null
+++ b/Sofa/framework/DefaultType/src/sofa/defaulttype/typeinfo/TypeInfo_SelectableItem.h
@@ -0,0 +1,50 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/helper/SelectableItem.h>
+#include <sofa/defaulttype/typeinfo/DataTypeInfo.h>
+
+namespace sofa::defaulttype
+{
+
+template<class TDataType>
+struct SelectableItemDataTypeInfo : defaulttype::DefaultDataTypeInfo<TDataType>
+{
+    enum { ValidInfo       = 1 /**< 1 if this type has valid infos*/ };
+    static const std::string name()
+    {
+        return "SelectableItem";
+    }
+    static const std::string GetTypeName()
+    {
+        return "SelectableItem";
+    }
+};
+
+
+template<class T>
+struct DataTypeInfo<T, std::enable_if_t<std::is_base_of_v<helper::BaseSelectableItem, T>>> :
+    SelectableItemDataTypeInfo<T>
+{
+};
+
+}
diff --git a/Sofa/framework/Helper/CMakeLists.txt b/Sofa/framework/Helper/CMakeLists.txt
index c35b5b41344..0fc3c9d2650 100644
--- a/Sofa/framework/Helper/CMakeLists.txt
+++ b/Sofa/framework/Helper/CMakeLists.txt
@@ -42,6 +42,7 @@ set(HEADER_FILES
     ${SRC_ROOT}/RandomGenerator.h
     ${SRC_ROOT}/SimpleTimer.h
     ${SRC_ROOT}/ScopedAdvancedTimer.h
+    ${SRC_ROOT}/SelectableItem.h
     ${SRC_ROOT}/SortedPermutation.h
     ${SRC_ROOT}/StringUtils.h
     ${SRC_ROOT}/TagFactory.h
diff --git a/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp b/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp
index 0386afd1323..5bb28b0964e 100644
--- a/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp
+++ b/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp
@@ -759,7 +759,8 @@ const std::map< std::string, Renamed, std::less<> > renamedComponents = {
     {"UniformConstraint", Renamed("v24.06","v25.06","UniformLagrangianConstraint")},
     {"UnilateralInteractionConstraint", Renamed("v24.06","v25.06","UnilateralLagrangianConstraint")},
     {"StiffSpringForceField", Renamed("v24.06","v25.06","SpringForceField")},
-    {"ParallelStiffSpringForceField", Renamed("v24.06","v25.06","ParallelSpringForceField")}
+    {"ParallelStiffSpringForceField", Renamed("v24.06","v25.06","ParallelSpringForceField")},
+    {"ShewchukPCGLinearSolver", Renamed("v24.12","v25.12","PCGLinearSolver")}
 
 };
 
@@ -799,7 +800,6 @@ const std::map< std::string, Dealiased, std::less<> > dealiasedComponents = {
     {"ConjugateGradient", Dealiased("v24.12","CGLinearSolver")},
     {"MINRESSolver", Dealiased("v24.12","MinResLinearSolver")},
     {"MinResSolver", Dealiased("v24.12","MinResLinearSolver")},
-    {"PCGLinearSolver", Dealiased("v24.12","ShewchukPCGLinearSolver")},
     {"JacobiLinearSolver", Dealiased("v24.12","JacobiPreconditioner")},
     {"JacobiSolver", Dealiased("v24.12","JacobiPreconditioner")},
     {"SSORLinearSolver", Dealiased("v24.12","SSORPreconditioner")},
diff --git a/Sofa/framework/Helper/src/sofa/helper/SelectableItem.h b/Sofa/framework/Helper/src/sofa/helper/SelectableItem.h
new file mode 100644
index 00000000000..87a32ecedd9
--- /dev/null
+++ b/Sofa/framework/Helper/src/sofa/helper/SelectableItem.h
@@ -0,0 +1,328 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <ostream>
+#include <string>
+#include <sofa/helper/StringUtils.h>
+#include <sofa/helper/logging/ComponentInfo.h>
+#include <sofa/helper/logging/Messaging.h>
+#include <sofa/type/fixed_array.h>
+#include <map>
+
+
+namespace sofa::helper
+{
+
+struct Item
+{
+    std::string_view key;
+    std::string_view description;
+};
+
+/**
+ * Type used to define the replacement key of a deprecated key
+ */
+struct DeprecatedItem
+{
+    std::string_view replacementKey;
+    std::string_view deprecationMessage;
+};
+
+
+#define MAKE_SELECTABLE_ITEMS_BEGIN(ClassName, ...) \
+    struct ClassName final : sofa::helper::SelectableItem<ClassName> \
+    { \
+        using sofa::helper::SelectableItem<ClassName>::SelectableItem; \
+        using sofa::helper::SelectableItem<ClassName>::operator=; \
+        using sofa::helper::SelectableItem<ClassName>::operator==; \
+        static constexpr std::array s_items { \
+            __VA_ARGS__ \
+        }; \
+        static_assert(std::is_same_v<decltype(s_items)::value_type, sofa::helper::Item>); \
+
+/**
+ * Helper macro to ease the creation of a class derived from SelectableItem
+ *
+ * @param ClassName Name to give to the class derived from SelectableItem
+ *
+ * Example:
+ * MAKE_SELECTABLE_ITEMS(ResolutionMethod,
+ *      sofa::helper::Item{"ProjectedGaussSeidel", "Projected Gauss-Seidel"},
+ *      sofa::helper::Item{"UnbuiltGaussSeidel", "Gauss-Seidel no matrix assembly"},
+ *      sofa::helper::Item{"NonsmoothNonlinearConjugateGradient", "Non-smooth non-linear conjugate gradient"}
+ *  );
+ */
+#define MAKE_SELECTABLE_ITEMS(ClassName, ...) \
+    MAKE_SELECTABLE_ITEMS_BEGIN(ClassName, __VA_ARGS__) }
+
+/**
+ * Version of the helper macro where a deprecation map is declared. Warning: It is not defined
+ */
+#define MAKE_SELECTABLE_ITEMS_WITH_DEPRECATION(ClassName, ...) \
+    MAKE_SELECTABLE_ITEMS_BEGIN(ClassName, __VA_ARGS__) \
+    static const std::map<std::string_view, sofa::helper::DeprecatedItem> s_deprecationMap; \
+    }
+
+struct BaseSelectableItem
+{
+    [[nodiscard]] virtual std::size_t getNumberOfItems() const { return 0; }
+    [[nodiscard]] virtual const Item* getItemsData() const { return nullptr; }
+    [[nodiscard]] virtual std::size_t getSelectedId() const { return 0; }
+    virtual void setSelectedId(std::size_t id) {}
+    virtual void print(std::ostream& out) const {}
+    virtual void readFromStream(std::istream& stream) {}
+};
+
+inline std::ostream& operator<<(std::ostream& out, const BaseSelectableItem& selectableItem)
+{
+    selectableItem.print(out);
+    return out;
+}
+
+inline std::istream& operator>>(std::istream& in, BaseSelectableItem& selectableItem)
+{
+    selectableItem.readFromStream(in);
+    return in;
+}
+
+template<class Derived, typename IdType>
+class SelectableItem;
+
+namespace logging
+{
+
+template<class Derived, typename IdType>
+bool notMuted(const SelectableItem<Derived, IdType>* );
+
+template<class Derived, typename IdType>
+ComponentInfo::SPtr getComponentInfo(const SelectableItem<Derived, IdType>* );
+
+}
+
+/**
+ * Selection of an item among a fixed list of items.
+ *
+ * The class is designed to have the list of items static, so the class can be
+ * used in a constexpr context. In particular, it can be used in a switch
+ * statement.
+ *
+ * The helper macro @MAKE_SELECTABLE_ITEMS can be used to make the code more concise.
+ *
+ * @tparam Derived The class derived from SelectableItem. The CRTP must be used.
+ * @tparam IdType A type for the indices
+ */
+template<class Derived, typename IdType = std::size_t>
+class SelectableItem : public BaseSelectableItem
+{
+public:
+
+    static constexpr std::size_t numberOfItems()
+    {
+        return std::tuple_size<decltype(Derived::s_items)>{};
+    }
+
+    static constexpr const auto &items()
+    {
+        return Derived::s_items;
+    }
+
+    using id_type = IdType;
+
+    constexpr SelectableItem() = default;
+    constexpr SelectableItem(const std::string_view key) : m_selected_id(findId(key))
+    {
+        if (m_selected_id >= numberOfItems())
+        {
+            keyError(key);
+        }
+    }
+
+    [[nodiscard]] std::string_view key() const
+    {
+        return Derived::s_items[m_selected_id].key;
+    }
+
+    [[nodiscard]] std::string_view description() const
+    {
+        return Derived::s_items[m_selected_id].description;
+    }
+
+    explicit operator std::string_view() const
+    {
+        return key();
+    }
+
+    constexpr operator id_type() const
+    {
+        return m_selected_id;
+    }
+
+    [[nodiscard]] bool operator==(const std::string_view key)
+    {
+        return key == Derived::s_items[m_selected_id].key;
+    }
+
+    [[nodiscard]] bool operator==(const SelectableItem& other) const
+    {
+        return m_selected_id == other.m_selected_id;
+    }
+
+    [[nodiscard]] bool operator!=(const SelectableItem& other) const
+    {
+        return m_selected_id != other.m_selected_id;
+    }
+
+    constexpr SelectableItem& operator=(const std::string_view key)
+    {
+        m_selected_id = findId(key);
+        if (m_selected_id >= numberOfItems())
+        {
+            keyError(key);
+        }
+        return *this;
+    }
+
+    [[nodiscard]] static const std::string& allKeysAsString()
+    {
+        static std::string allKeys = ::sofa::helper::join(Derived::s_items.begin(), Derived::s_items.end(),
+            [](const Item& item){ return item.key;}, ',');
+        return allKeys;
+    }
+
+    static const std::string& dataDescription()
+    {
+        static const std::string dataDescription =
+            sofa::helper::join(Derived::s_items.begin(), Derived::s_items.end(),
+            [](const Item& item)
+            {
+                return "- " + std::string{item.key} + ": " + std::string{item.description};
+            }, '\n');
+        return dataDescription;
+    }
+
+    void print(std::ostream& out) const final
+    {
+        out << key();
+    }
+
+    void readFromStream(std::istream& stream) final
+    {
+        std::string tmp;
+        std::getline(stream, tmp);
+        m_selected_id = findId(tmp);
+        if (m_selected_id >= numberOfItems())
+        {
+            keyError(tmp);
+        }
+    }
+
+protected:
+    id_type m_selected_id {};
+
+    static constexpr id_type findId(const std::string_view key)
+    {
+        return findId_impl(key, std::make_index_sequence<numberOfItems()>{});
+    }
+
+    template<id_type... Is>
+    static constexpr id_type findId_impl(const std::string_view key, std::index_sequence<Is...>)
+    {
+        id_type result = static_cast<id_type>(-1);
+        ((Derived::s_items[Is].key == key ? (result = Is, true) : false) || ...);
+        return result;
+    }
+
+    template <typename, typename = std::void_t<>>
+    struct has_deprecation_map : std::false_type {};
+
+    // Specialization when T has a static member s_foo
+    template <typename T>
+    struct has_deprecation_map<T, std::void_t<decltype(T::s_deprecationMap)>> : std::true_type {};
+
+    void keyError(const std::string_view key)
+    {
+        if constexpr (has_deprecation_map<Derived>::value)
+        {
+            static_assert(std::is_same_v<std::remove_cv_t<decltype(Derived::s_deprecationMap)>, std::map<std::string_view, DeprecatedItem>>);
+            const auto it = Derived::s_deprecationMap.find(key);
+            if (it != Derived::s_deprecationMap.end())
+            {
+                if (it->second.replacementKey == key)
+                {
+                    dmsg_fatal() << "Item '" << key << "' is deprecated, but new key is also '" << it->second.replacementKey << "'.";
+                }
+                else
+                {
+                    msg_warning() << "Item '" << key << "' is deprecated. New key to use instead is '" << it->second.replacementKey << "'. " << it->second.deprecationMessage;
+                    this->operator=(it->second.replacementKey);
+                    return;
+                }
+            }
+        }
+
+        static_assert(numberOfItems() > 0);
+        msg_error() << "Item '" << key << "' does not exist in the list ["
+                << this->allKeysAsString() << "]. Fall back to item '"
+                << items().front().key << "'";
+        m_selected_id = 0;
+    }
+
+private:
+    [[nodiscard]] std::size_t getNumberOfItems() const final
+    {
+        return numberOfItems();
+    }
+
+    [[nodiscard]] const Item* getItemsData() const final
+    {
+        return Derived::s_items.data();
+    }
+
+    [[nodiscard]] std::size_t getSelectedId() const final
+    {
+        return m_selected_id;
+    }
+
+    void setSelectedId(std::size_t id) final
+    {
+        m_selected_id = id;
+    }
+};
+
+namespace logging
+{
+
+template<class Derived, typename IdType>
+bool notMuted(const SelectableItem<Derived, IdType>* )
+{
+    return true;
+}
+
+template<class Derived, typename IdType>
+ComponentInfo::SPtr getComponentInfo(const SelectableItem<Derived, IdType>* )
+{
+    return std::make_shared<ComponentInfo>("SelectableItem");
+}
+}
+
+}
+
diff --git a/Sofa/framework/Helper/src/sofa/helper/system/DynamicLibrary.cpp b/Sofa/framework/Helper/src/sofa/helper/system/DynamicLibrary.cpp
index feb08471881..fa90b3195f2 100644
--- a/Sofa/framework/Helper/src/sofa/helper/system/DynamicLibrary.cpp
+++ b/Sofa/framework/Helper/src/sofa/helper/system/DynamicLibrary.cpp
@@ -111,7 +111,7 @@ void * DynamicLibrary::getSymbolAddress(Handle handle,
 # endif
     if(symbolAddress == nullptr)
         fetchLastError();
-#if not defined (WIN32)
+#if !defined (WIN32)
     else // checking that the symbol really comes from the provided library
     {
 
diff --git a/Sofa/framework/Helper/test/CMakeLists.txt b/Sofa/framework/Helper/test/CMakeLists.txt
index 4189c2da7e1..b6120ff71b9 100644
--- a/Sofa/framework/Helper/test/CMakeLists.txt
+++ b/Sofa/framework/Helper/test/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SOURCE_FILES
     KdTree_test.cpp
     NameDecoder_test.cpp
     OptionsGroup_test.cpp
+    SelectableItem_test.cpp
     StringUtils_test.cpp
     TagFactory_test.cpp
     Utils_test.cpp
diff --git a/Sofa/framework/Helper/test/OptionsGroup_test.cpp b/Sofa/framework/Helper/test/OptionsGroup_test.cpp
index 77c4f56e777..da90a1d0647 100644
--- a/Sofa/framework/Helper/test/OptionsGroup_test.cpp
+++ b/Sofa/framework/Helper/test/OptionsGroup_test.cpp
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #include <sofa/helper/OptionsGroup.h>
 #include <gtest/gtest.h>
+#include <sofa/defaulttype/typeinfo/DataTypeInfo.h>
+
 
 namespace sofa
 {
@@ -50,4 +52,9 @@ TEST(OptionsGroup, constructors)
     EXPECT_EQ(opt5.getSelectedItem(), "optionB");
 }
 
+TEST(OptionsGroup, DataTypeInfo)
+{
+    EXPECT_EQ(defaulttype::DataTypeInfo<OptionsGroup>::GetTypeName(), "OptionsGroup");
+}
+
 }
diff --git a/Sofa/framework/Helper/test/SelectableItem_test.cpp b/Sofa/framework/Helper/test/SelectableItem_test.cpp
new file mode 100644
index 00000000000..f9a023ed496
--- /dev/null
+++ b/Sofa/framework/Helper/test/SelectableItem_test.cpp
@@ -0,0 +1,214 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <sofa/helper/SelectableItem.h>
+#include <gtest/gtest.h>
+#include <sofa/testing/TestMessageHandler.h>
+#include <sofa/core/objectmodel/Data.h>
+
+
+struct TestSelectableItem final : sofa::helper::SelectableItem<TestSelectableItem>
+{
+    using Inherited = sofa::helper::SelectableItem<TestSelectableItem>;
+    using Inherited::SelectableItem;
+    using Inherited::operator=;
+
+    static constexpr std::array s_items {
+        sofa::helper::Item{"small", "small displacements"},
+        sofa::helper::Item{"large", "large displacements"},
+    };
+
+    inline static const std::map<std::string_view, sofa::helper::DeprecatedItem> s_deprecationMap {
+        {"deprecated_key", sofa::helper::DeprecatedItem{"large", "Now use 'large'"}}
+    };
+
+    [[nodiscard]] id_type selectedId() const
+    {
+        return m_selected_id;
+    }
+};
+
+MAKE_SELECTABLE_ITEMS(TestSelectableItemNoDeprecation,
+    sofa::helper::Item{"small", "small displacements"},
+    sofa::helper::Item{"large", "large displacements"});
+
+MAKE_SELECTABLE_ITEMS_WITH_DEPRECATION(TestSelectableItemWithDeprecation,
+    sofa::helper::Item{"small", "small displacements"},
+    sofa::helper::Item{"large", "large displacements"});
+
+const std::map<std::string_view, sofa::helper::DeprecatedItem> TestSelectableItemWithDeprecation::s_deprecationMap{
+    {"deprecated_key", sofa::helper::DeprecatedItem{"large", "Now use 'large'"}}
+};
+
+TEST(SelectableItem, numberOfItems)
+{
+    static constexpr auto size = TestSelectableItem::numberOfItems();
+    EXPECT_EQ(size, 2);
+}
+
+TEST(SelectableItem, allKeysAsString)
+{
+    static const auto& allKeys = TestSelectableItem::allKeysAsString();
+    EXPECT_EQ(allKeys, "small,large");
+}
+
+TEST(SelectableItem, operator_equal_string_view)
+{
+    static constexpr std::string_view key_small = "small";
+
+    TestSelectableItem foo;
+    foo = key_small;
+
+    EXPECT_EQ(foo.selectedId(), 0);
+}
+
+TEST(SelectableItem, operator_equal_string)
+{
+    static std::string key_large = "large";
+
+    TestSelectableItem foo;
+    foo = key_large;
+
+    EXPECT_EQ(foo.selectedId(), 1);
+}
+
+TEST(SelectableItem, description)
+{
+    TestSelectableItem foo;
+    foo = "large";
+
+    EXPECT_EQ(foo.description(), "large displacements");
+}
+
+TEST(SelectableItem, operator_equal_to)
+{
+    TestSelectableItem foo;
+    foo = "large";
+
+    const bool is_equal = (foo == "large");
+    EXPECT_TRUE(is_equal);
+}
+
+TEST(SelectableItem, convertion_string_view)
+{
+    TestSelectableItem foo;
+    foo = "large";
+
+    const auto key = static_cast<std::string_view>(foo);
+    EXPECT_EQ(key, "large");
+}
+
+TEST(SelectableItem, stream_out)
+{
+    TestSelectableItem foo;
+    foo = "large";
+
+    std::stringstream ss;
+    ss << foo;
+
+    EXPECT_EQ(ss.str(), "large");
+}
+
+TEST(SelectableItem, stream_in)
+{
+    TestSelectableItem foo;
+
+    std::istringstream input_stream("large");
+    input_stream >> foo;
+
+    EXPECT_EQ(foo.selectedId(), 1);
+}
+
+TEST(SelectableItem, constexpr_constructor)
+{
+    static constexpr TestSelectableItem foo("large");
+    EXPECT_EQ(foo.selectedId(), 1);
+
+    // static constexpr TestSelectableItem bar("bar"); //does not compile because "bar" does not exist in the list
+}
+
+TEST(SelectableItem, wrong_key)
+{
+    // required to be able to use EXPECT_MSG_NOEMIT and EXPECT_MSG_EMIT
+    sofa::helper::logging::MessageDispatcher::addHandler(sofa::testing::MainGtestMessageHandler::getInstance() ) ;
+
+    TestSelectableItem foo;
+
+    {
+        EXPECT_MSG_EMIT(Error);
+        foo = "foo";
+    }
+
+    EXPECT_EQ(foo.selectedId(), 0);
+}
+
+TEST(SelectableItem, data)
+{
+    const sofa::Data<TestSelectableItem> data;
+    EXPECT_EQ(data.getValueTypeString(), "SelectableItem");
+}
+
+TEST(SelectableItem, DataTypeInfo_BaseSelectableItems)
+{
+    EXPECT_EQ(sofa::defaulttype::DataTypeInfo<sofa::helper::BaseSelectableItem>::GetTypeName(), "SelectableItem");
+}
+
+TEST(SelectableItem, DataTypeInfoValidInfo_SelectableItem)
+{
+    EXPECT_TRUE(sofa::defaulttype::DataTypeInfo<TestSelectableItem>::ValidInfo);
+}
+
+TEST(SelectableItem, DataTypeInfo_SelectableItem)
+{
+    EXPECT_EQ(sofa::defaulttype::DataTypeInfo<TestSelectableItem>::GetTypeName(), "SelectableItem");
+}
+
+TEST(SelectableItem, BaseData_typeName)
+{
+    EXPECT_EQ(sofa::core::objectmodel::BaseData::typeName<TestSelectableItem>(), "SelectableItem");
+}
+
+TEST(SelectableItem, dynamic_cast_base)
+{
+    TestSelectableItem foo;
+    EXPECT_NE(nullptr, dynamic_cast<sofa::helper::BaseSelectableItem*>(&foo));
+}
+
+template<class T>
+void testDeprecatedKey()
+{
+    // required to be able to use EXPECT_MSG_NOEMIT and EXPECT_MSG_EMIT
+    sofa::helper::logging::MessageDispatcher::addHandler(sofa::testing::MainGtestMessageHandler::getInstance() ) ;
+
+    T foo;
+
+    EXPECT_MSG_EMIT(Warning);
+    foo = "deprecated_key";
+
+    const auto key = static_cast<std::string_view>(foo);
+    EXPECT_EQ(key, "large");
+}
+
+TEST(SelectableItem, deprecated_key)
+{
+    testDeprecatedKey<TestSelectableItem>();
+    testDeprecatedKey<TestSelectableItemWithDeprecation>();
+}
diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h
index 0ff0e0455a9..e3a6d3db154 100644
--- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h
+++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h
@@ -620,31 +620,25 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress
     /// write to an output stream
     inline friend std::ostream& operator << ( std::ostream& out, const CompressedRowSparseMatrixConstraint<TBlock, Policy>& sc)
     {
-        for (RowConstIterator rowIt = sc.begin(); rowIt !=  sc.end(); ++rowIt)
+        std::ostringstream ossrow;
+        std::size_t nbLines = 0;
+        for (RowConstIterator rowIt = sc.begin(); rowIt != sc.end(); ++rowIt)
         {
-            out << "Constraint ID : ";
-            out << rowIt.index();
-            const auto colToString = [](const ColConstIterator& colIt)
-            {
-                std::stringstream ss;
-                ss << "dof ID : " << colIt.index() << "  value : " << colIt.val();
-                return ss.str();
-            };
+            ossrow << rowIt.index() << " ";
 
-            ColConstIterator colIt = rowIt.begin();
-            const ColConstIterator colItEnd = rowIt.end();
-            if (colIt != colItEnd)
+            std::ostringstream ossline;
+            std::size_t n = 0;
+            for (ColConstIterator colIt = rowIt.begin(); colIt != rowIt.end(); ++colIt)
             {
-                out << "  " << colToString(colIt++);
-                while(colIt != colItEnd)
-                {
-                    out << "  " << colToString(colIt++);
-                }
+                ossline << colIt.index() << " " << colIt.val() << " ";
+                n++;
             }
-
-            out << "\n";
+            ossrow << n << " " << ossline.str();
+            nbLines++ ;
         }
 
+        out << nbLines << " " << ossrow.str();
+
         return out;
     }
 
@@ -653,16 +647,25 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress
     {
         sc.clear();
 
+        unsigned int nbLines;
         unsigned int c_id;
         unsigned int c_number;
         unsigned int c_dofIndex;
         TBlock c_value;
 
-        while (!(in.rdstate() & std::istream::eofbit))
+        if (in.rdstate() & std::istream::eofbit)
+        {
+            return in;
+        }
+
+        in >> nbLines;
+
+        unsigned int currentNbLines = 0;
+        while (currentNbLines < nbLines && !(in.rdstate() & std::istream::eofbit))
         {
             in >> c_id;
             in >> c_number;
-
+         
             auto c_it = sc.writeLine(c_id);
 
             for (unsigned int i = 0; i < c_number; i++)
@@ -671,11 +674,44 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress
                 in >> c_value;
                 c_it.addCol(c_dofIndex, c_value);
             }
+            currentNbLines++;
         }
 
+        assert(nbLines == currentNbLines);
+
         sc.compress();
+
         return in;
     }
+    
+    /// write into output stream (default is standard output)
+    void prettyPrint(std::ostream& out = std::cout) const
+    {
+        for (RowConstIterator rowIt = this->begin(); rowIt !=  this->end(); ++rowIt)
+        {
+            out << "Constraint ID : ";
+            out << rowIt.index();
+            const auto colToString = [](const ColConstIterator& colIt)
+            {
+                std::stringstream ss;
+                ss << "dof ID : " << colIt.index() << "  value : " << colIt.val();
+                return ss.str();
+            };
+
+            ColConstIterator colIt = rowIt.begin();
+            const ColConstIterator colItEnd = rowIt.end();
+            if (colIt != colItEnd)
+            {
+                out << "  " << colToString(colIt++);
+                while(colIt != colItEnd)
+                {
+                    out << "  " << colToString(colIt++);
+                }
+            }
+
+            out << "\n";
+        }
+    }
 
     static const char* Name()
     {
diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp
index bb5a6fe5a5a..9b164ce6bd1 100644
--- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp
+++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp
@@ -281,6 +281,51 @@ TYPED_TEST(SparseMatrixTest, CheckThatTheSizeOfASparseMatrixIsOneAfterTryingToWr
     EXPECT_EQ(1u, matrix.size());
 }
 
+
+//////////////////////////////////////////////////
+TYPED_TEST(SparseMatrixTest, CheckOutInSerialization)
+{
+    typedef TypeParam Matrix;
+    Matrix outMatrix;
+    auto lineIndices = TestHelpers::Populate(outMatrix, 3u, 3u);
+    EXPECT_EQ(3u, outMatrix.size());
+
+    std::ostringstream oss;
+    oss << outMatrix;
+
+    Matrix inMatrix;
+    std::istringstream iss(oss.str());
+    iss >> inMatrix;
+
+    EXPECT_TRUE(!lineIndices.empty());
+
+    for (const auto& lineIndex : lineIndices)
+    {
+        auto inRowIt = inMatrix.readLine(lineIndex);
+        auto outRowIt = outMatrix.readLine(lineIndex);
+
+        EXPECT_TRUE(inRowIt == outRowIt);
+
+        if (inRowIt != inMatrix.end() && outRowIt != outMatrix.end())
+        {
+            auto inColItEnd = inRowIt.end();
+            auto outColItEnd = outRowIt.end();
+
+            auto inColIt = inRowIt.begin();
+            auto outColIt = outRowIt.begin();
+
+            for (; 
+                inColIt != inColItEnd && outColIt != outColItEnd;
+                ++inColIt, ++outColIt)
+            {
+                EXPECT_TRUE(inColIt.index() == outColIt.index());
+                EXPECT_TRUE(inColIt.val() == outColIt.val());
+            }
+        }
+    }
+
+}
+
 //////////////////////////////////////////////////
 TYPED_TEST(SparseMatrixTest, CheckThatAMatrixIsConsideredEmptyAfterIsHasBeenCleared)
 {
@@ -2121,14 +2166,29 @@ TEST(CompressedRowSparseMatrixConstraint, ostream)
     std::stringstream ss;
     ss << A;
 
+    static const std::string expectedOutput =
+R"(4 0 2 1 0 0.360985 0 12 0.926981 0 0 2 1 9 0.451858 0 0 3 1 6 0.777417 0 0 4 3 5 0 0 0.474108 7 0 0.983937 0 9 0.238781 0 0 )";
+
+    EXPECT_EQ(ss.str(), expectedOutput);
+}
+
+TEST(CompressedRowSparseMatrixConstraint, prettyPrint)
+{
+    sofa::linearalgebra::CompressedRowSparseMatrixConstraint<sofa::type::Vec3> A;
+
+    generateMatrix(A, 5, 15, 0.1, 7);
+
     static const std::string expectedOutput =
 R"(Constraint ID : 0  dof ID : 1  value : 0 0.360985 0  dof ID : 12  value : 0.926981 0 0
 Constraint ID : 2  dof ID : 9  value : 0.451858 0 0
 Constraint ID : 3  dof ID : 6  value : 0.777417 0 0
 Constraint ID : 4  dof ID : 5  value : 0 0 0.474108  dof ID : 7  value : 0 0.983937 0  dof ID : 9  value : 0.238781 0 0
 )";
+    
+    std::ostringstream oss;
+    A.prettyPrint(oss);
 
-    EXPECT_EQ(ss.str(), expectedOutput);
+    EXPECT_EQ(oss.str(), expectedOutput);
 }
 
 } // namespace sofa
diff --git a/Sofa/framework/SimpleApi/test/SimpleApi_test.cpp b/Sofa/framework/SimpleApi/test/SimpleApi_test.cpp
index dc5ae56c599..f9c26dfccd7 100644
--- a/Sofa/framework/SimpleApi/test/SimpleApi_test.cpp
+++ b/Sofa/framework/SimpleApi/test/SimpleApi_test.cpp
@@ -90,3 +90,55 @@ TEST_F(SimpleApi_test, createParamString )
 {
     ASSERT_TRUE( testParamString() );
 }
+
+TEST(SimpleApi_test_solo, testIsSetWithDataLink)
+{
+    const Simulation::SPtr simu = createSimulation("DAG");
+    const Node::SPtr root = createRootNode(simu, "root");
+
+    // test not set
+    const auto obj1 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop1"}
+    });
+    auto* objdata1 = obj1->findData("printLog");
+    ASSERT_FALSE(objdata1->isSet());
+
+    // test set
+    const auto obj2 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop2"}, 
+        {"printLog", "false"}
+    });
+    auto* objdata2 = obj2->findData("printLog");
+    ASSERT_TRUE(objdata2->isSet());
+
+    // test set through a link of a already created object
+    const auto obj3 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop3"},
+        {"printLog", "@/loop2.printLog"}
+    });
+    auto* objdata3 = obj3->findData("printLog");
+    ASSERT_TRUE(objdata3->isSet());
+
+    // test set through a link of a not yet created object (deferred linking)
+    const auto obj4 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop4"},
+        {"printLog", "@/loop5.printLog"}
+    });
+    const auto obj5 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop5"},
+        {"printLog", "true"}
+    });
+    auto* objdata4 = obj4->findData("printLog");
+    ASSERT_TRUE(objdata4->isSet());
+    
+    // test link with a wrong path (or non existent parent)
+    const auto obj6 = createObject(root, "DefaultAnimationLoop", {
+        {"name", "loop6"},
+        {"printLog", "@/loop7.printLog"}
+    });
+
+    auto* objdata6 = obj6->findData("printLog");
+    ASSERT_TRUE(objdata6->isSet());
+    ASSERT_EQ(objdata6->getParent(), nullptr);
+
+}
diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp
index 9167f002344..a050b22b631 100644
--- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp
+++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp
@@ -288,7 +288,11 @@ void MechanicalOperations::computeDfV(core::MultiVecDerivId df, bool clear, bool
 }
 
 /// accumulate $ df += (m M + b B + k K) dx $ (given the latest propagated displacement)
-void MechanicalOperations::addMBKdx(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear, bool accumulate)
+void MechanicalOperations::addMBKdx(core::MultiVecDerivId df,
+                                    const MatricesFactors::M m,
+                                    const MatricesFactors::B b,
+                                    const MatricesFactors::K k,
+                                    const bool clear, const bool accumulate)
 {
     setDf(df);
     if (clear)
@@ -296,14 +300,18 @@ void MechanicalOperations::addMBKdx(core::MultiVecDerivId df, SReal m, SReal b,
         executeVisitor( MechanicalResetForceVisitor(&mparams, df, true) );
         //finish();
     }
-    mparams.setBFactor(b);
-    mparams.setKFactor(k);
-    mparams.setMFactor(m);
+    mparams.setBFactor(b.get());
+    mparams.setKFactor(k.get());
+    mparams.setMFactor(m.get());
     executeVisitor( MechanicalAddMBKdxVisitor(&mparams, df, accumulate) );
 }
 
 /// accumulate $ df += (m M + b B + k K) velocity $
-void MechanicalOperations::addMBKv(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear, bool accumulate)
+void MechanicalOperations::addMBKv(core::MultiVecDerivId df,
+                                   const MatricesFactors::M m,
+                                   const MatricesFactors::B b,
+                                   const MatricesFactors::K k,
+                                   const bool clear, const bool accumulate)
 {
     const core::ConstMultiVecDerivId dx = mparams.dx();
     mparams.setDx(mparams.v());
@@ -313,9 +321,9 @@ void MechanicalOperations::addMBKv(core::MultiVecDerivId df, SReal m, SReal b, S
         executeVisitor( MechanicalResetForceVisitor(&mparams, df, true) );
         //finish();
     }
-    mparams.setBFactor(b);
-    mparams.setKFactor(k);
-    mparams.setMFactor(m);
+    mparams.setBFactor(b.get());
+    mparams.setKFactor(k.get());
+    mparams.setMFactor(m.get());
     /* useV = true */
     executeVisitor( MechanicalAddMBKdxVisitor(&mparams, df, accumulate) );
     mparams.setDx(dx);
@@ -396,14 +404,15 @@ void MechanicalOperations::resetSystem(core::behavior::LinearSolver* linearSolve
     }
 }
 
-void MechanicalOperations::setSystemMBKMatrix(SReal mFact, SReal bFact,
-    SReal kFact, core::behavior::LinearSolver* linearSolver)
+void MechanicalOperations::setSystemMBKMatrix(
+    MatricesFactors::M m, MatricesFactors::B b, MatricesFactors::K k,
+    core::behavior::LinearSolver* linearSolver)
 {
     if (linearSolver)
     {
-        mparams.setMFactor(mFact);
-        mparams.setBFactor(bFact);
-        mparams.setKFactor(kFact);
+        mparams.setMFactor(m.get());
+        mparams.setBFactor(b.get());
+        mparams.setKFactor(k.get());
         mparams.setSupportOnlySymmetricMatrix(!linearSolver->supportNonSymmetricSystem());
         linearSolver->setSystemMBKMatrix(&mparams);
     }
@@ -518,7 +527,7 @@ void MechanicalOperations::m_setSystemMBKMatrix(SReal mFact, SReal bFact, SReal
         showMissingLinearSolverError();
         return;
     }
-    setSystemMBKMatrix(mFact, bFact, kFact, s);
+    setSystemMBKMatrix(MatricesFactors::M(mFact), MatricesFactors::B(bFact), MatricesFactors::K(kFact), s);
 }
 
 void MechanicalOperations::m_setSystemRHVector(core::MultiVecDerivId v)
@@ -622,6 +631,24 @@ void MechanicalOperations::printWithElapsedTime( core::ConstMultiVecId /*v*/, un
 {
 }
 
+void MechanicalOperations::addMBKdx(core::MultiVecDerivId df, SReal m, SReal b,
+    SReal k, bool clear, bool accumulate)
+{
+    addMBKdx(df, core::MatricesFactors::M(m), core::MatricesFactors::B(b), core::MatricesFactors::K(k), clear, accumulate);
+}
+
+void MechanicalOperations::addMBKv(core::MultiVecDerivId df, SReal m, SReal b,
+    SReal k, bool clear, bool accumulate)
+{
+    addMBKv(df, core::MatricesFactors::M(m), core::MatricesFactors::B(b), core::MatricesFactors::K(k), clear, accumulate);
+}
+
+void MechanicalOperations::setSystemMBKMatrix(SReal mFact, SReal bFact,
+    SReal kFact, core::behavior::LinearSolver* linearSolver)
+{
+    setSystemMBKMatrix(core::MatricesFactors::M(mFact), core::MatricesFactors::B(bFact), core::MatricesFactors::K(kFact), linearSolver);
+}
+
 void MechanicalOperations::showMissingLinearSolverError() const
 {
     if (!hasShownMissingLinearSolverMap[ctx])
diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h
index 6b56a17e7bf..b968b4d53b8 100644
--- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h
+++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.h
@@ -28,6 +28,8 @@
 #include <sofa/core/behavior/MultiMatrixAccessor.h>
 #include <sofa/simulation/VisitorExecuteFunc.h>
 #include <sofa/core/ConstraintParams.h>
+#include <sofa/core/MatricesFactors.h>
+
 
 namespace sofa::simulation::common
 {
@@ -81,9 +83,9 @@ class SOFA_SIMULATION_CORE_API MechanicalOperations
     /// Compute the current force delta (given the latest propagated velocity)
     void computeDfV(core::MultiVecDerivId df, bool clear = true, bool accumulate = true);
     /// accumulate $ df += (m M + b B + k K) dx $ (given the latest propagated displacement)
-    void addMBKdx(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear = true, bool accumulate = true);
+    void addMBKdx(core::MultiVecDerivId df, core::MatricesFactors::M m, core::MatricesFactors::B b, core::MatricesFactors::K k, bool clear = true, bool accumulate = true);
     /// accumulate $ df += (m M + b B + k K) velocity $
-    void addMBKv(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear = true, bool accumulate = true);
+    void addMBKv(core::MultiVecDerivId df, core::MatricesFactors::M m, core::MatricesFactors::B b, core::MatricesFactors::K k, bool clear = true, bool accumulate = true);
     /// Add dt*Gravity to the velocity
     void addSeparateGravity(SReal dt, core::MultiVecDerivId result = core::VecDerivId::velocity() );
 
@@ -101,7 +103,7 @@ class SOFA_SIMULATION_CORE_API MechanicalOperations
 /// @{
 
     void resetSystem(core::behavior::LinearSolver* linearSolver);
-    void setSystemMBKMatrix(SReal mFact, SReal bFact, SReal kFact, core::behavior::LinearSolver* linearSolver);
+    void setSystemMBKMatrix(core::MatricesFactors::M m, core::MatricesFactors::B b, core::MatricesFactors::K k, core::behavior::LinearSolver* linearSolver);
     void setSystemRHVector(core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver);
     void setSystemLHVector(core::MultiVecDerivId v, core::behavior::LinearSolver* linearSolver);
     void solveSystem(core::behavior::LinearSolver* linearSolver);
@@ -149,6 +151,14 @@ class SOFA_SIMULATION_CORE_API MechanicalOperations
 
     /// @}
 
+
+    SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKDX()
+    void addMBKdx(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear = true, bool accumulate = true);
+    SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKV()
+    void addMBKv(core::MultiVecDerivId df, SReal m, SReal b, SReal k, bool clear = true, bool accumulate = true);
+    SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_SETSYSTEMMBKMATRIX()
+    void setSystemMBKMatrix(SReal mFact, SReal bFact, SReal kFact, core::behavior::LinearSolver* linearSolver);
+
 protected:
     VisitorExecuteFunc executeVisitor;
 
diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/config.h.in b/Sofa/framework/Simulation/Core/src/sofa/simulation/config.h.in
index 51abc357090..0497c54b26d 100644
--- a/Sofa/framework/Simulation/Core/src/sofa/simulation/config.h.in
+++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/config.h.in
@@ -278,3 +278,27 @@ SOFA_ATTRIBUTE_DISABLED("v23.12", "v24.06", "rootdata/nodedata feature was never
     SOFA_ATTRIBUTE_DEPRECATED( \
     "v24.06", "v24.12", "Use print instead")
 #endif // SOFA_BUILD_SOFA_SIMULATION_CORE
+
+#ifdef SOFA_BUILD_SOFA_SIMULATION_CORE
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKDX()
+#else
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKDX() \
+    SOFA_ATTRIBUTE_DEPRECATED( \
+    "v24.12", "v25.06", "Use the other addMBKdx overload instead.")
+#endif // SOFA_BUILD_SOFA_SIMULATION_CORE
+
+#ifdef SOFA_BUILD_SOFA_SIMULATION_CORE
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKV()
+#else
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_ADDMBKV() \
+    SOFA_ATTRIBUTE_DEPRECATED( \
+    "v24.12", "v25.06", "Use the other addMBKv overload instead.")
+#endif // SOFA_BUILD_SOFA_SIMULATION_CORE
+
+#ifdef SOFA_BUILD_SOFA_SIMULATION_CORE
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_SETSYSTEMMBKMATRIX()
+#else
+#define SOFA_ATTRIBUTE_DEPRECATED_MECHANICALOPERATIONS_SETSYSTEMMBKMATRIX() \
+    SOFA_ATTRIBUTE_DEPRECATED( \
+    "v24.12", "v25.06", "Use the other setSystemMBKMatrix overload instead.")
+#endif // SOFA_BUILD_SOFA_SIMULATION_CORE
diff --git a/Sofa/framework/Type/CMakeLists.txt b/Sofa/framework/Type/CMakeLists.txt
index 533ecda0f47..f44dd533442 100644
--- a/Sofa/framework/Type/CMakeLists.txt
+++ b/Sofa/framework/Type/CMakeLists.txt
@@ -34,6 +34,7 @@ set(HEADER_FILES
     ${SOFATYPESRC_ROOT}/fwd.h
     ${SOFATYPESRC_ROOT}/isRigidType.h
     ${SOFATYPESRC_ROOT}/stable_vector.h
+    ${SOFATYPESRC_ROOT}/StrongType.h
     ${SOFATYPESRC_ROOT}/trait/Rebind.h
     ${SOFATYPESRC_ROOT}/trait/is_container.h
     ${SOFATYPESRC_ROOT}/trait/is_fixed_array.h
diff --git a/Sofa/framework/Type/src/sofa/type/RGBAColor.cpp b/Sofa/framework/Type/src/sofa/type/RGBAColor.cpp
index 611ed9d0687..79c4aebe25a 100644
--- a/Sofa/framework/Type/src/sofa/type/RGBAColor.cpp
+++ b/Sofa/framework/Type/src/sofa/type/RGBAColor.cpp
@@ -23,6 +23,7 @@
 
 #include <sstream>
 #include <locale>
+#include <map>
 
 #include <sofa/type/fixed_array.h>
 #include <sofa/type/fixed_array_algorithms.h>
@@ -176,6 +177,30 @@ static std::istream& trimInitialSpaces(std::istream& in)
     return in;
 }
 
+const std::map<std::string, RGBAColor> stringToColorMap {
+    {"white", g_white},
+    {"black", g_black},
+    {"red", g_red},
+    {"green", g_green},
+    {"blue", g_blue},
+    {"cyan", g_cyan},
+    {"magenta", g_magenta},
+    {"yellow", g_yellow},
+    {"gray", g_gray},
+    {"darkgray", g_darkgray},
+    {"lightgray", g_lightgray},
+    {"orange", g_orange},
+    {"purple", g_purple},
+    {"pink", g_pink},
+    {"brown", g_brown},
+    {"lime", g_lime},
+    {"teal", g_teal},
+    {"navy", g_navy},
+    {"olive", g_olive},
+    {"maroon", g_maroon},
+    {"silver", g_silver},
+    {"gold", g_gold}
+};
 
 SOFA_TYPE_API std::istream& operator>>(std::istream& i, RGBAColor& t)
 {
@@ -230,18 +255,14 @@ SOFA_TYPE_API std::istream& operator>>(std::istream& i, RGBAColor& t)
 
         /// if end of line is returned before encountering ' ' or 7... is it fine
         /// so we can clear the failure bitset.
-
-        /// Compare the resulting string with supported colors.
-        /// If you add more colors... please also add them in the test file.
-        if (str == "white")    { r = 1.0f; g = 1.0f; b = 1.0f; }
-        else if (str == "black")    { r = 0.0f; g = 0.0f; b = 0.0f; }
-        else if (str == "red")      { r = 1.0f; g = 0.0f; b = 0.0f; }
-        else if (str == "green")    { r = 0.0f; g = 1.0f; b = 0.0f; }
-        else if (str == "blue")     { r = 0.0f; g = 0.0f; b = 1.0f; }
-        else if (str == "cyan")     { r = 0.0f; g = 1.0f; b = 1.0f; }
-        else if (str == "magenta")  { r = 1.0f; g = 0.0f; b = 1.0f; }
-        else if (str == "yellow")   { r = 1.0f; g = 1.0f; b = 0.0f; }
-        else if (str == "gray")     { r = 0.5f; g = 0.5f; b = 0.5f; }
+        if (const auto it = stringToColorMap.find(str);
+            it != stringToColorMap.end())
+        {
+            r = it->second.r();
+            g = it->second.g();
+            b = it->second.b();
+            a = it->second.a();
+        }
         else {
             /// If we cannot parse the field we returns that with the fail bit.
             i.setstate(std::ios_base::failbit) ;
diff --git a/Sofa/framework/Type/src/sofa/type/RGBAColor.h b/Sofa/framework/Type/src/sofa/type/RGBAColor.h
index 13e2bbbf70a..44fefc4eacc 100644
--- a/Sofa/framework/Type/src/sofa/type/RGBAColor.h
+++ b/Sofa/framework/Type/src/sofa/type/RGBAColor.h
@@ -91,6 +91,17 @@ class SOFA_TYPE_API RGBAColor
     constexpr static const RGBAColor& gray();
     constexpr static const RGBAColor& darkgray();
     constexpr static const RGBAColor& lightgray();
+    constexpr static const RGBAColor& orange();
+    constexpr static const RGBAColor& purple();
+    constexpr static const RGBAColor& pink();
+    constexpr static const RGBAColor& brown();
+    constexpr static const RGBAColor& lime();
+    constexpr static const RGBAColor& teal();
+    constexpr static const RGBAColor& navy();
+    constexpr static const RGBAColor& olive();
+    constexpr static const RGBAColor& maroon();
+    constexpr static const RGBAColor& silver();
+    constexpr static const RGBAColor& gold();
 
     /// @brief enlight a color by a given factor.
     static RGBAColor lighten(const RGBAColor& in, const SReal factor);
@@ -226,17 +237,28 @@ constexpr RGBAColor operator/(const RGBAColor& l, const float div)
 }
 
 
-constexpr RGBAColor g_white     {1.0f,1.0f,1.0f,1.0f};
-constexpr RGBAColor g_black     {0.0f,0.0f,0.0f,1.0f};
-constexpr RGBAColor g_red       {1.0f,0.0f,0.0f,1.0f};
-constexpr RGBAColor g_green     {0.0f,1.0f,0.0f,1.0f};
-constexpr RGBAColor g_blue      {0.0f,0.0f,1.0f,1.0f};
-constexpr RGBAColor g_cyan      {0.0f,1.0f,1.0f,1.0f};
-constexpr RGBAColor g_magenta   {1.0f,0.0f,1.0f,1.0f};
-constexpr RGBAColor g_yellow    {1.0f,1.0f,0.0f,1.0f};
-constexpr RGBAColor g_gray      {0.5f,0.5f,0.5f,1.0f};
-constexpr RGBAColor g_darkgray  {0.25f,0.25f,0.25f,1.0f};
-constexpr RGBAColor g_lightgray {0.75f,0.75f,0.75f,1.0f};
+constexpr RGBAColor g_white     {1.0f, 1.0f, 1.0f, 1.0f};
+constexpr RGBAColor g_black     {0.0f, 0.0f, 0.0f, 1.0f};
+constexpr RGBAColor g_red       {1.0f, 0.0f, 0.0f, 1.0f};
+constexpr RGBAColor g_green     {0.0f, 1.0f, 0.0f, 1.0f};
+constexpr RGBAColor g_blue      {0.0f, 0.0f, 1.0f, 1.0f};
+constexpr RGBAColor g_cyan      {0.0f, 1.0f, 1.0f, 1.0f};
+constexpr RGBAColor g_magenta   {1.0f, 0.0f, 1.0f, 1.0f};
+constexpr RGBAColor g_yellow    {1.0f, 1.0f, 0.0f, 1.0f};
+constexpr RGBAColor g_gray      {0.5f, 0.5f, 0.5f, 1.0f};
+constexpr RGBAColor g_darkgray  {0.25f, 0.25f,0.25f, 1.0f};
+constexpr RGBAColor g_lightgray {0.75f, 0.75f,0.75f, 1.0f};
+constexpr RGBAColor g_orange    { 1.0f, 0.5f, 0.0f, 1.0f };
+constexpr RGBAColor g_purple    { 0.5f, 0.0f, 0.5f, 1.0f };
+constexpr RGBAColor g_pink      { 1.0f, 0.0f, 0.5f, 1.0f };
+constexpr RGBAColor g_brown     { 0.6f, 0.3f, 0.0f, 1.0f };
+constexpr RGBAColor g_lime      { 0.5f, 1.0f, 0.0f, 1.0f };
+constexpr RGBAColor g_teal      { 0.0f, 0.5f, 0.5f, 1.0f };
+constexpr RGBAColor g_navy      { 0.0f, 0.0f, 0.5f, 1.0f };
+constexpr RGBAColor g_olive     { 0.5f, 0.5f, 0.0f, 1.0f };
+constexpr RGBAColor g_maroon    { 0.5f, 0.0f, 0.0f, 1.0f };
+constexpr RGBAColor g_silver    { 0.75f, 0.75f, 0.75f, 1.0f };
+constexpr RGBAColor g_gold      { 1.0f, 0.84f, 0.0f, 1.0f };
 
 constexpr const RGBAColor& RGBAColor::white()    { return g_white;     }
 constexpr const RGBAColor& RGBAColor::black()    { return g_black;     }
@@ -249,5 +271,17 @@ constexpr const RGBAColor& RGBAColor::yellow()   { return g_yellow;    }
 constexpr const RGBAColor& RGBAColor::gray()     { return g_gray;      }
 constexpr const RGBAColor& RGBAColor::darkgray() { return g_darkgray;  }
 constexpr const RGBAColor& RGBAColor::lightgray(){ return g_lightgray; }
+constexpr const RGBAColor& RGBAColor::orange()   { return g_orange; }
+constexpr const RGBAColor& RGBAColor::purple()   { return g_purple; }
+constexpr const RGBAColor& RGBAColor::pink()     { return g_pink  ; }
+constexpr const RGBAColor& RGBAColor::brown()    { return g_brown ; }
+constexpr const RGBAColor& RGBAColor::lime()     { return g_lime  ; }
+constexpr const RGBAColor& RGBAColor::teal()     { return g_teal  ; }
+constexpr const RGBAColor& RGBAColor::navy()     { return g_navy  ; }
+constexpr const RGBAColor& RGBAColor::olive()    { return g_olive ; }
+constexpr const RGBAColor& RGBAColor::maroon()   { return g_maroon; }
+constexpr const RGBAColor& RGBAColor::silver()   { return g_silver; }
+constexpr const RGBAColor& RGBAColor::gold()     { return g_gold  ; }
+
 
 } // namespace sofa::type
diff --git a/Sofa/framework/Type/src/sofa/type/StrongType.h b/Sofa/framework/Type/src/sofa/type/StrongType.h
new file mode 100644
index 00000000000..eff9c92711f
--- /dev/null
+++ b/Sofa/framework/Type/src/sofa/type/StrongType.h
@@ -0,0 +1,302 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+
+namespace sofa::type
+{
+
+/**
+ *  A template class to create strong, type-safe aliases for existing types,
+ *  with optional extensible functionalities.
+ *
+ * The StrongType class is designed to add type safety by creating wrapper types
+ * around primitive types or existing classes. By defining unique, named types
+ * for different uses of the same underlying type, StrongType helps prevent
+ * mistakes caused by passing incorrect parameters of the same type.
+ *
+ * For example, instead of passing two `int` parameters to a function (e.g., for
+ * "width" and "height"), StrongType allows you to define distinct types for
+ * width and height:
+ *
+ * @code
+ * using Width = StrongType<int, struct WidthTag>;
+ * using Height = StrongType<int, struct HeightTag>;
+ *
+ * void resize(Width width, Height height);
+ *
+ * resize(Width{5}, Height{10});  // Correct
+ * resize(Height{10}, Width{5});  // Compilation error
+ * @endcode
+ *
+ * Strongly inspired by https://github.com/joboccara/NamedType
+ *
+ * @tparam TUnderlyingType The type to be wrapped by StrongType
+ * @tparam UniqueIdentifyingTag A unique tag type to create a distinct StrongType from other wrappers around the same `TUnderlyingType`.
+ * @tparam Functionalities A set of optional functionalities to extend StrongType
+ */
+template<class TUnderlyingType, class UniqueIdentifyingTag, template <typename> class... Functionalities>
+struct StrongType : Functionalities<StrongType<TUnderlyingType, UniqueIdentifyingTag, Functionalities...>>...
+{
+    using UnderlyingType = TUnderlyingType;
+
+    constexpr explicit StrongType(const UnderlyingType& v) : m_value(v) {}
+
+    [[nodiscard]] constexpr UnderlyingType& get() noexcept
+    {
+        return m_value;
+    }
+
+    [[nodiscard]] constexpr const std::remove_reference_t<UnderlyingType>& get() const
+    {
+        return m_value;
+    }
+
+private:
+    UnderlyingType m_value;
+};
+
+namespace functionality
+{
+
+template <typename T, template <typename> class Functionality>
+struct UnderlyingAccess
+{
+    [[nodiscard]] T& underlying()
+    {
+        return static_cast<T&>(*this);
+    }
+    [[nodiscard]] constexpr T const& underlying() const
+    {
+        return static_cast<T const&>(*this);
+    }
+};
+
+template <typename T>
+struct PreIncrementable : UnderlyingAccess<T, PreIncrementable>
+{
+    constexpr T& operator++()
+    {
+        ++static_cast<T&>(*this).get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct PostIncrementable : UnderlyingAccess<T, PostIncrementable>
+{
+    constexpr T operator++(int)
+    {
+        return T(this->underlying().get()++);
+    }
+};
+
+
+template <typename T>
+struct Incrementable
+    : PreIncrementable<T>
+    , PostIncrementable<T>
+{
+    using PostIncrementable<T>::operator++;
+    using PreIncrementable<T>::operator++;
+};
+
+template <typename T>
+struct PreDecrementable : UnderlyingAccess<T, PreDecrementable>
+{
+    constexpr T& operator--()
+    {
+        --this->underlying().get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct PostDecrementable : UnderlyingAccess<T, PostDecrementable>
+{
+    constexpr T operator--(int)
+    {
+        return T( this->underlying().get()-- );
+    }
+};
+
+
+template <typename T>
+struct Decrementable
+    : PreDecrementable<T>
+    , PostDecrementable<T>
+{
+    using PostDecrementable<T>::operator--;
+    using PreDecrementable<T>::operator--;
+};
+
+template <typename T>
+struct BinaryAddable : UnderlyingAccess<T, BinaryAddable>
+{
+    constexpr T operator+(T const& other) const
+    {
+        return T(this->underlying().get() + other.get());
+    }
+    T& operator+=(T const& other)
+    {
+        this->underlying().get() += other.get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct UnaryAddable : UnderlyingAccess<T, UnaryAddable>
+{
+    constexpr T operator+() const
+    {
+        return T(+this->underlying().get());
+    }
+};
+
+template <typename T>
+struct Addable
+    : BinaryAddable<T>
+    , UnaryAddable<T>
+{
+    using BinaryAddable<T>::operator+;
+    using UnaryAddable<T>::operator+;
+};
+
+template <typename T>
+struct BinarySubtractable : UnderlyingAccess<T, BinarySubtractable>
+{
+    [[nodiscard]] constexpr T operator-(T const& other) const
+    {
+        return T(this->underlying().get() - other.get());
+    }
+    constexpr T& operator-=(T const& other)
+    {
+        this->underlying().get() -= other.get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct UnarySubtractable : UnderlyingAccess<T, UnarySubtractable>
+{
+    [[nodiscard]] constexpr T operator-() const
+    {
+        return T(-this->underlying().get());
+    }
+};
+
+template <typename T>
+struct Subtractable
+    : BinarySubtractable<T>
+    , UnarySubtractable<T>
+{
+    using UnarySubtractable<T>::operator-;
+    using BinarySubtractable<T>::operator-;
+};
+
+
+template <typename T>
+struct Multiplicable : UnderlyingAccess<T, Multiplicable>
+{
+    [[nodiscard]] constexpr T operator*(T const& other) const
+    {
+        return T(this->underlying().get() * other.get());
+    }
+    constexpr T& operator*=(T const& other)
+    {
+        this->underlying().get() *= other.get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct Divisible : UnderlyingAccess<T, Divisible>
+{
+    [[nodiscard]] constexpr T operator/(T const& other) const
+    {
+        return T(this->underlying().get() / other.get());
+    }
+    constexpr T& operator/=(T const& other)
+    {
+        this->underlying().get() /= other.get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct Modulable : UnderlyingAccess<T, Modulable>
+{
+    [[nodiscard]] constexpr T operator%(T const& other) const
+    {
+        return T(this->underlying().get() % other.get());
+    }
+    constexpr T& operator%=(T const& other)
+    {
+        this->underlying().get() %= other.get();
+        return this->underlying();
+    }
+};
+
+template <typename T>
+struct Comparable : UnderlyingAccess<T, Comparable>
+{
+    [[nodiscard]] constexpr bool operator<(Comparable<T> const& other) const
+    {
+        return this->underlying().get() < other.underlying().get();
+    }
+    [[nodiscard]] constexpr bool operator>(Comparable<T> const& other) const
+    {
+        return other.underlying().get() < this->underlying().get();
+    }
+    [[nodiscard]] constexpr bool operator<=(Comparable<T> const& other) const
+    {
+        return !(other < *this);
+    }
+    [[nodiscard]] constexpr bool operator>=(Comparable<T> const& other) const
+    {
+        return !(*this < other);
+    }
+    [[nodiscard]] constexpr bool operator==(Comparable<T> const& other) const
+    {
+        return !(*this < other) && !(other < *this);
+    }
+    [[nodiscard]] constexpr bool operator!=(Comparable<T> const& other) const
+    {
+        return !(*this == other);
+    }
+};
+
+template <typename T>
+struct Arithmetic
+    : Incrementable<T>
+    , Decrementable<T>
+    , Addable<T>
+    , Subtractable<T>
+    , Multiplicable<T>
+    , Divisible<T>
+    , Modulable<T>
+    , Comparable<T>
+{
+};
+
+}
+
+}
diff --git a/Sofa/framework/Type/test/CMakeLists.txt b/Sofa/framework/Type/test/CMakeLists.txt
index 7211d45a336..84f70ddc9d4 100644
--- a/Sofa/framework/Type/test/CMakeLists.txt
+++ b/Sofa/framework/Type/test/CMakeLists.txt
@@ -3,11 +3,12 @@ cmake_minimum_required(VERSION 3.22)
 project(Sofa.Type_test)
 
 set(SOURCE_FILES
-    Color_test.cpp
     MatSym_test.cpp
     MatTypes_test.cpp
     Material_test.cpp
     Quater_test.cpp
+    RGBAColor_test.cpp
+    StrongType_test.cpp
     SVector_test.cpp
     VecTypes_test.cpp
     fixed_array_test.cpp
diff --git a/Sofa/framework/Type/test/Color_test.cpp b/Sofa/framework/Type/test/RGBAColor_test.cpp
similarity index 90%
rename from Sofa/framework/Type/test/Color_test.cpp
rename to Sofa/framework/Type/test/RGBAColor_test.cpp
index fb30479c4af..1f4d37dc5b1 100644
--- a/Sofa/framework/Type/test/Color_test.cpp
+++ b/Sofa/framework/Type/test/RGBAColor_test.cpp
@@ -65,6 +65,19 @@ void Color_Test::checkCreateFromString()
     EXPECT_EQ( RGBAColor::fromString("magenta"), RGBAColor(1.0,0.0,1.0,1.0) ) ;
     EXPECT_EQ( RGBAColor::fromString("yellow"), RGBAColor(1.0,1.0,0.0,1.0) ) ;
     EXPECT_EQ( RGBAColor::fromString("gray"), RGBAColor(0.5,0.5,0.5,1.0) ) ;
+    EXPECT_EQ( RGBAColor::fromString("darkgray"), RGBAColor(0.25f, 0.25f,0.25f, 1.0f ) );
+    EXPECT_EQ( RGBAColor::fromString("lightgray"), RGBAColor(0.75f, 0.75f,0.75f, 1.0f ) );
+    EXPECT_EQ( RGBAColor::fromString("orange"), RGBAColor(1.0f, 0.5f, 0.0f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("purple"), RGBAColor(0.5f, 0.0f, 0.5f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("pink"), RGBAColor(1.0f, 0.0f, 0.5f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("brown"), RGBAColor(0.6f, 0.3f, 0.0f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("lime"), RGBAColor(0.5f, 1.0f, 0.0f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("teal"), RGBAColor(0.0f, 0.5f, 0.5f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("navy"), RGBAColor(0.0f, 0.0f, 0.5f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("olive"), RGBAColor(0.5f, 0.5f, 0.0f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("maroon"), RGBAColor(0.5f, 0.0f, 0.0f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("silver"), RGBAColor(0.75f, 0.75f, 0.75f, 1.0f ));
+    EXPECT_EQ( RGBAColor::fromString("gold"), RGBAColor(1.0f, 0.84f, 0.0f, 1.0f ));
 
     RGBAColor color;
     EXPECT_TRUE( RGBAColor::read("white", color) ) ;
diff --git a/Sofa/framework/Type/test/StrongType_test.cpp b/Sofa/framework/Type/test/StrongType_test.cpp
new file mode 100644
index 00000000000..767813b0baa
--- /dev/null
+++ b/Sofa/framework/Type/test/StrongType_test.cpp
@@ -0,0 +1,83 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <gtest/gtest.h>
+#include <sofa/type/StrongType.h>
+
+namespace sofa
+{
+
+TEST(StrongType, constructor)
+{
+    using myint = sofa::type::StrongType<int, struct _myinttag>;
+
+    // myint a; //this is not possible as the default constructor is explicit
+    constexpr myint a ( 2 );
+}
+
+TEST(StrongType, preIncrementable)
+{
+    using myint = sofa::type::StrongType<int, struct _preincrementable, sofa::type::functionality::PreIncrementable>;
+    myint a ( 2 );
+    ++a;
+    EXPECT_EQ(a.get(), 3);
+}
+
+TEST(StrongType, postIncrementable)
+{
+    using myint = sofa::type::StrongType<int, struct _postincrementable, sofa::type::functionality::PostIncrementable>;
+    myint a ( 2 );
+    a++;
+    EXPECT_EQ(a.get(), 3);
+}
+
+TEST(StrongType, floatAddable)
+{
+    using myfloat = sofa::type::StrongType<float, struct _floataddable, sofa::type::functionality::Addable>;
+    static constexpr myfloat a ( 2 );
+    static constexpr myfloat b ( 3 );
+    EXPECT_EQ((a + b).get(), 5);
+}
+
+TEST(StrongType, floatBinarySubtractable)
+{
+    using myfloat = sofa::type::StrongType<float, struct _floatBinarySubstractable, sofa::type::functionality::BinarySubtractable>;
+    static constexpr myfloat a ( 2 );
+    static constexpr myfloat b ( 3 );
+    EXPECT_EQ((a - b).get(), -1);
+}
+
+TEST(StrongType, floatUnarySubtractable)
+{
+    using myfloat = sofa::type::StrongType<float, struct _floatUnarySubstractable, sofa::type::functionality::UnarySubtractable>;
+    static constexpr myfloat a ( 2 );
+    EXPECT_EQ((-a).get(), -2);
+}
+
+TEST(StrongType, floatMultiplicable)
+{
+    using myfloat = sofa::type::StrongType<float, struct _floatMultiplicable, sofa::type::functionality::Multiplicable>;
+    static constexpr myfloat a ( 2 );
+    static constexpr myfloat b ( 3 );
+    EXPECT_EQ((a * b).get(), 6);
+}
+
+}
diff --git a/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyBVHController.cpp b/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyBVHController.cpp
index a7b0672fd02..f735d6b58cb 100644
--- a/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyBVHController.cpp
+++ b/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyBVHController.cpp
@@ -166,8 +166,6 @@ void ArticulatedHierarchyBVHController::applyController(void)
     }
 }
 
-SOFA_DECL_CLASS(ArticulatedHierarchyBVHController)
-
 // Register in the Factory
 int ArticulatedHierarchyBVHControllerClass = core::RegisterObject("Implements a handler that controls the values of the articulations of an articulated hierarchy container using a .bvh file.")
         .add< ArticulatedHierarchyBVHController >()
diff --git a/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyController.cpp b/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyController.cpp
index 2542e59290e..89f8830e987 100644
--- a/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyController.cpp
+++ b/applications/plugins/ArticulatedSystemPlugin/src/ArticulatedSystemPlugin/ArticulatedHierarchyController.cpp
@@ -407,8 +407,6 @@ void ArticulatedHierarchyController::applyController(void)
     }
 }
 
-SOFA_DECL_CLASS(ArticulatedHierarchyController)
-
 // Register in the Factory
 int ArticulatedHierarchyControllerClass = core::RegisterObject("Implements an user interaction handler that controls the values of the articulations of an articulated hierarchy container.")
         .add< ArticulatedHierarchyController >()
diff --git a/applications/plugins/SofaCUDA/examples/raptor-cpu.scn b/applications/plugins/SofaCUDA/examples/raptor-cpu.scn
index 13de694e82d..a29d72c86d4 100644
--- a/applications/plugins/SofaCUDA/examples/raptor-cpu.scn
+++ b/applications/plugins/SofaCUDA/examples/raptor-cpu.scn
@@ -6,7 +6,7 @@
 <!--    <MeshOBJLoader name="surface" filename="mesh/raptor8k1.obj" />-->
 		<EulerImplicitSolver rayleighMass="0.01" rayleighStiffness="0.01" />
 		<CGLinearSolver iterations="25" tolerance="1e-6" threshold="1e-20"/>
-<!--	<ShewchukPCGLinearSolver iterations="25" tolerance="1e-3" preconditioner="@precond" />
+<!--	<PCGLinearSolver iterations="25" tolerance="1e-3" preconditioner="@precond" />
 		<JacobiPreconditioner template="CudaDiagonalMatrixf" name="precond" /> -->
 		<MeshTopology src="@volume" />
 		<MechanicalObject template="Vec3" />
diff --git a/applications/plugins/SofaCUDA/examples/raptor-cuda.scn b/applications/plugins/SofaCUDA/examples/raptor-cuda.scn
index 4b3c414284d..962f7f57077 100644
--- a/applications/plugins/SofaCUDA/examples/raptor-cuda.scn
+++ b/applications/plugins/SofaCUDA/examples/raptor-cuda.scn
@@ -13,7 +13,7 @@
 	<!--    <MeshOBJLoader name="surface" filename="mesh/raptor8k1.obj" />-->
 		<EulerImplicitSolver rayleighStiffness="0.1" rayleighMass="0.1" />
 		<CGLinearSolver iterations="25" tolerance="1e-6" threshold="1e-20"/>
-        <!--	<ShewchukPCGLinearSolver iterations="25" tolerance="1e-3" preconditioner="@precond" />
+        <!--	<PCGLinearSolver iterations="25" tolerance="1e-3" preconditioner="@precond" />
 		<JacobiPreconditioner template="CudaDiagonalMatrixf" name="precond" /> -->
 		<MeshTopology src="@volume" />
 		<MechanicalObject template="CudaVec3f" />
diff --git a/applications/plugins/SofaCUDA/examples/raptor.scn b/applications/plugins/SofaCUDA/examples/raptor.scn
index 0a6c0019fe9..69b658ffb1a 100644
--- a/applications/plugins/SofaCUDA/examples/raptor.scn
+++ b/applications/plugins/SofaCUDA/examples/raptor.scn
@@ -1,7 +1,7 @@
 <Node name="root" dt="0.04" showBehaviorModels="0" showCollisionModels="0" showMappings="0" showForceFields="1">
     <RequiredPlugin name="Sofa.Component.Engine.Transform"/> <!-- Needed to use components [IndexValueMapper] -->
     <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <!-- Needed to use components [MeshVTKLoader] -->
-    <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [ShewchukPCGLinearSolver] -->
+    <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [PCGLinearSolver] -->
     <RequiredPlugin name="Sofa.Component.LinearSolver.Preconditioner"/> <!-- Needed to use components [JacobiPreconditioner] -->
     <RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [IdentityMapping] -->
     <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
@@ -16,7 +16,7 @@
 <!--    <MeshOBJLoader name="surface" filename="mesh/raptor8k1.obj" />-->
     <EulerImplicitSolver rayleighStiffness="0.1" rayleighMass="0.1" />
 <!--    <CGLinearSolver iterations="30" tolerance="1e-10" threshold="1e-10"/>-->
-        <ShewchukPCGLinearSolver iterations="30" tolerance="1e-10" preconditioner="@precond" />
+        <PCGLinearSolver iterations="30" tolerance="1e-10" preconditioner="@precond" />
 	<JacobiPreconditioner template="CudaDiagonalMatrixf" name="precond" />
     <MeshTopology src="@volume" />
 	<MechanicalObject template="CudaVec3f" />
diff --git a/cmake/Modules/FindQGLViewer.cmake b/cmake/Modules/FindQGLViewer.cmake
index ae6fe0d2091..65eb618dfff 100644
--- a/cmake/Modules/FindQGLViewer.cmake
+++ b/cmake/Modules/FindQGLViewer.cmake
@@ -14,14 +14,14 @@ if(NOT TARGET QGLViewer)
 
   if(NOT QGLViewer_INCLUDE_DIR)
     find_path(QGLViewer_INCLUDE_DIR
-      NAMES qglviewer.h
-      PATH_SUFFIXES include/QGLViewer
+      NAMES QGLViewer/qglviewer.h
+      PATH_SUFFIXES include Headers
     )
   endif()
 
   if(NOT QGLViewer_LIBRARY)
   find_library(QGLViewer_LIBRARY
-    NAMES QGLViewer QGLViewer2 QGLViewer-qt5
+    NAMES QGLViewer QGLViewer2 QGLViewer-qt5 QGLViewer-qt6
     PATH_SUFFIXES lib
   )
   endif()
diff --git a/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn b/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn
index c11b2de4255..89b5871bb32 100644
--- a/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn
+++ b/examples/Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver iterations="1000" />
+    <PCGLinearSolver iterations="1000" />
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
 
 </Node>
\ No newline at end of file
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_AsyncSparseLDLSolver.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_AsyncSparseLDLSolver.scn
index 588a6fc78ed..082a192cb47 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_AsyncSparseLDLSolver.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_AsyncSparseLDLSolver.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <AsyncSparseLDLSolver name="preconditioner"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.45" method="large" />
 
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_BlockJacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_BlockJacobiPreconditioner.scn
index e657eb38a31..797e956cb24 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_BlockJacobiPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_BlockJacobiPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <BlockJacobiPreconditioner name="preconditioner"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
 
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn
index 2940aa82c53..b8a29eb500e 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_JacobiPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <JacobiPreconditioner name="preconditioner"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
 
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_NoPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_NoPreconditioner.scn
index dacb8fed69d..91a74ddd909 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_NoPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_NoPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" use_precond="false"/>
+    <PCGLinearSolver name="PCG" iterations="1000" use_precond="false"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
 
 </Node>
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn
index de1f8a7a434..fcb6eea9c34 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_PrecomputedWarpPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <PrecomputedWarpPreconditioner name="preconditioner" linearSolver="@initSolver" printLog="true"/>
     <SparseLDLSolver name="initSolver" template="CompressedRowSparseMatrix" />
     <TetrahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_SSORPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_SSORPreconditioner.scn
index 295669b9709..7f0307c18f3 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_SSORPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_SSORPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <SSORPreconditioner name="preconditioner"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
 
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn
index 24134226f6b..aa4c61c823e 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpPreconditioner.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner"/>
     <WarpPreconditioner name="preconditioner" linearSolver="@initSolver" printLog="true"/>
     <SparseLDLSolver name="initSolver" template="CompressedRowSparseMatrixMat3x3" />
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
diff --git a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn
index fb72a11ea40..5198f1f8e9e 100644
--- a/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn
+++ b/examples/Component/LinearSolver/Preconditioner/FEMBAR_PCG_WarpedAsyncSparseLDLSolver.scn
@@ -2,7 +2,7 @@
 
     <include href="../FEMBAR-common.xml"/>
 
-    <ShewchukPCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner" update_step="1"/>
+    <PCGLinearSolver name="PCG" iterations="1000" preconditioner="@preconditioner" update_step="1"/>
     <WarpPreconditioner name="preconditioner" linearSolver="@factorization" update_step="15"/>
     <AsyncSparseLDLSolver name="factorization" template="CompressedRowSparseMatrixMat3x3"/>
     <HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />
diff --git a/examples/Component/Mapping/NonLinear/AreaMapping.scn b/examples/Component/Mapping/NonLinear/AreaMapping.scn
index 15d707aae2c..1591dceb42e 100644
--- a/examples/Component/Mapping/NonLinear/AreaMapping.scn
+++ b/examples/Component/Mapping/NonLinear/AreaMapping.scn
@@ -48,7 +48,9 @@
         </Node>
 
         <Node name="Visu">
-            <OglModel name="Visual" color="green" />
+            <VisualStyle displayFlags="hideWireframe"/>
+            <DataDisplay name="Visual" triangleData="@../constraintSpace/areaDoFs.position"/>
+            <OglColorMap colorScheme="HSV" showLegend="true" legendTitle="Triangle area" min="@Visual.currentMin" max="@Visual.currentMax"/>
             <IdentityMapping input="@.." output="@Visual" />
         </Node>
     </Node>
@@ -73,7 +75,7 @@
         </Node>
 
         <Node name="Visu">
-            <OglModel name="Visual" color="red" />
+            <OglModel name="Visual" color="darkgray" />
             <IdentityMapping input="@.." output="@Visual" />
         </Node>
     </Node>
diff --git a/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn b/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn
index 0fe9a314781..557186571f7 100644
--- a/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn
+++ b/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn
@@ -8,7 +8,7 @@
     <RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
     <RequiredPlugin name="Sofa.Component.Engine.Transform"/> <!-- Needed to use components [IndexValueMapper] -->
     <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
-    <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [ShewchukPCGLinearSolver] -->
+    <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [PCGLinearSolver] -->
     <RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [BarycentricMapping IdentityMapping] -->
     <RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
     <RequiredPlugin name="Sofa.Component.MechanicalLoad"/> <!-- Needed to use components [ConstantForceField] -->
@@ -39,7 +39,7 @@
 
     <Node name="FEM">
         <EulerImplicitSolver rayleighMass="0.1" rayleighStiffness="0.1" _velocityDamping="0.5"/>
-        <ShewchukPCGLinearSolver preconditioner="@precond" update_step="15"/>
+        <PCGLinearSolver preconditioner="@precond" update_step="15"/>
         <SparseLDLSolver name="precond" template="CompressedRowSparseMatrixMat3x3"/>
 
         <TetrahedronSetTopologyContainer name="Container" position="@../Topology/HexaTop.position" tetrahedra="@../Topology/Container.tetrahedra"/>
diff --git a/examples/Component/SolidMechanics/Spring/TriangularBendingSprings.scn b/examples/Component/SolidMechanics/Spring/TriangularBendingSprings.scn
index 8a54f619b88..06cace33be4 100644
--- a/examples/Component/SolidMechanics/Spring/TriangularBendingSprings.scn
+++ b/examples/Component/SolidMechanics/Spring/TriangularBendingSprings.scn
@@ -4,7 +4,7 @@
         <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
         <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <!-- Needed to use components [MeshGmshLoader] -->
         <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [AsyncSparseLDLSolver] -->
-        <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [ShewchukPCGLinearSolver] -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [PCGLinearSolver] -->
         <RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [IdentityMapping] -->
         <RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [DiagonalMass] -->
         <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
@@ -20,7 +20,7 @@
 
     <Node name="SquareGravity">
         <EulerImplicitSolver name="cg_odesolver" printLog="false"  rayleighStiffness="0.1" rayleighMass="0.1" />
-        <ShewchukPCGLinearSolver preconditioner="@preconditioner"/>
+        <PCGLinearSolver preconditioner="@preconditioner"/>
         <AsyncSparseLDLSolver name="preconditioner" template="CompressedRowSparseMatrixMat3x3"/>
         <MeshGmshLoader name="loader" filename="mesh/square3.msh" createSubelements="true"/>
         <MechanicalObject src="@loader" scale="10" />
diff --git a/examples/Demos/fallingBeamLagrangianCollision.scn b/examples/Demos/fallingBeamLagrangianCollision.scn
index 9910730d396..3706a2e8791 100644
--- a/examples/Demos/fallingBeamLagrangianCollision.scn
+++ b/examples/Demos/fallingBeamLagrangianCollision.scn
@@ -10,7 +10,7 @@
         <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Correction"/> <!-- Needed to use components [LinearSolverConstraintCorrection] -->
         <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Solver"/> <!-- Needed to use components [GenericConstraintSolver] -->
         <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
-        <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [ShewchukPCGLinearSolver] -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [PCGLinearSolver] -->
         <RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [IdentityMapping] -->
         <RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
         <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
@@ -44,7 +44,7 @@
 
     <Node name="FEM">
         <EulerImplicitSolver firstOrder="false" rayleighMass="0.1" rayleighStiffness="0.1"/>
-        <ShewchukPCGLinearSolver preconditioner="@precond"/>
+        <PCGLinearSolver preconditioner="@precond"/>
         <SparseLDLSolver name="precond" template="CompressedRowSparseMatrixMat3x3" parallelInverseProduct="true" />
 
         <TetrahedronSetTopologyContainer name="Container" position="@../BEAMVOLUME/HexaTop.position"
diff --git a/flake.lock b/flake.lock
new file mode 100644
index 00000000000..26492ec359a
--- /dev/null
+++ b/flake.lock
@@ -0,0 +1,108 @@
+{
+  "nodes": {
+    "flake-parts": {
+      "inputs": {
+        "nixpkgs-lib": "nixpkgs-lib"
+      },
+      "locked": {
+        "lastModified": 1727826117,
+        "narHash": "sha256-K5ZLCyfO/Zj9mPFldf3iwS6oZStJcU4tSpiXTMYaaL0=",
+        "owner": "hercules-ci",
+        "repo": "flake-parts",
+        "rev": "3d04084d54bedc3d6b8b736c70ef449225c361b1",
+        "type": "github"
+      },
+      "original": {
+        "owner": "hercules-ci",
+        "repo": "flake-parts",
+        "type": "github"
+      }
+    },
+    "flake-utils": {
+      "locked": {
+        "lastModified": 1659877975,
+        "narHash": "sha256-zllb8aq3YO3h8B/U0/J1WBgAL8EX5yWf5pMj3G0NAmc=",
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "rev": "c0e246b9b83f637f4681389ecabcb2681b4f3af0",
+        "type": "github"
+      },
+      "original": {
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "type": "github"
+      }
+    },
+    "nixgl": {
+      "inputs": {
+        "flake-utils": "flake-utils",
+        "nixpkgs": "nixpkgs"
+      },
+      "locked": {
+        "lastModified": 1713543440,
+        "narHash": "sha256-lnzZQYG0+EXl/6NkGpyIz+FEOc/DSEG57AP1VsdeNrM=",
+        "owner": "nix-community",
+        "repo": "nixGL",
+        "rev": "310f8e49a149e4c9ea52f1adf70cdc768ec53f8a",
+        "type": "github"
+      },
+      "original": {
+        "owner": "nix-community",
+        "repo": "nixGL",
+        "type": "github"
+      }
+    },
+    "nixpkgs": {
+      "locked": {
+        "lastModified": 1660551188,
+        "narHash": "sha256-a1LARMMYQ8DPx1BgoI/UN4bXe12hhZkCNqdxNi6uS0g=",
+        "owner": "nixos",
+        "repo": "nixpkgs",
+        "rev": "441dc5d512153039f19ef198e662e4f3dbb9fd65",
+        "type": "github"
+      },
+      "original": {
+        "owner": "nixos",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-lib": {
+      "locked": {
+        "lastModified": 1727825735,
+        "narHash": "sha256-0xHYkMkeLVQAMa7gvkddbPqpxph+hDzdu1XdGPJR+Os=",
+        "type": "tarball",
+        "url": "https://github.com/NixOS/nixpkgs/archive/fb192fec7cc7a4c26d51779e9bab07ce6fa5597a.tar.gz"
+      },
+      "original": {
+        "type": "tarball",
+        "url": "https://github.com/NixOS/nixpkgs/archive/fb192fec7cc7a4c26d51779e9bab07ce6fa5597a.tar.gz"
+      }
+    },
+    "nixpkgs_2": {
+      "locked": {
+        "lastModified": 1728930054,
+        "narHash": "sha256-mCaSyViQyiLgZKVLpDcVacbPSNjvsBzfWnGaxvhTzs8=",
+        "owner": "nim65s",
+        "repo": "nixpkgs",
+        "rev": "8e906b3e2aa2274c9d0c555393b4801d3c0badee",
+        "type": "github"
+      },
+      "original": {
+        "owner": "nim65s",
+        "ref": "qt6-libqglviewer",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "root": {
+      "inputs": {
+        "flake-parts": "flake-parts",
+        "nixgl": "nixgl",
+        "nixpkgs": "nixpkgs_2"
+      }
+    }
+  },
+  "root": "root",
+  "version": 7
+}
diff --git a/flake.nix b/flake.nix
new file mode 100644
index 00000000000..054afca0134
--- /dev/null
+++ b/flake.nix
@@ -0,0 +1,37 @@
+{
+  description = "SOFA is an open-source framework for interactive physics simulation, with emphasis on biomechanical and robotic simulations";
+
+  inputs = {
+    flake-parts.url = "github:hercules-ci/flake-parts";
+    #nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable";
+    # ref. https://github.com/NixOS/nixpkgs/pull/348549
+    nixpkgs.url = "github:nim65s/nixpkgs/qt6-libqglviewer";
+    nixgl.url = "github:nix-community/nixGL";
+  };
+
+  outputs =
+    inputs@{ flake-parts, nixgl, nixpkgs, ... }:
+    flake-parts.lib.mkFlake { inherit inputs; } {
+      systems = nixpkgs.lib.systems.flakeExposed;
+      perSystem =
+        { pkgs, self', system, ... }:
+        {
+          _module.args.pkgs = import nixpkgs {
+            inherit system;
+            overlays = [ nixgl.overlay ];
+          };
+          apps.nixgl = {
+            type = "app";
+            program = pkgs.writeShellApplication {
+              name = "nixgl-sofa";
+              text = "${pkgs.lib.getExe pkgs.nixgl.auto.nixGLDefault} ${pkgs.lib.getExe self'.packages.sofa}";
+            };
+          };
+          devShells.default = pkgs.mkShell { inputsFrom = [ self'.packages.default ]; };
+          packages = {
+            default = self'.packages.sofa;
+            sofa = pkgs.callPackage ./package.nix { };
+          };
+        };
+    };
+}
diff --git a/package.nix b/package.nix
new file mode 100644
index 00000000000..1b879c51d14
--- /dev/null
+++ b/package.nix
@@ -0,0 +1,82 @@
+{
+  boost,
+  cmake,
+  cxxopts,
+  eigen,
+  #fetchFromGitHub,
+  glew,
+  gtest,
+  lib,
+  qt6Packages,
+  libGL,
+  metis,
+  stdenv,
+  tinyxml-2,
+  zlib,
+}:
+
+stdenv.mkDerivation {
+  pname = "sofa";
+  version = "24.06.00";
+
+  src = lib.fileset.toSource {
+    root = ./.;
+    fileset = lib.fileset.unions [
+      ./applications
+      ./Authors.txt
+      ./cmake
+      ./CHANGELOG.md
+      ./CMakeLists.txt
+      ./CMakePresets.json
+      ./examples
+      ./extlibs
+      ./LICENSE-LGPL.md
+      ./package.cmake
+      ./README.md
+      ./scripts
+      ./share
+      ./Sofa
+      ./tools
+    ];
+  };
+
+  propagatedNativeBuildInputs = [
+    cmake
+    qt6Packages.wrapQtAppsHook
+  ];
+  propagatedBuildInputs = [
+    boost
+    cxxopts
+    eigen
+    glew
+    gtest
+    qt6Packages.libqglviewer
+    qt6Packages.qtbase
+    libGL
+    metis
+    tinyxml-2
+    zlib
+  ];
+
+  cmakeFlags = [
+    (lib.cmakeBool "SOFA_ALLOW_FETCH_DEPENDENCIES" false)
+  ];
+
+  doCheck = true;
+
+  postFixup = lib.optionalString stdenv.hostPlatform.isDarwin ''
+    install_name_tool -change \
+      $out/lib/libSceneChecking.24.12.99.dylib \
+      $out/plugins/SceneChecking/lib/libSceneChecking.24.12.99.dylib \
+      $out/bin/.runSofa-24.12.99-wrapped
+  '';
+
+  meta = {
+    description = "SOFA is an open-source framework for interactive physics simulation, with emphasis on biomechanical and robotic simulations";
+    homepage = "https://github.com/sofa-framework/sofa";
+    license = lib.licenses.lgpl21Only;
+    maintainers = with lib.maintainers; [ nim65s ];
+    mainProgram = "runSofa";
+    platforms = lib.platforms.unix ++ lib.platforms.windows;
+  };
+}