[SOLVED] Extracting Dynamics from Sofa

[jlorenze]

Hello,

I am interesting in trying to extract information about the dynamics of a soft robot from Sofa, similar to what is discussed in this paper. In particular I am interested in extracting the Jacobians dF/dq and dF/dv in equations (7)/(8).

In my simulation I am using an EulerImplicitSolver with a SparseLDLSolver linear solver and Rayleigh damping with Rayleigh mass/stiffness both set to 0.1. I then use the savingMatrixToFile method from the SparseLDLSolver class to save the A matrix at a time step. From the EulerImplicitSolver.cpp file:
matrix = MechanicalMatrix(1+tr*h*f_rayleighMass.getValue(),-tr*h,-tr*h*(h+f_rayleighStiffness.getValue()));

Therefore, it seems like I save the matrix = (1+h*f_rayleighMass)*M – h*B – h*(h+f_rayleighStiffness)*K. Since I have not added any damping in the force fields the B matrix should be 0. Therefore since I know my matrix M (using uniformMass), I can solve for K.

However I would appreciate some clarity on what the matrix K in this case actually is. In the documentation, K is the Jacobian df/dx, but its not clear exactly what the function f is in my case. For example consider the soft robot case where I have gravity, a RestShapeSpringsForceField, a TetrahedronFEMForceField, and a cable constraint from the SoftRobot library.
1. Is K the accumulated Jacobian of gravity, the ForceFields, and the cable constraint with respect to x? Or is it of just the ForceFields?

Additionally, I was wondering if some clarification about the Rayleigh damping could be provided. Specifically:
2. In the addition to b:
mop.addMBKv(b, -f_rayleighMass.getValue(), 1, h+f_rayleighStiffness.getValue())
why are there Rayleigh terms being added in?
3. Why in the definition of the A matrix is it (1+h*f_rayleighMass) and not (1-h*f_rayleighMass)? Isn’t Rayleigh damping where you have a damping term being added to the left hand side that is -C dv where C = f_rayleighMass*M + f_rayleighStiffness*K?
4. I am confused as to why the third argument is 1 and not 0? In the case where the Rayleigh mass/stiffness are both set to 0, the current code seems to set

b = h*f0 + h*B*v + h^2*K*v

but this seems contrary to the documentation.

[Hugo]

Hi @jlorenze

1. With a constraint resolution using the Lagrange multiplier method, f corresponds to ForceFields (including the gravity). It means any internal or external forces. Constraints would be in the H^T lambda part.

2. It comes from the fact that using Rayleigh damping assumes that df/dv (the damping matrix, here noted A) equals:

(your question made me notice a typo in the doc, thx!)

3. this was the typo error. You can also see the full Taylor expansion here.

4. It’s 1.0 because we want to be able to take into account physical damping (B, i.e. 1.0*B ). We actually compute here: b = f0 + 1.0*B*v + h*K*v
Then only, the line b.teq(h*tr); multiplies the entire term with h.

Do not hesitate to point out any other inconsistency in the doc. You can even make PR to make proposal of evolution in the markdown files.
I hope this helps.
Best

Hugo

[jlorenze]

Thanks!

[Author]

Posts