[SOLVED] Measuring contact pressure

[twxu]

Hi,

Has anyone had any experience with measuring contact pressure in SOFA? i.e. when object A contacts object B, measure the pressure experienced by object B at the contact area.

I have been using the Monitor object so far to measure forces on individual nodes (which has been working very well), but is there some way to convert the force measurements to pressure measurements?

For example, could the pressure on a single triangular mesh element be calculated from the force measured at the 3 triangle corners? As far as I know, the Monitor and ExtraMonitor objects do not have an option to directly measure pressure.

Any guidance is appreciated, thanks!

[Hugo]

Hi @twxu

Which contact responsive method are you using:
– penalty method
– or Lagrange-multiplier based resolution ?

Depending on the method, I would point you out two different solution.
(If you are not sure about the method, could you just share your scene file?)

Best

Hugo

[twxu]

Hi @Hugo,

Thanks for the response.

I am not 100% sure, but I believe I am using the Lagrange-multiplier based method.

My scene and object files are here: https://github.com/thomaswxu/SOFA
I hope it is not too annoying to use; you will have to change some file paths/names for the objects to be loaded correctly.
(i.e. change “fea_variant” variable in scene file, then change corresponding file paths in “filenames.py”)

Please feel free to contact me if it doesn’t work or if you have any questions!

[Hugo]

Hey

From your scene singleFinger_wObject_XML.scn, it appears that you indeed use Lagrange multipliers: use of the FreeMotionAnimationLoop + GenericConstraintSolver.

Recovering contact forces can be done in the GenericConstraintSolver:
– activate the “computeConstaintForces” boolean
– this fills the nodal data field “constraintForces” that you can use to compute your pressure

I hope this helps.
Best

Hugo

[twxu]

Hi @Hugo,

Thanks for your suggestion! I have activated the “computeConstraintForces” boolean, and can see the “constraintForces” data under GenericConstraintSolver in the SOFA GUI.

But, I’m not sure how to interpret this data. It appears to be a single row vector; I see in the SOFA documentation that it represents intensities of constraint forces, but I could not find further information.

Could you provide more details on how to calculate contact pressure from the “constraintForces” data? I’m not sure how to proceed with the pressure calculation currently.

[Hugo]

Hi @twxu

“constraintForces” is a vector which size equals the number of constraints (maybe contacts in your case). If you have 3 contacts you would have :

“Contact1_ForceInX Contact1_ForceInY Contact1_ForceInZ Contact2_ForceInX Contact2_ForceInY Contact2_ForceInZ Contact3_ForceInX Contact3_ForceInY Contact3_ForceInZ”

Best,

Hugo

[Hugo]

To compute the pressure, use directly the P = F/S relationship.
Best,

Hugo

[twxu]

Hi @Hugo,

Thank you for the clarification; that makes sense.

How would you recommend computing the area term (S)? Is there a way in SOFA to calculate the area of individual mesh elements?
Or, would it be better to sum all the forces and divide by the entire area?

[Hugo]

Hi @twxu

Good question.
To get such geometrical information, your C++ code must define a pointer a SetGeometryAlgorithms. For instance, you can find such computation on triangles in the MeshMatrixMass.h code :

sofa::component::topology::TriangleSetGeometryAlgorithms<GeometricalTypes>* triangleGeo;

and its usage in MeshMatrixMass.inl :

triangleGeo->computeRestTriangleArea(triangleAdded[i])

I hope this helps.
Best

Hugo

[nicklj]

Hi @Hugo

Sorry for continuing to ask questions on this “SOLVED” thread as I think my question is highly relevant to this one.

According to your reply in this post, “constraintForces” is a vector which size equals the number of constraints, and should have the form “Contact1_ForceInX Contact1_ForceInY Contact1_ForceInZ…”.

I’ve checked the source code of “GenericConstraintSolver.cpp” and “ConstraintSolverImpl.h”. It seems this “constraintForces” is nothing but read out the lambda value for Lagrange multipliers, as: “current_cp->getF()” which is actually “current_cp->f.ptr()”.

