Merge branch 'master' into add_fixedarray_accessor

.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

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;
 

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

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 {}
 
 
 

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

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;

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)

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/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/ShewchukPCGLinearSolver.cpp → Sofa/Component/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp

@@ -19,21 +19,21 @@
 *                                                                             *
 * Contact information: [email protected]                             *
 ******************************************************************************/
-#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>
 
 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