diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.h
index 236ab760b61..2da102806d3 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.h
@@ -71,7 +71,11 @@ class BarycentricMapping : public LinearMapping<TIn, TOut>
     void init() override;
     void reinit() override;
     void apply(const core::MechanicalParams *mparams, Data< typename Out::VecCoord >& out, const Data< typename In::VecCoord >& in) override;
-    void applyJ(const core::MechanicalParams *mparams, Data< typename Out::VecDeriv >& out, const Data< typename In::VecDeriv >& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
     void applyJT(const core::MechanicalParams *mparams, Data< typename In::VecDeriv >& out, const Data< typename Out::VecDeriv >& in) override;
     void applyJT(const core::ConstraintParams *cparams, Data< typename In::MatrixDeriv >& out, const Data< typename Out::MatrixDeriv >& in) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.inl
index ef6ae436d07..384fb6232ca 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BarycentricMapping.inl
@@ -303,14 +303,19 @@ void BarycentricMapping<TIn, TOut>::apply(const core::MechanicalParams * mparams
 }
 
 template <class TIn, class TOut>
-void BarycentricMapping<TIn, TOut>::applyJ (const core::MechanicalParams * mparams, Data< typename Out::VecDeriv >& _out, const Data< typename In::VecDeriv >& in)
+void BarycentricMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
     if (d_mapper != nullptr)
     {
-        auto outWriteAccessor = sofa::helper::getWriteAccessor(_out);
-        d_mapper->applyJ(outWriteAccessor.wref(), in.getValue());
+        auto outWriteAccessor = sofa::helper::getWriteAccessor(tangentOut);
+        d_mapper->applyJ(outWriteAccessor.wref(), tangentIn.getValue());
     }
 }
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.h
index 95abbb15dc1..4bef072f30b 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.h
@@ -99,7 +99,11 @@ class BeamLinearMapping : public LinearMapping<TIn, TOut>
 
     void apply(const core::MechanicalParams *mparams, Data< typename Out::VecCoord >& out, const Data< typename In::VecCoord >& in) override;
 
-    void applyJ(const core::MechanicalParams *mparams, Data< typename Out::VecDeriv >& out, const Data< typename In::VecDeriv >& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT(const core::MechanicalParams *mparams, Data< typename In::VecDeriv >& out, const Data< typename Out::VecDeriv >& in) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.inl
index 99eba171599..654a6022a60 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/BeamLinearMapping.inl
@@ -134,10 +134,15 @@ void BeamLinearMapping<TIn, TOut>::apply(const core::MechanicalParams * /*mparam
 }
 
 template <class TIn, class TOut>
-void BeamLinearMapping<TIn, TOut>::applyJ(const core::MechanicalParams * /*mparams*/, Data< typename Out::VecDeriv >& _out, const Data< typename In::VecDeriv >& _in)
+void BeamLinearMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut,
+    const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    helper::WriteAccessor< Data< typename Out::VecDeriv > > out = _out;
-    helper::ReadAccessor< Data< typename In::VecDeriv > > in = _in;
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    helper::WriteAccessor< Data< typename Out::VecDeriv > > out = tangentOut;
+    helper::ReadAccessor< Data< typename In::VecDeriv > > in = tangentIn;
 
     out.resize(points.size());
     for(unsigned int i=0; i<points.size(); i++)
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.h
index 92079a44b9c..adb370f7f4e 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.h
@@ -78,8 +78,11 @@ class CenterOfMassMapping : public LinearMapping<TIn, TOut>
     void apply( const sofa::core::MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in) override;
     //void apply( typename Out::VecCoord& out, const typename In::VecCoord& in );
 
-    void applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) override;
-    //void applyJ( typename Out::VecDeriv& out, const typename In::VecDeriv& in );
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT( const sofa::core::MechanicalParams* mparams, InDataVecDeriv& out, const OutDataVecDeriv& in) override;
     //void applyJT( typename In::VecDeriv& out, const typename Out::VecDeriv& in );
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.inl
index 15dcd77f4a5..d7adac4981d 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMapping.inl
@@ -80,14 +80,16 @@ void CenterOfMassMapping<TIn, TOut>::apply( const sofa::core::MechanicalParams*
     outData.endEdit();
 }
 
-
 template <class TIn, class TOut>
-void CenterOfMassMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& outData, const InDataVecDeriv& inData)
+void CenterOfMassMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut,
+    const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
-    OutVecDeriv& childForces = *outData.beginEdit();
-    const InVecDeriv& parentForces = inData.getValue();
+    OutVecDeriv& childForces = *tangentOut.beginEdit();
+    const InVecDeriv& parentForces = tangentIn.getValue();
 
     if(!masses || totalMass==0.0)
     {
@@ -106,10 +108,9 @@ void CenterOfMassMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalParams*
 
     childForces[0] = outF / totalMass;
 
-    outData.endEdit();
+    tangentOut.endEdit();
 }
 
-
 template <class TIn, class TOut>
 void CenterOfMassMapping<TIn, TOut>::applyJT( const sofa::core::MechanicalParams* mparams, InDataVecDeriv& outData, const OutDataVecDeriv& inData)
 {
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.h
index e4ecbfc49d3..e7f00c5f924 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.h
@@ -84,10 +84,13 @@ class CenterOfMassMulti2Mapping : public LinearMulti2Mapping<TIn1, TIn2, TOut>
         const core::MechanicalParams* mparams, const type::vector<OutDataVecCoord*>& dataVecOutPos,
         const type::vector<const In1DataVecCoord*>& dataVecIn1Pos ,
         const type::vector<const In2DataVecCoord*>& dataVecIn2Pos) override;
-    void applyJ(
-        const core::MechanicalParams* mparams, const type::vector< OutDataVecDeriv*>& dataVecOutVel,
-        const type::vector<const In1DataVecDeriv*>& dataVecIn1Vel,
-        const type::vector<const In2DataVecDeriv*>& dataVecIn2Vel) override;
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In1>*>& dataVecTangentIn1,
+        const type::vector<const DataVecDeriv_t<In2>*>& dataVecTangentIn2,
+        const type::vector<const DataVecCoord_t<In1>*>& dataVecPositionIn1,
+        const type::vector<const DataVecCoord_t<In2>*>& dataVecPositionIn2) override;
     void applyJT(
         const core::MechanicalParams* mparams, const type::vector< In1DataVecDeriv*>& dataVecOut1Force,
         const type::vector< In2DataVecDeriv*>& dataVecOut2Force,
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.inl
index 164c5d8b7e8..c161cabae72 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMulti2Mapping.inl
@@ -94,31 +94,34 @@ void CenterOfMassMulti2Mapping< TIn1, TIn2, TOut >::apply(
     for(unsigned int i=0; i<dataVecOutPos.size(); i++)
         dataVecOutPos[i]->endEdit();
 }
-
 template <class TIn1, class TIn2, class TOut>
-void CenterOfMassMulti2Mapping< TIn1, TIn2, TOut >::applyJ(
-        const core::MechanicalParams* mparams, const type::vector< OutDataVecDeriv*>& dataVecOutVel,
-        const type::vector<const In1DataVecDeriv*>& dataVecIn1Vel,
-        const type::vector<const In2DataVecDeriv*>& dataVecIn2Vel)
-//applyJ(const type::vector< OutVecDeriv*>& outDeriv, const type::vector<const In1VecDeriv*>& inDeriv1, const type::vector<const In2VecDeriv*>& inDeriv2)
+void CenterOfMassMulti2Mapping<TIn1, TIn2, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    const type::vector<DataVecDeriv_t<Out>*>& dataVecTangentOut,
+    const type::vector<const DataVecDeriv_t<In1>*>& dataVecTangentIn1,
+    const type::vector<const DataVecDeriv_t<In2>*>& dataVecTangentIn2,
+    const type::vector<const DataVecCoord_t<In1>*>& dataVecPositionIn1,
+    const type::vector<const DataVecCoord_t<In2>*>& dataVecPositionIn2)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(dataVecPositionIn1);
+    SOFA_UNUSED(dataVecPositionIn2);
 
-    assert( dataVecOutVel.size() == 1 );
+    assert( dataVecTangentOut.size() == 1 );
     typedef typename type::vector<In1Deriv>::iterator                     iter_deriv1;
     typedef typename type::vector<In2Deriv>::iterator                     iter_deriv2;
 
     //Not optimized at all...
     type::vector<OutVecDeriv*> outDeriv;
-    for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-        outDeriv.push_back(dataVecOutVel[i]->beginEdit());
+    for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+        outDeriv.push_back(dataVecTangentOut[i]->beginEdit());
 
     type::vector<const In1VecDeriv*> inDeriv1;
-    for(unsigned int i=0; i<dataVecIn1Vel.size(); i++)
-        inDeriv1.push_back(&dataVecIn1Vel[i]->getValue());
+    for(unsigned int i=0; i<dataVecTangentIn1.size(); i++)
+        inDeriv1.push_back(&dataVecTangentIn1[i]->getValue());
     type::vector<const In2VecDeriv*> inDeriv2;
-    for(unsigned int i=0; i<dataVecIn2Vel.size(); i++)
-        inDeriv2.push_back(&dataVecIn2Vel[i]->getValue());
+    for(unsigned int i=0; i<dataVecTangentIn2.size(); i++)
+        inDeriv2.push_back(&dataVecTangentIn2[i]->getValue());
 
     SReal px=0,py=0,pz=0;
 