As a result, I think your reply in this thread seems different from that you posted before, as in:

Getting forces generated by a collision

which lambda also needs to be multiplied by the constraint values in the collision model to retrieve the contact force, which seems more reasonable to me.

However, I’ve tried both, and in any way, I still cannot successfully extract the contact force correctly. The values seem to be more than an order lower than that I simulated in Abaqus, or measured in the experiment.

By the way, I used FreeMotionAnimationLoop and GenericConstraintSolver in my simulation.

Much appreciate if you may help to explain a bit more on this. I’ve been troubled by this problem for months…

[nicklj]

I think I may have a better understanding on the program now. lambda is the force in local contact space, and the first component is always the normal force. The “constraint” values in the “mechanicalObject” of the collision model, is actaully a transformation from local space to global coordinate.

Besides, the force obtained in constraint problem is actually dt*lambda. So in order to obtain the force in terms of unit Newton, we need to divided by dt.

Please kindly help to correct if I have any misunderstanding. Thanks.

[btt]

Hi @hugo, @nicklj,
Sorry but my comment is somehow disappeared eventhough I got reported from email that comment has been successfully summitted. I am interested in this topic as well. However, I have been struggling in either extraction of the vector of lambda or computing contact pressure as guided in this

Could anyone give me a little more detail regarding this procedure? or any other clarification would be appreciated.
Furthermore, I have the following questions:
1. let say constraint number 0 1 2 impact on node 55 then the total resulting force, so-called Fs[55], generated from collision (in x y z) on node 55 will be the summation of those calculated from constraint 0, 1, 2 right?
2. Secondly, the term F in the equation p=F/S is the summation of all contact forces Fs[i], where i are the impacted nodes, generated in the contact area, is it right?
Thank you for any help
best regards

[nicklj]

@btt
For question 1. Just search the “constraint” in the collision MechanicalObject, and find all the node 55 related rows, for each row, multiples the lambda value to get the force vector in global coordinate, then add all these values together to get the total reaction force related to node 55.

[btt]

Thank you @nicklj for the explanation, I did exactly what you described above however calculated values for reaction forces seem unreliable (too big). For example, the constraint from the collision model for node 8 is:

9 1 8 0.969057 -0.104965 -0.223409
10 1 8 0.0651265 0.981735 -0.17876  
11 1 8 0.238092 0.158679 0.958193 
...
15 1 8 0.971997 0.0886831 -0.217617  
16 1 8 -0.0778631 0.995286 0.0578188  
17 1 8 0.221718 -0.0392554 0.97432

Lambda values for row number 9 to 11 are:
[11293731.263678096], [-1448629.6310082236], [4136659.638746046]
whereas, those for row number 15 to 17 are zeros. Then the force for node 8 is
F = (constraint[9]*lamda[9]+constraint[10]*lamda[10]+constraint[10]*lamda[10])=[1.1835 -0.1951 0.1700] *10^7
The unit system i used is: mm, tonne, N
Even when I multiply this value with time step dt = 0.01 as indicated in sofa documentation, the final values are still unbelievably large. Do you have any comments on this issue? Or any other help will also be appreciated.
best

[btt]

Hi @hugo, @nicklj,
I am interested in this topic as well, I also want to extract the force generated by collision as pressure on contact area. I followed all related post in forum especially this one. From Sofa documentation, I understand the principle of Lagrange-multiplier based method is to create forces on two objects potentially in contact to create constraint (in this case is a non-holonomic constraint)between them. However, I have been struggling in either extraction of the vector of lambda or computing contact pressure as guided in this. Could anyone give me a little more detail regarding this procedure? or any other clarification would be appreciated.
Furthermore, I have the following questions:
1. let say constraint number 0 1 2 impact on node 55

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

then the total resulting force, so called Fs[55], generated from collision (in x y z) on node 55 will be the summation of those calculated from constraint 0, 1, 2 right?
2. Secondly, the term F in the equation p=F/S is the summation of all contact forces Fs[i], where i are the impacted nodes, generated in contact area, is it right?
Thank you for any help
best regards

[Author]

Posts