diff --git a/SofaKernel/modules/SofaBaseMechanics/SofaBaseMechanics_test/DiagonalMass_test.cpp b/SofaKernel/modules/SofaBaseMechanics/SofaBaseMechanics_test/DiagonalMass_test.cpp
index c1cb3078d4d..d916f2c9033 100644
--- a/SofaKernel/modules/SofaBaseMechanics/SofaBaseMechanics_test/DiagonalMass_test.cpp
+++ b/SofaKernel/modules/SofaBaseMechanics/SofaBaseMechanics_test/DiagonalMass_test.cpp
@@ -830,30 +830,57 @@ class DiagonalMass_test : public BaseTest
         ASSERT_NE(modifier, nullptr);
 
         const VecMass& vMasses = mass->d_vertexMass.getValue();
-        static const Real refValue = Real(1.0 / 4.0);  // 0.125
+        static const Real refValue = Real(1.0 / 4.0);  // 0.25
         static const Real initMass = mass->getTotalMass();
 
-        // check value at init
-        EXPECT_EQ(vMasses.size(), 9);
+        // check value at init. Grid of 4 quads, should be first and third line: 0.25, 0.5, 0.25, second line: 0.5, 1.0, 0.5
+        // 1/4 ---- 1/2 ---- 1/4
+        //  |   m:1  |   m:1  |
+        //  |  id:2  |  id:3  |
+        // 1/2 ----  1  ---- 1/2
+        //  |   m:1  |   m:1  |
+        //  |  id:0  |  id:1  |
+        // 1/4 ---- 1/2 ---- 1/4
+        EXPECT_EQ(vMasses.size(), 9);  
         EXPECT_NEAR(vMasses[0], refValue, 1e-4);
         EXPECT_NEAR(vMasses[1], refValue * 2, 1e-4);
+        EXPECT_NEAR(vMasses[4], refValue * 4, 1e-4);
         EXPECT_NEAR(initMass, 4, 1e-4);
-
-        sofa::type::vector<sofa::Index> ids = { 0 };
+        Real lastV = vMasses[8];
+       
+        
         // remove quad id: 0
+        // 1/4 ---- 1/2 ---- 1/4
+        //  |   m:1  |   m:1  |
+        //  |  id:2  |  id:0  |
+        // 1/4 ---- 3/4 ---- 1/2
+        //           |   m:1  |
+        //           |  id:1  |
+        //          1/4 ---- 1/4
+        sofa::type::vector<sofa::Index> ids = { 0 };
         modifier->removeQuads(ids, true, true);
         EXPECT_EQ(vMasses.size(), 8);
-        EXPECT_NEAR(vMasses[0], refValue, 1e-4); // check update of Mass when removing tetra
-        EXPECT_NEAR(vMasses[1], refValue * 2, 1e-4);
-        EXPECT_NEAR(mass->getTotalMass(), initMass - refValue, 1e-4);
-        const Real lastV = vMasses[7];
+        EXPECT_NEAR(vMasses[0], lastV, 1e-4); // check update of Mass when removing quad, vMasses[0] is now mapped to point position 8.
+        EXPECT_NEAR(vMasses[1], refValue, 1e-4); // one neighboord quad removed. 
+        EXPECT_NEAR(vMasses[4], refValue * 3, 1e-4); // one neighboord quad removed. 
+
+        EXPECT_NEAR(mass->getTotalMass(), initMass - 1, 1e-4);
+        lastV = vMasses[7];
 
         // remove quad id: 0
+        // 1/4 ---- 1/4
+        //  |   m:1  |
+        //  |  id:0  |
+        // 1/4 ---- 2/4 ---- 1/4
+        //           |   m:1  |
+        //           |  id:1  |
+        //          1/4 ---- 1/4
         modifier->removeQuads(ids, true, true);
         EXPECT_EQ(vMasses.size(), 7);
-        EXPECT_NEAR(vMasses[0], lastV, 1e-4); // check swap value
-        EXPECT_NEAR(vMasses[1], refValue * 2, 1e-4);
-        EXPECT_NEAR(mass->getTotalMass(), initMass - 2*refValue, 1e-4);
+        EXPECT_NEAR(vMasses[0], lastV - refValue, 1e-4); // check swap value
+        EXPECT_NEAR(vMasses[1], refValue, 1e-4);
+        EXPECT_NEAR(vMasses[4], refValue * 2, 1e-4); // one neighboord quad removed. 
+        EXPECT_NEAR(mass->getTotalMass(), initMass - 2, 1e-4);
 
         ids.push_back(1);
         // remove quad id: 0, 1
diff --git a/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.h b/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.h
index b23002c5aba..6f15cc08cac 100644
--- a/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.h
+++ b/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.h
@@ -82,6 +82,7 @@ class DiagonalMass : public core::behavior::Mass<DataTypes>
     typedef core::topology::BaseMeshTopology::Edge Edge;
     typedef core::topology::BaseMeshTopology::EdgeID EdgeID;
     typedef core::topology::BaseMeshTopology::Quad Quad;
+    typedef core::topology::BaseMeshTopology::QuadID QuadID;
     typedef core::topology::BaseMeshTopology::Triangle Triangle;
     typedef core::topology::BaseMeshTopology::TriangleID TriangleID;
     typedef core::topology::BaseMeshTopology::Tetrahedron Tetrahedron;
@@ -89,83 +90,6 @@ class DiagonalMass : public core::behavior::Mass<DataTypes>
     typedef core::topology::BaseMeshTopology::Hexahedron Hexahedron;
     typedef core::topology::BaseMeshTopology::HexahedronID HexahedronID;
 
