[FEM] Components override the template API instead of the generic one (…

…#4982)

* Components override the template API instead of the generic one

* Fix compilation

* Fix masses

---------

Co-authored-by: erik pernod <[email protected]>

Sofa/Component/Mass/src/sofa/component/mass/DiagonalMass.h

@@ -310,7 +310,7 @@ class DiagonalMass : public core::behavior::Mass<DataTypes>
     void addGravityToV(const core::MechanicalParams* mparams, DataVecDeriv& d_v) override;
 
     /// Add Mass contribution to global Matrix assembling
-    void addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
+    void addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset) override;
     void buildMassMatrix(sofa::core::behavior::MassMatrixAccumulator* matrices) override;
     void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override {}
     void buildDampingMatrix(core::behavior::DampingMatrix* /* matrices */) override {}

Sofa/Component/Mass/src/sofa/component/mass/DiagonalMass.inl

@@ -616,16 +616,14 @@ DiagonalMass<DataTypes, GeometricalTypes>::getMomentum ( const core::MechanicalP
 }
 
 template <class DataTypes, class GeometricalTypes>
-void DiagonalMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
+void DiagonalMass<DataTypes, GeometricalTypes>::addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset)
 {
     const MassVector &masses= d_vertexMass.getValue();
     static constexpr auto N = Deriv::total_size;
     AddMToMatrixFunctor<Deriv,MassType, linearalgebra::BaseMatrix> calc;
-    sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
-    const Real mFactor = Real(sofa::core::mechanicalparams::mFactorIncludingRayleighDamping(mparams, this->rayleighMass.getValue()));
     for (unsigned int i = 0; i < masses.size(); i++)
     {
-        calc(r.matrix, masses[i], r.offset + N * i, mFactor);
+        calc(mat, masses[i], offset + N * i, mFact);
     }
 }
 

Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.h

@@ -219,7 +219,7 @@ class MeshMatrixMass : public core::behavior::Mass<DataTypes>
 
 
     /// Add Mass contribution to global Matrix assembling
-    void addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
+    void addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset) override;
     void buildMassMatrix(sofa::core::behavior::MassMatrixAccumulator* matrices) override;
     void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override {}
     void buildDampingMatrix(core::behavior::DampingMatrix* /* matrices */) override {}

Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl

@@ -2190,7 +2190,7 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::addGravityToV(const core::Mech
 
 
 template <class DataTypes, class GeometricalTypes>
-void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
+void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset)
 {
     const auto &vertexMass= d_vertexMass.getValue();
     const auto &edgeMass= d_edgeMass.getValue();
@@ -2200,9 +2200,6 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::Mecha
 
     static constexpr auto N = Deriv::total_size;
     AddMToMatrixFunctor<Deriv,MassType, sofa::linearalgebra::BaseMatrix> calc;
-    sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
-    sofa::linearalgebra::BaseMatrix* mat = r.matrix;
-    const Real mFactor = Real(sofa::core::mechanicalparams::mFactorIncludingRayleighDamping(mparams, this->rayleighMass.getValue()));
 
     if((mat->colSize()) != (linearalgebra::BaseMatrix::Index)(l_topology->getNbPoints()*N) || (mat->rowSize()) != (linearalgebra::BaseMatrix::Index)(l_topology->getNbPoints()*N))
     {
@@ -2217,7 +2214,7 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::Mecha
         unsigned int i {};
         for (const auto& v : vertexMass)
         {
-            calc(r.matrix, v * m_massLumpingCoeff, r.offset + N * i, mFactor);
+            calc(mat, v * m_massLumpingCoeff, offset + N * i, mFact);
             massTotal += v * m_massLumpingCoeff;
             ++i;
         }
@@ -2239,7 +2236,7 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::Mecha
         unsigned int i {};
         for (const auto& v : vertexMass)
         {
-            calc(r.matrix, v, r.offset + N * i, mFactor);
+            calc(mat, v, offset + N * i, mFact);
             massTotal += v;
             ++i;
         }
@@ -2250,8 +2247,8 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::addMToMatrix(const core::Mecha
             v0 = edges[j][0];
             v1 = edges[j][1];
 
-            calc(r.matrix, edgeMass[j], r.offset + N*v0, r.offset + N*v1, mFactor);
-            calc(r.matrix, edgeMass[j], r.offset + N*v1, r.offset + N*v0, mFactor);
+            calc(mat, edgeMass[j], offset + N*v0, offset + N*v1, mFact);
+            calc(mat, edgeMass[j], offset + N*v1, offset + N*v0, mFact);
 
             massTotal += 2 * edgeMass[j];
         }

Sofa/Component/Mass/src/sofa/component/mass/UniformMass.h

@@ -142,7 +142,7 @@ class UniformMass : public core::behavior::Mass<DataTypes>
 
     void addGravityToV(const core::MechanicalParams* mparams, DataVecDeriv& d_v) override;
 
-    void addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override; /// Add Mass contribution to global Matrix assembling
+    void addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset) override; /// Add Mass contribution to global Matrix assembling
     void buildMassMatrix(sofa::core::behavior::MassMatrixAccumulator* matrices) override;
     void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override {}
     void buildDampingMatrix(core::behavior::DampingMatrix* /* matrices */) override {}

Sofa/Component/Mass/src/sofa/component/mass/UniformMass.inl

@@ -555,22 +555,18 @@ UniformMass<DataTypes>::getMomentum ( const core::MechanicalParams* params,
 
 /// Add Mass contribution to global Matrix assembling
 template <class DataTypes>
-void UniformMass<DataTypes>::addMToMatrix (const MechanicalParams *mparams,
-                                                     const MultiMatrixAccessor* matrix)
+void UniformMass<DataTypes>::addMToMatrix (sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset)
 {
     const MassType& m = d_vertexMass.getValue();
 
     static constexpr auto N = Deriv::total_size;
 
     AddMToMatrixFunctor<Deriv,MassType, linearalgebra::BaseMatrix> calc;
-    MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(mstate);
-
-    const Real mFactor = Real(sofa::core::mechanicalparams::mFactorIncludingRayleighDamping(mparams, this->rayleighMass.getValue()));
 
     const ReadAccessor<Data<SetIndexArray > > indices = d_indices;
     for (auto id : *indices)
     {
-        calc ( r.matrix, m, int(r.offset + N * id), mFactor);
+        calc ( mat, m, int(offset + N * id), mFact);
     }
 }
 

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BeamFEMForceField.h

@@ -172,7 +172,7 @@ class BeamFEMForceField : public BaseLinearElasticityFEMForceField<DataTypes>
     virtual void reinitBeam(Index i);
     void addForce(const MechanicalParams* mparams, DataVecDeriv &  dataF, const DataVecCoord &  dataX , const DataVecDeriv & dataV ) override;
     void addDForce(const MechanicalParams* mparams, DataVecDeriv&   datadF , const DataVecDeriv&   datadX ) override;
-    void addKToMatrix(const MechanicalParams* mparams, const MultiMatrixAccessor* matrix ) override;
+    void addKToMatrix(sofa::linearalgebra::BaseMatrix * matrix, SReal kFact, unsigned int &offset) override;
     void buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix) override;
     void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final;
     SReal getPotentialEnergy(const MechanicalParams* mparams, const DataVecCoord&  x) const override;

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BeamFEMForceField.inl

@@ -480,115 +480,104 @@ void BeamFEMForceField<DataTypes>::applyStiffnessLarge(VecDeriv& df, const VecDe
 }
 
 template<class DataTypes>
-void BeamFEMForceField<DataTypes>::addKToMatrix(const sofa::core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix )
+void BeamFEMForceField<DataTypes>::addKToMatrix(sofa::linearalgebra::BaseMatrix * matrix, SReal kFact, unsigned int &offset)
 {
-    sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
-    Real k = (Real)sofa::core::mechanicalparams::kFactorIncludingRayleighDamping(mparams, this->rayleighStiffness.getValue());
-    linearalgebra::BaseMatrix* mat = r.matrix;
-
     if (!m_indexedElements)
         return;
 
-    if (r)
+    if (m_partialListSegment)
     {
-        const unsigned int &offset = r.offset;
 
-        if (m_partialListSegment)
+        for (unsigned int i : d_listSegment.getValue())
         {
+            const auto& [a, b] = (*m_indexedElements)[i].array();
 
-            for (unsigned int i : d_listSegment.getValue())
-            {
-                const auto& [a, b] = (*m_indexedElements)[i].array();
-
-                type::Quat<SReal>& q = beamQuat(i);
-                q.normalize();
-                Transformation R,Rt;
-                q.toMatrix(R);
-                Rt.transpose(R);
-                const StiffnessMatrix& K0 = d_beamsData.getValue()[i]._k_loc;
-                StiffnessMatrix K;
-                for (int x1=0; x1<12; x1+=3)
-                    for (int y1=0; y1<12; y1+=3)
+            type::Quat<SReal>& q = beamQuat(i);
+            q.normalize();
+            Transformation R,Rt;
+            q.toMatrix(R);
+            Rt.transpose(R);
+            const StiffnessMatrix& K0 = d_beamsData.getValue()[i]._k_loc;
+            StiffnessMatrix K;
+            for (int x1=0; x1<12; x1+=3)
+                for (int y1=0; y1<12; y1+=3)
+                {
+                    type::Mat<3,3,Real> m;
+                    K0.getsub(x1,y1, m);
+                    m = R*m*Rt;
+                    K.setsub(x1,y1, m);
+                }
+            int index[12];
+            for (int x1=0; x1<6; x1++)
+                index[x1] = offset+a*6+x1;
+            for (int x1=0; x1<6; x1++)
+                index[6+x1] = offset+b*6+x1;
+            for (int x1=0; x1<12; ++x1)
+                for (int y1=0; y1<12; ++y1)
+                    matrix->add(index[x1], index[y1], - K(x1,y1)*kFact);
+
+        }
+
+    }
+    else
+    {
+        unsigned int i {};
+        for(auto it = m_indexedElements->begin() ; it != m_indexedElements->end() ; ++it, ++i)
+        {
+            const auto& [a, b] = it->array();
+
+            type::Quat<SReal>& q = beamQuat(i);
+            q.normalize();
+            Transformation R,Rt;
+            q.toMatrix(R);
+            Rt.transpose(R);
+            const StiffnessMatrix& K0 = d_beamsData.getValue()[i]._k_loc;
+            StiffnessMatrix K;
+            const bool exploitSymmetry = d_useSymmetricAssembly.getValue();
+
+            if (exploitSymmetry) {
+                for (int x1=0; x1<12; x1+=3) {
+                    for (int y1=x1; y1<12; y1+=3)
                     {
                         type::Mat<3,3,Real> m;
                         K0.getsub(x1,y1, m);
                         m = R*m*Rt;
-                        K.setsub(x1,y1, m);
-                    }
-                int index[12];
-                for (int x1=0; x1<6; x1++)
-                    index[x1] = offset+a*6+x1;
-                for (int x1=0; x1<6; x1++)
-                    index[6+x1] = offset+b*6+x1;
-                for (int x1=0; x1<12; ++x1)
-                    for (int y1=0; y1<12; ++y1)
-                        mat->add(index[x1], index[y1], - K(x1,y1)*k);
-
-            }
-
-        }
-        else
-        {
-            unsigned int i {};
-            for(auto it = m_indexedElements->begin() ; it != m_indexedElements->end() ; ++it, ++i)
-            {
-                const auto& [a, b] = it->array();
-
-                type::Quat<SReal>& q = beamQuat(i);
-                q.normalize();
-                Transformation R,Rt;
-                q.toMatrix(R);
-                Rt.transpose(R);
-                const StiffnessMatrix& K0 = d_beamsData.getValue()[i]._k_loc;
-                StiffnessMatrix K;
-                const bool exploitSymmetry = d_useSymmetricAssembly.getValue();
-
-                if (exploitSymmetry) {
-                    for (int x1=0; x1<12; x1+=3) {
-                        for (int y1=x1; y1<12; y1+=3)
-                        {
-                            type::Mat<3,3,Real> m;
-                            K0.getsub(x1,y1, m);
-                            m = R*m*Rt;
 
+                        for (int i=0; i<3; i++)
+                            for (int j=0; j<3; j++) {
+                                K.elems[i+x1][j+y1] += m[i][j];
+                                K.elems[j+y1][i+x1] += m[i][j];
+                            }
+                        if (x1 == y1)
                             for (int i=0; i<3; i++)
-                                for (int j=0; j<3; j++) {
-                                    K.elems[i+x1][j+y1] += m[i][j];
-                                    K.elems[j+y1][i+x1] += m[i][j];
-                                }
-                            if (x1 == y1)
-                                for (int i=0; i<3; i++)
-                                    for (int j=0; j<3; j++)
-                                        K.elems[i+x1][j+y1] *= SReal(0.5);
-
-                        }
+                                for (int j=0; j<3; j++)
+                                    K.elems[i+x1][j+y1] *= SReal(0.5);
+
                     }
-                } else  {
-                    for (int x1=0; x1<12; x1+=3) {
-                        for (int y1=0; y1<12; y1+=3)
-                        {
-                            type::Mat<3,3,Real> m;
-                            K0.getsub(x1,y1, m);
-                            m = R*m*Rt;
-                            K.setsub(x1,y1, m);
-                        }
+                }
+            } else  {
+                for (int x1=0; x1<12; x1+=3) {
+                    for (int y1=0; y1<12; y1+=3)
+                    {
+                        type::Mat<3,3,Real> m;
+                        K0.getsub(x1,y1, m);
+                        m = R*m*Rt;
+                        K.setsub(x1,y1, m);
                     }
                 }
+            }
 
-                int index[12];
-                for (int x1=0; x1<6; x1++)
-                    index[x1] = offset+a*6+x1;
-                for (int x1=0; x1<6; x1++)
-                    index[6+x1] = offset+b*6+x1;
-                for (int x1=0; x1<12; ++x1)
-                    for (int y1=0; y1<12; ++y1)
-                        mat->add(index[x1], index[y1], - K(x1,y1)*k);
+            int index[12];
+            for (int x1=0; x1<6; x1++)
+                index[x1] = offset+a*6+x1;
+            for (int x1=0; x1<6; x1++)
+                index[6+x1] = offset+b*6+x1;
+            for (int x1=0; x1<12; ++x1)
+                for (int y1=0; y1<12; ++y1)
+                    matrix->add(index[x1], index[y1], - K(x1,y1)*kFact);
 
-            }
         }
-
     }
-
 }
 
 template <class DataTypes>

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedralFEMForceField.h

@@ -154,7 +154,7 @@ class HexahedralFEMForceField : virtual public BaseLinearElasticityFEMForceField
         return 0.0;
     }
 
-    void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
+    void addKToMatrix(sofa::linearalgebra::BaseMatrix * matrix, SReal kFact, unsigned int &offset) override;
     void buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix) override;
 
     void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final;

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedralFEMForceField.inl

@@ -584,15 +584,13 @@ void HexahedralFEMForceField<DataTypes>::accumulateForcePolar(WDataRefVecDeriv&
 /////////////////////////////////////////////////
 
 template<class DataTypes>
-void HexahedralFEMForceField<DataTypes>::addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
+void HexahedralFEMForceField<DataTypes>::addKToMatrix(sofa::linearalgebra::BaseMatrix * matrix, SReal kFact, unsigned int &offset)
 {
     // Build Matrix Block for this ForceField
     int i,j,n1, n2;
 
     Index node1, node2;
 
-    sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
-    const Real kFactor = (Real)sofa::core::mechanicalparams::kFactorIncludingRayleighDamping(mparams, this->rayleighStiffness.getValue());
     const type::vector<typename HexahedralFEMForceField<DataTypes>::HexahedronInformation>& hexahedronInf = d_hexahedronInfo.getValue();
 
     for(Size e=0 ; e<this->l_topology->getNbHexahedra() ; ++e)
@@ -613,7 +611,7 @@ void HexahedralFEMForceField<DataTypes>::addKToMatrix(const core::MechanicalPara
                         Coord(Ke[3*n1+2][3*n2+0],Ke[3*n1+2][3*n2+1],Ke[3*n1+2][3*n2+2])) ) * hexahedronInf[e].rotation;
                 for(i=0; i<3; i++)
                     for (j=0; j<3; j++)
-                        r.matrix->add(r.offset+3*node1+i, r.offset+3*node2+j, - tmp[i][j]*kFactor);
+                        matrix->add(offset+3*node1+i, offset+3*node2+j, - tmp[i][j]*kFact);
             }
         }
     }

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/HexahedralFEMForceFieldAndMass.h