@@ -155,8 +158,8 @@ void CenterOfMassMulti2Mapping< TIn1, TIn2, TOut >::applyJ(
     OutDataTypes::set((*outVecDeriv)[0], px,py,pz);
 
     //Really Not optimized at all...
-    for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-        dataVecOutVel[i]->endEdit();
+    for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+        dataVecTangentOut[i]->endEdit();
 }
 
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.h
index e6143072db6..87105a78b9c 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.h
@@ -73,8 +73,11 @@ class CenterOfMassMultiMapping : public LinearMultiMapping<TIn, TOut>
     void apply(const core::MechanicalParams* mparams, const type::vector<OutDataVecCoord*>& dataVecOutPos, const type::vector<const InDataVecCoord*>& dataVecInPos) override;
     //virtual void apply(const type::vector<OutVecCoord*>& outPos, const vecConstInVecCoord& inPos);
 
-    void applyJ(const core::MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel) override;
-    //virtual void applyJ(const type::vector<OutVecDeriv*>& outDeriv, const type::vector<const  InVecDeriv*>& inDeriv);
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+        const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn) override;
 
     void applyJT(const core::MechanicalParams* mparams, const type::vector<InDataVecDeriv*>& dataVecOutForce, const type::vector<const OutDataVecDeriv*>& dataVecInForce) override;
     //virtual void applyJT(const type::vector< InVecDeriv*>& outDeriv, const type::vector<const OutVecDeriv*>& inDeriv);
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.inl
index 0ad2447ea9c..7ba9e251187 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/CenterOfMassMultiMapping.inl
@@ -66,22 +66,26 @@ void CenterOfMassMultiMapping< TIn, TOut >::apply(const core::MechanicalParams*
         dataVecOutPos[i]->endEdit();
 }
 
-
 template <class TIn, class TOut>
-void CenterOfMassMultiMapping< TIn, TOut >::applyJ(const core::MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel)
+void CenterOfMassMultiMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    const type::vector<DataVecDeriv_t<Out>*>& dataVecTangentOut,
+    const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+    const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(dataVecPositionIn);
 
     typedef typename InVecDeriv::iterator iter_deriv;
 
     //Not optimized at all...
     type::vector<OutVecDeriv*> outDeriv;
-    for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-        outDeriv.push_back(dataVecOutVel[i]->beginEdit());
+    for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+        outDeriv.push_back(dataVecTangentOut[i]->beginEdit());
 
     type::vector<const InVecDeriv*> inDeriv;
-    for(unsigned int i=0; i<dataVecInVel.size(); i++)
-        inDeriv.push_back(&dataVecInVel[i]->getValue());
+    for(unsigned int i=0; i<dataVecTangentIn.size(); i++)
+        inDeriv.push_back(&dataVecTangentIn[i]->getValue());
 
     assert( outDeriv.size() == 1 );
 
@@ -98,8 +102,8 @@ void CenterOfMassMultiMapping< TIn, TOut >::applyJ(const core::MechanicalParams*
     OutDataTypes::set((*outVecDeriv)[0], x,y,z);
 
     //Really Not optimized at all...
-    for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-        dataVecOutVel[i]->endEdit();
+    for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+        dataVecTangentOut[i]->endEdit();
 }
 
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.h
index 2daa51fb3ea..a7d08063880 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.h
@@ -127,10 +127,13 @@ class DeformableOnRigidFrameMapping : public LinearMulti2Mapping<TIn, TInRoot, T
 
     //ApplyJ
     void applyJ( OutVecDeriv& out, const InVecDeriv& in, const InRootVecDeriv* inroot );
-    void applyJ(
-        const core::MechanicalParams* /* mparams */, const type::vector< OutDataVecDeriv*>& dataVecOutVel,
-        const type::vector<const InDataVecDeriv*>& dataVecInVel,
-        const type::vector<const InRootDataVecDeriv*>& dataVecInRootVel) override;
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<TIn>*>& dataVecTangentIn1,
+        const type::vector<const DataVecDeriv_t<TInRoot>*>& dataVecTangentIn2,
+        const type::vector<const DataVecCoord_t<TIn>*>& dataVecPositionIn1,
+        const type::vector<const DataVecCoord_t<TInRoot>*>& dataVecPositionIn2) override;
 
     //ApplyJT Force
     void applyJT( InVecDeriv& out, const OutVecDeriv& in, InRootVecDeriv* outroot );
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.inl
index 195104551be..557943aeb00 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DeformableOnRigidFrameMapping.inl
@@ -281,26 +281,33 @@ void DeformableOnRigidFrameMapping<TIn, TInRoot, TOut>::applyJ( typename Out::Ve
 }
 
 template <class TIn, class TInRoot, class TOut>
-void DeformableOnRigidFrameMapping<TIn, TInRoot, TOut>::applyJ(
-    const core::MechanicalParams* /* mparams */, const type::vector< OutDataVecDeriv*>& dataVecOutVel,
-    const type::vector<const InDataVecDeriv*>& dataVecInVel,
-    const type::vector<const InRootDataVecDeriv*>& dataVecInRootVel)
+void DeformableOnRigidFrameMapping<TIn, TInRoot, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    const type::vector<DataVecDeriv_t<Out>*>& dataVecTangentOut,
+    const type::vector<const DataVecDeriv_t<TIn>*>& dataVecTangentIn1,
+    const type::vector<const DataVecDeriv_t<TInRoot>*>& dataVecTangentIn2,
+    const type::vector<const DataVecCoord_t<TIn>*>& dataVecPositionIn1,
+    const type::vector<const DataVecCoord_t<TInRoot>*>& dataVecPositionIn2)
 {
-    if(dataVecOutVel.empty() || dataVecInVel.empty())
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(dataVecPositionIn1);
+    SOFA_UNUSED(dataVecPositionIn2);
+
+    if(dataVecTangentOut.empty() || dataVecTangentIn1.empty())
         return;
 
     const InRootVecDeriv* inroot = nullptr;
 
     //We need only one input In model and input Root model (if present)
-    OutVecDeriv& out = *dataVecOutVel[0]->beginEdit();
-    const InVecDeriv& in = dataVecInVel[0]->getValue();
+    OutVecDeriv& out = *dataVecTangentOut[0]->beginEdit();
+    const InVecDeriv& in = dataVecTangentIn1[0]->getValue();
 
-    if (!dataVecInRootVel.empty())
-        inroot = &dataVecInRootVel[0]->getValue();
+    if (!dataVecTangentIn2.empty())
+        inroot = &dataVecTangentIn2[0]->getValue();
 
     applyJ(out,in, inroot);
 
-    dataVecOutVel[0]->endEdit();
+    dataVecTangentOut[0]->endEdit();
 }
 
 template <class TIn, class TInRoot, class TOut>
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.h
index 78ec90b91c7..9bc51291f5d 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.h
@@ -46,7 +46,11 @@ class DistanceToPlaneMapping : public LinearMapping<TIn, defaulttype::Vec1Types>
 
     void apply(const core::MechanicalParams *mparams, Data<VecCoord_t<TOut>>& out, const Data<VecCoord_t<TIn>>& in) override;
 
-    void applyJ(const core::MechanicalParams *mparams, Data<VecDeriv_t<TOut>>& out, const Data<VecDeriv_t<TIn>>& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<TOut>& tangentOut,
+        const DataVecDeriv_t<TIn>& tangentIn,
+        const DataVecCoord_t<TIn>& positionIn) override;
 
     void applyJT(const core::MechanicalParams *mparams, Data<VecDeriv_t<TIn>>& out, const Data<VecDeriv_t<TOut>>& in) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.inl
index 2e1cabe3043..5f54d0d5e61 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/DistanceToPlaneMapping.inl
@@ -100,17 +100,19 @@ void DistanceToPlaneMapping<TIn>::apply(const core::MechanicalParams *mparams, D
         writeOut[i] = type::dot(TIn::getCPos(readIn[i]),planeNormal) - planeDistanceToOrigin;
     }
 }
-
 template <class TIn>
-void DistanceToPlaneMapping<TIn>::applyJ(const core::MechanicalParams *mparams, Data<VecDeriv_t<TOut>>& out, const Data<VecDeriv_t<TIn>>& in)
+void DistanceToPlaneMapping<TIn>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<TOut>& tangentOut,
+    const DataVecDeriv_t<TIn>& tangentIn, const DataVecCoord_t<TIn>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
     if (this-> d_componentState.getValue() != sofa::core::objectmodel::ComponentState::Valid)
         return;
 
-    auto writeOut = helper::getWriteAccessor(out);
-    const auto readIn = helper::getReadAccessor(in);
+    auto writeOut = helper::getWriteAccessor(tangentOut);
+    const auto readIn = helper::getReadAccessor(tangentIn);
     const auto planeNormal = d_planeNormal.getValue();
 
     for ( unsigned i = 0; i<readIn.size(); i++ )
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.h
index 5b926849ced..c593850df33 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.h
@@ -104,7 +104,11 @@ class IdentityMapping : public LinearMapping<TIn, TOut>
 
     void apply(const core::MechanicalParams *mparams, Data<VecCoord>& out, const Data<InVecCoord>& in) override;
 