-    class DMassPointEngine : public topology::TopologyDataHandler<Point,MassVector>
-    {
-    public:
-        typedef typename DiagonalMass<DataTypes,TMassType>::MassVector MassVector;
-        DMassPointEngine(DiagonalMass<DataTypes,TMassType>* _dm, sofa::component::topology::PointData<MassVector>* _data)
-            : topology::TopologyDataHandler<Point,MassVector>(_data), dm(_dm)
-        {}
-
-        void applyCreateFunction(PointID pointIndex, TMassType& m, const Point&, const sofa::type::vector< PointID > &,
-                                 const sofa::type::vector< double > &);
-
-        using topology::TopologyDataHandler<Point,MassVector>::ApplyTopologyChange;
-
-        ///////////////////////// Functions on Points //////////////////////////////////////
-        /// Apply removing points.
-        void applyPointDestruction(const sofa::type::vector<PointID> & /*indices*/);
-        /// Callback to remove points.
-        virtual void ApplyTopologyChange(const core::topology::PointsRemoved* /*event*/);
-
-        ///////////////////////// Functions on Edges //////////////////////////////////////
-        /// Apply adding edges elements.
-        void applyEdgeCreation(const sofa::type::vector< EdgeID >& /*indices*/,
-                               const sofa::type::vector< Edge >& /*elems*/,
-                               const sofa::type::vector< sofa::type::vector< EdgeID > >& /*ancestors*/,
-                               const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
-        /// Apply removing edges elements.
-        void applyEdgeDestruction(const sofa::type::vector<EdgeID> & /*indices*/);
-
-        /// Callback to add edges elements.
-        virtual void ApplyTopologyChange(const core::topology::EdgesAdded* /*event*/);
-        /// Callback to remove edges elements.
-        virtual void ApplyTopologyChange(const core::topology::EdgesRemoved* /*event*/);
-
-        ///////////////////////// Functions on Triangles //////////////////////////////////////
-        /// Apply adding triangles elements.
-        void applyTriangleCreation(const sofa::type::vector< TriangleID >& /*indices*/,
-                                   const sofa::type::vector< Triangle >& /*elems*/,
-                                   const sofa::type::vector< sofa::type::vector< TriangleID > >& /*ancestors*/,
-                                   const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
-        /// Apply removing triangles elements.
-        void applyTriangleDestruction(const sofa::type::vector<TriangleID> & /*indices*/);
-
-        /// Callback to add triangles elements.
-        virtual void ApplyTopologyChange(const core::topology::TrianglesAdded* /*event*/);
-        /// Callback to remove triangles elements.
-        virtual void ApplyTopologyChange(const core::topology::TrianglesRemoved* /*event*/);
-
-        ///////////////////////// Functions on Tetrahedron //////////////////////////////////////
-        /// Apply adding tetrahedron elements.
-        void applyTetrahedronCreation(const sofa::type::vector< TetrahedronID >& /*indices*/,
-                                      const sofa::type::vector< Tetrahedron >& /*elems*/,
-                                      const sofa::type::vector< sofa::type::vector< TetrahedronID > >& /*ancestors*/,
-                                      const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
-        /// Apply removing tetrahedron elements.
-        void applyTetrahedronDestruction(const sofa::type::vector<TetrahedronID> & /*indices*/);
-
-        /// Callback to add tetrahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::TetrahedraAdded* /*event*/);
-        /// Callback to remove tetrahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::TetrahedraRemoved* /*event*/);
-
-        ///////////////////////// Functions on Hexahedron //////////////////////////////////////
-        /// Apply adding hexahedron elements.
-        void applyHexahedronCreation(const sofa::type::vector< HexahedronID >& /*indices*/,
-                                     const sofa::type::vector< Hexahedron >& /*elems*/,
-                                     const sofa::type::vector< sofa::type::vector< HexahedronID > >& /*ancestors*/,
-                                     const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
-        /// Apply removing hexahedron elements.
-        void applyHexahedronDestruction(const sofa::type::vector<HexahedronID> & /*indices*/);
-        /// Callback to add hexahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::HexahedraAdded* /*event*/);
-        /// Callback to remove hexahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::HexahedraRemoved* /*event*/);
-
-    protected:
-        DiagonalMass<DataTypes,TMassType>* dm;
-    };
     /// the mass density used to compute the mass from a mesh topology and geometry
     Data< Real > d_massDensity;
 
@@ -181,8 +105,6 @@ class DiagonalMass : public core::behavior::Mass<DataTypes>
     Data< float > d_showAxisSize; ///< factor length of the axis displayed (only used for rigids)
     core::objectmodel::DataFileName d_fileMass; ///< an Xsp3.0 file to specify the mass parameters
 
-    DMassPointEngine* m_pointEngine;
-
     /// value defining the initialization process of the mass (0 : totalMass, 1 : massDensity, 2 : vertexMass)
     int m_initializationProcess;
 
@@ -247,6 +169,88 @@ class DiagonalMass : public core::behavior::Mass<DataTypes>
     /// Compute the vertexMass using input density and return the corresponding full mass.
     Real computeVertexMass(Real density);
 
+    /** Method to initialize @sa MassVector when a new Point is created.
+    * Will be set as creation callback in the PointData @sa d_vertexMass
+    */
+    void applyPointCreation(PointID pointIndex, MassType& m, const Point&,
+        const sofa::type::vector< PointID >&,
+        const sofa::type::vector< double >&);
+
+    /** Method to update @sa d_vertexMass when a Point is removed.
+    * Will be set as destruction callback in the PointData @sa d_vertexMass
+    */
+    void applyPointDestruction(Index id, MassType& VertexMass);
+
+
+    /** Method to update @sa d_vertexMass when a new Edge is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when EDGESADDED event is fired.    
+    */
+    void applyEdgeCreation(const sofa::type::vector< EdgeID >& /*indices*/,
+        const sofa::type::vector< Edge >& /*elems*/,
+        const sofa::type::vector< sofa::type::vector< EdgeID > >& /*ancestors*/,
+        const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
+
+    /** Method to update @sa d_vertexMass when a Edge is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when EDGESREMOVED event is fired.
+    */
+    void applyEdgeDestruction(const sofa::type::vector<EdgeID>& /*indices*/);
+
+
+    /** Method to update @sa d_vertexMass when a new Triangle is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TRIANGLESADDED event is fired.
+    */
+    void applyTriangleCreation(const sofa::type::vector< TriangleID >& /*indices*/,
+        const sofa::type::vector< Triangle >& /*elems*/,
+        const sofa::type::vector< sofa::type::vector< TriangleID > >& /*ancestors*/,
+        const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
+
+    /** Method to update @sa d_vertexMass when a Triangle is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TRIANGLESREMOVED event is fired.
+    */
+    void applyTriangleDestruction(const sofa::type::vector<TriangleID>& /*indices*/);
+
+
+    /** Method to update @sa d_vertexMass when a new Quad is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when QUADSADDED event is fired.
+    */
+    void applyQuadCreation(const sofa::type::vector< QuadID >& /*indices*/,
+        const sofa::type::vector< Quad >& /*elems*/,
+        const sofa::type::vector< sofa::type::vector< QuadID > >& /*ancestors*/,
+        const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
+
+    /** Method to update @sa d_vertexMass when a Quad is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when QUADSREMOVED event is fired.
+    */
+    void applyQuadDestruction(const sofa::type::vector<QuadID>& /*indices*/);
+    
+
+    /** Method to update @sa d_vertexMass when a new Tetrahedron is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TETRAHEDRAADDED event is fired.
+    */
+    void applyTetrahedronCreation(const sofa::type::vector< TetrahedronID >& /*indices*/,
+        const sofa::type::vector< Tetrahedron >& /*elems*/,
+        const sofa::type::vector< sofa::type::vector< TetrahedronID > >& /*ancestors*/,
+        const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
+
+    /** Method to update @sa d_vertexMass when a Tetrahedron is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TETRAHEDRAREMOVED event is fired.
+    */
+    void applyTetrahedronDestruction(const sofa::type::vector<TetrahedronID>& /*indices*/);
+
+
+    /** Method to update @sa d_vertexMass when a new Hexahedron is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when HEXAHEDRAADDED event is fired.
+    */
+    void applyHexahedronCreation(const sofa::type::vector< HexahedronID >& /*indices*/,
+        const sofa::type::vector< Hexahedron >& /*elems*/,
+        const sofa::type::vector< sofa::type::vector< HexahedronID > >& /*ancestors*/,
+        const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/);
+    
+    /** Method to update @sa d_vertexMass when a Hexahedron is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when HEXAHEDRAREMOVED event is fired.
+    */
+    void applyHexahedronDestruction(const sofa::type::vector<HexahedronID>& /*indices*/);
+
 public:
 
     SReal getTotalMass() const { return d_totalMass.getValue(); }
diff --git a/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.inl b/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.inl
index f93dbc8b6bb..cef5af2efb2 100644
--- a/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.inl
+++ b/SofaKernel/modules/SofaBaseMechanics/src/SofaBaseMechanics/DiagonalMass.inl
@@ -49,7 +49,6 @@ DiagonalMass<DataTypes, MassType>::DiagonalMass()
     , d_showCenterOfGravity( initData(&d_showCenterOfGravity, false, "showGravityCenter", "Display the center of gravity of the system" ) )
     , d_showAxisSize( initData(&d_showAxisSize, 1.0f, "showAxisSizeFactor", "Factor length of the axis displayed (only used for rigids)" ) )
     , d_fileMass( initData(&d_fileMass,  "filename", "Xsp3.0 file to specify the mass parameters" ) )
-    , m_pointEngine(nullptr)
     , l_topology(initLink("topology", "link to the topology container"))
     , m_massTopologyType(TopologyElementType::UNKNOWN)
     , m_topology(nullptr)
@@ -60,84 +59,50 @@ DiagonalMass<DataTypes, MassType>::DiagonalMass()
 template <class DataTypes, class MassType>
 DiagonalMass<DataTypes, MassType>::~DiagonalMass()
 {
-    if (m_pointEngine)
-        delete m_pointEngine;
+
 }
 
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyCreateFunction(PointID, MassType &m, const Point &, const sofa::type::vector<PointID> &, const sofa::type::vector<double> &)
+void DiagonalMass<DataTypes,MassType>::applyPointCreation(PointID, MassType &m, const Point &, const sofa::type::vector<PointID> &, const sofa::type::vector<double> &)
 {
     m=0;
 }
 
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyPointDestruction(const sofa::type::vector<PointID> & pointsRemoved)
+void DiagonalMass<DataTypes,MassType>::applyPointDestruction(Index id, MassType& VertexMass)
 {
-    helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-    helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
-
-    Index last = masses.size() - 1;
-    for (std::size_t i = 0; i < pointsRemoved.size(); ++i)
-    {
-        Index id = pointsRemoved[i];
-
-        if (id >= masses.size())
-        {
-            msg_error("DMassPointEngine") << "Point to remove is out of bounds from mass vector: " << id << " out of size: " << masses.size();
-            continue;
-        }
-        
-        // Remove the mass of the removed points from totalMass
-        totalMass -= masses[id];
-
-        // swap value before resize
-        masses[id] = masses[last];
-        --last;
-    }
-
-    masses.resize(last + 1);
+    SOFA_UNUSED(id);
+    helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
+    totalMass -= VertexMass;
+    this->cleanTracker();
 }
 
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::PointsRemoved* e)
-{
-    if(!dm->d_computeMassOnRest.getValue())
-        msg_warning("DiagonalMassPointEngine") << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-    const auto& pointsRemoved = e->getArray();
-    applyPointDestruction(pointsRemoved);
-    dm->cleanTracker();
-    dm->printMass();
-
-    msg_info(dm) << "ApplyTopologyChange: EdgesRemoved";
-    msg_info(dm) << "Size of vertexMass: "<< dm->d_vertexMass.getValue().size();
-}
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyEdgeCreation(const sofa::type::vector< EdgeID >& edgeAdded,
+void DiagonalMass<DataTypes,MassType>::applyEdgeCreation(const sofa::type::vector< EdgeID >& edgeAdded,
         const sofa::type::vector< Edge >& /*elems*/,
         const sofa::type::vector< sofa::type::vector< EdgeID > >& /*ancestors*/,
         const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::EDGE)
+    if (this->getMassTopologyType() == TopologyElementType::EDGE)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<edgeAdded.size(); ++i)
         {
             /// get the edge to be added
-            const Edge &e=dm->m_topology->getEdge(edgeAdded[i]);
+            const Edge &e=this->m_topology->getEdge(edgeAdded[i]);
             // compute its mass based on the mass density and the edge length
-            if(dm->edgeGeo)
+            if(this->edgeGeo)
             {
-                mass=(md*dm->edgeGeo->computeRestEdgeLength(edgeAdded[i]))/(typename DataTypes::Real(2.0));
+                mass=(md*this->edgeGeo->computeRestEdgeLength(edgeAdded[i]))/(typename DataTypes::Real(2.0));
             }
             // added mass on its two vertices
             masses[e[0]]+=mass;
@@ -145,29 +110,32 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyEdgeCreation(const
 
             totalMass += 2.0*mass;
         }
+
+        this->cleanTracker();
+        printMass();
     }
 }
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyEdgeDestruction(const sofa::type::vector<EdgeID> & edgeRemoved)
+void DiagonalMass<DataTypes,MassType>::applyEdgeDestruction(const sofa::type::vector<EdgeID> & edgeRemoved)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::EDGE)
+    if (this->getMassTopologyType() == TopologyElementType::EDGE)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<edgeRemoved.size(); ++i)
         {
             /// get the edge to be added
-            const Edge &e=dm->m_topology->getEdge(edgeRemoved[i]);
+            const Edge &e= this->m_topology->getEdge(edgeRemoved[i]);
             // compute its mass based on the mass density and the edge length
-            if(dm->edgeGeo)
+            if(this->edgeGeo)
             {
-                mass=(md*dm->edgeGeo->computeRestEdgeLength(edgeRemoved[i]))/(typename DataTypes::Real (2.0));
+                mass=(md*this->edgeGeo->computeRestEdgeLength(edgeRemoved[i]))/(typename DataTypes::Real (2.0));
             }
             // removed mass on its two vertices
             masses[e[0]]-=mass;
@@ -175,73 +143,36 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyEdgeDestruction(co
 
             totalMass -= 2.0*mass;
         }
-    }
-}
-
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::EdgesAdded* e)
-{
-    const auto& edgeIndex = e->getIndexArray();
-    const sofa::type::vector< Edge >& edges = e->getArray();
-    const auto& ancestors = e->ancestorsList;
-    const sofa::type::vector< sofa::type::vector< double > >& coeffs = e->coefs;
-
-    if(dm->edgeGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning("DiagonalMassPointEngine") << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        if(dm->f_printLog.getValue())
-            msg_info("DiagonalMassPointEngine")<<"ApplyTopologyChange: EdgesAdded";
-
-        applyEdgeCreation(edgeIndex, edges, ancestors, coeffs);
-        dm->cleanTracker();
-        dm->printMass();
-    }
-}
 
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::EdgesRemoved* e)
-{
-    const auto& edgeRemoved = e->getArray();
-
-    if(dm->edgeGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        msg_info(dm) << "ApplyTopologyChange: EdgesRemoved";
-
-        applyEdgeDestruction(edgeRemoved);
-        dm->cleanTracker();
-        dm->printMass();
+        this->cleanTracker();
+        printMass();
     }
 }
 
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTriangleCreation(const sofa::type::vector< TriangleID >& triangleAdded,
+void DiagonalMass<DataTypes,MassType>::applyTriangleCreation(const sofa::type::vector< TriangleID >& triangleAdded,
         const sofa::type::vector< Triangle >& /*elems*/,
         const sofa::type::vector< sofa::type::vector< TriangleID > >& /*ancestors*/,
         const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<triangleAdded.size(); ++i)
         {
             /// get the triangle to be added
-            const Triangle &t=dm->m_topology->getTriangle(triangleAdded[i]);
+            const Triangle &t= this->m_topology->getTriangle(triangleAdded[i]);
             // compute its mass based on the mass density and the triangle area
-            if(dm->triangleGeo)
+            if(this->triangleGeo)
             {
-                mass=(md*dm->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(3.0));
+                mass=(md*this->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(3.0));
             }
             // added mass on its three vertices
             masses[t[0]]+=mass;
@@ -250,30 +181,33 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTriangleCreation(c
 
             totalMass+= 3.0*mass;
         }
+
+        this->cleanTracker();
+        printMass();
     }
 }
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTriangleDestruction(const sofa::type::vector<TriangleID > & triangleRemoved)
+void DiagonalMass<DataTypes,MassType>::applyTriangleDestruction(const sofa::type::vector<TriangleID > & triangleRemoved)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<triangleRemoved.size(); ++i)
         {
             /// get the triangle to be added
-            const Triangle &t=dm->m_topology->getTriangle(triangleRemoved[i]);
+            const Triangle &t= this->m_topology->getTriangle(triangleRemoved[i]);
 
             /// compute its mass based on the mass density and the triangle area
-            if(dm->triangleGeo)
+            if(this->triangleGeo)
             {
-                mass=(md*dm->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(3.0));
+                mass=(md*this->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(3.0));
             }
 
             /// removed  mass on its three vertices
@@ -283,72 +217,113 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTriangleDestructio
 
             totalMass -= 3.0 * mass;
         }
+
+        this->cleanTracker();
+        printMass();
     }
 }
 
+
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::TrianglesAdded* e)
+void DiagonalMass<DataTypes, MassType>::applyQuadCreation(const sofa::type::vector< QuadID >& quadAdded,
+    const sofa::type::vector< Quad >& /*elems*/,
+    const sofa::type::vector< sofa::type::vector< QuadID > >& /*ancestors*/,
+    const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/)
 {
-    const auto& triangleAdded = e->getIndexArray();
-    const sofa::type::vector< Triangle >& elems = e->getElementArray();
-    const auto& ancestors = e->ancestorsList;
-    const sofa::type::vector< sofa::type::vector< double > >& coefs = e->coefs;
-
-    if(dm->triangleGeo)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
+
+        typename DataTypes::Real md = getMassDensity();
+        typename DataTypes::Real mass = typename DataTypes::Real(0);
+        unsigned int i;
+
+        for (i = 0; i < quadAdded.size(); ++i)
+        {
+            /// get the quad to be added
+            const Quad& q = this->m_topology->getQuad(quadAdded[i]);
+            // compute its mass based on the mass density and the quad area
+            if (this->quadGeo)
+            {
+                mass = (md * this->quadGeo->computeRestQuadArea(quadAdded[i])) / (typename DataTypes::Real(4.0));
+            }
+            // added mass on its three vertices
+            masses[q[0]] += mass;
+            masses[q[1]] += mass;
+            masses[q[2]] += mass;
+            masses[q[3]] += mass;
 
-        msg_info(dm) << "ApplyTopologyChange: TrianglesAdded";
+            totalMass += 4.0 * mass;
+        }
 
-        applyTriangleCreation(triangleAdded,elems,ancestors,coefs);
-        dm->cleanTracker();
-        dm->printMass();
+        this->cleanTracker();
+        printMass();
     }
 }
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::TrianglesRemoved* e)
+void DiagonalMass<DataTypes, MassType>::applyQuadDestruction(const sofa::type::vector<QuadID >& quadRemoved)
 {
-    const auto& triangleRemoved = e->getArray();
-
-    if(dm->triangleGeo)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
+
+        typename DataTypes::Real md = getMassDensity();
+        typename DataTypes::Real mass = typename DataTypes::Real(0);
+        unsigned int i;
+
+        for (i = 0; i < quadRemoved.size(); ++i)
+        {
+            /// get the quad to be added
+            const Quad& q = this->m_topology->getQuad(quadRemoved[i]);
+
+            /// compute its mass based on the mass density and the quad area
+            if (this->quadGeo)
+            {
+                mass = (md * this->quadGeo->computeRestQuadArea(quadRemoved[i])) / (typename DataTypes::Real(4.0));
+            }
 
-        msg_info(dm) << "ApplyTopologyChange: TrianglesRemoved";
+            /// removed  mass on its three vertices
+            masses[q[0]] -= mass;
+            masses[q[1]] -= mass;
+            masses[q[2]] -= mass;
+            masses[q[3]] -= mass;
+
+            totalMass -= 4.0 * mass;
+        }
 
-        applyTriangleDestruction(triangleRemoved);
-        dm->cleanTracker();
-        dm->printMass();
+        this->cleanTracker();
+        printMass();
     }
 }
 
