MeshMatrixMass

This component belongs to the category of Masses. In the dynamic equation (see Physics integration page), the mass density results from the first derivative in time of the momentum term. The MeshMatrixMass computes the integral of this mass density over the volume of the object geometry. To do so and for any given topology (triangles, quads, tetrahedra or hexahedra), the MeshMatrixMass integrates the mass density inside each elements and sums the mass matrix in the system matrix .

Volume integration

As detailed in the Physics integration page, the left hand side part of the linear momentum conservation equals . To integrate over the domain, its weak form will result in the mass matrix:

where are the test functions, which are basis functions ensuring the existence of a solution. Since no exact integration can be performed on a random domain , the MeshMatrixMass relies on the Finite Element Method (FEM) and accumulates the result of the integral over each finite element (triangles, quads, tetrahedra or hexahedra):

The FEM relies on simple geometries in which any field can be interpolated using shape functions (see FEM at a glance). Note that the same basis functions are chosen for both the test and the shape functions. The interpolation of the acceleration term thus gives:

By change of variables, the computation of the mass matrix results in solving the following integration of the shape functions in each element:

Case of a linear tetrahedron

In the case of a linear tetrahedron, the shape functions are:

By replacing the shape functions, we therefore obtain:

We can note that the matrix is symmetric. The integration in the reference (or parent) space can be numerically computed using a Gauss quadrature (or Gauss point integration). The resulting mass matrix is:

API

Depending on the type of LinearSolver used:

• for iterative solvers, the result of the multiplication between the mass matrix and an approximated solution is computed by the function:

template <class DataTypes, class MassType>
void MeshMatrixMass<DataTypes, MassType>::addMDx(const core::MechanicalParams*, DataVecDeriv& vres, const DataVecDeriv& vdx, SReal factor)
{
    const MassVector &vertexMass= d_vertexMassInfo.getValue();
    const MassVector &edgeMass= d_edgeMassInfo.getValue();

    helper::WriteAccessor< DataVecDeriv > res = vres;
    helper::ReadAccessor< DataVecDeriv > dx = vdx;

    size_t v0,v1,nbEdges=_topology->getNbEdges();

    for (unsigned int i=0; i<dx.size(); i++)
    {
        res[i] += dx[i] * vertexMass[i] * (Real)factor;
    }

    for (unsigned int j=0; j<nbEdges; ++j)
    {
        v0=_topology->getEdge(j)[0];
        v1=_topology->getEdge(j)[1];

        res[v0] += dx[v1] * edgeMass[j] * (Real)factor;
        res[v1] += dx[v0] * edgeMass[j] * (Real)factor;
    }
}

• for direct solvers, the mass matrix is built by the function:

template <class DataTypes, class MassType>
void MeshMatrixMass<DataTypes, MassType>::addMToMatrix(const core::MechanicalParams *mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
{
    const MassVector &vertexMass= d_vertexMassInfo.getValue();
    const MassVector &edgeMass= d_edgeMassInfo.getValue();

    size_t v0,v1,nbEdges=_topology->getNbEdges();

    const int N = defaulttype::DataTypeInfo<Deriv>::size();
    AddMToMatrixFunctor<Deriv,MassType> calc;
    sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
    sofa::defaulttype::BaseMatrix* mat = r.matrix;
    Real mFactor = (Real)mparams->mFactorIncludingRayleighDamping(this->rayleighMass.getValue());


    for (size_t i = 0; i < vertexMass.size(); i++)
    {
        calc(r.matrix, vertexMass[i], r.offset + N*i, mFactor);
    }

    for (size_t j = 0; j < nbEdges; ++j)
    {
        v0 = _topology->getEdge(j)[0];
        v1 = _topology->getEdge(j)[1];

        calc(r.matrix, edgeMass[j], r.offset + N*v0, r.offset + N*v1, mFactor);
        calc(r.matrix, edgeMass[j], r.offset + N*v1, r.offset + N*v0, mFactor);
    }
}

Data

The MeshMatrixMass can be initialized using two different input data:

• totalMass: corresponding to the total mass of the object, which will be distributed over its volume taking into account the geometry
• massDensity: corresponding to the mass density used for the integration detailed above

Note that using the optional data lumping, it is possible to simply the mass matrix by making it diagonal. This is called mass lumping and it consists in summing all mass values of a line on the diagonal. The DiagonalMass is an optimized version of this mass lumping approach. In case of a linear tetrahedron, if the data lumping is true, the (lumped) mass matrix becomes:

Use lumping with caution since it is a numerical approximation, thus decreasing the accuracy of the integration.

Usage

The MeshMatrixMass requires a MechanicalObject to store the degrees of freedom associated to the nodes, as well as a Topology. An integration scheme and a solver are also necessary to solve the linear system at each time step.

Several topologies are handled by the MeshMatrixMass, namely: triangles, quads, tetrahedra or hexahedra. Only the beam model (edge topology) is not handled by this component.

Example

This component is used as follows in XML format:

<MeshMatrixMass massDensity="1000" />

or using SofaPython3:

node.addObject('MeshMatrixMass', massDensity='1000')

An example scene involving a MeshMatrixMass is available in examples/Components/mass/MeshMatrixMass.scn