[SOLVED] Damping forces in SOFA

[jjcasmar]

Hello,

I have been looking at the DiagonalVelocityDampingForceField class to understand damping models in SOFA and I have some questions.

First of all, lets remind what simulation in computer graphics tries to solve.

$\left(M – h\frac{\partial f}{\partial v} – h^2\frac{\partial f}{\partial x} \right) \Delta v = h \left( f_0 + h \frac{\partial f}{\partial x} \right)$

$\frac{\partial f}{\partial v} = B$
$\frac{\partial f}{\partial x} = K$

Usually f can be derived from a potential ($f = -\frac{\partial \phi}{\partial x}$) but this is not necessary.

In computer graphics Ralyeigh damping model is usually used. This model says that $B = \alpha M + \beta K$ and completely ignores $\frac{\partial f}{\partial v}$.

My main question is: how sofa supports damping? Does it only allows rayleigh model?
Reading the DiagonalVelocityDampingForceField I see that it creates a force in each node that depends on the DOF velocity.

                unsigned nbDampingCoeff = dampingCoefficients.getValue().size();

                if( nbDampingCoeff )
                {
                    sofa::helper::WriteAccessor<DataVecDeriv> f(_f);
                    const VecDeriv& v = _v.getValue();

                    for(unsigned i=0; i<v.size();i++)
                        for(unsigned j=0; j<Deriv::total_size; j++)
                            if( i<nbDampingCoeff )
                                f[i][j] -= v[i][j]*dampingCoefficients.getValue()[i][j];
                            else
                                f[i][j] -= v[i][j]*dampingCoefficients.getValue().back()[j];

This is clear to me, as the velocitiy damping usually is $f = -b \cdot v\[i\]$. My main questions are:
why are you not computing the derivatives? It is clear the K is zero, but B isn’t.
how does rayligh damping model affects this forcefield? I have seen that the rayleigh stiffness parameter ($\beta$) is defined in the BaseForceField class. I am not sure if I am going to get rayleigh damping for every ForceField that I define.

I’m not sure how the pure virtual functions that I must define for each forcefield work.

                virtual void addForce (const core::MechanicalParams*, DataVecDeriv&, const DataVecCoord&, const DataVecDeriv&);

                virtual void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& d_df , const DataVecDeriv& d_dx);

                virtual void addKToMatrix(sofa::defaulttype::BaseMatrix *, SReal, unsigned int &){}

                virtual void addBToMatrix(sofa::defaulttype::BaseMatrix * mat, SReal bFact, unsigned int& offset);

I understand that addForce computes the forces but I’m not sure about the other :s
addForce: Computes $K \cdot x + B \cdot v$
addDForce: Am I taking the derivative in x or in v? Or is it just in x for the RHS n the first equation in this post?
addKToMatrix:
addBToMatrix: This two should add K and B to the system matrix (the LHS in the first equation in this post), but how do they work? Why are these methods empty in the DiagonalVelocitityDampingForceField? The base implementation in ForceField class for the addKToMatrix takes into account the rayleigh damping value

    if (r)
        addKToMatrix(r.matrix, mparams->kFactorIncludingRayleighDamping(rayleighStiffness.getValue()), r.offset);

and this is confusing me.

Can you please give some light in these topics?

Thank you!

[Hugo]

First, you are well describing the two different type of damping.
This might help other users!

In SOFA, the Rayleigh damping can be directly (and easily defined) by defining the Rayleigh mass or Rayleigh stiffness parameter in the integration scheme component (e.g. EulerImplicitSolver). Both parameters are set to zero by default.

To implement your own damping in a ForceField, you can follow what is done in addForce() in DiagonalVelocityDampingForceField.inl to computing damping based on the derivative of your function f.
$B = \frac{\partial f}{\partial v}$
However, in the case of DiagonalVelocityDampingForceField, it is an “explicit” damping. Assuming your are solving a system Ax = b, the damping is included in b since it uses the addForce function.

To implement your own implicit damping (based on an unknown velocity v(t+1)), you need to implement your damping in the addDForce function (as you can see on the commented code in DiagonalVelocityDampingForceField). This addDForce() function will include the damping in the A matrix of our linear system.

My explanation is valid if you are using an iterative linear solver (e.g. CG). If you are using a direct solver, replace addDForce by addBToMatrix.

This topic remains complex to explain and understand on a forum. For deeper explanation about this, we can propose you a specific training.

Best regards,

Hugo

[jjcasmar]

I think I have understood it.

Depending on the equation system solver method, the addDForce (for krylov methods) or the addBToMatrix (for direct methods) is called. Inside those functions is where I must defined how K (df/dx) and B (df/dv) are computed, isnt it?

[Hugo]

Hi Juan,

You have perfectly understood the principle!
Enjoy damping with SOFA 😉

Cheers,

Hugo

[Author]

Posts