+
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTetrahedronCreation(const sofa::type::vector< TetrahedronID >& tetrahedronAdded,
+void DiagonalMass<DataTypes,MassType>::applyTetrahedronCreation(const sofa::type::vector< TetrahedronID >& tetrahedronAdded,
         const sofa::type::vector< Tetrahedron >& /*elems*/,
         const sofa::type::vector< sofa::type::vector< TetrahedronID > >& /*ancestors*/,
         const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<tetrahedronAdded.size(); ++i)
         {
             /// get the tetrahedron to be added
-            const Tetrahedron &t=dm->m_topology->getTetrahedron(tetrahedronAdded[i]);
+            const Tetrahedron &t= this->m_topology->getTetrahedron(tetrahedronAdded[i]);
 
             /// compute its mass based on the mass density and the tetrahedron volume
-            if(dm->tetraGeo)
+            if(this->tetraGeo)
             {
-                mass=(md*dm->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(4.0));
+                mass=(md*this->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(4.0));
             }
 
             /// added  mass on its four vertices
@@ -360,29 +335,31 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTetrahedronCreatio
             totalMass += 4.0*mass;
         }
 
+        this->cleanTracker();
+        printMass();
     }
 }
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTetrahedronDestruction(const sofa::type::vector<TetrahedronID> & tetrahedronRemoved)
+void DiagonalMass<DataTypes,MassType>::applyTetrahedronDestruction(const sofa::type::vector<TetrahedronID> & tetrahedronRemoved)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<tetrahedronRemoved.size(); ++i)
         {
             /// get the tetrahedron to be added
-            const Tetrahedron &t=dm->m_topology->getTetrahedron(tetrahedronRemoved[i]);
-            if(dm->tetraGeo)
+            const Tetrahedron &t= this->m_topology->getTetrahedron(tetrahedronRemoved[i]);
+            if(this->tetraGeo)
             {
                 // compute its mass based on the mass density and the tetrahedron volume
-                mass=(md*dm->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(4.0));
+                mass=(md*this->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(4.0));
             }
             // removed  mass on its four vertices
             masses[t[0]]-=mass;
@@ -393,71 +370,35 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyTetrahedronDestruc
             totalMass -= 4.0*mass;
         }
 
+        this->cleanTracker();
+        printMass();
     }
 }
 
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::TetrahedraAdded* e)
-{
-    const auto& tetrahedronAdded = e->getIndexArray();
-    const sofa::type::vector< Tetrahedron >& elems = e->getElementArray();
-    const auto& ancestors = e->ancestorsList;
-    const sofa::type::vector< sofa::type::vector< double > >& coefs = e->coefs;
-
-    if(dm->tetraGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        msg_info(dm) << "ApplyTopologyChange: TetrahedraAdded";
-
-        applyTetrahedronCreation(tetrahedronAdded, elems, ancestors, coefs);
-        dm->cleanTracker();
-        dm->printMass();
-    }
-}
-
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::TetrahedraRemoved* e)
-{
-    const auto& tetrahedronRemoved = e->getArray();
-
-    if(dm->tetraGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        msg_info(dm) << "ApplyTopologyChange: TetrahedraRemoved";
-
-        applyTetrahedronDestruction(tetrahedronRemoved);
-        dm->cleanTracker();
-        dm->printMass();
-    }
-}
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyHexahedronCreation(const sofa::type::vector< HexahedronID >& hexahedronAdded,
+void DiagonalMass<DataTypes,MassType>::applyHexahedronCreation(const sofa::type::vector< HexahedronID >& hexahedronAdded,
         const sofa::type::vector< Hexahedron >& /*elems*/,
         const sofa::type::vector< sofa::type::vector< HexahedronID > >& /*ancestors*/,
         const sofa::type::vector< sofa::type::vector< double > >& /*coefs*/)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=typename DataTypes::Real(0);
         unsigned int i;
 
         for (i=0; i<hexahedronAdded.size(); ++i)
         {
             /// get the tetrahedron to be added
-            const Hexahedron &t=dm->m_topology->getHexahedron(hexahedronAdded[i]);
+            const Hexahedron &t=this->m_topology->getHexahedron(hexahedronAdded[i]);
             // compute its mass based on the mass density and the tetrahedron volume
-            if(dm->hexaGeo)
+            if(this->hexaGeo)
             {
-                mass=(md*dm->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(8.0));
+                mass=(md*this->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(8.0));
             }
             // added  mass on its four vertices
             for (unsigned int j=0; j<8; ++j)
@@ -466,29 +407,31 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyHexahedronCreation
             totalMass += 8.0*mass;
         }
 
+        this->cleanTracker();
+        printMass();
     }
 }
 
 template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyHexahedronDestruction(const sofa::type::vector<HexahedronID> & hexahedronRemoved)
+void DiagonalMass<DataTypes,MassType>::applyHexahedronDestruction(const sofa::type::vector<HexahedronID> & hexahedronRemoved)
 {
-    if (dm->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor<Data<MassVector> > masses(dm->d_vertexMass);
-        helper::WriteAccessor<Data<Real> > totalMass(dm->d_totalMass);
+        helper::WriteAccessor<Data<MassVector> > masses(d_vertexMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
-        typename DataTypes::Real md=dm->getMassDensity();
+        typename DataTypes::Real md=getMassDensity();
         typename DataTypes::Real mass=(typename DataTypes::Real) 0;
         unsigned int i;
 
         for (i=0; i<hexahedronRemoved.size(); ++i)
         {
             /// get the tetrahedron to be added
-            const Hexahedron &t=dm->m_topology->getHexahedron(hexahedronRemoved[i]);
-            if(dm->hexaGeo)
+            const Hexahedron &t=this->m_topology->getHexahedron(hexahedronRemoved[i]);
+            if(this->hexaGeo)
             {
                 // compute its mass based on the mass density and the tetrahedron volume
-                mass=(md*dm->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(8.0));
+                mass=(md*this->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(8.0));
             }
             // removed  mass on its four vertices
             for (unsigned int j=0; j<8; ++j)
@@ -497,47 +440,10 @@ void DiagonalMass<DataTypes,MassType>::DMassPointEngine::applyHexahedronDestruct
             totalMass -= 8.0*mass;
         }
 
+        this->cleanTracker();
+        printMass();
     }
 }
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::HexahedraAdded* e)
-{
-    const auto& hexahedronAdded = e->getIndexArray();
-    const sofa::type::vector< Hexahedron >& elems = e->getElementArray();
-    const auto& ancestors = e->ancestorsList;
-    const sofa::type::vector< sofa::type::vector< double > >& coefs = e->coefs;
-
-    if(dm->hexaGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        msg_info(dm) << "ApplyTopologyChange: HexahedraAdded";
-
-        applyHexahedronCreation(hexahedronAdded,elems,ancestors,coefs);
-        dm->cleanTracker();
-        dm->printMass();
-    }
-}
-
-template <class DataTypes, class MassType>
-void DiagonalMass<DataTypes,MassType>::DMassPointEngine::ApplyTopologyChange(const core::topology::HexahedraRemoved* e)
-{
-    const auto& hexahedronRemoved = e->getArray();
-
-    if(dm->hexaGeo)
-    {
-        if(!dm->d_computeMassOnRest.getValue())
-            msg_warning(dm) << "ApplyTopologyChange: option computeMassOnRest should be true to have consistent topological change";
-
-        msg_info(dm) << "ApplyTopologyChange: HexahedraRemoved";
-
-        applyHexahedronDestruction(hexahedronRemoved);
-        dm->cleanTracker();
-        dm->printMass();
-    }
-}
-
 
 
 template <class DataTypes, class MassType>
@@ -686,18 +592,79 @@ template <class DataTypes, class MassType>
 void DiagonalMass<DataTypes, MassType>::initTopologyHandlers()
 {
     // add the functions to handle topology changes.
-    m_pointEngine = new DMassPointEngine(this, &d_vertexMass);
-    d_vertexMass.createTopologyHandler(m_topology, m_pointEngine);
-    if (edgeGeo)
+    d_vertexMass.createTopologyHandler(m_topology);
+    d_vertexMass.setCreationCallback([this](Index pointIndex, MassType& m,
+        const core::topology::BaseMeshTopology::Point& point,
+        const sofa::type::vector< Index >& ancestors,
+        const sofa::type::vector< double >& coefs)
+    {
+        applyPointCreation(pointIndex, m, point, ancestors, coefs);
+    });
+    d_vertexMass.setDestructionCallback([this](Index pointIndex, MassType& m)
+    {
+        applyPointDestruction(pointIndex, m);
+    });
+
+    if (edgeGeo) 
+    {
         d_vertexMass.linkToEdgeDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::EDGESADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::EdgesAdded* edgeAdd = static_cast<const core::topology::EdgesAdded*>(eventTopo);
+            applyEdgeCreation(edgeAdd->getIndexArray(), edgeAdd->getElementArray(), edgeAdd->ancestorsList, edgeAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::EDGESREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::EdgesRemoved* edgeRemove = static_cast<const core::topology::EdgesRemoved*>(eventTopo);
+            applyEdgeDestruction(edgeRemove->getArray());
+        });
+    }
     if (triangleGeo)
+    {
         d_vertexMass.linkToTriangleDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TrianglesAdded* tAdd = static_cast<const core::topology::TrianglesAdded*>(eventTopo);
+            applyTriangleCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TrianglesRemoved* tRemove = static_cast<const core::topology::TrianglesRemoved*>(eventTopo);
+            applyTriangleDestruction(tRemove->getArray());
+        });
+    }
     if (quadGeo)
+    {
         d_vertexMass.linkToQuadDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::QuadsAdded* qAdd = static_cast<const core::topology::QuadsAdded*>(eventTopo);
+            applyQuadCreation(qAdd->getIndexArray(), qAdd->getElementArray(), qAdd->ancestorsList, qAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::QuadsRemoved* qRemove = static_cast<const core::topology::QuadsRemoved*>(eventTopo);
+            applyQuadDestruction(qRemove->getArray());
+        });
+    }
     if (tetraGeo)
+    {
         d_vertexMass.linkToTetrahedronDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TetrahedraAdded* tAdd = static_cast<const core::topology::TetrahedraAdded*>(eventTopo);
+            applyTetrahedronCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TetrahedraRemoved* tRemove = static_cast<const core::topology::TetrahedraRemoved*>(eventTopo);
+            applyTetrahedronDestruction(tRemove->getArray());
+        });
+    }
     if (hexaGeo)
+    {
         d_vertexMass.linkToHexahedronDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::HexahedraAdded* hAdd = static_cast<const core::topology::HexahedraAdded*>(eventTopo);
+            applyHexahedronCreation(hAdd->getIndexArray(), hAdd->getElementArray(), hAdd->ancestorsList, hAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::HexahedraRemoved* hRemove = static_cast<const core::topology::HexahedraRemoved*>(eventTopo);
+            applyHexahedronDestruction(hRemove->getArray());
+        });
+    }
 }
 
 template <class DataTypes, class MassType>
