[Mapping.NonLinear] A base class for (almost) all non linear mappings

Sofa/Component/Mapping/NonLinear/CMakeLists.txt

@@ -8,6 +8,8 @@ set(HEADER_FILES
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/init.h
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.h
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.inl
+    ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/BaseNonLinearMapping.h
+    ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/BaseNonLinearMapping.inl
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.h
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.inl
     ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMapping.h

Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h

@@ -22,7 +22,7 @@
 #pragma once
 
 #include <sofa/component/mapping/nonlinear/config.h>
-#include <sofa/core/Mapping.h>
+#include <sofa/component/mapping/nonlinear/BaseNonLinearMapping.h>
 #include <sofa/component/mapping/nonlinear/NonLinearMappingData.h>
 #include <sofa/core/topology/BaseMeshTopology.h>
 #include <sofa/linearalgebra/EigenSparseMatrix.h>
@@ -32,17 +32,16 @@ namespace sofa::component::mapping::nonlinear
 {
 
 template <class TIn, class TOut>
-class AreaMapping : public core::Mapping<TIn, TOut>, public StabilizedNonLinearMappingData
+class AreaMapping : public BaseNonLinearMapping<TIn, TOut, true>
 {
 public:
-    SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut));
+    SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, true));
 
     using In = TIn;
     using Out = TOut;
 
     using Real = Real_t<Out>;
 
-    typedef linearalgebra::EigenSparseMatrix<TIn,TOut> SparseMatrixEigen;
     static constexpr auto Nin = In::deriv_total_size;
 
     SingleLink<AreaMapping<TIn, TOut>, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology;
@@ -52,44 +51,25 @@ class AreaMapping : public core::Mapping<TIn, TOut>, public StabilizedNonLinearM
     void init() override;
 
     void apply(const core::MechanicalParams* mparams, DataVecCoord_t<Out>& out, const DataVecCoord_t<In>& in) override;
-    void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& out, const DataVecDeriv_t<In>& in) override;
-    void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t<In>& out, const DataVecDeriv_t<Out>& in) override;
-    void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t<In>& out, const DataMatrixDeriv_t<Out>& in) override;
-    void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) override;
-
-    void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) override;
-    const linearalgebra::BaseMatrix* getK() override;
     void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override;
 
-    const type::vector<sofa::linearalgebra::BaseMatrix*>* getJs() override;
-
 protected:
     AreaMapping();
 
-    using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix<TIn,TIn>;
-
-    SparseMatrixEigen jacobian; ///< Jacobian of the mapping
-    type::vector<linearalgebra::BaseMatrix*> baseMatrices; ///< Jacobian of the mapping, in a vector
-    typename AreaMapping::SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix
-
-    /**
-     * @brief Represents an entry in the Jacobian matrix.
-     *
-     * The JacobianEntry struct is used to store information about an entry in the
-     * Jacobian matrix, specifically the vertex identifier and the corresponding
-     * Jacobian value. It also provides a comparison operator for sorting entries
-     * by vertex ID.
-     */
-    struct JacobianEntry
-    {
-        sofa::Index vertexId;
-        typename In::Coord jacobianValue;
-        bool operator<(const JacobianEntry& other) const { return vertexId < other.vertexId;}
-    };
+    void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor,
+                            Data<VecDeriv_t<In> >& parentForce,
+                            const Data<VecDeriv_t<In> >& parentDisplacement,
+                            const Data<VecDeriv_t<Out> >& childForce) override;
 
-    const VecCoord_t<In>* m_vertices{nullptr};
+    using typename Inherit1::SparseKMatrixEigen;
 
+    void doUpdateK(
+        const core::MechanicalParams* mparams, const Data<VecDeriv_t<Out> >& childForce,
+        SparseKMatrixEigen& matrix) override;
+
+    const VecCoord_t<In>* m_vertices{nullptr};
 
+    using JacobianEntry = typename Inherit1::JacobianEntry;
 };
 
 #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_AREAMAPPING_CPP)

Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl

@@ -21,8 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/component/mapping/nonlinear/AreaMapping.h>
+#include <sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl>
 #include <sofa/core/BaseLocalMappingMatrix.h>
