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 Ax=b. 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 A:

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 x of your system. This force will contribute to the righthand side b:
Ax = b = F_{ext}
In any forcefield, the vector b 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
Looking at continuum mechanics, the linear system Ax=b becomes:
$$M ddot{x} = K(x)$$
where K may depend on the state x. K can either be linear or nonlinear regarding x. In case K is nonlinear, the resulting K(x) must be recomputed at each time step. The choice of the space integration used for the physical law will determine how the matrix K is filled. 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 global K matrix.
The contribution of the physical law in the linear system will depend on the integration scheme:

with explicit scheme, we have K(x) = K(x(t)) where x(t) is the known current state. The initial equation becomes:
Mdv = dt ⋅ K(x(t))
The physical law therefore only contributes to the righthand side b of the linear system through the function:addForce() // corresponding to the term : dt K(x(t))

with the implicit scheme, we have
$$K(x) = K(x(t+dt)) = K(x(t)) + frac{dK(x(t+dt))}{dx} dx$$
where x(t+dt) is the unknown current state and the initial equation becomes:
$$M dv = dt cdot (K(x(t)) + frac{dK(x(t+dt))}{dx})$$$$left(M – dt cdot frac{dK(x(t+dt))}{dx}right) dv= dt cdot K(x(t))$$
Therefore, the implicit scheme will also implement the function:addForce() // corresponding to the term : dt K(x(t))
and the part of the lefthand side depending on the unknown state x(t+dt) is implemented in the function:
addKToMatrix() // corresponding to the term :  dt dK(x(t+dt))/dx
Last modified: 14 November 2017