TD8 : Système masse-ressort

L'objectif de ce TD est d'implémenter un système masse-ressort pour représenter, par exemple, un tissu. N'hésitez pas à vous référer au cours Siggraph d'Andrew Witkin et David Baraff.

1. Téléchargez le nouveau squelette de code suivant : sia_td8.zip. Analysez l'architecture du programme, en particulier les nouvelles classes dans les fichiers integration.h et particules.h.

2. Complétez les méthodes addForces pour les 3 classes GravityForce, DragForce et SpringForce.

3. Complétez ensuite les méthodes getDimensions, get/setState et getDerivative de la classe ParticleSystem. La méthode segment des vecteurs d'Eigen vous simplifiera la tâche.

4. Complétez enfin la fonction explicitEulerStep qui réalisera une étape du schéma d'intégration d'Euler. Implémentez ParticleSystem::step et n'oubliez pas d'appeler ParticleSystem::updateGL pour mettre à jour l'affichage. Vous pouvez alors tester votre simulateur physique. Pour interagir avec la grille, vous pouvez ajouter une nouvelle force d'accroche temporaire liée à la position de la souris avec CTRL + clic.

5. Expérimentez avec différents pas de temps en modifiant l'appel à _psys.step dans Viewer::updateScene(). Jouez également avec les différents paramètres des ressort.

6. Pour augmenter la stabilité du tissu et éviter les cisaillement, ajoutez des ressorts le long de la diagonale des cellules de la grille. Quel coefficient de rigidité $k_s$ faut-il choisir ?

7. Implémentez la méthode du point milieu. Observez son effet sur la stabilité et les performances en faisant varier le pas de temps.

8. Ajoutez un plan pour représenter le sol de la scène. Je vous propose de créer une classe abstraite Collider qui représentera tous les objets (dont la trajectoire n'est pas simulée) pouvant provoquer des collisions, et une classe PlaneCollider qui en dérivera. Vous pouvez le visualiser ce plan en 3D avec une instance de Quad. Complétez la classe PlaneCollider pour détecter et réponde aux collisions entre les masses et ce plan.

9. Chargez le maillage 'cloth.obj' représentant un vêtement. Utilisez-le pour instancier un nouveau réseau de masses-ressorts remplaçant la grille. Plutôt que les ressorts, affichez le maillage du vêtement avec le modèle d'éclairage de Blinn-Phong. Réalisez le couplage entre le modèle physique et le modèle visuel en mettant à jour les positions du maillage après la simulation.

  

10. Chargez le maillage du mannequin 'girl.obj' et affichez-le. Essayez de gérer les collisions entre le vêtement et ce maillage. Exploitez la BVH pour accélérer les tests de collision.

Bonus : Il existe une méthode d'intégration numérique stable quelle que soit la rigidité du système : la méthode d'Euler implicite (ou backward Euler). Pour en comprendre le principe, lisez tout d'abord les sections 3 et 4 des notes de Witkin et Baraff.

Vous allez désormais implémentez un tel solveur ; à chaque itération, vous devrez :

  1. obtenir le vecteur d'état $(\vx{x_0},\vx{v_0})$ ainsi que les matrice des masses $M$, des forces $\vx{f}_0$ et de leurs dérivées (Jacobiens) $J_\vx{x} = \partial \vx{f} / \partial \vx{x}$ et $J_\vx{v} = \partial \vx{f} / \partial \vx{v}$ depuis le système de particules,
  2. calculer la matrice $A = M-J_\vx{x} \Delta t^2 - J_\vx{v}\Delta t$ et le vecteur $\vx{b} = (\vx{f}_0 + J_\vx{x} \vx{v}_0 \Delta t)\Delta t$,
  3. résoudre le système linéaire $A\Delta \vx{v} = \vx{b}$ en utilisant la fonction solve() fournie ci-dessous pour déterminer $\Delta \vx{v}$,
  4. mettre à jour l'état du système $\vx{v}_1 = \vx{v}_0 + \Delta \vx{v}$ et $\vx{x}_1 = \vx{x}_0 + \vx{v}_1 \Delta t$.

Quitte à être moins générique, modifiez et ajoutez des méthodes à ODESystem pour calculer et transmettre les valeurs nécessaires au solveur. Vous devrez notamment calculer les Jacobiens pour les ressorts, en vous aidant de cette page.

#include <Eigen/Sparse>
#include <Eigen/IterativeLinearSolvers>

VectorXd solve(const MatrixXd &A, const VectorXd &b) {
    SparseMatrix<double> spA = A.sparseView();
    ConjugateGradient< SparseMatrix<double> > solver;
    return solver.compute(spA).solve(b);
}

Notez que cette méthode d'intégration n'est garantie d'être stable que si la matrice des Jacobiens est définie négative, ce qui n'est pas forcément le cas lorsque les ressorts sont compressés. Un hack simple pour éviter ce problème est de clamper la distance $1 - r / ||\vx{x}_ij||$ à zéro.