-#include <sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h>
 
 namespace sofa::component::mapping::nonlinear
 {
@@ -123,10 +123,7 @@ void AreaMapping<TIn, TOut>::init()
     typename core::behavior::MechanicalState<In>::ReadVecCoord pos = this->getFromModel()->readPositions();
 
     this->getToModel()->resize( nbTriangles );
-    jacobian.resizeBlocks(nbTriangles, pos.size());
-
-    baseMatrices.resize( 1 );
-    baseMatrices[0] = &jacobian;
+    this->m_jacobian.resizeBlocks(nbTriangles, pos.size());
 
     Inherit1::init();
 
@@ -140,14 +137,16 @@ template <class TIn, class TOut>
 void AreaMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams,
     DataVecCoord_t<Out>& out, const DataVecCoord_t<In>& in)
 {
+    SOFA_UNUSED( mparams );
+
     helper::WriteOnlyAccessor< Data<VecCoord_t<Out>> > _out = out;
     helper::ReadAccessor< Data<VecCoord_t<In> > > _in = in;
 
     m_vertices = _in.operator->();
 
     const auto& triangles = l_topology->getTriangles();
 
-    jacobian.clear();
+    this->m_jacobian.clear();
 
     for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
     {
@@ -177,132 +176,79 @@ void AreaMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams,
         //insertion in increasing column order
         std::sort(jacobianEntries.begin(), jacobianEntries.end());
 
-        jacobian.beginRow(triangleId);
+        this->m_jacobian.beginRow(triangleId);
         for (const auto& [vertexId, jacobianValue] : jacobianEntries)
         {
             for (unsigned d = 0; d < In::spatial_dimensions; ++d)
             {
-                jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]);
+                this->m_jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]);
             }
         }
     }
 
-    jacobian.compress();
-}
-
-template <class TIn, class TOut>
-void AreaMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mparams,
-    DataVecDeriv_t<Out>& out, const DataVecDeriv_t<In>& in)
-{
-    if( jacobian.rowSize() )
-    {
-        auto dOutWa = sofa::helper::getWriteOnlyAccessor(out);
-        auto dInRa = sofa::helper::getReadAccessor(in);
-        jacobian.mult(dOutWa.wref(),dInRa.ref());
-    }
-}
-
-template <class TIn, class TOut>
-void AreaMapping<TIn, TOut>::applyJT(const core::MechanicalParams* mparams,
-    DataVecDeriv_t<In>& out, const DataVecDeriv_t<Out>& in)
-{
-    if( jacobian.rowSize() )
-    {
-        auto dOutRa = sofa::helper::getReadAccessor(in);
-        auto dInWa = sofa::helper::getWriteOnlyAccessor(out);
-        jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref());
-    }
+    this->m_jacobian.compress();
 }
 
 template <class TIn, class TOut>
-void AreaMapping<TIn, TOut>::applyJT(const core::ConstraintParams* cparams,
-    DataMatrixDeriv_t<In>& out, const DataMatrixDeriv_t<Out>& in)
+void AreaMapping<TIn, TOut>::matrixFreeApplyDJT(
+    const core::MechanicalParams* mparams, Real kFactor,
+    Data<VecDeriv_t<In>>& parentForce,
+    const Data<VecDeriv_t<In>>& parentDisplacement,
+    const Data<VecDeriv_t<Out>>& childForce)
 {
-    SOFA_UNUSED(cparams);
-    auto childMatRa  = sofa::helper::getReadAccessor(in);
-    auto parentMatWa = sofa::helper::getWriteAccessor(out);
-    addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref());
-}
+    SOFA_UNUSED(mparams);
+    const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();
 
-template <class TIn, class TOut>
-void AreaMapping<TIn, TOut>::applyDJT(const core::MechanicalParams* mparams,
-    core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId)
-{
-    if (!m_vertices)
-    {
-        return;
-    }
+    helper::WriteAccessor parentForceAccessor(parentForce);
+    helper::ReadAccessor parentDisplacementAccessor(parentDisplacement);
+    helper::ReadAccessor childForceAccessor(childForce);
 
-    const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
-    if( !geometricStiffness )
+    const auto& triangles = l_topology->getTriangles();
+    for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
     {
-        return;
-    }
-
-    helper::WriteAccessor<Data<VecDeriv_t<In> > > parentForceAccessor(*parentForceId[this->fromModel.get()].write());
-    helper::ReadAccessor<Data<VecDeriv_t<In> > > parentDisplacementAccessor(*mparams->readDx(this->fromModel.get()));
-    const SReal kFactor = mparams->kFactor();
-    helper::ReadAccessor<Data<VecDeriv_t<Out> > > childForceAccessor(mparams->readF(this->toModel.get()));
+        const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];
 