@@ -1270,7 +1237,7 @@ template <class DataTypes, class MassType>
 bool DiagonalMass<DataTypes, MassType>::checkTotalMass()
 {
     //Check for negative or null value, if wrongly set use the default value totalMass = 1.0
-    if(d_totalMass.getValue() <= 0.0)
+    if(d_totalMass.getValue() < 0.0)
     {
         msg_warning() << "totalMass data can not have a negative value.\n"
                       << "To remove this warning, you need to set a strictly positive value to the totalMass data";
@@ -1312,7 +1279,7 @@ bool DiagonalMass<DataTypes, MassType>::checkVertexMass()
         //Check that the vertexMass vector has only strictly positive values
         for(size_t i=0; i<vertexMass.size(); i++)
         {
-            if(vertexMass[i]<=0)
+            if(vertexMass[i]<0)
             {
                 msg_warning() << "Negative value of vertexMass vector: vertexMass[" << i << "] = " << vertexMass[i];
                 return false;
@@ -1359,7 +1326,7 @@ bool DiagonalMass<DataTypes, MassType>::checkMassDensity()
     const Real &massDensity = d_massDensity.getValue();
 
     //Check that the massDensity is strictly positive
-    if(massDensity <= 0.0)
+    if(massDensity < 0.0)
     {
         msg_warning() << "Negative value of massDensity: massDensity = " << massDensity;
         return false;
diff --git a/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.h b/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.h
index 4339458f13f..61e4a9714ae 100644
--- a/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.h
+++ b/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.h
@@ -237,187 +237,146 @@ class MeshMatrixMass : public core::behavior::Mass<DataTypes>
 
 
 protected:
+    /** Method to initialize @sa MassType when a new Point is created to compute mass coefficient matrix.
+    * Will be set as creation callback in the PointData @sa d_vertexMass
+    */
+    void applyVertexMassCreation(Index pointIndex, MassType& VertexMass,
+        const core::topology::BaseMeshTopology::Point& point,
+        const sofa::type::vector< Index >&,
+        const sofa::type::vector< double >&);
+
+    /** Method to update @sa d_vertexMass when a Point is removed.
+    * Will be set as destruction callback in the PointData @sa d_vertexMass
+    */
+    void applyVertexMassDestruction(Index, MassType&);
+
+
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a new Triangle is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TRIANGLESADDED event is fired.
+    */
+    void applyVertexMassTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-    class VertexMassHandler : public topology::TopologyDataHandler<core::topology::BaseMeshTopology::Point,MassVector>
-    {
-    public:
-        VertexMassHandler(MeshMatrixMass<DataTypes,TMassType>* _m, topology::PointData<type::vector<TMassType> >* _data) : topology::TopologyDataHandler<core::topology::BaseMeshTopology::Point,type::vector<TMassType> >(_data), m(_m) {}
-
-        /// Mass initialization Creation Functions:
-        /// Vertex mass coefficient matrix creation function
-        void applyCreateFunction(Index pointIndex, TMassType & VertexMass,
-                const sofa::type::vector< Index > &,
-                const sofa::type::vector< double >&);
-
-        /// Apply removing 
-        void applyDestroyFunction(Index, TMassType&);
-
-
-        ///////////////////////// Functions on Triangles //////////////////////////////////////
-
-        /// Mass coefficient Creation/Destruction functions for Triangular Mesh:
-        /// Vertex coefficient of mass matrix creation function to handle creation of new triangles
-        void applyTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
-
-        /// Vertex coefficient of mass matrix destruction function to handle creation of new triangles
-        void applyTriangleDestruction(const sofa::type::vector<Index> & triangleRemoved);
-
-        using topology::TopologyDataHandler<core::topology::BaseMeshTopology::Point,MassVector>::ApplyTopologyChange;
-        /// Callback to add triangles elements.
-        void ApplyTopologyChange(const core::topology::TrianglesAdded* topoEvent);
-        /// Callback to remove triangles elements.
-        void ApplyTopologyChange(const core::topology::TrianglesRemoved* topoEvent);
-
-
-        ///////////////////////// Functions on Quads //////////////////////////////////////
-
-        /// Mass coefficient Creation/Destruction functions for Quad Mesh:
-        /// Vertex coefficient of mass matrix creation function to handle creation of new quads
-        void applyQuadCreation(const sofa::type::vector< Index >& quadAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
-
-        /// Vertex coefficient of mass matrix destruction function to handle creation of new quads
-        void applyQuadDestruction(const sofa::type::vector<Index> & quadRemoved);
-
-        /// Callback to add quads elements.
-        void ApplyTopologyChange(const core::topology::QuadsAdded* topoEvent);
-        /// Callback to remove quads elements.
-        void ApplyTopologyChange(const core::topology::QuadsRemoved* topoEvent);
-
-
-        ///////////////////////// Functions on Tetrahedron //////////////////////////////////////
-
-        /// Mass coefficient Creation/Destruction functions for Tetrahedral Mesh:
-        /// Vertex coefficient of mass matrix creation function to handle creation of new tetrahedra
-        void applyTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
-
-        /// Vertex coefficient of mass matrix destruction function to handle creation of new tetrahedra
-        void applyTetrahedronDestruction(const sofa::type::vector<Index> & tetrahedronRemoved);
-
-        /// Callback to add tetrahedron elements.
-        void ApplyTopologyChange(const core::topology::TetrahedraAdded* topoEvent);
-        /// Callback to remove tetrahedron elements.
-        void ApplyTopologyChange(const core::topology::TetrahedraRemoved* topoEvent);
-
-
-        ///////////////////////// Functions on Hexahedron //////////////////////////////////////
-
-        /// Mass coefficient Creation/Destruction functions for Hexahedral Mesh:
-        /// Vertex coefficient of mass matrix creation function to handle creation of new hexahedra
-        void applyHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
-
-        /// Vertex coefficient of mass matrix destruction function to handle creation of new hexahedra
-        void applyHexahedronDestruction(const sofa::type::vector<Index> & hexahedronRemoved);
-
-        /// Callback to add hexahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::HexahedraAdded* topoEvent);
-         /// Callback to remove hexahedron elements.
-        virtual void ApplyTopologyChange(const core::topology::HexahedraRemoved* topoEvent);
-
-    protected:
-        MeshMatrixMass<DataTypes,TMassType>* m;
-    };
-    VertexMassHandler* m_vertexMassHandler;
-
-    class EdgeMassHandler : public topology::TopologyDataHandler<core::topology::BaseMeshTopology::Edge,MassVector>
-    {
-    public:
-        EdgeMassHandler(MeshMatrixMass<DataTypes,TMassType>* _m, topology::EdgeData<type::vector<TMassType> >* _data) : topology::TopologyDataHandler<core::topology::BaseMeshTopology::Edge,type::vector<TMassType> >(_data), m(_m) {}
-
-        /// Edge mass coefficient matrix creation function
-        void applyCreateFunction(Index edgeIndex, MassType & EdgeMass,
-                const core::topology::BaseMeshTopology::Edge&,
-                const sofa::type::vector< Index > &,
-                const sofa::type::vector< double >&);
-
-        /// Apply removing 
-        void applyDestroyFunction(Index, MassType&);
-
-        using topology::TopologyDataHandler<core::topology::BaseMeshTopology::Edge,MassVector>::ApplyTopologyChange;
-
-        ///////////////////////// Functions on Triangles //////////////////////////////////////
-
-        /// Edge coefficient of mass matrix creation function to handle creation of new triangles
-        void applyTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
-
-        /// Edge coefficient of mass matrix destruction function to handle creation of new triangles
-        void applyTriangleDestruction(const sofa::type::vector<Index> & triangleRemoved);
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a Triangle is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TRIANGLESREMOVED event is fired.
+    */
+    void applyVertexMassTriangleDestruction(const sofa::type::vector<Index>& triangleRemoved);
 
-        /// Callback to add triangles elements.
-        void ApplyTopologyChange(const core::topology::TrianglesAdded* topoEvent);
-        /// Callback to remove triangles elements.
-        void ApplyTopologyChange(const core::topology::TrianglesRemoved* topoEvent);
 
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a new Quad is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when QUADSADDED event is fired.
+    */
+    void applyVertexMassQuadCreation(const sofa::type::vector< Index >& quadAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        ///////////////////////// Functions on Quads //////////////////////////////////////
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a Quad is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when QUADSREMOVED event is fired.
+    */
+    void applyVertexMassQuadDestruction(const sofa::type::vector<Index>& quadRemoved);
 
-        /// Edge coefficient of mass matrix creation function to handle creation of new quads
-        void applyQuadCreation(const sofa::type::vector< Index >& quadAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        /// Edge coefficient of mass matrix destruction function to handle creation of new quads
-        void applyQuadDestruction(const sofa::type::vector<Index> & quadRemoved);
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a new Tetrahedron is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TETRAHEDRAADDED event is fired.
+    */
+    void applyVertexMassTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        /// Callback to add quads elements.
-        void ApplyTopologyChange(const core::topology::QuadsAdded* topoEvent);
-        /// Callback to remove quads elements.
-        void ApplyTopologyChange(const core::topology::QuadsRemoved* topoEvent);
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a Tetrahedron is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when TETRAHEDRAREMOVED event is fired.
+    */
+    void applyVertexMassTetrahedronDestruction(const sofa::type::vector<Index>& tetrahedronRemoved);
+
+    
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a new Hexahedron is created.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when HEXAHEDRAADDED event is fired.
+    */
+    void applyVertexMassHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a Hexahedron is removed.
+    * Will be set as callback in the PointData @sa d_vertexMass to update the mass vector when HEXAHEDRAREMOVED event is fired.
+    */
+    void applyVertexMassHexahedronDestruction(const sofa::type::vector<Index>& hexahedronRemoved);
+   
+
+
+    /** Method to initialize @sa MassType when a new Edge is created to compute mass coefficient matrix.
+    * Will be set as creation callback in the EdgeData @sa d_edgeMass
+    */
+    void applyEdgeMassCreation(Index edgeIndex, MassType& EdgeMass,
+        const core::topology::BaseMeshTopology::Edge&,
+        const sofa::type::vector< Index >&,
+        const sofa::type::vector< double >&);
+
+    /** Method to update @sa d_edgeMass when a Edge is removed.
+    * Will be set as destruction callback in the EdgeData @sa d_edgeMass
+    */
+    void applyEdgeMassDestruction(Index, MassType&);
+
+    
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a new Triangle is created.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when TRIANGLESADDED event is fired.
+    */
+    void applyEdgeMassTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        ///////////////////////// Functions on Tetrahedron //////////////////////////////////////
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a Triangle is removed.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when TRIANGLESREMOVED event is fired.
+    */
+    void applyEdgeMassTriangleDestruction(const sofa::type::vector<Index>& triangleRemoved);
 
-        /// Edge coefficient of mass matrix creation function to handle creation of new tetrahedra
-        void applyTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        /// Edge coefficient of mass matrix destruction function to handle creation of new tetrahedra
-        void applyTetrahedronDestruction(const sofa::type::vector<Index> & tetrahedronRemoved);
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a new Quad is created.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when QUADSADDED event is fired.
+    */
+    void applyEdgeMassQuadCreation(const sofa::type::vector< Index >& quadAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        /// Callback to add tetrahedron elements.
-        void ApplyTopologyChange(const core::topology::TetrahedraAdded* topoEvent);
-        /// Callback to remove tetrahedron elements.
-        void ApplyTopologyChange(const core::topology::TetrahedraRemoved* topoEvent);
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a Quad is removed.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when QUADSREMOVED event is fired.
+    */
+    void applyEdgeMassQuadDestruction(const sofa::type::vector<Index>& quadRemoved);
 
 
-        ///////////////////////// Functions on Hexahedron //////////////////////////////////////
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a new Tetrahedron is created.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when TETRAHEDRAADDED event is fired.
+    */
+    void applyEdgeMassTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-        /// Edge coefficient of mass matrix creation function to handle creation of new hexahedra
-        void applyHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
-                const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
-                const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
-                const sofa::type::vector< sofa::type::vector< double > >& coefs);
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a Tetrahedron is removed.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when TETRAHEDRAREMOVED event is fired.
+    */
+    void applyEdgeMassTetrahedronDestruction(const sofa::type::vector<Index>& tetrahedronRemoved);
 
-        /// Edge coefficient of mass matrix destruction function to handle creation of new hexahedra
-        void applyHexahedronDestruction(const sofa::type::vector<Index> & /*indices*/);
 
-        /// Callback to add hexahedron elements.
-        void ApplyTopologyChange(const core::topology::HexahedraAdded* topoEvent);
-         /// Callback to remove hexahedron elements.
-        void ApplyTopologyChange(const core::topology::HexahedraRemoved* topoEvent);
+    /** Method to update @sa d_edgeMass using mass matrix coefficient when a new Hexahedron is created.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when HEXAHEDRAADDED event is fired.
+    */
+    void applyEdgeMassHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
+        const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
+        const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
+        const sofa::type::vector< sofa::type::vector< double > >& coefs);
 
-    protected:
-        MeshMatrixMass<DataTypes,TMassType>* m;
-    };
+    /** Method to update @sa d_vertexMass using mass matrix coefficient when a Hexahedron is removed.
+    * Will be set as callback in the EdgeData @sa d_edgeMass to update the mass vector when HEXAHEDRAREMOVED event is fired.
+    */
+    void applyEdgeMassHexahedronDestruction(const sofa::type::vector<Index>& /*indices*/);
 
-    EdgeMassHandler* m_edgeMassHandler;
 
     sofa::core::topology::BaseMeshTopology* m_topology;
 };
diff --git a/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.inl b/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.inl
index 01c7cd4f96a..c308a0a9097 100644
--- a/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.inl
+++ b/modules/SofaMiscForceField/src/SofaMiscForceField/MeshMatrixMass.inl
@@ -60,8 +60,6 @@ MeshMatrixMass<DataTypes, MassType>::MeshMatrixMass()
     , f_graph( initData(&f_graph,"graph","Graph of the controlled potential") )
     , l_topology(initLink("topology", "link to the topology container"))
     , m_massTopologyType(TopologyElementType::UNKNOWN)
-    , m_vertexMassHandler(nullptr)
-    , m_edgeMassHandler(nullptr)
     , m_topology(nullptr)
 {
     f_graph.setWidget("graph");
@@ -72,44 +70,44 @@ MeshMatrixMass<DataTypes, MassType>::MeshMatrixMass()
 template <class DataTypes, class MassType>
 MeshMatrixMass<DataTypes, MassType>::~MeshMatrixMass()
 {
-    if (m_vertexMassHandler) delete m_vertexMassHandler;
-    if (m_edgeMassHandler) delete m_edgeMassHandler;
+
 }
 
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyCreateFunction(Index, MassType & VertexMass,
-        const sofa::type::vector< Index > &,
-        const sofa::type::vector< double >&)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassCreation(Index, MassType & VertexMass,
+    const core::topology::BaseMeshTopology::Point&,
+    const sofa::type::vector< Index > &,
+    const sofa::type::vector< double >&)
 {
     VertexMass = 0;
 }
 
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyCreateFunction(Index, MassType & EdgeMass,
-        const core::topology::BaseMeshTopology::Edge&,
-        const sofa::type::vector< Index > &,
-        const sofa::type::vector< double >&)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassCreation(Index, MassType & EdgeMass,
+    const core::topology::BaseMeshTopology::Edge&,
+    const sofa::type::vector< Index > &,
+    const sofa::type::vector< double >&)
 {
     EdgeMass = 0;
 }
 
 
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyDestroyFunction(Index id, MassType& VertexMass)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassDestruction(Index id, MassType& VertexMass)
 {
     SOFA_UNUSED(id);
-    helper::WriteAccessor<Data<Real> > totalMass(this->m->d_totalMass);
+    helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
     totalMass -= VertexMass;
 }
 
 
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyDestroyFunction(Index id, MassType& EdgeMass)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassDestruction(Index id, MassType& EdgeMass)
 {
     SOFA_UNUSED(id);
-    if(!this->m->isLumped())
+    if(!isLumped())
     {
-        helper::WriteAccessor<Data<Real> > totalMass(this->m->d_totalMass);
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
         totalMass -= EdgeMass;
     }
 }
@@ -123,44 +121,42 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyDestroyFunction(
 
 /// Creation fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);        
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses (d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);        
 
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbTri = MMM->m_topology->getNbTriangles();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbTri = this->m_topology->getNbTriangles();
         if (nbMass < nbTri)
-            MMM->addMassDensity(triangleAdded, ancestors, coefs);
+            addMassDensity(triangleAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<triangleAdded.size(); ++i)
         {
             // Get the triangle to be added
-            const core::topology::BaseMeshTopology::Triangle &t = MMM->m_topology->getTriangle(triangleAdded[i]);
+            const core::topology::BaseMeshTopology::Triangle &t = this->m_topology->getTriangle(triangleAdded[i]);
 
             // Compute rest mass of conserne triangle = density * triangle surface.
-            if(MMM->triangleGeo)
+            if(this->triangleGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[triangleAdded[i]] * MMM->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(6.0));
+                    mass=(densityM[triangleAdded[i]] * this->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(6.0));
                 }
                 else
                 {
-                    mass=(densityM[triangleAdded[i]] * MMM->triangleGeo->computeTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(6.0));
+                    mass=(densityM[triangleAdded[i]] * this->triangleGeo->computeTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(6.0));
                 }
             }
 
@@ -169,13 +165,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleCreati
                 VertexMasses[ t[j] ] += mass;
 
             // update total mass: 
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 3.0 * mass;
             }
             else
             {
-                totalMass += 3.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass += 3.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -183,44 +179,42 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleCreati
 
 /// Creation fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassTriangleCreation(const sofa::type::vector< Index >& triangleAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Triangle >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbTri = MMM->m_topology->getNbTriangles();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbTri = this->m_topology->getNbTriangles();
         if (nbMass < nbTri)
-            MMM->addMassDensity(triangleAdded, ancestors, coefs);
+            addMassDensity(triangleAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<triangleAdded.size(); ++i)
         {
             // Get the edgesInTriangle to be added
-            const core::topology::BaseMeshTopology::EdgesInTriangle &te = MMM->m_topology->getEdgesInTriangle(triangleAdded[i]);
+            const core::topology::BaseMeshTopology::EdgesInTriangle &te = this->m_topology->getEdgesInTriangle(triangleAdded[i]);
 
             // Compute rest mass of conserne triangle = density * triangle surface.
-            if(MMM->triangleGeo)
+            if(this->triangleGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[triangleAdded[i]] * MMM->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(12.0));
+                    mass=(densityM[triangleAdded[i]] * this->triangleGeo->computeRestTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(12.0));
                 }
                 else
                 {
-                    mass=(densityM[triangleAdded[i]] * MMM->triangleGeo->computeTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(12.0));
+                    mass=(densityM[triangleAdded[i]] * this->triangleGeo->computeTriangleArea(triangleAdded[i]))/(typename DataTypes::Real(12.0));
                 }
             }
 
@@ -229,7 +223,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleCreation
                 EdgeMasses[ te[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 3.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -239,33 +233,32 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleCreation
 
 /// Destruction fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleDestruction(const sofa::type::vector< Index >& triangleRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassTriangleDestruction(const sofa::type::vector< Index >& triangleRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses (d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<triangleRemoved.size(); ++i)
         {
             // Get the triangle to be removed
-            const core::topology::BaseMeshTopology::Triangle &t = MMM->m_topology->getTriangle(triangleRemoved[i]);
+            const core::topology::BaseMeshTopology::Triangle &t = this->m_topology->getTriangle(triangleRemoved[i]);
 
             // Compute rest mass of conserne triangle = density * triangle surface.
-            if(MMM->triangleGeo)
+            if(this->triangleGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[triangleRemoved[i]] * MMM->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(6.0));
+                    mass=(densityM[triangleRemoved[i]] * this->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(6.0));
                 }
                 else
                 {
-                    mass=(densityM[triangleRemoved[i]] * MMM->triangleGeo->computeTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(6.0));
+                    mass=(densityM[triangleRemoved[i]] * this->triangleGeo->computeTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(6.0));
                 }
             }
 
@@ -274,13 +267,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleDestru
                 VertexMasses[ t[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 3.0 * mass;
             }
             else
             {
-                totalMass -= 3.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass -= 3.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -288,33 +281,32 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTriangleDestru
 
 /// Destruction fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleDestruction(const sofa::type::vector< Index >& triangleRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassTriangleDestruction(const sofa::type::vector< Index >& triangleRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TRIANGLE)
+    if (this->getMassTopologyType() == TopologyElementType::TRIANGLE)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses (d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<triangleRemoved.size(); ++i)
         {
             // Get the triangle to be removed
-            const core::topology::BaseMeshTopology::EdgesInTriangle &te = MMM->m_topology->getEdgesInTriangle(triangleRemoved[i]);
+            const core::topology::BaseMeshTopology::EdgesInTriangle &te = this->m_topology->getEdgesInTriangle(triangleRemoved[i]);
 
             // Compute rest mass of conserne triangle = density * triangle surface.
-            if(MMM->triangleGeo)
+            if(this->triangleGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[triangleRemoved[i]] * MMM->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(12.0));
+                    mass=(densityM[triangleRemoved[i]] * this->triangleGeo->computeRestTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(12.0));
                 }
                 else
                 {
-                    mass=(densityM[triangleRemoved[i]] * MMM->triangleGeo->computeTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(12.0));
+                    mass=(densityM[triangleRemoved[i]] * this->triangleGeo->computeTriangleArea(triangleRemoved[i]))/(typename DataTypes::Real(12.0));
                 }
             }
 
@@ -323,7 +315,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleDestruct
                 EdgeMasses[ te[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 3.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -331,50 +323,6 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTriangleDestruct
     }
 }
 
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::TrianglesAdded* topoEvent)
-{
-    const auto &triangleAdded = topoEvent->getIndexArray();
-    const sofa::type::vector<core::topology::BaseMeshTopology::Triangle> &elems = topoEvent->getElementArray();
-    const auto & ancestors = topoEvent->ancestorsList;
-    const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-    applyTriangleCreation(triangleAdded, elems, ancestors, coefs);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::TrianglesRemoved* e)
-{
-    const auto &triangleRemoved = e->getArray();
-
-    applyTriangleDestruction(triangleRemoved);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::TrianglesAdded* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &triangleAdded = topoEvent->getIndexArray();
-        const sofa::type::vector<core::topology::BaseMeshTopology::Triangle> &elems = topoEvent->getElementArray();
-        const auto & ancestors = topoEvent->ancestorsList;
-        const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-        applyTriangleCreation(triangleAdded, elems, ancestors, coefs);
-    }
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::TrianglesRemoved* e)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &triangleRemoved = e->getArray();
-
-        applyTriangleDestruction(triangleRemoved);
-    }
-}
-
 // }
 
 // ---------------------------------------------------
@@ -384,44 +332,42 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(c
 
 /// Creation fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadCreation(const sofa::type::vector< Index >& quadAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassQuadCreation(const sofa::type::vector< Index >& quadAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::QUAD)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbQ = MMM->m_topology->getNbQuads();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbQ = this->m_topology->getNbQuads();
         if (nbMass < nbQ)
-            MMM->addMassDensity(quadAdded, ancestors, coefs);
+            addMassDensity(quadAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<quadAdded.size(); ++i)
         {
             // Get the quad to be added
-            const core::topology::BaseMeshTopology::Quad &q = MMM->m_topology->getQuad(quadAdded[i]);
+            const core::topology::BaseMeshTopology::Quad &q = this->m_topology->getQuad(quadAdded[i]);
 
             // Compute rest mass of conserne quad = density * quad surface.
-            if(MMM->quadGeo)
+            if(this->quadGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[quadAdded[i]] * MMM->quadGeo->computeRestQuadArea(quadAdded[i]))/(typename DataTypes::Real(8.0));
+                    mass=(densityM[quadAdded[i]] * this->quadGeo->computeRestQuadArea(quadAdded[i]))/(typename DataTypes::Real(8.0));
                 }
                 else
                 {
-                    mass=(densityM[quadAdded[i]] * MMM->quadGeo->computeQuadArea(quadAdded[i]))/(typename DataTypes::Real(8.0));
+                    mass=(densityM[quadAdded[i]] * this->quadGeo->computeQuadArea(quadAdded[i]))/(typename DataTypes::Real(8.0));
                 }
             }
 
@@ -430,13 +376,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadCreation(c
                 VertexMasses[ q[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 4.0 * mass;
             }
             else
             {
-                totalMass += 4.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass += 4.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -444,44 +390,43 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadCreation(c
 
 /// Creation fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadCreation(const sofa::type::vector< Index >& quadAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassQuadCreation(const sofa::type::vector< Index >& quadAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Quad >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
 
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::QUAD)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
         
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbQ = MMM->m_topology->getNbQuads();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbQ = this->m_topology->getNbQuads();
         if (nbMass < nbQ)
-            MMM->addMassDensity(quadAdded, ancestors, coefs);
+            addMassDensity(quadAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<quadAdded.size(); ++i)
         {
             // Get the EdgesInQuad to be added
-            const core::topology::BaseMeshTopology::EdgesInQuad &qe = MMM->m_topology->getEdgesInQuad(quadAdded[i]);
+            const core::topology::BaseMeshTopology::EdgesInQuad &qe = this->m_topology->getEdgesInQuad(quadAdded[i]);
 
             // Compute rest mass of conserne quad = density * quad surface.
-            if(MMM->quadGeo)
+            if(this->quadGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[quadAdded[i]] * MMM->quadGeo->computeRestQuadArea(quadAdded[i]))/(typename DataTypes::Real(16.0));
+                    mass=(densityM[quadAdded[i]] * this->quadGeo->computeRestQuadArea(quadAdded[i]))/(typename DataTypes::Real(16.0));
                 }
                 else
                 {
-                    mass=(densityM[quadAdded[i]] * MMM->quadGeo->computeQuadArea(quadAdded[i]))/(typename DataTypes::Real(16.0));
+                    mass=(densityM[quadAdded[i]] * this->quadGeo->computeQuadArea(quadAdded[i]))/(typename DataTypes::Real(16.0));
                 }
             }
 
@@ -490,7 +435,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadCreation(con
                 EdgeMasses[ qe[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 4.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -500,33 +445,32 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadCreation(con
 
 /// Destruction fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadDestruction(const sofa::type::vector< Index >& quadRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassQuadDestruction(const sofa::type::vector< Index >& quadRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::QUAD)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<quadRemoved.size(); ++i)
         {
             // Get the quad to be removed
-            const core::topology::BaseMeshTopology::Quad &q = MMM->m_topology->getQuad(quadRemoved[i]);
+            const core::topology::BaseMeshTopology::Quad &q = this->m_topology->getQuad(quadRemoved[i]);
 
             // Compute rest mass of conserne quad = density * quad surface.
-            if(MMM->quadGeo)
+            if(this->quadGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[quadRemoved[i]] * MMM->quadGeo->computeRestQuadArea(quadRemoved[i]))/(typename DataTypes::Real(8.0));
+                    mass=(densityM[quadRemoved[i]] * this->quadGeo->computeRestQuadArea(quadRemoved[i]))/(typename DataTypes::Real(8.0));
                 }
                 else
                 {
-                    mass=(densityM[quadRemoved[i]] * MMM->quadGeo->computeQuadArea(quadRemoved[i]))/(typename DataTypes::Real(8.0));
+                    mass=(densityM[quadRemoved[i]] * this->quadGeo->computeQuadArea(quadRemoved[i]))/(typename DataTypes::Real(8.0));
                 }
             }
 
@@ -535,13 +479,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadDestructio
                 VertexMasses[ q[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 4.0 * mass;
             }
             else
             {
-                totalMass -= 4.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass -= 4.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -549,33 +493,32 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyQuadDestructio
 
 /// Destruction fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadDestruction(const sofa::type::vector< Index >& quadRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassQuadDestruction(const sofa::type::vector< Index >& quadRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::QUAD)
+    if (this->getMassTopologyType() == TopologyElementType::QUAD)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<quadRemoved.size(); ++i)
         {
             // Get the EdgesInQuad to be removed
-            const core::topology::BaseMeshTopology::EdgesInQuad &qe = MMM->m_topology->getEdgesInQuad(quadRemoved[i]);
+            const core::topology::BaseMeshTopology::EdgesInQuad &qe = this->m_topology->getEdgesInQuad(quadRemoved[i]);
 
             // Compute rest mass of conserne quad = density * quad surface.
-            if(MMM->quadGeo)
+            if(this->quadGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[quadRemoved[i]] * MMM->quadGeo->computeRestQuadArea(quadRemoved[i]))/(typename DataTypes::Real(16.0));
+                    mass=(densityM[quadRemoved[i]] * this->quadGeo->computeRestQuadArea(quadRemoved[i]))/(typename DataTypes::Real(16.0));
                 }
                 else
                 {
-                    mass=(densityM[quadRemoved[i]] * MMM->quadGeo->computeQuadArea(quadRemoved[i]))/(typename DataTypes::Real(16.0));
+                    mass=(densityM[quadRemoved[i]] * this->quadGeo->computeQuadArea(quadRemoved[i]))/(typename DataTypes::Real(16.0));
                 }
             }
 
@@ -584,7 +527,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadDestruction(
                 EdgeMasses[ qe[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 4.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -592,50 +535,6 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyQuadDestruction(
     }
 }
 
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::QuadsAdded* topoEvent)
-{
-    const auto &quadAdded = topoEvent->getIndexArray();
-    const sofa::type::vector<core::topology::BaseMeshTopology::Quad> &elems = topoEvent->getElementArray();
-    const auto & ancestors = topoEvent->ancestorsList;
-    const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-    applyQuadCreation(quadAdded, elems, ancestors, coefs);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::QuadsRemoved* topoEvent)
-{
-    const auto &quadRemoved = topoEvent->getArray();
-
-    applyQuadDestruction(quadRemoved);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::QuadsAdded* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &quadAdded = topoEvent->getIndexArray();
-        const sofa::type::vector<core::topology::BaseMeshTopology::Quad> &elems = topoEvent->getElementArray();
-        const auto& ancestors = topoEvent->ancestorsList;
-        const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-        applyQuadCreation(quadAdded, elems, ancestors, coefs);
-    }
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::QuadsRemoved* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &quadRemoved = topoEvent->getArray();
-
-        applyQuadDestruction(quadRemoved);
-    }
-}
-
 // }
 
 
@@ -647,44 +546,42 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(c
 
 /// Creation fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbT = MMM->m_topology->getNbTetrahedra();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbT = this->m_topology->getNbTetrahedra();
         if (nbMass < nbT)
-            MMM->addMassDensity(tetrahedronAdded, ancestors, coefs);
+            addMassDensity(tetrahedronAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<tetrahedronAdded.size(); ++i)
         {
             // Get the tetrahedron to be added
-            const core::topology::BaseMeshTopology::Tetrahedron &t = MMM->m_topology->getTetrahedron(tetrahedronAdded[i]);
+            const core::topology::BaseMeshTopology::Tetrahedron &t = this->m_topology->getTetrahedron(tetrahedronAdded[i]);
 
             // Compute rest mass of conserne tetrahedron = density * tetrahedron volume.
-            if(MMM->tetraGeo)
+            if(this->tetraGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[tetrahedronAdded[i]] * MMM->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(10.0));
+                    mass=(densityM[tetrahedronAdded[i]] * this->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(10.0));
                 }
                 else
                 {
-                    mass=(densityM[tetrahedronAdded[i]] * MMM->tetraGeo->computeTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(10.0));
+                    mass=(densityM[tetrahedronAdded[i]] * this->tetraGeo->computeTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(10.0));
                 }
             }
 
@@ -693,13 +590,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronCre
                 VertexMasses[ t[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 4.0 * mass;
             }
             else
             {
-                totalMass += 4.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass += 4.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -707,44 +604,42 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronCre
 
 /// Creation fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassTetrahedronCreation(const sofa::type::vector< Index >& tetrahedronAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Tetrahedron >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
         
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbT = MMM->m_topology->getNbTetrahedra();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbT = this->m_topology->getNbTetrahedra();
         if (nbMass < nbT)
-            MMM->addMassDensity(tetrahedronAdded, ancestors, coefs);
+            addMassDensity(tetrahedronAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<tetrahedronAdded.size(); ++i)
         {
             // Get the edgesInTetrahedron to be added
-            const core::topology::BaseMeshTopology::EdgesInTetrahedron &te = MMM->m_topology->getEdgesInTetrahedron(tetrahedronAdded[i]);
+            const core::topology::BaseMeshTopology::EdgesInTetrahedron &te = this->m_topology->getEdgesInTetrahedron(tetrahedronAdded[i]);
 
             // Compute rest mass of conserne triangle = density * tetrahedron volume.
-            if(MMM->tetraGeo)
+            if(this->tetraGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[tetrahedronAdded[i]] * MMM->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[tetrahedronAdded[i]] * this->tetraGeo->computeRestTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(20.0));
                 }
                 else
                 {
-                    mass=(densityM[tetrahedronAdded[i]] * MMM->tetraGeo->computeTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[tetrahedronAdded[i]] * this->tetraGeo->computeTetrahedronVolume(tetrahedronAdded[i]))/(typename DataTypes::Real(20.0));
                 }
             }
 
@@ -753,7 +648,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronCreat
                 EdgeMasses[ te[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 6.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -763,33 +658,32 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronCreat
 
 /// Destruction fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronDestruction(const sofa::type::vector< Index >& tetrahedronRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassTetrahedronDestruction(const sofa::type::vector< Index >& tetrahedronRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<tetrahedronRemoved.size(); ++i)
         {
             // Get the tetrahedron to be removed
-            const core::topology::BaseMeshTopology::Tetrahedron &t = MMM->m_topology->getTetrahedron(tetrahedronRemoved[i]);
+            const core::topology::BaseMeshTopology::Tetrahedron &t = this->m_topology->getTetrahedron(tetrahedronRemoved[i]);
 
             // Compute rest mass of conserne tetrahedron = density * tetrahedron volume.
-            if(MMM->tetraGeo)
+            if(this->tetraGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[tetrahedronRemoved[i]] * MMM->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(10.0));
+                    mass=(densityM[tetrahedronRemoved[i]] * this->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(10.0));
                 }
                 else
                 {
-                    mass=(densityM[tetrahedronRemoved[i]] * MMM->tetraGeo->computeTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(10.0));
+                    mass=(densityM[tetrahedronRemoved[i]] * this->tetraGeo->computeTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(10.0));
                 }
             }
 
@@ -798,13 +692,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronDes
                 VertexMasses[ t[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 4.0 * mass;
             }
             else
             {
-                totalMass -= 4.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass -= 4.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -812,33 +706,32 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyTetrahedronDes
 
 /// Destruction fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronDestruction(const sofa::type::vector< Index >& tetrahedronRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassTetrahedronDestruction(const sofa::type::vector< Index >& tetrahedronRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::TETRAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<tetrahedronRemoved.size(); ++i)
         {
             // Get the edgesInTetrahedron to be removed
-            const core::topology::BaseMeshTopology::EdgesInTetrahedron &te = MMM->m_topology->getEdgesInTetrahedron(tetrahedronRemoved[i]);
+            const core::topology::BaseMeshTopology::EdgesInTetrahedron &te = this->m_topology->getEdgesInTetrahedron(tetrahedronRemoved[i]);
 
             // Compute rest mass of conserne triangle = density * tetrahedron volume.
-            if(MMM->tetraGeo)
+            if(this->tetraGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[tetrahedronRemoved[i]] * MMM->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[tetrahedronRemoved[i]] * this->tetraGeo->computeRestTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(20.0));
                 }
                 else
                 {
-                    mass=(densityM[tetrahedronRemoved[i]] * MMM->tetraGeo->computeTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[tetrahedronRemoved[i]] * this->tetraGeo->computeTetrahedronVolume(tetrahedronRemoved[i]))/(typename DataTypes::Real(20.0));
                 }
             }
 
@@ -847,7 +740,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronDestr
                 EdgeMasses[ te[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 6.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -855,50 +748,6 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyTetrahedronDestr
     }
 }
 
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::TetrahedraAdded* topoEvent)
-{
-    const auto &tetraAdded = topoEvent->getIndexArray();
-    const sofa::type::vector<core::topology::BaseMeshTopology::Tetrahedron> &elems = topoEvent->getElementArray();
-    const auto & ancestors = topoEvent->ancestorsList;
-    const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-    applyTetrahedronCreation(tetraAdded, elems, ancestors, coefs);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::TetrahedraRemoved* topoEvent)
-{
-    const auto &tetraRemoved = topoEvent->getArray();
-
-    applyTetrahedronDestruction(tetraRemoved);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::TetrahedraAdded* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &tetraAdded = topoEvent->getIndexArray();
-        const sofa::type::vector<core::topology::BaseMeshTopology::Tetrahedron> &elems = topoEvent->getElementArray();
-        const auto& ancestors = topoEvent->ancestorsList;
-        const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-        applyTetrahedronCreation(tetraAdded, elems, ancestors, coefs);
-    }
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::TetrahedraRemoved* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto& tetraRemoved = topoEvent->getArray();
-
-        applyTetrahedronDestruction(tetraRemoved);
-    }
-}
-
 // }
 
 
@@ -909,44 +758,42 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(c
 
 /// Creation fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
         
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbT = MMM->m_topology->getNbHexahedra();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbT = this->m_topology->getNbHexahedra();
         if (nbMass < nbT)
-            MMM->addMassDensity(hexahedronAdded, ancestors, coefs);
+            addMassDensity(hexahedronAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<hexahedronAdded.size(); ++i)
         {
             // Get the hexahedron to be added
-            const core::topology::BaseMeshTopology::Hexahedron &h = MMM->m_topology->getHexahedron(hexahedronAdded[i]);
+            const core::topology::BaseMeshTopology::Hexahedron &h = this->m_topology->getHexahedron(hexahedronAdded[i]);
 
             // Compute rest mass of conserne hexahedron = density * hexahedron volume.
-            if(MMM->hexaGeo)
+            if(this->hexaGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[hexahedronAdded[i]] * MMM->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[hexahedronAdded[i]] * this->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(20.0));
                 }
                 else
                 {
-                    mass=(densityM[hexahedronAdded[i]] * MMM->hexaGeo->computeHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[hexahedronAdded[i]] * this->hexaGeo->computeHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(20.0));
                 }
             }
 
@@ -955,13 +802,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronCrea
                 VertexMasses[ h[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 8.0 * mass;
             }
             else
             {
-                totalMass += 8.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass += 8.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -969,44 +816,42 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronCrea
 
 /// Creation fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassHexahedronCreation(const sofa::type::vector< Index >& hexahedronAdded,
         const sofa::type::vector< core::topology::BaseMeshTopology::Hexahedron >& elems,
         const sofa::type::vector< sofa::type::vector< Index > >& ancestors,
         const sofa::type::vector< sofa::type::vector< double > >& coefs)
 {
     SOFA_UNUSED(elems);
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // update mass density vector
-        sofa::Size nbMass = MMM->getMassDensity().size();
-        sofa::Size nbT = MMM->m_topology->getNbHexahedra();
+        sofa::Size nbMass = getMassDensity().size();
+        sofa::Size nbT = this->m_topology->getNbHexahedra();
         if (nbMass < nbT)
-            MMM->addMassDensity(hexahedronAdded, ancestors, coefs);
+            addMassDensity(hexahedronAdded, ancestors, coefs);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<hexahedronAdded.size(); ++i)
         {
             // Get the EdgesInHexahedron to be added
-            const core::topology::BaseMeshTopology::EdgesInHexahedron &he = MMM->m_topology->getEdgesInHexahedron(hexahedronAdded[i]);
+            const core::topology::BaseMeshTopology::EdgesInHexahedron &he = this->m_topology->getEdgesInHexahedron(hexahedronAdded[i]);
 
             // Compute rest mass of conserne hexahedron = density * hexahedron volume.
-            if(MMM->hexaGeo)
+            if(this->hexaGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[hexahedronAdded[i]] * MMM->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(40.0));
+                    mass=(densityM[hexahedronAdded[i]] * this->hexaGeo->computeRestHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(40.0));
                 }
                 else
                 {
-                    mass=(densityM[hexahedronAdded[i]] * MMM->hexaGeo->computeHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(40.0));
+                    mass=(densityM[hexahedronAdded[i]] * this->hexaGeo->computeHexahedronVolume(hexahedronAdded[i]))/(typename DataTypes::Real(40.0));
                 }
             }
 
@@ -1015,7 +860,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronCreati
                 EdgeMasses[ he[j] ] += mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass += 12.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -1025,33 +870,32 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronCreati
 
 /// Destruction fonction for mass stored on vertices
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronDestruction(const sofa::type::vector< Index >& hexahedronRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyVertexMassHexahedronDestruction(const sofa::type::vector< Index >& hexahedronRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( MMM->d_vertexMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > VertexMasses ( d_vertexMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<hexahedronRemoved.size(); ++i)
         {
             // Get the hexahedron to be removed
-            const core::topology::BaseMeshTopology::Hexahedron &h = MMM->m_topology->getHexahedron(hexahedronRemoved[i]);
+            const core::topology::BaseMeshTopology::Hexahedron &h = this->m_topology->getHexahedron(hexahedronRemoved[i]);
 
             // Compute rest mass of conserne hexahedron = density * hexahedron volume.
-            if(MMM->hexaGeo)
+            if(this->hexaGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[hexahedronRemoved[i]] * MMM->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[hexahedronRemoved[i]] * this->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(20.0));
                 }
                 else
                 {
-                    mass=(densityM[hexahedronRemoved[i]] * MMM->hexaGeo->computeHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(20.0));
+                    mass=(densityM[hexahedronRemoved[i]] * this->hexaGeo->computeHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(20.0));
                 }
             }
 
@@ -1060,13 +904,13 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronDest
                 VertexMasses[ h[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 8.0 * mass;
             }
             else
             {
-                totalMass -= 8.0 * mass * MMM->m_massLumpingCoeff;
+                totalMass -= 8.0 * mass * m_massLumpingCoeff;
             }
         }
     }
@@ -1074,33 +918,32 @@ void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::applyHexahedronDest
 
 /// Destruction fonction for mass stored on edges
 template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronDestruction(const sofa::type::vector< Index >& hexahedronRemoved)
+void MeshMatrixMass<DataTypes, MassType>::applyEdgeMassHexahedronDestruction(const sofa::type::vector< Index >& hexahedronRemoved)
 {
-    MeshMatrixMass<DataTypes, MassType> *MMM = this->m;
-    if (MMM && MMM->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
+    if (this->getMassTopologyType() == TopologyElementType::HEXAHEDRON)
     {
-        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( MMM->d_edgeMass );
-        helper::WriteAccessor<Data<Real> > totalMass(MMM->d_totalMass);
+        helper::WriteAccessor< Data< type::vector<MassType> > > EdgeMasses ( d_edgeMass );
+        helper::WriteAccessor<Data<Real> > totalMass(d_totalMass);
 
         // Initialisation
-        const type::vector<Real> densityM = MMM->getMassDensity();
+        const type::vector<Real> densityM = getMassDensity();
         typename DataTypes::Real mass = typename DataTypes::Real(0);
 
         for (unsigned int i = 0; i<hexahedronRemoved.size(); ++i)
         {
             // Get the EdgesInHexahedron to be removed
-            const core::topology::BaseMeshTopology::EdgesInHexahedron &he = MMM->m_topology->getEdgesInHexahedron(hexahedronRemoved[i]);
+            const core::topology::BaseMeshTopology::EdgesInHexahedron &he = this->m_topology->getEdgesInHexahedron(hexahedronRemoved[i]);
 
             // Compute rest mass of conserne hexahedron = density * hexahedron volume.
-            if(MMM->hexaGeo)
+            if(this->hexaGeo)
             {
-                if(MMM->d_computeMassOnRest.getValue())
+                if(d_computeMassOnRest.getValue())
                 {
-                    mass=(densityM[hexahedronRemoved[i]] * MMM->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(40.0));
+                    mass=(densityM[hexahedronRemoved[i]] * this->hexaGeo->computeRestHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(40.0));
                 }
                 else
                 {
-                    mass=(densityM[hexahedronRemoved[i]] * MMM->hexaGeo->computeHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(40.0));
+                    mass=(densityM[hexahedronRemoved[i]] * this->hexaGeo->computeHexahedronVolume(hexahedronRemoved[i]))/(typename DataTypes::Real(40.0));
                 }
             }
 
@@ -1109,7 +952,7 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronDestru
                 EdgeMasses[ he[j] ] -= mass;
 
             // update total mass
-            if(!this->m->isLumped())
+            if(!this->isLumped())
             {
                 totalMass -= 12.0 * mass * 2.0; // x 2 because mass is actually splitted over half-edges
             }
@@ -1117,49 +960,6 @@ void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::applyHexahedronDestru
     }
 }
 
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::HexahedraAdded* topoEvent)
-{
-    const auto &hexaAdded = topoEvent->getIndexArray();
-    const sofa::type::vector<core::topology::BaseMeshTopology::Hexahedron> &elems = topoEvent->getElementArray();
-    const auto & ancestors = topoEvent->ancestorsList;
-    const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-    applyHexahedronCreation(hexaAdded, elems, ancestors, coefs);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::VertexMassHandler::ApplyTopologyChange(const core::topology::HexahedraRemoved* topoEvent)
-{
-    const auto &hexaRemoved = topoEvent->getArray();
-
-    applyHexahedronDestruction(hexaRemoved);
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::HexahedraAdded* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &hexaAdded = topoEvent->getIndexArray();
-        const sofa::type::vector<core::topology::BaseMeshTopology::Hexahedron> &elems = topoEvent->getElementArray();
-        const auto & ancestors = topoEvent->ancestorsList;
-        const sofa::type::vector<sofa::type::vector<double> > & coefs = topoEvent->coefs;
-
-        applyHexahedronCreation(hexaAdded, elems, ancestors, coefs);
-    }
-}
-
-template< class DataTypes, class MassType>
-void MeshMatrixMass<DataTypes, MassType>::EdgeMassHandler::ApplyTopologyChange(const core::topology::HexahedraRemoved* topoEvent)
-{
-    if(!this->m->isLumped())
-    {
-        const auto &hexaRemoved = topoEvent->getArray();
-
-        applyHexahedronDestruction(hexaRemoved);
-    }
-}
 
 // }
 
@@ -1298,13 +1098,32 @@ template <class DataTypes, class MassType>
 void MeshMatrixMass<DataTypes, MassType>::initTopologyHandlers(sofa::core::topology::TopologyElementType topologyType)
 {
     // add the functions to handle topology changes for Vertex informations
-    m_vertexMassHandler = new VertexMassHandler(this, &d_vertexMass);
-    d_vertexMass.createTopologyHandler(m_topology, m_vertexMassHandler);
+    d_vertexMass.createTopologyHandler(m_topology);
+    d_vertexMass.setCreationCallback([this](Index pointIndex, MassType& m,
+        const core::topology::BaseMeshTopology::Point& point,
+        const sofa::type::vector< Index >& ancestors,
+        const sofa::type::vector< double >& coefs)
+    {
+        applyVertexMassCreation(pointIndex, m, point, ancestors, coefs);
+    });
+    d_vertexMass.setDestructionCallback([this](Index pointIndex, MassType& m)
+    {
+        applyVertexMassDestruction(pointIndex, m);
+    });
 
     // add the functions to handle topology changes for Edge informations
-    m_edgeMassHandler = new EdgeMassHandler(this, &d_edgeMass);
-    d_edgeMass.createTopologyHandler(m_topology, m_edgeMassHandler);
-
+    d_edgeMass.createTopologyHandler(m_topology);
+    d_edgeMass.setCreationCallback([this](Index edgeIndex, MassType& EdgeMass,
+        const core::topology::BaseMeshTopology::Edge& edge,
+        const sofa::type::vector< Index >& ancestors,
+        const sofa::type::vector< double >& coefs)
+    {
+        applyEdgeMassCreation(edgeIndex, EdgeMass, edge, ancestors, coefs);
+    });
+    d_edgeMass.setDestructionCallback([this](Index edgeIndex, MassType& m)
+    {
+        applyEdgeMassDestruction(edgeIndex, m);
+    });
 
     // register engines to the corresponding toplogy containers depending on current topology type
     bool hasTriangles = false;
@@ -1312,13 +1131,54 @@ void MeshMatrixMass<DataTypes, MassType>::initTopologyHandlers(sofa::core::topol
     if (topologyType == TopologyElementType::HEXAHEDRON)
     {
         d_vertexMass.linkToHexahedronDataArray();
-        d_edgeMass.linkToHexahedronDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::HexahedraAdded* hAdd = static_cast<const core::topology::HexahedraAdded*>(eventTopo);
+            applyVertexMassHexahedronCreation(hAdd->getIndexArray(), hAdd->getElementArray(), hAdd->ancestorsList, hAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::HexahedraRemoved* hRemove = static_cast<const core::topology::HexahedraRemoved*>(eventTopo);
+            applyVertexMassHexahedronDestruction(hRemove->getArray());
+        });
+
+        if (!this->isLumped())
+        {
+            d_edgeMass.linkToHexahedronDataArray();
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::HexahedraAdded* hAdd = static_cast<const core::topology::HexahedraAdded*>(eventTopo);
+                applyEdgeMassHexahedronCreation(hAdd->getIndexArray(), hAdd->getElementArray(), hAdd->ancestorsList, hAdd->coefs);
+            });
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::HEXAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::HexahedraRemoved* hRemove = static_cast<const core::topology::HexahedraRemoved*>(eventTopo);
+                applyEdgeMassHexahedronDestruction(hRemove->getArray());
+            });
+        }
+
         hasQuads = true; // hexahedron imply quads
     }
     else if (topologyType == TopologyElementType::TETRAHEDRON)
     {
         d_vertexMass.linkToTetrahedronDataArray();
-        d_edgeMass.linkToTetrahedronDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TetrahedraAdded* tAdd = static_cast<const core::topology::TetrahedraAdded*>(eventTopo);
+            applyVertexMassTetrahedronCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TetrahedraRemoved* tRemove = static_cast<const core::topology::TetrahedraRemoved*>(eventTopo);
+            applyVertexMassTetrahedronDestruction(tRemove->getArray());
+        });
+
+        if (!this->isLumped())
+        {
+            d_edgeMass.linkToTetrahedronDataArray();
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAADDED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::TetrahedraAdded* tAdd = static_cast<const core::topology::TetrahedraAdded*>(eventTopo);
+                applyEdgeMassTetrahedronCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+            });
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TETRAHEDRAREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::TetrahedraRemoved* tRemove = static_cast<const core::topology::TetrahedraRemoved*>(eventTopo);
+                applyEdgeMassTetrahedronDestruction(tRemove->getArray());
+            });
+        }
 
         hasTriangles = true; // Tetrahedron imply triangles
     }
@@ -1326,17 +1186,54 @@ void MeshMatrixMass<DataTypes, MassType>::initTopologyHandlers(sofa::core::topol
     if (topologyType == TopologyElementType::QUAD || hasQuads)
     {
         d_vertexMass.linkToQuadDataArray();
-        d_edgeMass.linkToQuadDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::QuadsAdded* qAdd = static_cast<const core::topology::QuadsAdded*>(eventTopo);
+            applyVertexMassQuadCreation(qAdd->getIndexArray(), qAdd->getElementArray(), qAdd->ancestorsList, qAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::QuadsRemoved* qRemove = static_cast<const core::topology::QuadsRemoved*>(eventTopo);
+            applyVertexMassQuadDestruction(qRemove->getArray());
+        });
+
+        if (!this->isLumped())
+        {
+            d_edgeMass.linkToQuadDataArray();
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSADDED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::QuadsAdded* qAdd = static_cast<const core::topology::QuadsAdded*>(eventTopo);
+                applyEdgeMassQuadCreation(qAdd->getIndexArray(), qAdd->getElementArray(), qAdd->ancestorsList, qAdd->coefs);
+            });
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::QUADSREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::QuadsRemoved* qRemove = static_cast<const core::topology::QuadsRemoved*>(eventTopo);
+                applyEdgeMassQuadDestruction(qRemove->getArray());
+            });
+        }
     }
 
     if (topologyType == TopologyElementType::TRIANGLE || hasTriangles)
     {
         d_vertexMass.linkToTriangleDataArray();
-        d_edgeMass.linkToTriangleDataArray();
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESADDED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TrianglesAdded* tAdd = static_cast<const core::topology::TrianglesAdded*>(eventTopo);
+            applyVertexMassTriangleCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+        });
+        d_vertexMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+            const core::topology::TrianglesRemoved* tRemove = static_cast<const core::topology::TrianglesRemoved*>(eventTopo);
+            applyVertexMassTriangleDestruction(tRemove->getArray());
+        });
+
+        if (!this->isLumped())
+        {
+            d_edgeMass.linkToTriangleDataArray();
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESADDED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::TrianglesAdded* tAdd = static_cast<const core::topology::TrianglesAdded*>(eventTopo);
+                applyEdgeMassTriangleCreation(tAdd->getIndexArray(), tAdd->getElementArray(), tAdd->ancestorsList, tAdd->coefs);
+            });
+            d_edgeMass.addTopologyEventCallBack(sofa::core::topology::TopologyChangeType::TRIANGLESREMOVED, [this](const core::topology::TopologyChange* eventTopo) {
+                const core::topology::TrianglesRemoved* tRemove = static_cast<const core::topology::TrianglesRemoved*>(eventTopo);
+                applyEdgeMassTriangleDestruction(tRemove->getArray());
+            });
+        }
     }