-    void applyJ(const core::MechanicalParams *mparams, Data<VecDeriv>& out, const Data<InVecDeriv>& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT(const core::MechanicalParams *mparams, Data<InVecDeriv>& out, const Data<VecDeriv>& in) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.inl
index ecfeb57aff8..2ca7ad0cdbb 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMapping.inl
@@ -73,12 +73,18 @@ void IdentityMapping<TIn, TOut>::apply(const core::MechanicalParams * /*mparams*
 }
 
 template <class TIn, class TOut>
-void IdentityMapping<TIn, TOut>::applyJ(const core::MechanicalParams * /*mparams*/, Data<VecDeriv>& dOut, const Data<InVecDeriv>& dIn)
+void IdentityMapping<TIn, TOut>::doApplyJacobianVectorProduct(const core::MechanicalParams* mparams,
+                                                              DataVecDeriv_t<Out>& tangentOut,
+                                                              const DataVecDeriv_t<In>& tangentIn,
+                                                              const DataVecCoord_t<In>& positionIn)
 {
-    helper::WriteOnlyAccessor< Data<VecDeriv> > out = dOut;
-    helper::ReadAccessor< Data<InVecDeriv> > in = dIn;
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
-    for( size_t i=0 ; i<out.size() ; ++i)
+    helper::WriteOnlyAccessor< Data<VecDeriv> > out = tangentOut;
+    helper::ReadAccessor< Data<InVecDeriv> > in = tangentIn;
+
+    for (size_t i = 0; i < out.size(); ++i)
     {
         core::eq(out[i], in[i]);
     }
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 c287432d43d..612130cb164 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
@@ -68,7 +68,11 @@ class IdentityMultiMapping : public LinearMultiMapping<TIn, TOut>
     void init() override;
 
     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;
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+        const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn) override;
     void applyJT(const core::MechanicalParams* mparams, const type::vector<InDataVecDeriv*>& dataVecOutForce, const type::vector<const OutDataVecDeriv*>& dataVecInForce) override;
     void applyDJT(const core::MechanicalParams* /*mparams*/, core::MultiVecDerivId /*enforce*/, core::ConstMultiVecDerivId /*outForce*/) override {}
     void applyJT( const core::ConstraintParams* cparams, const type::vector< InDataMatrixDeriv* >& dataMatOutConst, const type::vector< const OutDataMatrixDeriv* >& dataMatInConst ) override;
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.inl
index 0cad33c7bc8..176c15bb855 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/IdentityMultiMapping.inl
@@ -115,16 +115,21 @@ void IdentityMultiMapping<TIn, TOut>::apply(const core::MechanicalParams* mparam
 }
 
 template <class TIn, class TOut>
-void IdentityMultiMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel)
+void IdentityMultiMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    const type::vector<DataVecDeriv_t<Out>*>& dataVecTangentOut,
+    const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+    const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(dataVecPositionIn);
 
-    OutVecDeriv& out = *(dataVecOutVel[0]->beginEdit());
+    OutVecDeriv& out = *(dataVecTangentOut[0]->beginEdit());
 
     unsigned offset = 0;
-    for(unsigned i=0; i<dataVecInVel.size(); i++ )
+    for(unsigned i=0; i<dataVecTangentIn.size(); i++ )
     {
-        const InVecDeriv& in = dataVecInVel[i]->getValue();
+        const InVecDeriv& in = dataVecTangentIn[i]->getValue();
 
         for(unsigned int j=0; j<in.size(); j++)
         {
@@ -133,7 +138,7 @@ void IdentityMultiMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mpara
         offset += in.size();
     }
 
-    dataVecOutVel[0]->endEdit();
+    dataVecTangentOut[0]->endEdit();
 }
 
 template <class TIn, class TOut>
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.h
index aefaaec6ff2..d12fb7e3876 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.h
@@ -87,8 +87,11 @@ class LineSetSkinningMapping : public LinearMapping<TIn, TOut>
     void apply( const sofa::core::MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in) override;
     //void apply( typename Out::VecCoord& out, const typename In::VecCoord& in );
 
-    void applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) override;
-    //void applyJ( typename Out::VecDeriv& out, const typename In::VecDeriv& in );
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT( const sofa::core::MechanicalParams* mparams, InDataVecDeriv& out, const OutDataVecDeriv& in) override;
     //void applyJT( typename In::VecDeriv& out, const typename Out::VecDeriv& in );
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.inl
index bb1802830e5..d88abb5d49f 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/LineSetSkinningMapping.inl
@@ -245,13 +245,18 @@ void LineSetSkinningMapping<TIn, TOut>::apply( const sofa::core::MechanicalParam
 }
 
 template <class TIn, class TOut>
-void LineSetSkinningMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& outData, const InDataVecDeriv& inData)
+void LineSetSkinningMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
     const InVecCoord& xfrom = this->fromModel->read(core::vec_id::read_access::position)->getValue();
-    OutVecDeriv& out = *outData.beginEdit();
-    const InVecDeriv& in = inData.getValue();
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     for (unsigned int verticeIndex=0; verticeIndex<out.size(); verticeIndex++)
     {
         out[verticeIndex] = typename Out::Deriv();
@@ -262,7 +267,7 @@ void LineSetSkinningMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalPara
             out[verticeIndex] += (getVCenter(in[m_topology->getLine(iline.lineIndex)[0]]) - IP.cross(getVOrientation(in[m_topology->getLine(iline.lineIndex)[0]]))) * iline.weight;
         }
     }
-    outData.endEdit();
+    tangentOut.endEdit();
 }
 
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.h
index 10aea9ea318..c5cdab8e5b3 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.h
@@ -70,7 +70,11 @@ class Mesh2PointMechanicalMapping : public LinearMapping<TIn, TOut>
 
     void apply(const core::MechanicalParams *mparams, Data<OutVecCoord>& out, const Data<InVecCoord>& in) override;
 
-    void applyJ(const core::MechanicalParams *mparams, Data<OutVecDeriv>& out, const Data<InVecDeriv>& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<TOut>& tangentOut,
+        const DataVecDeriv_t<TIn>& tangentIn,
+        const DataVecCoord_t<TIn>& positionIn) override;
 
     void applyJT(const core::MechanicalParams *mparams, Data<InVecDeriv>& out, const Data<OutVecDeriv>& in) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.inl
index c7b854648da..429a97a8f5a 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/Mesh2PointMechanicalMapping.inl
@@ -192,16 +192,20 @@ void Mesh2PointMechanicalMapping<TIn, TOut>::apply(const core::MechanicalParams
         }
     }
 }
-
 template <class TIn, class TOut>
-void Mesh2PointMechanicalMapping<TIn, TOut>::applyJ(const core::MechanicalParams * /*mparams*/, Data<OutVecDeriv>& dOut, const Data<InVecDeriv>& dIn)
+void Mesh2PointMechanicalMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<TOut>& tangentOut,
+    const DataVecDeriv_t<TIn>& tangentIn, const DataVecCoord_t<TIn>& positionIn)
 {
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
     if (!l_topologicalMapping) return;
 
     using sofa::InvalidID;
 
-    helper::WriteAccessor< Data<OutVecDeriv> > out = dOut;
-    helper::ReadAccessor< Data<InVecDeriv> > in = dIn;
+    helper::WriteAccessor< Data<OutVecDeriv> > out = tangentOut;
+    helper::ReadAccessor< Data<InVecDeriv> > in = tangentIn;
 
     const auto& pointMap = l_topologicalMapping->getPointsMappedFromPoint();
     const auto& edgeMap = l_topologicalMapping->getPointsMappedFromEdge();
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.h
index 20c8aeabc4c..fb3eea83a0a 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.h
@@ -75,7 +75,11 @@ class SimpleTesselatedTetraMechanicalMapping : public LinearMapping<TIn, TOut>
 
     void apply(const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn) override;
 
-    void applyJ(const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<TOut>& tangentOut,
+        const DataVecDeriv_t<TIn>& tangentIn,
+        const DataVecCoord_t<TIn>& positionIn) override;
 
     void applyJT(const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.inl
index a972cb86ff6..32842953d15 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SimpleTesselatedTetraMechanicalMapping.inl
@@ -80,8 +80,12 @@ void SimpleTesselatedTetraMechanicalMapping<TIn, TOut>::apply ( const core::Mech
 }
 
 template <class TIn, class TOut>
-void SimpleTesselatedTetraMechanicalMapping<TIn, TOut>::applyJ( const core::MechanicalParams* /* mparams */, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SimpleTesselatedTetraMechanicalMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<TOut>& tangentOut,
+    const DataVecDeriv_t<TIn>& tangentIn, const DataVecCoord_t<TIn>& positionIn)
 {
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
     if (!topoMap) return;
     const auto& pointMap = topoMap->getPointMappedFromPoint();
@@ -89,8 +93,8 @@ void SimpleTesselatedTetraMechanicalMapping<TIn, TOut>::applyJ( const core::Mech
     if (pointMap.empty() && edgeMap.empty()) return;
     const core::topology::BaseMeshTopology::SeqEdges& edges = inputTopo->getEdges();
 
-    helper::ReadAccessor<InDataVecDeriv> in = dIn;
-    helper::WriteAccessor<OutDataVecDeriv> out = dOut;
+    helper::ReadAccessor<InDataVecDeriv> in = tangentIn;
+    helper::WriteAccessor<OutDataVecDeriv> out = tangentOut;
 
     out.resize(outputTopo->getNbPoints());
     for(Size i = 0; i < pointMap.size(); ++i)
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.h
index 2b03b9459f7..a687e327ac4 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.h
@@ -106,8 +106,11 @@ class SkinningMapping : public LinearMapping<TIn, TOut>
     void apply( const sofa::core::MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in) override;
     //void apply( typename Out::VecCoord& out, const typename In::VecCoord& in );
 
-    void applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) override;
-    //void applyJ( typename Out::VecDeriv& out, const typename In::VecDeriv& in );
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT( const sofa::core::MechanicalParams* mparams, InDataVecDeriv& out, const OutDataVecDeriv& in) override;
     //void applyJT( typename In::VecDeriv& out, const typename Out::VecDeriv& in );
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.inl
index e0790d829b0..39116da20c3 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SkinningMapping.inl
@@ -233,15 +233,17 @@ void SkinningMapping<TIn, TOut>::apply( const sofa::core::MechanicalParams* mpar
     }
     outData.endEdit();
 }
-
-
 template <class TIn, class TOut>
-void SkinningMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& outData, const InDataVecDeriv& inData)
+void SkinningMapping<TIn, TOut>::doApplyJacobianVectorProduct(const core::MechanicalParams* mparams,
+                                                              DataVecDeriv_t<Out>& tangentOut,
+                                                              const DataVecDeriv_t<In>& tangentIn,
+                                                              const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
-    OutVecDeriv& out = *outData.beginWriteOnly();
-    const InVecDeriv& in = inData.getValue();
+    OutVecDeriv& out = *tangentOut.beginWriteOnly();
+    const InVecDeriv& in = tangentIn.getValue();
 
     unsigned int nbref=d_nbRef.getValue()[0];
     sofa::helper::ReadAccessor<Data<type::vector<sofa::type::SVector<InReal> > > > m_weights  (d_weight );
@@ -261,7 +263,7 @@ void SkinningMapping<TIn, TOut>::applyJ( const sofa::core::MechanicalParams* mpa
         }
     }
 
-    outData.endEdit();
+    tangentOut.endEdit();
 }
 
 template <class TIn, class TOut>
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.h
index 21858156950..227b1794df2 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.h
@@ -110,7 +110,11 @@ class SubsetMapping : public LinearMapping<TIn, TOut>
 
     void apply ( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn ) override;
 
