[FEM] Fix polar rotation of non-square elements

[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.