1 August 2019 at 17 h 32 min #14054
I was trying to find a way to read the forces generated by a collision of two objects. As I understood from this thread, externalForce and force vectors are used to access forces exerted on the node. I used keyboardControl.scn from SofaPython plugin to test it. The issue for me is that force vector of the Raptor has only zeros all the time, and externalForce is just empty. Could you explain to me this behavior? Is this the right way to access forces exerted on individual nodes of the mesh? Thank you in advance.
Vlad2 August 2019 at 10 h 18 min #14056
- Zykl.io UG
A few years back, I encountered the same problem as you described (force vector containing only zeros, externalForce vector just empty) when trying to extract contact forces from a SOFA scene.
I did manage to obtain contact forces between objects by directly tapping into the “raw” results that were produced by the SOFA solver component(s), but this naturally involved modifications to the source code of some SOFA core classes.
I can see if I could easily isolate the changes I made for this at least, if this helps.
With best regards,
Fabian7 August 2019 at 8 h 53 min #14069
Thank you very much for the tip! Is it sort of an expected behavior? If not, do you think it’s worth to open an issue on github?
Vlad8 August 2019 at 14 h 22 min #14073
I totally missed the last line of your answer. If you could share the changes you did that would help me a lot 🙂 Thanks!
Vlad13 August 2019 at 18 h 18 min #1410314 August 2019 at 16 h 14 min #14118
The example I am referring to (keyboardControl.scn) uses FreeMotionAnimationLoop. I found out that it is possible read forces if I change LCPConstraintSolver with GenericConstraintSolver (even though externalForce is empty anyway), but using LCPConstraintSolver force vector is just filled with zeros.
Vlad2 September 2019 at 15 h 37 min #14171
- SOFA Consortium
Using the LCPConstraintSolver means that you rely on the Lagrange Multiplier method to solve the constraints using a Gauss Seidel. The lambda (Lagrange Multiplier) resulting from the resolution actually corresponds to the collision forces. But recovering them is a bit tricky. Fasten your seatbelt.
In the MechanicalObject of your object, there is a data field named “constraint”. It should look like:
0 1 55 -0.172748 0.984895 0.0118137 1 1 55 -0.0161293 -0.0148211 0.99976 2 1 55 0.984834 0.172516 0.018446 … 9 2 58 0.487691 -0.305337 0.146266 59 0.333767 -0.208967 0.100102 10 2 58 0.336271 0.407063 -0.271458 59 0.230138 0.278586 -0.185781 11 2 58 0.0393245 0.305838 0.50733 59 0.026913 0.20931 0.347208 …
This constraint field can be octained from the MechanicalObject using:
const MatrixDeriv& constraints = l_mState->read(core::ConstMatrixDerivId::constraintJacobian())->getValue();
This corresponds to: , the matrix containing the constraint directions.
The first digit of each line corresponds to the constraint ID. Therefore one line corresponds to one constraint. Here the data “constraint” has 3 constraints per contact because it corresponds to friction contact.
To get the forces, you need to multiply each line (each constraint) with the associated lambda multiplier λ, resulting from the Gauss Seidel resolution. The vector of λ can be obtained from the ConstraintProblem of the simulation, as follows:
sofa::component::constraintset::ConstraintProblem* cp = l_constraintSolver->getConstraintProblem(); sofa::component::linearsolver::FullVector<double> lambdaVector = cp->f;
Let’s compute the force associated to the constraint number 9:
9 2 58 0.487691 -0.305337 0.146266 59 0.333767 -0.208967 0.100102
To make the multiplication, you need to get the ID saved in the column 2 (
constraint = 2here). This digits gives you the number of nodes in the collision mesh are impacted by this constraint. For each node impacted, you have a block corresponding to the node ID and its force in x,y,z [ID fx fy fz].
The constraint number 9 impacts 2 nodes: the node 48 and node 49. To compute the force resulting from the collision applied on each of these node, you need to compute:
F_collis += Deriv( 0.487691 -0.305337 0.146266) * lambda F_collis += Deriv( 0.333767 -0.208967 0.100102) * lambda
I hope this is clear enough. If not, I was thinking about creating a specific component to recover the collision forces. I could implement it.
Hugo16 September 2019 at 8 h 55 min #14239
You must be logged in to reply to this topic.