-    void applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn ) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT ( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn ) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.inl
index 5f6b68e08c2..212a5932703 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/SubsetMapping.inl
@@ -209,12 +209,18 @@ void SubsetMapping<TIn, TOut>::apply ( const core::MechanicalParams* /*mparams*/
 }
 
 template <class TIn, class TOut>
-void SubsetMapping<TIn, TOut>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SubsetMapping<TIn, TOut>::doApplyJacobianVectorProduct(const core::MechanicalParams* mparams,
+                                                            DataVecDeriv_t<Out>& tangentOut,
+                                                            const DataVecDeriv_t<In>& tangentIn,
+                                                            const DataVecCoord_t<In>& positionIn)
 {
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
     const IndexArray& indices = d_indices.getValue();
 
-    const InVecDeriv& in = dIn.getValue();
-    OutVecDeriv& out = *dOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
+    OutVecDeriv& out = *tangentOut.beginEdit();
     const std::size_t fromSize = in.size();
 
     out.resize(indices.size());
@@ -224,7 +230,7 @@ void SubsetMapping<TIn, TOut>::applyJ( const core::MechanicalParams* /*mparams*/
             out[i] = in[ indices[i] ];
     }
 
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <class TIn, class 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 8fc5781b9a3..892d6b2ec9e 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
@@ -54,7 +54,11 @@ class SubsetMultiMapping : public LinearMultiMapping<TIn, TOut>
 
     // usual Mapping API
     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;
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+        const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn) 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 {}
 
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 b788abdd744..ffd572c4f7a 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
@@ -160,10 +160,16 @@ void SubsetMultiMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams,
 }
 
 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)
+void SubsetMultiMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    const type::vector<DataVecDeriv_t<Out>*>& dataVecTangentOut,
+    const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+    const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn)
 {
     SOFA_UNUSED(mparams);
-    apply_impl( dataVecOutVel, dataVecInVel, d_indexPairs.getValue());
+    SOFA_UNUSED(dataVecPositionIn);
+
+    apply_impl( dataVecTangentOut, dataVecTangentIn, d_indexPairs.getValue());
 }
 
 template <class TIn, class TOut>
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.h
index 4677d59287a..2bac27a709b 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.h
@@ -80,7 +80,11 @@ class TubularMapping : public LinearMapping<TIn, TOut>
 
     void apply ( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn ) override;
 
-    void applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn ) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT ( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn ) override;
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.inl b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.inl
index cf8dd51108a..391db4f96a3 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.inl
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/TubularMapping.inl
@@ -124,14 +124,21 @@ void TubularMapping<TIn, TOut>::apply ( const core::MechanicalParams* /* mparams
 
 
 template <class TIn, class TOut>
-void TubularMapping<TIn, TOut>::applyJ( const core::MechanicalParams* /* mparams */, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void TubularMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn)
 {
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
     if (this->d_componentState.getValue() != sofa::core::objectmodel::ComponentState::Valid)
         return;
 
     // Propagation of velocities from the input DOFs to the output DOFs
-    const InVecDeriv& in = dIn.getValue();
-    OutVecDeriv& out = *dOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
+    OutVecDeriv& out = *tangentOut.beginEdit();
 
     unsigned int N = d_nbPointsOnEachCircle.getValue();
 
@@ -154,7 +161,7 @@ void TubularMapping<TIn, TOut>::applyJ( const core::MechanicalParams* /* mparams
         }
     }
 
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 
diff --git a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/VoidMapping.h b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/VoidMapping.h
index 8a52296593b..e8003b0a9b1 100644
--- a/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/VoidMapping.h
+++ b/Sofa/Component/Mapping/Linear/src/sofa/component/mapping/linear/VoidMapping.h
@@ -83,8 +83,16 @@ class VoidMapping : public LinearBaseMapping
     {
     }
 
-    void applyJ(const core::MechanicalParams* /* mparams */, core::MultiVecDerivId /* outVel */, core::ConstMultiVecDerivId /* inVel */) override
+    void applyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        core::MultiVecDerivId tangentOutId,
+        core::ConstMultiVecDerivId tangentInId,
+        core::ConstMultiVecCoordId positionInId) override
     {
+        SOFA_UNUSED(mparams);
+        SOFA_UNUSED(tangentOutId);
+        SOFA_UNUSED(tangentInId);
+        SOFA_UNUSED(positionInId);
     }
 
     void applyJT(const core::MechanicalParams* /* mparams */, core::MultiVecDerivId /* inForce */, core::ConstMultiVecDerivId /* outForce */) override
diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h
index 847c8628cdf..759fc69bcbc 100644
--- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h
+++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h
@@ -47,7 +47,11 @@ class BaseNonLinearMapping : public core::Mapping<TIn, TOut>, public NonLinearMa
 
     void init() override;
 
-    void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& out, const DataVecDeriv_t<In>& in) final;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) final;
     void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t<In>& out, const DataVecDeriv_t<Out>& in) final;
     void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t<In>& out, const DataMatrixDeriv_t<Out>& in) final;
     void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) final;
diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl
index 742d3303e87..c0fe99674e7 100644
--- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl
+++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl
@@ -38,15 +38,17 @@ init()
 }
 
 template <class TIn, class TOut, bool HasStabilizedGeometricStiffness>
-void BaseNonLinearMapping<TIn, TOut, HasStabilizedGeometricStiffness>::applyJ(
-    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& out,
-    const DataVecDeriv_t<In>& in)
+void BaseNonLinearMapping<TIn, TOut, HasStabilizedGeometricStiffness>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams,
+    DataVecDeriv_t<Out>& tangentOut,
+    const DataVecDeriv_t<In>& tangentIn,
+    const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED( mparams );
     if( m_jacobian.rowSize() )
     {
-        auto dOutWa = sofa::helper::getWriteOnlyAccessor(out);
-        auto dInRa = sofa::helper::getReadAccessor(in);
+        auto dOutWa = sofa::helper::getWriteOnlyAccessor(tangentOut);
+        auto dInRa = sofa::helper::getReadAccessor(tangentIn);
         m_jacobian.mult(dOutWa.wref(),dInRa.ref());
     }
 }
diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMultiMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMultiMapping.h
index 763b6a752b7..7208287e3e1 100644
--- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMultiMapping.h
+++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMultiMapping.h
@@ -114,24 +114,29 @@ class DistanceMultiMapping : public core::MultiMapping<TIn, TOut>, public Stabil
 
     }
 
-    void applyJ(const core::MechanicalParams *mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel) override
+    virtual void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+        const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn) override
     {
         SOFA_UNUSED(mparams);
+        SOFA_UNUSED(dataVecPositionIn);
 
         //Not optimized at all...
         type::vector<OutVecDeriv*> vecOutVel;
-        for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-            vecOutVel.push_back(dataVecOutVel[i]->beginEdit());
+        for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+            vecOutVel.push_back(dataVecTangentOut[i]->beginEdit());
 
         type::vector<const InVecDeriv*> vecInVel;
-        for(unsigned int i=0; i<dataVecInVel.size(); i++)
-            vecInVel.push_back(&dataVecInVel[i]->getValue());
+        for(unsigned int i=0; i<dataVecTangentIn.size(); i++)
+            vecInVel.push_back(&dataVecTangentIn[i]->getValue());
 
         this->applyJ(vecOutVel, vecInVel);
 
         //Really Not optimized at all...
-        for(unsigned int i=0; i<dataVecOutVel.size(); i++)
-            dataVecOutVel[i]->endEdit();
+        for(unsigned int i=0; i<dataVecTangentOut.size(); i++)
+            dataVecTangentOut[i]->endEdit();
 
     }
 
@@ -157,7 +162,6 @@ class DistanceMultiMapping : public core::MultiMapping<TIn, TOut>, public Stabil
     }
 
     using Inherit::apply;
-    using Inherit::applyJ;
     using Inherit::applyJT;
 
     virtual void apply(const type::vector<OutVecCoord*>& outPos, const vecConstInVecCoord& inPos);
diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.h
index 25a384772fb..fc9a7c2bf5f 100644
--- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.h
+++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.h
@@ -116,7 +116,11 @@ class RigidMapping : public core::Mapping<TIn, TOut>, public NonLinearMappingDat
 
     void apply(const core::MechanicalParams *mparams, Data<OutVecCoord>& out, const Data<InVecCoord>& in) override;
 
-    void applyJ(const core::MechanicalParams *mparams, Data<OutVecDeriv>& out, const Data<InVecDeriv>& in) override;
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        OutDataVecDeriv& tangentOut,
+        const InDataVecDeriv& tangentIn,
+        const InDataVecCoord& positionIn) override;
 
     void applyJT(const core::MechanicalParams *mparams, Data<InVecDeriv>& out, const Data<OutVecDeriv>& in) override;
 
diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.inl
index 6f6b9829a31..dd5e0c04306 100644
--- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.inl
+++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/RigidMapping.inl
@@ -327,10 +327,17 @@ void RigidMapping<TIn, TOut>::apply(const core::MechanicalParams * /*mparams*/,
 }
 
 template <class TIn, class TOut>
-void RigidMapping<TIn, TOut>::applyJ(const core::MechanicalParams * /*mparams*/, Data<OutVecDeriv>& dOut, const Data<InVecDeriv>& dIn)
+void RigidMapping<TIn, TOut>::doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        OutDataVecDeriv& tangentOut,
+        const InDataVecDeriv& tangentIn,
+        const InDataVecCoord& positionIn)
 {
-    helper::WriteOnlyAccessor< Data<OutVecDeriv> > out = dOut;
-    helper::ReadAccessor< Data<InVecDeriv> > in = dIn;
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    helper::WriteOnlyAccessor< Data<OutVecDeriv> > out = tangentOut;
+    helper::ReadAccessor< Data<InVecDeriv> > in = tangentIn;
 
     const OutVecCoord& pts = this->getPoints();
     out.resize(pts.size());
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 87a9f23eaca..8f9459a6bf2 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
@@ -400,7 +400,7 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
         /// Updated to parentNew
         positionAccessorIn.wref() = parentNew;
         mapping->apply(&mparams, core::vec_id::write_access::position, core::vec_id::write_access::position);
-        mapping->applyJ(&mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity);
+        mapping->applyJacobianVectorProduct(&mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity, core::vec_id::write_access::position);
 
         bool succeed = true;
 
@@ -447,7 +447,7 @@ struct Mapping_test: public BaseSimulationTest, NumericTest<typename _Mapping::I
     void computeVelocityOutFromVelocityIn(core::MechanicalParams mparams, VecDeriv_t<Out>& velocityOut, const VecDeriv_t<In>& velocityIn)
     {
         inDofs->writeVelocities().wref() = velocityIn;
-        mapping->applyJ( &mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity );
+        mapping->applyJacobianVectorProduct( &mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity, core::vec_id::read_access::position );
         velocityOut = outDofs->readVelocities().ref();
     }
 
diff --git a/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MultiMappingTestCreation.h b/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MultiMappingTestCreation.h
index 6b1cbd357e8..f57f9681641 100644
--- a/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MultiMappingTestCreation.h
+++ b/Sofa/Component/Mapping/Testing/src/sofa/component/mapping/testing/MultiMappingTestCreation.h
@@ -177,7 +177,7 @@ struct MultiMapping_test : public BaseSimulationTest, NumericTest<typename _Mult
 
         /// apply the mapping
         mapping->apply(&mparams, core::vec_id::write_access::position, core::vec_id::write_access::position);
-        mapping->applyJ(&mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity);
+        mapping->applyJacobianVectorProduct( &mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity, core::vec_id::read_access::position );
 
         /// test apply: check if the child positions are the expected ones
         bool succeed=true;
@@ -249,7 +249,7 @@ struct MultiMapping_test : public BaseSimulationTest, NumericTest<typename _Mult
             sofa::testing::copyToData( vin, vp[p] );
         }
         mparams.setDx(core::vec_id::read_access::velocity);
-        mapping->applyJ( &mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity );
+        mapping->applyJacobianVectorProduct( &mparams, core::vec_id::write_access::velocity, core::vec_id::write_access::velocity, core::vec_id::read_access::position );
         ReadOutVecDeriv vout = outDofs->readVelocities();
         sofa::testing::copyFromData( vc, vout);
         //        cout<<"child velocity vc = " << vc << endl;
diff --git a/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h b/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h
index 73359ad108d..9a54f3c4ee7 100644
--- a/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h
+++ b/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h
@@ -94,8 +94,11 @@ class HexahedronCompositeFEMMapping : public BasicMapping
     void apply( const sofa::core::MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in) override;
     //void apply( typename Out::VecCoord& out, const typename In::VecCoord& in );
 
-    void applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) override;
-    //void applyJ( typename Out::VecDeriv& out, const typename In::VecDeriv& in );
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        DataVecDeriv_t<Out>& tangentOut,
+        const DataVecDeriv_t<In>& tangentIn,
+        const DataVecCoord_t<In>& positionIn) override;
 
     void applyJT( const sofa::core::MechanicalParams* mparams, InDataVecDeriv& out, const OutDataVecDeriv& in) override;
     //void applyJT( typename In::VecDeriv& out, const typename Out::VecDeriv& in );
diff --git a/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.inl b/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.inl
index 324b53b2281..13367319071 100644
--- a/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.inl
+++ b/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.inl
@@ -258,15 +258,16 @@ void HexahedronCompositeFEMMapping<BasicMapping>::apply( const sofa::core::Mecha
 
 
 }
-
-
 template <class BasicMapping>
-void HexahedronCompositeFEMMapping<BasicMapping>::applyJ( const sofa::core::MechanicalParams* mparams, OutDataVecDeriv& outData, const InDataVecDeriv& inData)
+void HexahedronCompositeFEMMapping<BasicMapping>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut,
+    const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
     SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
 
-    OutVecDeriv& out = *outData.beginEdit();
-    const InVecDeriv& in = inData.getValue();
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
 
     // les deplacements des noeuds grossiers
     type::vector< sofa::type::Vec< 24, SReal >  > coarseDisplacements( _sparseGrid->getNbHexahedra() );
@@ -313,7 +314,7 @@ void HexahedronCompositeFEMMapping<BasicMapping>::applyJ( const sofa::core::Mech
             out[i] += (fineDisplacements[ finehexa[w] ]  * _finestBarycentricCoord[i].second[w] );
         }
     }
-    outData.endEdit();
+    tangentOut.endEdit();
 }
 
 
diff --git a/Sofa/framework/Core/src/sofa/core/BaseMapping.cpp b/Sofa/framework/Core/src/sofa/core/BaseMapping.cpp
index b70661fdcec..b9933e02217 100644
--- a/Sofa/framework/Core/src/sofa/core/BaseMapping.cpp
+++ b/Sofa/framework/Core/src/sofa/core/BaseMapping.cpp
@@ -19,12 +19,14 @@
 *                                                                             *
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
+#include <sofa/core/BaseLocalMappingMatrix.h>
 #include <sofa/core/BaseMapping.h>
-#include <sofa/core/objectmodel/BaseNode.h>
+#include <sofa/core/BaseMatrixAccumulatorComponent.h>
 #include <sofa/core/BaseState.h>
-#include <sofa/core/BaseLocalMappingMatrix.h>
+#include <sofa/core/objectmodel/BaseNode.h>
 #include <sofa/linearalgebra/BaseMatrix.h>
-#include <sofa/core/BaseMatrixAccumulatorComponent.h>
+
+#include "MechanicalParams.h"
 
 namespace sofa::core
 {
@@ -44,6 +46,13 @@ BaseMapping::BaseMapping()
 BaseMapping::~BaseMapping()
 {}
 
+void BaseMapping::applyJ(
+    const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel)
+{
+    this->applyJacobianVectorProduct(
+        mparams, outVel, inVel, mparams ? mparams->x() : sofa::core::vec_id::read_access::position);
+}
+
 bool BaseMapping::setFrom(BaseState*  )
 {
     dmsg_error() << "Calling a virtual method not implemented." ;
diff --git a/Sofa/framework/Core/src/sofa/core/BaseMapping.h b/Sofa/framework/Core/src/sofa/core/BaseMapping.h
index 0a2e897bb36..2494c26f49c 100644
--- a/Sofa/framework/Core/src/sofa/core/BaseMapping.h
+++ b/Sofa/framework/Core/src/sofa/core/BaseMapping.h
@@ -20,6 +20,7 @@
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
 #pragma once
+#include <sofa/core/config.h>
 #include <sofa/core/behavior/MechanicalState.h>
 
 namespace sofa::core
@@ -65,8 +66,33 @@ class SOFA_CORE_API BaseMapping : public virtual objectmodel::BaseComponent
 
     /// Apply the transformation from the input model to the output model (like apply displacement from BehaviorModel to VisualModel)
     virtual void apply (const MechanicalParams* mparams = mechanicalparams::defaultInstance(), MultiVecCoordId outPos = vec_id::write_access::position, ConstMultiVecCoordId inPos = vec_id::read_access::position ) = 0;
+
     /// Compute output velocity based on input velocity, using the linearized transformation (tangent operator). Also used to propagate small displacements.
-    virtual void applyJ(const MechanicalParams* mparams = mechanicalparams::defaultInstance(), MultiVecDerivId outVel = vec_id::write_access::velocity, ConstMultiVecDerivId inVel = vec_id::read_access::velocity ) = 0;
+    SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
+    virtual void applyJ(const MechanicalParams* mparams = mechanicalparams::defaultInstance(), MultiVecDerivId outVel = vec_id::write_access::velocity, ConstMultiVecDerivId inVel = vec_id::read_access::velocity ) final;
+
+    /**
+     * @brief Compute the product of the Jacobian matrix with a vector (JVP).
+     *
+     * This function computes $v_{out} = J(x_{in}) \cdot v_{in}$. It propagates velocities, small
+     * displacements, or generalized tangents through non-linear mappings.
+     * The Jacobian matrix $J$ maps derivatives (velocities) from the input space to the output space.
+     * Since the Jacobian often depends on position, this method requires access to the
+     * coordinates of the input model ($x_{in}$).
+     *
+     * @param mparams The current simulation parameters.
+     * @param tangentOutId The identifier for the output velocity/tangent vector ($\dot{q}_{out}$).
+     *                     This is where the result must be written.
+     * @param tangentInId The identifier for the input velocity/tangent vector ($\dot{q}_{in}$),
+     *                    which serves as the multiplication operand.
+     * @param positionInId The identifier for the coordinates of the input model ($x_{in}$).
+     *                     The Jacobian $J$ must be evaluated at this specific configuration.
+     */
+    virtual void applyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        MultiVecDerivId tangentOutId,
+        ConstMultiVecDerivId tangentInId,
+        ConstMultiVecCoordId positionInId) = 0;
 
     /// Accessor to the input model of this mapping
     virtual type::vector<BaseState*> getFrom() = 0;
diff --git a/Sofa/framework/Core/src/sofa/core/Mapping.h b/Sofa/framework/Core/src/sofa/core/Mapping.h
index ef555eb9da0..54ef0d3e949 100644
--- a/Sofa/framework/Core/src/sofa/core/Mapping.h
+++ b/Sofa/framework/Core/src/sofa/core/Mapping.h
@@ -110,14 +110,20 @@ class Mapping : public BaseMapping
     /// This method must be reimplemented by all mappings.
     virtual void apply( const MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in)= 0;
 
-    /// ApplyJ ///
-    /// Apply the mapping to derived (velocity, displacement) vectors.
-    /// $ out = J in $
-    /// where J is the tangent operator (the linear approximation) of the mapping
-    void applyJ(const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel ) override;
-
-    /// This method must be reimplemented by all mappings.
-    virtual void applyJ( const MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) = 0;
+    void applyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        MultiVecDerivId tangentOutId,
+        ConstMultiVecDerivId tangentInId,
+        ConstMultiVecCoordId positionInId) final;
+
+    virtual void doApplyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        OutDataVecDeriv& tangentOut,
+        const InDataVecDeriv& tangentIn,
+        const InDataVecCoord& positionIn) = 0;
+
+    SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
+    virtual void applyJ( const MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) final {}
 
     /// ApplyJT (Force)///
     /// Apply the reverse mapping to force vectors.