-
-    // PointData need to be linked to Edge container in any topology. d_edgeMass as EdgeData is automatically register to Edge container
-    d_vertexMass.linkToEdgeDataArray();
 }
 
 
@@ -1499,11 +1396,11 @@ void MeshMatrixMass<DataTypes, MassType>::computeMass()
 
     // set vertex tensor to 0
     for (Index i = 0; i<ndof; ++i)
-        m_vertexMassHandler->applyCreateFunction(i, my_vertexMassInfo[i], emptyAncestor, emptyCoefficient);
+        applyVertexMassCreation(i, my_vertexMassInfo[i], i, emptyAncestor, emptyCoefficient);
 
     // set edge tensor to 0
     for (Index i = 0; i<nbEdges; ++i)
-        m_edgeMassHandler->applyCreateFunction(i, my_edgeMassInfo[i], edges[i], emptyAncestor, emptyCoefficient);
+        applyEdgeMassCreation(i, my_edgeMassInfo[i], edges[i], emptyAncestor, emptyCoefficient);
 
     // Create mass matrix depending on current Topology:
     if (m_topology->getNbHexahedra()>0 && hexaGeo)  // Hexahedron topology
@@ -1518,10 +1415,10 @@ void MeshMatrixMass<DataTypes, MassType>::computeMass()
         m_massLumpingCoeff = 2.5;
         if (!isLumped())
         {
-            m_edgeMassHandler->applyHexahedronCreation(hexahedraAdded, m_topology->getHexahedra(), emptyAncestors, emptyCoefficients);
+            applyEdgeMassHexahedronCreation(hexahedraAdded, m_topology->getHexahedra(), emptyAncestors, emptyCoefficients);
         }
         
