[SOLVED] Damping forces in SOFA

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!