diff --git a/Sofa/framework/Core/src/sofa/core/Mapping.inl b/Sofa/framework/Core/src/sofa/core/Mapping.inl
index 5639eda600f..3bb0153fdd2 100644
--- a/Sofa/framework/Core/src/sofa/core/Mapping.inl
+++ b/Sofa/framework/Core/src/sofa/core/Mapping.inl
@@ -24,6 +24,7 @@
 #include <sofa/core/State.h>
 #include <sofa/core/behavior/BaseMechanicalState.h>
 #include <sofa/core/Mapping.h>
+#include <sofa/core/MechanicalParams.h>
 #include <iostream>
 
 namespace sofa::core
@@ -106,7 +107,7 @@ void Mapping<In,Out>::init()
     }
 
     apply(mechanicalparams::defaultInstance(), vec_id::write_access::position, vec_id::read_access::position);
-    applyJ(mechanicalparams::defaultInstance(), vec_id::write_access::velocity, vec_id::read_access::velocity);
+    applyJacobianVectorProduct(mechanicalparams::defaultInstance(), vec_id::write_access::velocity, vec_id::read_access::velocity, vec_id::read_access::position);
     if (f_applyRestPosition.getValue())
         apply(mechanicalparams::defaultInstance(), vec_id::write_access::restPosition, vec_id::read_access::restPosition);
 }
@@ -146,23 +147,28 @@ void Mapping<In,Out>::apply(const MechanicalParams* mparams, MultiVecCoordId out
             this->apply(mparams, *out, *in);
         }
     }
-}// Mapping::apply
+}  // Mapping::apply
 
-template <class In, class Out>
-void Mapping<In,Out>::applyJ(const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel)
+template <class TIn, class TOut>
+void Mapping<TIn, TOut>::applyJacobianVectorProduct(const MechanicalParams* mparams,
+                                                    MultiVecDerivId tangentOutId,
+                                                    ConstMultiVecDerivId tangentInId,
+                                                    ConstMultiVecCoordId positionInId)
 {
     State<In>* from = this->fromModel.get();
-    State<Out>*  to = this->toModel.get();
-    if(from && to)
+    State<Out>* to = this->toModel.get();
+    if (from && to)
     {
-        OutDataVecDeriv* out = outVel[to].write();
-        const InDataVecDeriv* in = inVel[from].read();
-        if(out && in)
+        OutDataVecDeriv* tangentOut = tangentOutId[to].write();
+        const InDataVecDeriv* tangentIn = tangentInId[from].read();
+        const InDataVecCoord* positionIn = positionInId[from].read();
+
+        if (tangentOut && tangentIn && positionIn)
         {
-                this->applyJ(mparams, *out, *in);
+            this->doApplyJacobianVectorProduct(mparams, *tangentOut, *tangentIn, *positionIn);
         }
     }
-}// Mapping::applyJ
+}
 
 template <class In, class Out>
 void Mapping<In,Out>::applyJT(const MechanicalParams *mparams, MultiVecDerivId inForce, ConstMultiVecDerivId outForce)
diff --git a/Sofa/framework/Core/src/sofa/core/Multi2Mapping.h b/Sofa/framework/Core/src/sofa/core/Multi2Mapping.h
index c0ead569710..e85b723bf0f 100644
--- a/Sofa/framework/Core/src/sofa/core/Multi2Mapping.h
+++ b/Sofa/framework/Core/src/sofa/core/Multi2Mapping.h
@@ -108,21 +108,30 @@ class Multi2Mapping : public BaseMapping
         const type::vector<const DataVecCoord_t<In1>*>& dataVecIn1Pos ,
         const type::vector<const DataVecCoord_t<In2>*>& dataVecIn2Pos) = 0;
 
-    /// ApplyJ ///
-    /// This method computes
-    /// $ out = J in $
-    /// where J is the tangent operator (the linear approximation) of the mapping
-    void applyJ (const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel ) override;
+    void applyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        MultiVecDerivId tangentOutId,
+        ConstMultiVecDerivId tangentInId,
+        ConstMultiVecCoordId positionInId) final;
+
+    virtual void doApplyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In1>*>& dataVecTangentIn1,
+        const type::vector<const DataVecDeriv_t<In2>*>& dataVecTangentIn2,
+        const type::vector<const DataVecCoord_t<In1>*>& dataVecPositionIn1,
+        const type::vector<const DataVecCoord_t<In2>*>& dataVecPositionIn2) = 0;
 
     /// This method must be reimplemented by all mappings.
     /// InDeriv and OutDeriv by default contains VecIds of type V_DERIV.
     /// The size of InDeriv vector is the same as the number of fromModels.
     /// The size of OutDeriv vector is the same as the number of OutModels.
+    SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
     virtual void applyJ(
         const MechanicalParams*,
         const type::vector< DataVecDeriv_t<Out>*>& dataVecOutVel,
         const type::vector<const DataVecDeriv_t<In1>*>& dataVecIn1Vel,
-        const type::vector<const DataVecDeriv_t<In2>*>& dataVecIn2Vel)
+        const type::vector<const DataVecDeriv_t<In2>*>& dataVecIn2Vel) final
     {
         //Not optimized at all...
         type::vector<VecDeriv_t<Out>*> vecOutVel;
@@ -141,11 +150,13 @@ class Multi2Mapping : public BaseMapping
         for(unsigned int i=0; i<dataVecOutVel.size(); i++)
             dataVecOutVel[i]->endEdit();
     }
+
     /// Compat Method
     /// @deprecated
+    SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
     virtual void applyJ(const type::vector< VecDeriv_t<Out>*>& /* outDeriv */,
             const type::vector<const VecDeriv_t<In1>*>& /* inDeriv1 */,
-            const type::vector<const VecDeriv_t<In2>*>& /* inDeriv2 */) {}
+            const type::vector<const VecDeriv_t<In2>*>& /* inDeriv2 */) final {}
 
     /// ApplyJT (Force)///
     /// Apply the mapping to Force vectors.
diff --git a/Sofa/framework/Core/src/sofa/core/Multi2Mapping.inl b/Sofa/framework/Core/src/sofa/core/Multi2Mapping.inl
index e95a1b2f195..f381830cc00 100644
--- a/Sofa/framework/Core/src/sofa/core/Multi2Mapping.inl
+++ b/Sofa/framework/Core/src/sofa/core/Multi2Mapping.inl
@@ -151,16 +151,28 @@ void Multi2Mapping<In1,In2,Out>::apply (const MechanicalParams* mparams, MultiVe
     this->apply(mparams, vecOutPos, vecIn1Pos, vecIn2Pos);
 }
 
-template < class In1, class In2,class Out>
-void Multi2Mapping<In1,In2,Out>::applyJ (const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel )
+template <class TIn1, class TIn2, class TOut>
+void Multi2Mapping<TIn1, TIn2, TOut>::applyJacobianVectorProduct(const MechanicalParams* mparams,
+                                                                 MultiVecDerivId tangentOutId,
+                                                                 ConstMultiVecDerivId tangentInId,
+                                                                 ConstMultiVecCoordId positionInId)
 {
-    type::vector<DataVecDeriv_t<Out>*> vecOutVel;
-    getVecOutDeriv(outVel, vecOutVel);
-    type::vector<const DataVecDeriv_t<In1>*> vecIn1Vel;
-    getConstVecIn1Deriv(inVel, vecIn1Vel);
-    type::vector<const DataVecDeriv_t<In2>*> vecIn2Vel;
-    getConstVecIn2Deriv(inVel, vecIn2Vel);
-    this->applyJ(mparams, vecOutVel, vecIn1Vel, vecIn2Vel);
+    type::vector<DataVecDeriv_t<Out>*> vecTangentOut;
+    getVecOutDeriv(tangentOutId, vecTangentOut);
+
+    type::vector<const DataVecDeriv_t<In1>*> vecTangentIn1;
+    getConstVecIn1Deriv(tangentInId, vecTangentIn1);
+
+    type::vector<const DataVecDeriv_t<In2>*> vecTangentIn2;
+    getConstVecIn2Deriv(tangentInId, vecTangentIn2);
+
+    type::vector<const DataVecCoord_t<In1>*> vecPositionIn1;
+    getConstVecIn1Coord(positionInId, vecPositionIn1);
+
+    type::vector<const DataVecCoord_t<In2>*> vecPositionIn2;
+    getConstVecIn2Coord(positionInId, vecPositionIn2);
+
+    this->doApplyJacobianVectorProduct(mparams, vecTangentOut, vecTangentIn1, vecTangentIn2, vecPositionIn1, vecPositionIn2);
 }
 
 template < class In1, class In2,class Out>
