[Core] Introduction of strong types (sofa-framework#5104)

* introduce strong types

* A strong type for matrices factors

* fix

* a more detailed description of StrongType

* More tests

* fix

* comments for MatricesFactors

* deprecate the previous methods

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
     //

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

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);

Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/StaticSolver.cpp

@@ -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.

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);

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

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

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: [email protected]                             *
+******************************************************************************/
+#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>;
+};
+
+}

Sofa/framework/Simulation/Core/src/sofa/simulation/MechanicalOperations.cpp

@@ -288,22 +288,30 @@ 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)
     {
         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])

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;
 

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

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_specialization_of.h