-        m_vertexMassHandler->applyHexahedronCreation(hexahedraAdded, m_topology->getHexahedra(), emptyAncestors, emptyCoefficients);
+        applyVertexMassHexahedronCreation(hexahedraAdded, m_topology->getHexahedra(), emptyAncestors, emptyCoefficients);
     }
     else if (m_topology->getNbTetrahedra()>0 && tetraGeo)  // Tetrahedron topology
     {
@@ -1536,10 +1433,10 @@ void MeshMatrixMass<DataTypes, MassType>::computeMass()
         m_massLumpingCoeff = 2.5;
         if (!isLumped())
         {
-            m_edgeMassHandler->applyTetrahedronCreation(tetrahedraAdded, m_topology->getTetrahedra(), emptyAncestors, emptyCoefficients);
+            applyEdgeMassTetrahedronCreation(tetrahedraAdded, m_topology->getTetrahedra(), emptyAncestors, emptyCoefficients);
         }
 
-        m_vertexMassHandler->applyTetrahedronCreation(tetrahedraAdded, m_topology->getTetrahedra(), emptyAncestors, emptyCoefficients);
+        applyVertexMassTetrahedronCreation(tetrahedraAdded, m_topology->getTetrahedra(), emptyAncestors, emptyCoefficients);
     }
     else if (m_topology->getNbQuads()>0 && quadGeo)  // Quad topology
     {
@@ -1554,10 +1451,10 @@ void MeshMatrixMass<DataTypes, MassType>::computeMass()
         m_massLumpingCoeff = 2.0;
         if (!isLumped())
         {
-            m_edgeMassHandler->applyQuadCreation(quadsAdded, m_topology->getQuads(), emptyAncestors, emptyCoefficients);
+            applyEdgeMassQuadCreation(quadsAdded, m_topology->getQuads(), emptyAncestors, emptyCoefficients);
         }
 
-        m_vertexMassHandler->applyQuadCreation(quadsAdded, m_topology->getQuads(), emptyAncestors, emptyCoefficients);
+        applyVertexMassQuadCreation(quadsAdded, m_topology->getQuads(), emptyAncestors, emptyCoefficients);
     }
     else if (m_topology->getNbTriangles()>0 && triangleGeo) // Triangle topology
     {
@@ -1572,10 +1469,10 @@ void MeshMatrixMass<DataTypes, MassType>::computeMass()
         m_massLumpingCoeff = 2.0;
         if (!isLumped())
         {
-            m_edgeMassHandler->applyTriangleCreation(trianglesAdded, m_topology->getTriangles(), emptyAncestors, emptyCoefficients);
+            applyEdgeMassTriangleCreation(trianglesAdded, m_topology->getTriangles(), emptyAncestors, emptyCoefficients);
         }
 
-        m_vertexMassHandler->applyTriangleCreation(trianglesAdded, m_topology->getTriangles(), emptyAncestors, emptyCoefficients);
+        applyVertexMassTriangleCreation(trianglesAdded, m_topology->getTriangles(), emptyAncestors, emptyCoefficients);
     }
 
 }