@@ -222,7 +234,7 @@ void Multi2Mapping<In1, In2, Out>::init()
     }
 
     apply(mechanicalparams::defaultInstance() , vec_id::write_access::position, vec_id::read_access::position);
-    applyJ(mechanicalparams::defaultInstance() , vec_id::write_access::velocity, vec_id::read_access::velocity);
+    applyJacobianVectorProduct(mechanicalparams::defaultInstance(), vec_id::write_access::velocity, vec_id::read_access::velocity, vec_id::read_access::position);
     if (f_applyRestPosition.getValue())
         apply(mechanicalparams::defaultInstance(), vec_id::write_access::restPosition, vec_id::read_access::restPosition);
 }
diff --git a/Sofa/framework/Core/src/sofa/core/MultiMapping.h b/Sofa/framework/Core/src/sofa/core/MultiMapping.h
index bc65f24d32f..0ba60d2c247 100644
--- a/Sofa/framework/Core/src/sofa/core/MultiMapping.h
+++ b/Sofa/framework/Core/src/sofa/core/MultiMapping.h
@@ -113,17 +113,20 @@ class MultiMapping : public BaseMapping
     /// The size of OutPos vector is the same as the number of OutModels.
     virtual void apply(const MechanicalParams* mparams, const type::vector<OutDataVecCoord*>& dataVecOutPos, const type::vector<const InDataVecCoord*>& dataVecInPos) = 0;
 
-    /// ApplyJ ///
-    /// Apply the mapping to derived (velocity, displacement) vectors.
-    /// \f$ out = J in \f$
-    /// where \f$ J \f$ is the tangent operator (the linear approximation) of the mapping
-    void applyJ (const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel ) override;
-
-    /// This method must be reimplemented by all mappings.
-    /// InDeriv and OutDeriv by default contains VecIds of type V_DERIV.
-    /// The size of InDeriv vector is the same as the number of fromModels.
-    /// The size of OutDeriv vector is the same as the number of OutModels.
-    virtual void applyJ(const MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel) = 0;
+    void applyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        MultiVecDerivId tangentOutId,
+        ConstMultiVecDerivId tangentInId,
+        ConstMultiVecCoordId positionInId) final;
+
+    virtual void doApplyJacobianVectorProduct(
+        const MechanicalParams* mparams,
+        const type::vector< DataVecDeriv_t<Out>*>& dataVecTangentOut,
+        const type::vector<const DataVecDeriv_t<In>*>& dataVecTangentIn,
+        const type::vector<const DataVecCoord_t<In>*>& dataVecPositionIn) = 0;
+
+    SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
+    virtual void applyJ(const MechanicalParams* mparams, const type::vector<OutDataVecDeriv*>& dataVecOutVel, const type::vector<const InDataVecDeriv*>& dataVecInVel) final {}
 
     /// ApplyJT (Force)///
     /// Apply the reverse mapping to force vectors.
diff --git a/Sofa/framework/Core/src/sofa/core/MultiMapping.inl b/Sofa/framework/Core/src/sofa/core/MultiMapping.inl
index b0c460635ba..58951254d0b 100644
--- a/Sofa/framework/Core/src/sofa/core/MultiMapping.inl
+++ b/Sofa/framework/Core/src/sofa/core/MultiMapping.inl
@@ -129,7 +129,7 @@ void MultiMapping<In,Out>::init()
     }
 
     apply(mechanicalparams::defaultInstance() , vec_id::write_access::position, vec_id::read_access::position);
-    applyJ(mechanicalparams::defaultInstance() , vec_id::write_access::velocity, vec_id::read_access::velocity);
+    applyJacobianVectorProduct(mechanicalparams::defaultInstance(), vec_id::write_access::velocity, vec_id::read_access::velocity, vec_id::read_access::position);
     if (f_applyRestPosition.getValue())
         apply(mechanicalparams::defaultInstance(), vec_id::write_access::restPosition, vec_id::read_access::restPosition);
 }
@@ -142,17 +142,25 @@ void MultiMapping<In,Out>::apply(const MechanicalParams* mparams, MultiVecCoordI
     type::vector<const InDataVecCoord*> vecInPos;
     getConstVecInCoord(inPos, vecInPos);
     this->apply(mparams, vecOutPos, vecInPos);
-}// MultiMapping::apply
+}  // MultiMapping::apply
 
-template <class In, class Out>
-void MultiMapping<In,Out>::applyJ(const MechanicalParams* mparams, MultiVecDerivId outVel, ConstMultiVecDerivId inVel)
+template <class TIn, class TOut>
+void MultiMapping<TIn, TOut>::applyJacobianVectorProduct(const MechanicalParams* mparams,
+                                                         MultiVecDerivId tangentOutId,
+                                                         ConstMultiVecDerivId tangentInId,
+                                                         ConstMultiVecCoordId positionInId)
 {
-    type::vector<OutDataVecDeriv*> vecOutVel;
-    getVecOutDeriv(outVel, vecOutVel);
-    type::vector<const InDataVecDeriv*> vecInVel;
-    getConstVecInDeriv(inVel, vecInVel);
-    this->applyJ(mparams, vecOutVel, vecInVel);
-}// MultiMapping::applyJ
+    type::vector<OutDataVecDeriv*> vecTangentOut;
+    getVecOutDeriv(tangentOutId, vecTangentOut);
+
+    type::vector<const InDataVecDeriv*> vecTangentIn;
+    getConstVecInDeriv(tangentInId, vecTangentIn);
+
+    type::vector<const InDataVecCoord*> vecPositionIn;
+    getConstVecInCoord(positionInId, vecPositionIn);
+
+    this->doApplyJacobianVectorProduct(mparams, vecTangentOut, vecTangentIn, vecPositionIn);
+}
 
 template <class In, class Out>
 void MultiMapping<In,Out>::applyJT(const MechanicalParams* mparams, MultiVecDerivId inForce, ConstMultiVecDerivId outForce)
diff --git a/Sofa/framework/Core/src/sofa/core/config.h.in b/Sofa/framework/Core/src/sofa/core/config.h.in
index 3ed70d58b2d..3b66d7c206d 100644
--- a/Sofa/framework/Core/src/sofa/core/config.h.in
+++ b/Sofa/framework/Core/src/sofa/core/config.h.in
@@ -132,3 +132,10 @@ SOFA_ATTRIBUTE_DEPRECATED("v26.06", "v29.06", "Use toBaseComponent instead.")
 #define SOFA_CORE_DEPRECATED_REMOVE_CONTEXT() \
 SOFA_ATTRIBUTE_DISABLED("v26.12", "v27.12", "Use BaseContext instead of Context.")
 #endif
+
+#ifdef SOFA_BUILD_SOFA_CORE
+#define SOFA_ATTRIBUTE_DEPRECATED__APPLYJ()
+#else
+#define SOFA_ATTRIBUTE_DEPRECATED__APPLYJ() \
+    SOFA_ATTRIBUTE_DEPRECATED("v26.12", "v27.06", "Use applyJacobianVectorProduct instead.")
+#endif
diff --git a/Sofa/framework/Simulation/Core/test/DifferentialOperations_test.cpp b/Sofa/framework/Simulation/Core/test/DifferentialOperations_test.cpp
index 9b09e41689a..3c5e319f27d 100644
--- a/Sofa/framework/Simulation/Core/test/DifferentialOperations_test.cpp
+++ b/Sofa/framework/Simulation/Core/test/DifferentialOperations_test.cpp
@@ -65,10 +65,17 @@ class ScaleMapping : public core::Mapping<DataTypes1D, DataTypes1D>
         }
     }
 
-    void applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn) override
+    void doApplyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        OutDataVecDeriv& tangentOut,
+        const InDataVecDeriv& tangentIn,
+        const InDataVecCoord& positionIn) override
     {
-        const auto in = sofa::helper::getReadAccessor(dIn);
-        auto out = sofa::helper::getWriteOnlyAccessor(dOut);
+        SOFA_UNUSED(mparams);
+        SOFA_UNUSED(positionIn);
+
+        const auto in = sofa::helper::getReadAccessor(tangentIn);
+        auto out = sofa::helper::getWriteOnlyAccessor(tangentOut);
         out.resize(in.size());
 
         for (size_t i = 0; i < in.size(); ++i)
diff --git a/Sofa/framework/Simulation/Core/test/Visitor_test.cpp b/Sofa/framework/Simulation/Core/test/Visitor_test.cpp
index 3dd42371924..d64788df2a4 100644
--- a/Sofa/framework/Simulation/Core/test/Visitor_test.cpp
+++ b/Sofa/framework/Simulation/Core/test/Visitor_test.cpp
@@ -202,11 +202,16 @@ class TestMapping : public core::BaseMapping
         SOFA_UNUSED(outPos);
         SOFA_UNUSED(inPos);
     }
