- This topic has 9 replies, 5 voices, and was last updated 3 years, 5 months ago by
.
-
Posts
-
Hello everyone,
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.
Regards,
VladHello, Vlad!
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,
FabianHello, @faichele!
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?
Regards,
VladHi, @faichele!
I totally missed the last line of your answer. If you could share the changes you did that would help me a lot 🙂 Thanks!
Regards,
VladHi @hugo,
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.
Regards,
VladHi @bobiko
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.100102To make the multiplication, you need to get the ID saved in the column 2 (
constraint[8][1] = 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[58] += Deriv( 0.487691 -0.305337 0.146266) * lambda[9] F_collis[59] += Deriv( 0.333767 -0.208967 0.100102) * lambda[9]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.
Best wishes,Hugo
Hi, @hugo!
Thank you for clarification, at least, now the behavior is known to me. I just have to think it through now.
Grtz,
Vladhi @hugo,
I have tried to extract the force (or more specifically total force in the contact area) caused by the collision between two collision models using Lagrange Multipliers. As far as I understand from your guidance and from SOFA code of ‘constraintsolver’ and ‘mechanicalObject’ components, the constraint forces in ‘constraintsolver’ are lambda and the constraint values provided in ‘mechanicalmodel’ actually belong to matrix H in this equation: A*deltav=dt*H_transpose*lambda. So if it is true, I suppose that constraint values are results of the partial derivation calculation of holonomic/non-holonomic constraints with respect to xi. Then, from your above response, what do you mean by calculating derivation of H as shown below:
F_collis[node] += Deriv( Hx Hy Hz) * lambda[i]
Deriv( Hx Hy Hz) is derivation of H with respect to time t?
Furthermore, could you suggest to me a way do this calculation in SOFA code or maybe in Sofapython? thank you for your support
best regards,
Nhanhi @hugo, @nhnhan
Interested as well. So far, I can export either lambda in constraint solver and constraint forces, let call it f, in the collision model by using sofapython. I presume the contact force at a node will be equal to dt*lambda*f. But I also feel confused with the @hugo comment as well:F_collis[58] += Deriv( 0.487691 -0.305337 0.146266) * lambda[9] F_collis[59] += Deriv( 0.333767 -0.208967 0.100102) * lambda[9]Could you please clarify this in more detail? It is very interesting to me but I have been struggled for weeks to really understand that. Thank you
best
- You must be logged in to reply to this topic.