@@ -1629,7 +1526,7 @@ template <class DataTypes, class MassType>
 bool MeshMatrixMass<DataTypes, MassType>::checkTotalMass()
 {
     //Check for negative or null value, if wrongly set use the default value totalMass = 1.0
-    if(d_totalMass.getValue() <= 0.0)
+    if(d_totalMass.getValue() < 0.0)
     {
         msg_warning(this) << "totalMass data can not have a negative value.\n"
                           << "To remove this warning, you need to set a strictly positive value to the totalMass data";
@@ -1669,7 +1566,7 @@ bool MeshMatrixMass<DataTypes, MassType>::checkVertexMass()
         //Check that the vertexMass vector has only strictly positive values
         for(size_t i=0; i<vertexMass.size(); i++)
         {
-            if(vertexMass[i]<=0)
+            if(vertexMass[i]<0)
             {
                 msg_warning() << "Negative value of vertexMass vector: vertexMass[" << i << "] = " << vertexMass[i];
                 return false;
@@ -1726,7 +1623,7 @@ bool MeshMatrixMass<DataTypes, MassType>::checkEdgeMass()
         //Check that the vertexMass vector has only strictly positive values
         for(size_t i=0; i<edgeMass.size(); i++)
         {
-            if(edgeMass[i]<=0)
+            if(edgeMass[i]<0)
             {
                 msg_warning() << "Negative value of edgeMass vector: edgeMass[" << i << "] = " << edgeMass[i];
                 return false;
@@ -1833,7 +1730,7 @@ bool MeshMatrixMass<DataTypes, MassType>::checkMassDensity()
     if(massDensity.size() == 1)
     {
         //Check that the massDensity is strictly positive
-        if(density <= 0.0)
+        if(density < 0.0)
         {
             msg_warning() << "Negative value of massDensity: massDensity = " << density;
             return false;
@@ -1857,7 +1754,7 @@ bool MeshMatrixMass<DataTypes, MassType>::checkMassDensity()
         //Check that the massDensity has only strictly positive values
         for(size_t i=0; i<massDensity.size(); i++)
         {
-            if(massDensity[i]<=0)
+            if(massDensity[i]<0)
             {
                 msg_warning() << "Negative value of massDensity vector: massDensity[" << i << "] = " << massDensity[i];
                 return false;