Physics in SOFA
In the previous doc pages, the integration scheme describes how to compute the configuration at the next time step from the current state. The linear solver explains the step performed to compute the solution of the linear system . This section explains how the physics contributes to this system.
Mass
The Mass of the system will contribute to the lefthand side within the matrix :

with direct solvers, the mass is included in the matrix using the function:
addMToMatrix()

with iterative solvers, the mass is taken into account through the function:
addMDx()
There is different way of integrating the mass in the system, described below.
UniformMass
This mass is independent of your modeling choices. The total mass is divided by the number of nodes, the mass matrix is diagonal and each value on the diagonal has the same value:
for ( unsigned int i=0; i<indices.size(); i++ )
res[indices[i]] += dx[indices[i]] * m; // m is constant
MeshMatrixMass
This mass integrates the mass density within the topology based on the Finite Element Method (FEM). Thus, the mass is spread on the nodes and edges of the topology. A large element will therefore have a higher associated mass due to its volume:
for (unsigned int i=0; i<dx.size(); i++)
{
res[i] += dx[i] * vertexMass[i] * (Real)factor;
massTotal += vertexMass[i] * (Real)factor;
}
for (unsigned int j=0; j<nbEdges; ++j)
{
tempMass = edgeMass[j] * (Real)factor;
v0=_topology>getEdge(j)[0];
v1=_topology>getEdge(j)[1];
res[v0] += dx[v1] * tempMass;
res[v1] += dx[v0] * tempMass;
massTotal += 2*edgeMass[j] * (Real)factor;
}
DiagonalMass
The diagonal mass is a simplification of the MeshMatrixMass approach. It integrates the mass density within the topology based on the Finite Element Method (FEM). However, the mass is moved only on integration points (only mass on the nodes, not on edges anymore): this is a numerical lumping. The mass matrix thus becomes diagonal:
for (size_t i=0; i<n; i++)
{
_res[i] += _dx[i] * masses[i]; // mass different on each node, depending on the topology
}
External forces
Forces can be applied on your physical object. Usually forces are sorted into external and internal forces. Let’s consider a simple external force independent from the state of your system (i.e. purely explicit force). This force will contribute to the righthand side :
In any forcefield, the vector is filled in the function:
addForce()
Taking the example of the ConstantForceField:
addForce(const core::MechanicalParams* params, DataVecDeriv& force, const DataVecCoord& position, const DataVecDeriv&)
{
sofa::helper::WriteAccessor< core::objectmodel::Data< VecDeriv > > _force = force;
for (unsigned int i=0; i<position.getValue().size(); i++)
_force[i] += constantForce; // constant value filling the b vector
}
Physical laws (internal forces)
Looking at continuum mechanics, the linear system arises from the dynamic equation whichis a timevarying partial differential equation which, after discretization, can be numerically solved as an ordinary differential equation. The system can then be written:
where is the degrees of freedom, the mass matrix, a scalar function of x—yields the material internal energy, and a function of and , describing other forces (constraint forces, internal damping, etc.) acting on our system. The derivative of the internal energy will lead to the computation of internal forces notes .
The contribution of the physical law in the linear system will depend on the integration scheme. For this explanation, we will rely on the explicit (or forward) Euler and the implicit (or backward) Euler scheme:
Explicit case
In case of an explicit integration, the above system becomes:
The physical law therefore only contributes to the righthand side of the linear system through the function:
addForce() // corresponding to the term : dt f(x(t))
SOFA is mainly based on the Finite Element Method to integrate in space the physical law, i.e. the contribution of each element of the mesh will be added to the numerical system (here the righthand side vector ).
Implicit case
In case of an implicit integration, the above system becomes:
In this equation, we can notice the same explicit contribution . Just like in the explicit case, this part is implemented in the same function addForce(). We can also notice the appearance of the stiffness matrix : . The stiffness matrix is a symetric matrix, can either be linear or nonlinear regarding .
The term resulting from the Taylor expansion with can be decomposed into one explicit term contributing to the righthand side vector and one implicit term contributing to the system matrix on the lefthand side.

for direct solvers, the addKToMatrix() function computes the implicit part of the stiffness which will build a matrix . The linear system matrix must be stored and built, since it will be inversed later on to solve the system.
addKToMatrix() // adding contribution in the A matrix :  dt^2 df(x(t+dt))/dx
The addDForce() function implements the explicit part of the stiffness by adding the result of the matrixvector multiplication directly in the righthand side vector :
addDForce() // adding contribution in b :  dt^2 df(x(t+dt))/dx

for iterative solvers, since the system matrix does not required to be built, both explicit and implicit contributions of the stiffness will be computed by the addDForce(): for the explicit part it compultes the matrixvector multiplication between and the known first derivate of the degrees of freedom, while for the implicit part it computes the matrixvector multiplication between and the unknown :
addDForce() // multiplying any vector by :  dt^2 df(x(t+dt))/dx
Last modified: 5 October 2018