[SOLVED] Many ‘nan’ values in EulerImplicitSolver log for TetrahedronFEMForceField

[pzingg]

I am trying to build my first real physics simulation using SofaPython3. I am trying to emulate the AttachConstraint.scn example in the SOFA repository. Instead of using the RegularGridTopology objects in that example, I create a Gmsh .msh file with hexahedron topology for each of three attached elements in the scene and load the meshes. This also requires me to add HexahedronSetTopologyContainer and HexahedronSetGeometryAlgorithms objects to the scene. Otherwise I tried to get as close to the original example as possible.

The only other issue I had in trying to replicate the scene is that there is no BoxConstraint object in SofaPython3, so I used a FixedConstraint instead to fix the first of three elements. As in the “TwoWay” example in the original AtttachConstraint.scn file, I build AttachConstraints between the first and second elements and between the second and third elements. The forces are set using TetrahedronFEMForceFields in the three elements as in the original.

I also have the “printLog” attribute set on the EulerImplicitSolver.

Now for the real problem:

The scene looks good (and “checkTopology” is set to “1” in the loader and no error messages are produced). But when I step or animate, I can see in the console that only the vertices in the second of my three elements have initial forces set, the ones in the first and third elements are all -nan values. Needless to say, that causes the simulation to fail.

How would I go about finding out what’s wrong with my scene? Could it be that the loaded mesh is somehow not right topologically? Or do I need some extra mapping between the hexahedral mesh and the TetrahedronFEMForceField?

I don’t understand why one of the three elements seems to have valid gravity force values on it, but the first and third ones are all -nan. The mesh files are created the same way. Also, the three meshes all load fine in Gmsh, but maybe SOFA is stricter. Is there a simple way to debug the initial setup of the force field that might show me what I’m doing wrong?

[pzingg]

Found out that my meshes were fine, but two of the elements were initialized with a poissonRatio of 0.5.

By using gdb and a debug version of SOFA, I found out that in the TetrahedronFEMForceField::computeMaterialStiffness() function that sets up the material stiffnesses for the force field, there is a calculation with a denominator of

((1+poissonRatio)*(1-2*poissonRatio))

This results in setting nan values for the DoF stiffnesses, which propogates to nan values for the forces, etc.

I am now using Poisson ratio values of 0.4 and 0.45, not 0.5, and everything is working as expected.

[Hugo]

Hi @pzingg

Thanks for writing us the solution to your problem.
The Poisson ratio must be <0.5 (strictly) for any elastic material. The value of 0.5 is a theoretically value corresponding to perfect imcompressibility of the material.

Did you have a warning from the TetrahedronFEMForceField at the initialization of the scene about your poissonRatio=0.5 ? If no, it would be nice to PR such a check. Would you mind doing it?

Best wishes,

Hugo

[Author]

Posts