[FEM] Compute and draw von Mises stress

[alxbilger]

This PR introduces the computation and visualization of von Mises stress, based on the FEM. It consists in the local least-squares projection of stress values from quadrature points to nodes. The projection is local: the value at the node depends on the considered element, i.e. a global node has multiple values depending on its neighbour element. The result is a discontinuous stress field. It's on purpose! A discontinuity would indicate a problem of discretization. Without any problem, the field should appear continuous.

The theory

In FEM, the Cauchy stress is approximated in any point of an element by:

$$ \sigma_\text{approx}(\xi) = \sum_i N_i(\xi) \sigma_i $$

where $\sigma_i$ are the values at nodes, $N_i$ are the shape functions.

$\sigma_i$ are unknowns and we want to compute them. For that, we introduce the error function:

$$ \rho(\sigma_0, \sigma_1, ...) = \frac{1}{2} \int_\Omega( \sum_i N_i(\xi) \sigma_i - \sigma )^2 d\Omega $$

We want this error to be minimal. Its minimum is when its derivative is zero:

$$ \int_\Omega N_i(\xi) \left( \sum_i N_i(\xi) \sigma_i - \sigma \right) d\Omega = 0 $$

We rearrange:

$$ \sum_j \left( \int_\Omega N_i N_j d\Omega \right) \sigma_j = \int_\Omega N_i \sigma d\Omega $$

that we can write as $M \sigma_j = b$. It's a linear system to solve. The solution gives the stress values at the nodes.

The $M$ matrix is the Gram matrix (mentioned in the code). It's constant over time, so it's pre-computed. In order to solve the linear system, we store the inverse of this matrix. The algorithm simply consists in assembling $b$ and multiply it by $M^{-1}$.

Note: for the future, one can think about smoothing the field to make it continuous. There are two ways:

1. Average the stress values around each node. It can be weighted by the volume of the element.
2. Instead of solving a local system, we assemble a global system to project the stress values in a global manner. It means assembling the global Gram matrix, which is sparse. It would require a linear solver, or a pre-computation of its decomposition for $M^{-1} * b$.

sofa_video_r60_0005_compressed.mp4

For future PRs:

• DrawElementColoredMesh is not implemented for all types of elements.
• Only LinearSmallStrainFEMForceField derives from CauchyStressEvaluator. The next candidate to support stress evaluation is naturally CorotationalFEMForceField. And HyperelasticMaterial once [FEM.HyperElastic] Integration of hyperelasticity from Elasticity plugin to SOFA #6206 is merged.

sofa_video_r60_0007.mp4

---

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

---

Reviewers will merge this pull-request only if

• it builds with SUCCESS for all platforms on the CI.
• it does not generate new warnings.
• it does not generate new unit test failures.
• it does not generate new scene test failures.
• it does not break API compatibility.
• it is more than 1 week old (or has fast-merge label).