[th-skam]

In short

Using the method polar to estimate the rotation $R$, triangular and tetrahedral elements produce spurious rotations and/or stretches. This is coming from failing tests in VnV for the CorotationalFEMForceField.

I added unit tests which are the examples below.

Longer

The PolarDecomposition.h extracts the rotation per element by approximating the deformation gradient $F$ as $H = Σ p_j q_jᵀ$, with $p_j$ and $q_j$ the current and rest node offsets from the centroid of the element. It then takes $R$ as the rotation of the polar decomposition of $H$, as if $H = F = R U$.

Generally, polar(F) is not equal to polar(H). They relate from $H = F M$ ($M = Σ q_j q_jᵀ$), in the case of affine node movement. Otherwise it's a least-squares best fit. So polar(H) = polar(F) only when $M$ is proportional to the identity, which is the case for squares and cubes but not for triangles and tetrahedra, and a pure stretch produces a rotation that does not exist. Our code does not account for this.

Examples as unit tests:

• If you rotate a square rigidly or only stretch it, SOFA computes a correct rotation. It does so for hexahedra as well.

• If you rotate a triangle, the polar decomposition will compute a correct $R$ but will also see stretches ($U$ =/= identity).

• If you only stretch a triangle (and that's the biggest error), the decomposition will return a rotation i.e. $R$ =/= identity.

[damienmarchal]

Thanks for the quality of the explanation.

[th-skam]

Originally, I treated this PR as a fix but it shouldn't be. We have more than one method for extracting the rotation matrix in SOFA, each one with its performances and limitations. I think it's better to keep them all.

I've made a PR with doc bringing it all together here:
sofa-framework/doc#230

The added unit tests (some of which are from the refs in the documentation) demonstrate the limitations of each method by expecting failures in some cases and success in others. For example triangle will fail extracting correct rotation from a stretched irregular triangle, or it will compute a different one based on the node numbering of the mesh. I've summarized them in this pdf:
rotation_methods_table.pdf
Blue cross is an expected failure.

[th-skam]

After what we talked about in the last dev meeting, I will need a bit more time to work on handling embedded elements (and document them) properly.
So I'd prefer to simplify this PR and keep it contained to full-dimensional ones.