-    void applyJ(const core::MechanicalParams* mparams, core::MultiVecDerivId outVel, core::ConstMultiVecDerivId inVel) override
+    void applyJacobianVectorProduct(
+        const core::MechanicalParams* mparams,
+        core::MultiVecDerivId tangentOutId,
+        core::ConstMultiVecDerivId tangentInId,
+        core::ConstMultiVecCoordId positionInId) override
     {
         SOFA_UNUSED(mparams);
-        SOFA_UNUSED(outVel);
-        SOFA_UNUSED(inVel);
+        SOFA_UNUSED(tangentOutId);
+        SOFA_UNUSED(tangentInId);
+        SOFA_UNUSED(positionInId);
     }
     type::vector<core::BaseState*> getFrom() override { return {} ;}
     type::vector<core::BaseState*> getTo() override { return {} ;}
diff --git a/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.h b/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.h
index cb9d55336e4..84000a3ec2c 100644
--- a/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.h
+++ b/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.h
@@ -67,7 +67,11 @@ class BeamLinearMapping_mt : public sofa::component::mapping::linear::BeamLinear
 
     virtual void apply(const sofa::core::MechanicalParams *mparams /* PARAMS FIRST */, sofa::Data< typename Out::VecCoord >& out, const sofa::Data< typename In::VecCoord >& in) override;
 
-    virtual void applyJ(const sofa::core::MechanicalParams *mparams /* PARAMS FIRST */, sofa::Data< typename Out::VecDeriv >& out, const sofa::Data< typename In::VecDeriv >& in) override;
+    virtual void doApplyJacobianVectorProduct(
+                const sofa::core::MechanicalParams* mparams,
+                sofa::DataVecDeriv_t<Out>& tangentOut,
+                const sofa::DataVecDeriv_t<In>& tangentIn,
+                const sofa::DataVecCoord_t<In>& positionIn) override;
 
     virtual void applyJT(const sofa::core::MechanicalParams *mparams /* PARAMS FIRST */, sofa::Data< typename In::VecDeriv >& out, const sofa::Data< typename Out::VecDeriv >& in) override;
 
diff --git a/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.inl b/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.inl
index 7cecb99ac36..eaa7ad6c62a 100644
--- a/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.inl
+++ b/applications/plugins/MultiThreading/src/MultiThreading/component/mapping/linear/BeamLinearMapping_mt.inl
@@ -155,21 +155,21 @@ namespace multithreading::component::mapping::linear
         //task_pool.purge_memory();
         
     }
-    
-    
-    
     template <class TIn, class TOut>
-    //void AdaptiveBeamMapping< TIn, TOut>::applyJ( typename Out::VecDeriv& out, const typename In::VecDeriv& in )
-    void BeamLinearMapping_mt< TIn, TOut>::applyJ(const sofa::core::MechanicalParams * params /* PARAMS FIRST */, sofa::Data< typename Out::VecDeriv >& _out, const sofa::Data< typename In::VecDeriv >& _in)
+    void BeamLinearMapping_mt<TIn, TOut>::doApplyJacobianVectorProduct(
+        const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t<Out>& tangentOut,
+        const sofa::DataVecDeriv_t<In>& tangentIn, const sofa::DataVecCoord_t<In>& positionIn)
     {
-        
+        SOFA_UNUSED(mparams);
+        SOFA_UNUSED(positionIn);
+
         //boost::pool<> task_pool(sizeof(BeamLinearMapping_mt< TIn, TOut>::applyJTask));
         unsigned int numPoints = this->points.size();
         
         if ( numPoints >  2*mGrainSize.getValue()  )
         {
-            sofa::helper::WriteAccessor< sofa::Data< typename Out::VecDeriv > > out = _out;
-            sofa::helper::ReadAccessor< sofa::Data< typename In::VecDeriv > > in = _in;
+            sofa::helper::WriteAccessor< sofa::Data< typename Out::VecDeriv > > out = tangentOut;
+            sofa::helper::ReadAccessor< sofa::Data< typename In::VecDeriv > > in = tangentIn;
             
             //const InVecDeriv& in= dIn.getValue();
             //VecDeriv& out = *dOut.beginEdit();
@@ -238,7 +238,7 @@ namespace multithreading::component::mapping::linear
         else
         {
             
-            sofa::component::mapping::linear::BeamLinearMapping<TIn,TOut>::applyJ( params, _out, _in );
+            sofa::component::mapping::linear::BeamLinearMapping<TIn,TOut>::doApplyJacobianVectorProduct( mparams, tangentOut, tangentIn, positionIn );
             
         }
         
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.h b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.h
index 0cf89ee8ebd..1fb806a4e75 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.h
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.h
@@ -31,7 +31,7 @@ template <>
 inline void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-inline void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+inline void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 inline void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
@@ -46,7 +46,7 @@ template <>
 inline void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-inline void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+inline void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 inline void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.inl b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.inl
index 4cd6a2f6198..0a81282b593 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.inl
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaIdentityMapping.inl
@@ -56,13 +56,17 @@ void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::appl
 }
 
 template <>
-void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void IdentityMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     out.fastResize(in.size());
     gpu::cuda::MechanicalObjectCudaVec3f_vAssign(out.size(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
@@ -117,13 +121,17 @@ void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::ap
 }
 
 template <>
-void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void IdentityMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     out.fastResize(in.size());
     gpu::cuda::MechanicalObjectCudaVec3f1_vAssign(out.size(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.h b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.h
index 683948d7c96..5490b5b2b7f 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.h
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.h
@@ -86,7 +86,7 @@ template <>
 void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
@@ -105,7 +105,7 @@ template <>
 void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
@@ -122,7 +122,7 @@ template <>
 void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
@@ -139,7 +139,7 @@ template <>
 void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.inl b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.inl
index 81f6ed32456..86fbd8b7935 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.inl
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/linear/CudaSubsetMapping.inl
@@ -81,14 +81,18 @@ void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::apply(
 }
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     const IndexArray& map = this->d_indices.getValue();
     out.fastResize(map.size());
     SubsetMappingCuda3f_applyJ(map.size(), map.deviceRead(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
@@ -131,14 +135,18 @@ void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::appl
 }
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     const IndexArray& map = this->d_indices.getValue();
     out.fastResize(map.size());
     SubsetMappingCuda3f1_applyJ(map.size(), map.deviceRead(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
@@ -179,14 +187,18 @@ void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::apply
 }
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SubsetMapping<gpu::cuda::CudaVec3f1Types, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     const IndexArray& map = this->d_indices.getValue();
     out.fastResize(map.size());
     SubsetMappingCuda3f1_3f_applyJ(map.size(), map.deviceRead(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
@@ -227,14 +239,18 @@ void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::apply
 }
 
 template <>
-void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void SubsetMapping<gpu::cuda::CudaVec3fTypes, gpu::cuda::CudaVec3f1Types>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
     const IndexArray& map = this->d_indices.getValue();
     out.fastResize(map.size());
     SubsetMappingCuda3f_3f1_applyJ(map.size(), map.deviceRead(), out.deviceWrite(), in.deviceRead());
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.h b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.h
index e3caec37ef6..351880b2d41 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.h
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.h
@@ -39,7 +39,7 @@ template <>
 void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
@@ -56,7 +56,7 @@ template <>
 void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::apply( const core::MechanicalParams* mparams, OutDataVecCoord& dOut, const InDataVecCoord& dIn );
 
 template <>
-void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::applyJ( const core::MechanicalParams* mparams, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn );
+void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::doApplyJacobianVectorProduct( const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn);
 
 template <>
 void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::applyJT( const core::MechanicalParams* mparams, InDataVecDeriv& dOut, const OutDataVecDeriv& dIn );
diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.inl b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.inl
index 4c312da1758..896fd586fda 100644
--- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.inl
+++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/mapping/nonlinear/CudaRigidMapping.inl
@@ -66,10 +66,14 @@ void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::apply
 }
 
 template <>
-void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
 
     const auto& points = this->d_points.getValue();
     gpu::cuda::CudaVec3fTypes::Deriv v, omega;
@@ -79,7 +83,7 @@ void RigidMapping<gpu::cuda::CudaRigid3fTypes, gpu::cuda::CudaVec3fTypes>::apply
 
     RigidMappingCuda3f_applyJ(points.size(), v, omega, out.deviceWrite(), m_rotatedPoints.deviceRead());
 
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>
@@ -128,10 +132,14 @@ void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::apply( co
 }
 
 template <>
-void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::applyJ( const core::MechanicalParams* /*mparams*/, OutDataVecDeriv& dOut, const InDataVecDeriv& dIn )
+void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::doApplyJacobianVectorProduct(
+    const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& tangentOut, const DataVecDeriv_t<In>& tangentIn, const DataVecCoord_t<In>& positionIn)
 {
-    OutVecDeriv& out = *dOut.beginEdit();
-    const InVecDeriv& in = dIn.getValue();
+    SOFA_UNUSED(mparams);
+    SOFA_UNUSED(positionIn);
+
+    OutVecDeriv& out = *tangentOut.beginEdit();
+    const InVecDeriv& in = tangentIn.getValue();
 
     const auto& points = this->d_points.getValue();
     gpu::cuda::CudaVec3fTypes::Deriv v, omega;
@@ -141,7 +149,7 @@ void RigidMapping<defaulttype::Rigid3Types, gpu::cuda::CudaVec3Types>::applyJ( c
 
     RigidMappingCuda3f_applyJ(points.size(), v, omega, out.deviceWrite(), m_rotatedPoints.deviceRead());
 
-    dOut.endEdit();
+    tangentOut.endEdit();
 }
 
 template <>