-    if( K.compressedMatrix.nonZeros() )
-    {
-        K.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), (typename In::Real)kFactor );
-    }
-    else
-    {
-        const auto& triangles = l_topology->getTriangles();
-        for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
+        if( childForceTri[0] < 0 || geometricStiffness==1 )
         {
-            const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];
-
-            if( childForceTri[0] < 0 || geometricStiffness==1 )
-            {
-                const auto& triangle = triangles[triangleId];
+            const auto& triangle = triangles[triangleId];
 
-                const sofa::type::fixed_array<Coord_t<In>, 3> v{
-                    (*m_vertices)[triangle[0]],
-                    (*m_vertices)[triangle[1]],
-                    (*m_vertices)[triangle[2]]
-                };
+            const sofa::type::fixed_array<Coord_t<In>, 3> v{
+                (*m_vertices)[triangle[0]],
+                (*m_vertices)[triangle[1]],
+                (*m_vertices)[triangle[2]]
+            };
 
-                //it's a 3x3 matrix, where each entry is a 3x3 matrix
-                const auto d2Area_d2x = computeSecondDerivativeArea(v);
+            //it's a 3x3 matrix, where each entry is a 3x3 matrix
+            const auto d2Area_d2x = computeSecondDerivativeArea(v);
 
-                for (unsigned int i = 0; i < 3; ++i)
+            for (unsigned int i = 0; i < 3; ++i)
+            {
+                for (unsigned int j = 0; j < 3; ++j)
                 {
-                    for (unsigned int j = 0; j < 3; ++j)
-                    {
-                        parentForceAccessor[triangle[i]] +=
-                            kFactor
-                            * d2Area_d2x[i][j]
-                            * parentDisplacementAccessor[triangle[j]]
-                            * childForceTri[0];
-                    }
+                    parentForceAccessor[triangle[i]] +=
+                        kFactor
+                        * d2Area_d2x[i][j]
+                        * parentDisplacementAccessor[triangle[j]]
+                        * childForceTri[0];
                 }
             }
         }
     }
 }
 
 template <class TIn, class TOut>
-void AreaMapping<TIn, TOut>::updateK(const core::MechanicalParams* mparams,
-    core::ConstMultiVecDerivId childForceId)
+void AreaMapping<TIn, TOut>::doUpdateK(const core::MechanicalParams* mparams,
+    const Data<VecDeriv_t<Out>>& childForce, SparseKMatrixEigen& matrix)
 {
     SOFA_UNUSED(mparams);
-    const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
-    if( !geometricStiffness ) { K.resize(0,0); return; }
+    const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();
 
-    helper::ReadAccessor<Data<VecDeriv_t<Out> > > childForce( *childForceId[this->toModel.get()].read() );
-
-    {
-        unsigned int kSize = this->fromModel->getSize();
-        K.resizeBlocks(kSize, kSize);
-    }
+    const helper::ReadAccessor childForceAccessor(childForce);
 
     const auto& triangles = l_topology->getTriangles();
     for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
     {
-        const Deriv_t<Out>& childForceTri = childForce[triangleId];
+        const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];
 
         if( childForceTri[0] < 0 || geometricStiffness==1 )
         {
@@ -321,26 +267,18 @@ void AreaMapping<TIn, TOut>::updateK(const core::MechanicalParams* mparams,
             {
                 for (unsigned int j = 0; j < 3; ++j)
                 {
-                    K.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]);
+                    matrix.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]);
                 }
             }
         }
     }
-
-    K.compress();
-}
-
-template <class TIn, class TOut>
-const linearalgebra::BaseMatrix* AreaMapping<TIn, TOut>::getK()
-{
-    return &K;
 }
 
 template <class TIn, class TOut>
 void AreaMapping<TIn, TOut>::buildGeometricStiffnessMatrix(
     sofa::core::GeometricStiffnessMatrix* matrices)
 {
-    const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
+    const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();
     if( !geometricStiffness )
     {
         return;
@@ -379,11 +317,5 @@ void AreaMapping<TIn, TOut>::buildGeometricStiffnessMatrix(
 
 }
 
-template <class TIn, class TOut>
-const type::vector<sofa::linearalgebra::BaseMatrix*>* AreaMapping<TIn, TOut>::
-getJs()
-{
-    return &baseMatrices;
-}
 
 }