@@ -72,14 +72,14 @@ class HexahedralFEMForceFieldAndMass : virtual public sofa::core::behavior::Mass
      void addMDx(const core::MechanicalParams* mparams, DataVecDeriv& f, const DataVecDeriv& dx, SReal factor) override;
 
     ///// WARNING this method only add diagonal elements in the given matrix !
-    void addMToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
+    void addMToMatrix(sofa::linearalgebra::BaseMatrix * mat, SReal mFact, unsigned int &offset) override;
 
     bool isDiagonal() const override { return d_useLumpedMass.getValue(); }
 
     using HexahedralFEMForceFieldT::addKToMatrix;
     using MassT::addKToMatrix;
     ///// WARNING this method only add diagonal elements in the given matrix !
-    void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
+    void addKToMatrix(sofa::linearalgebra::BaseMatrix * matrix, SReal kFact, unsigned int &offset) override;
 
     ///// WARNING this method only add diagonal elements in the given matrix !
     void addMBKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
@@ -88,6 +88,8 @@ class HexahedralFEMForceFieldAndMass : virtual public sofa::core::behavior::Mass
 
      void addForce(const core::MechanicalParams* mparams, DataVecDeriv& f, const DataVecCoord& x, const DataVecDeriv& v) override;
 
+    using HexahedralFEMForceFieldT::getPotentialEnergy;
+    using MassT::getPotentialEnergy;
     SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord&  /* x */) const override
     {
         msg_warning() << "Method getPotentialEnergy not implemented yet.";