Implémentation d'un moteur physique avec des contraintes rigides

Pierre Bénet

Dans un jeu vidéo, les simulations physiques jouent un rôle très particulier : grâce à elles, les mouvements à l'écran ont l'air plausibles, voire réalistes. Par exemple, ces simulations calculent comment la pomme tombera de l'arbre (au lieu de la voir rester pendue).

Cet article présente la simulation dynamique de solides rigides soumis à des contraintes (rotule, pivot…). Nous trouverons les équations générales qui régissent le mouvement de solides rigides sous contrainte. Nous utiliserons une méthode d'intégration précise (RK4) et nous présenterons un algorithme spécialement adapté à notre problème pour résoudre les systèmes linéaires. Nous aborderons aussi le calcul de cinématique inverse, utile en robotique.

Cet article fait appel à des notions de niveau licence ou classe préparatoire : calcul matriciel par blocs, calcul de dérivées partielles sur des vecteurs, des connaissances en physique comme le travail des forces et les lois de Newton.

Commentez

- Titre : Implémentation d'un moteur physique avec des contraintes rigides

- Auteur : Pierre Bénet

- Parution : 6 février 2017

- Mise à jour : 8 février 2017

- Public visé :

- Licence : Les sources présentées sur cette page sont libres de droits et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une œuvre intellectuelle protégée par les droits d'auteur. Copyright ® 2016 Pierre Bénet. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- I. Introduction

- II. Équations de la dynamique

- II-A. Mouvement d'un solide rigide

- II-B. Pour simplifier les équations

- II-C. Principe de D'Alembert

- II-D. Exemple : la liaison rotule

- III. Intégration

- III-A. Réduction au premier ordre

- III-B. Méthode d'Euler

- III-C. Dégénérescence des rotations

- III-D. Méthode de Runge Kutta d'ordre 4

- III-E. Méthodes d'Adams-Bashforth

- III-F. Précision et convergence

- IV. Correction des erreurs numériques

- IV-A. Contraintes arborescentes

- IV-B. Force fictive

- IV-C. Newton-Raphson

- IV-D. Impulsion

- IV-E. Cinématique inverse

- V. Résolution du système linéaire

- V-A. Propriété du problème

- V-B. Algorithme des corps articulés

- V-C. Élargissement de l'algorithme

- V-D. Pivot de Gauss creux

- V-E. Implémentation du pivot de Gauss creux

- VI. Dynamique inverse

- VII. Structures de données

- VII-A. Articulation

- VII-B. Solide

- VII-C. solveur

- VIII. Sources

- IX. Remerciements

I - Introduction

La simulation de corps rigides dans un moteur physique peut être usuellement découpée en plusieurs tâches, qui s'effectuent à tour de rôle dans une boucle itérative. Il y a, par exemple, la détection de collision, ensuite viennent la gestion des chocs (le calcul des impulsions) et enfin l'intégration du mouvement.

Il existe différentes approches pour réaliser un moteur physique. Il est important de connaître les technologies récentes pour avoir un moteur physique efficace. Voici un article qui regroupe ces différentes technologies :Interactive Simulation of Rigid Body Dynamics in Computer Graphics.

Dans notre cas, nous nous intéresserons à la gestion des contraintes rigides dans une telle simulation. Par exemple, les liaisons rotule et pivot sont des contraintes. Dans les jeux vidéo, cela sert à animer un corps : notamment, on peut utiliser une contrainte rotule entre le bras d'un personnage et le buste. Une roue de voiture ne se détache pas de la voiture, mais peut tourner, c'est une liaison pivot.

Liaisons rotule, pivot et pivot glissant (source : ODE user-guide).

Les contraintes réduisent la mobilité des corps entre eux. Ce sont des équations des positions et vitesses des corps et aussi du temps, qui peuvent être des contraintes d'égalité ou d'inégalité. Ici, nous ne traiterons que les contraintes d'égalité qui dépendent des positions, ce sont les contraintes holonomes (rotules, pivots, glissières…). Par contre, une roue de voiture qui roule sur le sol sans déraper est une contrainte non holonome (elle est connue sous le nom de « contrainte de roulement sans glissement ») : la vitesse instantanée du point de contact entre la roue et le sol doit être nulle, la contrainte n'est pas en position, mais en vitesse.

Toute la difficulté du problème est de trouver les forces qui agissent entre les solides au niveau des contraintes. Ces forces sont interdépendantes, ce qui nous oblige à prendre en compte leurs interactions. Ici, l'exemple que nous utiliserons est un bonhomme de bois composé de 14 liaisons rotules parfaites, soumis à sa propre inertie en apesanteur.

Corps articulé.

Corps soumis à sa propre inertie, les segments sont représentés en jaune ; les forces sont représentées en vert.

Dans notre cas, nous attacherons une importance à l'efficacité et la précision de la simulation des contraintes. Ces critères sont importants dans l'ingénierie, par exemple le prototypage de mécanisme, la simulation industrielle et la robotique. Au contraire, dans le cas d'un jeu, la précision est un critère moins important : la performance domine les choix d'implémentation.

Après avoir défini les équations de la dynamique avec contraintes, les méthodes d'intégration, les méthodes de correction (pour les contraintes bouclées, notamment) et la manière de traiter les impulsions dans un système avec contraintes, voilà ce que nous pourrons simuler :

Animation en format "mp4" non supportée.
Contenu uniquement consultable en ligne à cette adresse :
https://pierre-benet.developpez.com/moteur_physique/simulation2.mp4

Simulation d'un corps rigide en apesanteur : le bonhomme a attrapé son pied, il rebondit sur une barrière invisible (rebond adhérent avec un coefficient de restitution de 75 %). N'ayant pas de détection de collision, le bonhomme s'interpénètre et ne ressemble plus à grand-chose.

II - Équations de la dynamique

Le but ici est de calculer l'accélération en fonction des autres variables de notre système.

II-A - Mouvement d'un solide rigide

Un solide rigide est un solide indéformable, c'est-à-dire que la distance entre deux points appartenant à un même solide ne varie pas. La dynamique du système devient alors beaucoup plus simple : il suffit de définir la position d'un point du solide et l'orientation du solide pour que la position de tous les points du solide soit définie.

Solide rigide ou indéformable

Les lois de Newton nous donnent l'équation de la conservation de quantité de mouvement et du moment cinétique pour un solide rigide.

et

avecla masse du solide, la position du centre de gravité du solide etles forces appliquées au solide, la matrice d'inertie du solide qui dépend de l'orientation (une matrice orthogonale),la vitesse de rotation du solide et les couples appliqués au solide.

II-B - Pour simplifier les équations

Plus généralement, nous aurons une équation de la forme :

, avec

Il faut passer par cette formulation, car la vitesse de rotation (dans) n'est pas directement la dérivée de la matrice de rotation (dans), il y a seulement une relation linaireentre les deux. Ici, est la vitesse etla position généralisée. est la masse, une matricediagonale par bloc de trois dans notre cas et est inversible. Réécrivons l'équation de la dynamique pour faire apparaître les variables :

est simplement composé des vitesses de rotation et vitesses linéaires. Remarquons que, contrairement à, nous n'avons pas encore défini la position. Dans tous nos calculs, nous allons nous passer de la variable position, puis nous la définirons au moment de l'intégration.

Nous avons aussicontraintes :

.

Dans notre cas, nous ne considérerons que les forces de contraintes (qui maintiennent le pantin en forme, par exemple) et pas les forces extérieures (comme la gravité). Il faut donc calculer ces forces issues des contraintes. Afin d'avoir une expression simple et générale (nous pourrons ainsi modéliser des contraintes exotiques facilement, par exemple définir comme contrainte x=cos(y) sur les coordonnées d'un point du solide).

II-C - Principe de D'Alembert

Nous allons utiliser le principe de D'Alembert, qui dit que les forces de contraintes ne travaillent pas lors d'un déplacement virtuel. Un déplacement virtuelest un déplacement instantané et infinitésimal du système de telle sorte qu'il vérifie toujours ses contraintes. Pour une contrainte holonome, on doit donc avoir. En développant au premier ordre, on obtient une équation linéaire que doit vérifier tout déplacement virtuel :, c'est-à-dire . La matrice est appelée jacobienne de la contrainte. Elle possède lignes (autant que de contraintes) et colonnes (autant que de variables), de sorte qu'elle puisse être multipliée à droite par un déplacement et à gauche par un vecteur deéléments que nous nommerons plus tard multiplicateur de Lagrange.

De plus, le travail des forces de contraintes doit être nul pour tout déplacement virtuel, soit .

Alors, l'ensemble des forces vérifiant cette propriété peut s'écrire sous la forme. En effet, une telle force vérifie toujours, pour tout déplacement virtuel, , car tout déplacement virtuel vérifie. De plus, une analyse dimensionnelle nous montre que cet ensemble est exactement celui que l'on cherche : .

Étant donné que la matrice est plus compliquée à calculer que la jacobienne, nous allons chercher la force de contrainte de la forme : cette manière de procéder est équivalente, car la matriceest inversible.

appelé multiplicateur de Lagrange, est un vecteur deéléments, un élément par contrainte. C'est en quelque sorte un coefficient de proportionnalité vectoriel. Nous pouvons aussi retrouver nos résultats en appliquant les équations de Lagrange. Cependant, les changements de variables entre vitesse et position deviennent plus difficiles.

Nous avons l'espace dans lequel les forces se trouvent, il nous faut encore les calculer. Les inconnues sont l'accélération (variables) et le multiplicateur de Lagrange (variables). Il faut donc considérer les équations de la dynamique et les équations des contraintes pour résoudre le système (équations). Nous dérivons l'équation des contraintes deux fois par rapport au temps afin de faire apparaître l'accélération :

Nous substituons alors l'équation de la dynamique dans l'équation de la contrainte pour faire disparaître l'accélération, nous obtenons alors le système d'équations suivant :

La résolution s'effectuera en calculant d'aborddans la deuxième équation, puis en injectant la valeur detrouvée dans l'équation de la dynamique pour calculer l'accélération. Le résultat est une fonction qui nous donne l'accélération à partir de la position et la vitesse :

II-D - Exemple : la liaison rotule

Schéma de deux solides liés par une liaison rotule.

Liaison rotule (source : ODE user-guide).

En mécanique, une liaison entre deux solides est définie par le nombre de degrés de liberté de la liaison. Un solide possède 6 degrés de liberté : 3 degrés de liberté en rotation et 3 en translation. La liaison rotule supprime les 3 degrés de translation. La contrainte doit être composée du même nombre d'équations que le nombre de degrés de liberté enlevés, soit 3 équations pour la rotule.

La liaison rotule est physiquement réalisée avec une coque sphérique dans laquelle une boule est prisonnière, comme le montre la figure ci-dessus. La description la plus simple de cette contrainte est que le centre de la boule doit coïncider avec le centre de la coque sphérique. D'une manière mathématique, on définit donc un point appartenant à un solide et un autre appartenant à l'autre solide et on pose comme contrainte que la coïncidence de ces deux points donne une liaison rotule.

La liaison entre les solidesets'écrit comme suit :

Cette équation revient à dire que le pointappartenant au solidedoit coïncider avec le pointappartenant au solide. Les vecteurs etsont constants ; ce sont les coordonnées des points d'attache par rapport aux centres de gravité des solides, tandis que varient au cours du temps : ce sont nos positions cartésiennes et nos matrices de rotations faisant partie des coordonnées généralisées .

Les variables de vitesses sont les vitesses en translations et les vitesses de rotation des deux solides, c'est-à-dire. Maintenant que les variables de vitesse sont définies nous pouvons calculer la jacobienne de cette contrainte. Pour simplifier le calcul, nous pouvons remarquer que cette jacobienne est donnée par car : il suffit de dériver la contrainte par rapport au temps ; ainsi, le gradient de la contrainte se calcule plus facilement, sans faire apparaître la matrice. En fait, lorsque nous dérivons la contrainte, nous multiplions par la matricesans nous en apercevoir :,

L'application crochet est la matrice de l'application associée au produit vectoriel, car elle vérifie :

La dérivée de la contrainte numérotéeest alors :

En utilisant les propriétés du produit vectoriel, les termes non nuls de la jacobienne sont :

La jacobienne est indexée par matrices de 3 par 3. Cela simplifie et accélère les calculs. est l'indice du bloc de 3 lignes correspondant à la contrainte numéro. est l'indice de la variable de vitesse linéaire du solide, est celui de sa variable de rotation.

On a aussi le terme.

Dans ce cas particulier, le multiplicateur de Lagrange est exactement égal à la force de contrainte entre les deux solides.

III - Intégration

Animation en format "mp4" non supportée.
Contenu uniquement consultable en ligne à cette adresse :
https://pierre-benet.developpez.com/moteur_physique/simulation1.mp4

Simulation d'un corps en apesanteur soumis à sa propre inertie : les segments sont en jaune, les forces sont représentées en vert.

L'intégration de nos équations différentielles consiste à résoudre de façon approchée l'équation différentielle. À partir des conditions initiales (vitesses et position des solides), on calcule de façon approchée à différents instantsles nouvelles positions et vitesses. On utilise le ou les pas de temps précédents pour calculer le pas de temps suivant. En général, ces instants t sont équidistants, car on intègre avec un pas de temps constant.

La résolution numérique d'équations différentielles est une discipline de l'analyse numérique. Il est important qu'une méthode numérique donne une solution proche de la réalité avec peu de calculs. Pour comparer différentes méthodes selon ce critère, on introduit l'ordre d'une méthode d'intégration : plus l'ordre est élevé, plus la méthode converge vite vers la valeur réelle de la solution.

Mathématiquement, on dit qu'une méthode d'intégration est d'ordresi, pour toute solution de l'équation, l'erreur d'intégration de la méthode en fonction du pas de temps est plus petite qu'une puissance dedu pas de temps :.

III-A - Réduction au premier ordre

Dans notre cas, l'équation est du second degré, avec un changement de variables. Or, la plupart des méthodes de calcul d'équations différentielles sont développées pour des équations du premier ordre. Cependant, nous pouvons réduire le degré de l'équation différentielle facilement pour trouver la méthode adéquate à appliquer :

etdeviennent. Nous avons donc bien une équation de la forme .

III-B - Méthode d'Euler

La méthode d'Euler est une manière très simple (voire basique) de résoudre des équations différentielles du premier ordre.

Grâce à la transformation précédente, nous pouvons alors appliquer cette méthode à notre problème de dynamique.

III-C - Dégénérescence des rotations

Cependant, l'incrémentation de la matrice de rotationgénère des erreurs et, rapidement, celle-ci n'est plus orthogonale. En particulier, l'erreur due à la rotation est linéaire par rapport au pas de temps, même si une méthode d'intégration avec un ordre plus élevé est utilisée : par exemple, avec Runge Kutta 4, une méthode d'ordre 4, l'erreur reste d'ordre 1.

Ainsi, il faut normaliser cette matrice de rotation à chaque incrémentation pour implémenter des méthodes d'intégration d'ordre plus élevé. Sinon, la méthode sera d'ordre linéaire, car la dégénérescence de la rotation sera prépondérante.

On fait donc intervenir une nouvelle variable secondaire, un quaternion, pour normaliser la matrice dans le schéma d'intégration. Le quaternion est intégré, comme toutes les autres variables, et normalisé :

Finalement, la matrice de rotation est calculée à partir du quaternion.

III-D - Méthode de Runge Kutta d'ordre 4

La méthode de Runge Kutta 4 (RK4) est une méthode d'ordre 4 pour résoudre les équations différentielles ordinaires. Il existe des méthodes de Runge Kutta pour tout ordre, mais la plus connue est celle d'ordre 4 : elle est simple, robuste et offre un bon compromis entre précision et calculs. Considérons le problème suivant :

,

La méthode de RK4 est donnée par l'équation :

avec les coefficients

,,,

III-E - Méthodes d'Adams-Bashforth

Il existe aussi des méthodes d'Adams-Bashforth pour tout ordre, mais nous présenterons ici la méthode d'ordre 4, qui offre un pendant à RK4.

Contrairement aux méthodes de Runge Kutta, où le nombre d'appels de la fonctionest égal au degré de la méthode, ces méthodes utilisent un seul appel par pas, quel que soit l'ordre de la méthode, après l'initialisation. En effet, cette méthode réutilise les évaluations précédentes de la fonction : il suffit donc de stocker les trois derniers appels, le seul à effectuer pour chaque itération est donc, qui sera stocké pour l'itération suivante. C'est une méthode multipas, c'est-à-dire qui utilise plusieurs pas précédents de la solution pour calculer le pas suivant.

Cette méthode apparaît donc plus efficace que la méthode de Runge Kutta. Cependant, la raison pour laquelle elle est peu connue est que le démarrage de la méthode nécessite de calculer les trois premiers pas à l'ordre 4. Ainsi, il faut, par exemple, effectuer trois itérations de la méthode RK4 pour démarrer la méthode, ce qui rend l'implémentation un peu lourde. De plus, il est difficile de changer le pas de temps d'une telle méthode en cours de route.

Un autre avantage de la méthode de Runge Kutta sur la méthode d'Adams-Bashforth est son rayon de convergence, plus de quatre fois plus élevé. En d'autres termes, le pas de temps minimal est plus de quatre fois plus élevé pour la méthode de Runge Kutta que pour la méthode d'Adams-Bashforth.

De plus, la méthode d'Adams-Bashforth n'est pas quatre fois moins chère que la méthode de Runge Kutta (au vu du nombre d'évaluations de la fonction réduit), mais seulement 1,67 fois : ses propriétés de convergence sont légèrement moins bonnes, le pas de temps doit être plus petit. En effet, pour obtenir la même précision entre les méthodes, il faut que le pas de temps de la méthode d'Adams-Bashforth soit 2,4 fois plus petit que le pas de temps de la méthode RK4. Sachant qu'un pas de la méthode de Runge Kutta est quatre fois plus cher (quatre calculs de la dérivée) qu'un pas de la méthode d'Adams-Bashforth, on trouve le facteur 1,67.

III-F - Précision et convergence

Les méthodes d'Euler, de Runge Kutta d'ordre 2 (RK2) et 4 (RK4) ont été implémentées afin de vérifier leur convergence. Elles ont respectivement une convergence d'ordre 1, d'ordre 2 et d'ordre 4. Voilà le résultat de quelques simulations.

L'erreur est calculée à partir de l'énergie initiale et de l'énergie finale du système comme. Cela nous procure une bonne idée de l'erreur d'intégration du système. Nous pouvons voir que la méthode RK4 fait diminuer l'énergie et que les méthodes d'Euler et RK2 la font augmenter. Cette propriété donne un caractère stable à la méthode RK4, un peu comme les méthodes implicites.

Erreur d'intégration lorsque le bonhomme fait un tour.

pas

20

40

80

160

320

640

RK4

RK2

Euler

Nous pouvons voir que, rapidement, la méthode d'Euler est inutilisable (nan).

Beaucoup de moteurs physiques utilisent la méthode « leapfrog » pour intégrer le mouvement, une variation de la méthode d'Euler. Cette méthode d'ordre 2 symplectique (qui conserve l'énergie) fonctionne pour les équations de la forme et est aussi simple à implémenter que la méthode d'Euler. En effet, il suffit de calculer la vitesse au pas de temps suivant à partir de la position précédente, puis de calculer la position au pas de temps suivant avec la vitesse qui vient d'être calculée (et non pas celle au pas de temps précédent).

Seulement, l'équation du mouvement est rarement de cette forme (elle est de cette forme pour la balistique, la mécanique faisant intervenir les forces à potentiel comme l'interaction gravitationnelle). Cependant, on peut quand même l'utiliser dans le cas général, mais la méthode retombe alors au premier ordre.

Si l'on cherche seulement une solution plausible, d'après le tableau ci-dessus, la méthode RK4 apparaît être la plus efficace. Il suffit de 20 pas pour résoudre notre problème, ce qui revient à 80 évaluations de la fonction dérivée. Avec RK2, il faut 80 pas, ce qui revient à 160 appels de la fonction dérivée. Euler nécessite 320 appels de la fonction dérivée. Seulement, dans les jeux vidéo, il y a des frottements dans le système qui assurent que le système ne diverge pas et le nombre d'images par seconde à afficher à l'écran oblige à avoir un pas de temps petit, là où la méthode d'Euler fonctionne. Ainsi, la plupart des moteurs physiques utilisent une méthode d'Euler, malgré ses problèmes de convergence : Comparaison de différents moteurs physiques.

Par contre, les logiciels de CAO et d'ingénierie utilisent plutôt la méthode RK4, qui offre une précision très élevée pour un coût faible, sans hypothèse restrictive sur les situations à simuler pour garantir une convergence décente.

IV - Correction des erreurs numériques

Au fur et à mesure de l'intégration de notre équation différentielle, des erreurs numériques s'ajoutent. Après quelques itérations, notre contrainte n'est plus satisfaite. Il faut cependant corriger ces contraintes pour pouvoir effectuer de longues simulations, sinon les erreurs deviennent trop importantes et le système diverge : la solution calculée n'a plus beaucoup de sens. Puisque nos contraintes possèdent certaines propriétés, il est possible de réparer ces erreurs numériques assez facilement. D'un autre côté, s'il faut à chaque calcul normaliser les quaternions, il n'est pas nécessaire de corriger les contraintes à chaque itération pour avoir une méthode à convergence d'ordre élevé.

Dans le reste de cette section, nous supposons donc qu'il y a une petite erreur de contrainte :

.

On dit qu'une suiteconverge linéairement verss'il existetel que.

On dit de même qu'une suite est convergente d'ordresi

.

En particulier, la méthode de Newton-Raphson est convergente d'ordre 2 (ses itérés successifs forment une suite qui converge avec un ordre 2), elle est dite quadratique.

IV-A - Contraintes arborescentes

Pour corriger ces erreurs numériques, il y a plusieurs méthodes. Dans un premier temps, on peut facilement corriger une erreur sur un système où les contraintes ne sont pas bouclées, c'est-à-dire qu'il n'y a pas de chaîne fermée de contraintes. Dans ce cas, les contraintes forment des liaisons comme dans un graphe et ce graphe est arborescent. On peut définir une racine et corriger les contraintes une par une.

Pour éviter d'avoir trop d'erreurs à corriger, dans ce cas particulier de contraintes bouclées, il est aussi possible d'utiliser des coordonnées réduites (ou coordonnées généralisées). Par exemple, pour une liaison pivot, on définit la position d'un objet par rapport à l'autre par seulement un angle. Au lieu d'avoir cinq inconnues en force, on a une inconnue en déplacement. Par cette définition, il ne peut pas y avoir d'erreurs sur les contraintes. En coordonnées réduites, il convient d'utiliser les équations de Lagrange pour trouver les équations qui décrivent le système. C'est un peu une formulation duale du problème : à la place de résoudre les forces (les degrés de liberté enlevés), on résout directement le déplacement (les degrés de liberté restants). Cependant, considérer les équations de Lagrange donne dans ce cas des matrices pleines et non creuses comme dans notre cas. Les calculs deviennent alors plus longs et ralentissent les performances du moteur physique.

Dès qu'il y a des contraintes bouclées, cette réécriture du problème ne suffit pas.

IV-B - Force fictive

Nous pouvons aussi utiliser des forces fictives qui agiront comme des ressorts entre les liaisons. Pour cela, il faut bien choisir la rigidité de notre ressort et son amortissement. Ces derniers devront être réglés en fonction du pas de temps utilisé. Cependant, cette méthode n'est pas très précise et, si elle a l'air simple, le choix des ressorts et amortissements est difficile à réaliser. Cette méthode est connue sous le nom de penalty force method.

La convergence d'une telle méthode est linéaire, c'est-à-dire que l'erreur diminue à chaque itération d'un facteur constant (0,5, dans le cas où le système est bien réglé).

IV-C - Newton-Raphson

Une autre méthode est la méthode de Newton-Raphson, qui consiste à suivre le schéma jusqu'à ce que, avec le système à corriger et l'itération d'arrêtle système corrigé.

La convergence d'un tel schéma est très rapide, elle est dite quadratique. À chaque itération, l'erreur suivante sera équivalente au carré de l'erreur précédente.

Nous allons utiliser cette dernière méthode, car, si elle peut paraître difficile à résoudre (inversion d'un système linéaire), nous allons montrer que le système linéaire est le même que celui introduit précédemment.

Dans notre cas, pour retrouver notre matrice, nous allons poser : .

Nous voulons résoudre le système, àéquations etinconnues, avec. Nous ne pouvons pas directement résoudre ce système, puisqu'il y a plusieurs solutions possibles. Nous allons donc essayer de réduire le nombre d'inconnues.

Tout d'abord, nous pouvons remarquer que les forces (multiplicateur de Lagrange) que nous appliquons à notre système ont pour rôle de conserver cette contrainte : on peut alors se dire que ces forces feront de très bonnes inconnues pour notre système. Nous allons donc essayer de faire intervenir de telles forces à la place de nos inconnues en déplacement.

Pour cela, nous allons considérer le problème sous un autre angle : il s'agit de minimiser, qui représente une énergie de déformation, sous la contrainte, qui est la correction que nous voulons obtenir.

Nous allons donc trouver un extremum à la fonction en annulant son gradient. Ce dernier est donné par, l'extremum correspond donc à l'équation .

De la même manière, en substituant cette équation dans notre contrainte, on obtient une solution qui ressemble beaucoup à ce que nous avons déjà vu :

IV-D - Impulsion

Pour la correction en vitesse, la correction est directe, car l'équation est linéaire. Nous pouvons y trouver un parallèle avec la théorie de l'impulsion, qui autorise des discontinuités de vitesses. Une impulsion est une force infiniment grande appliquée pendant un temps infiniment court, de sorte qu'il en résulte une différence de vitesse finie. Nous pourrons d'ailleurs grâce à cette correction simuler des impulsions.

La correction en vitesse permet effectivement de simuler des chocs. Par exemple, supposons que nous voulons qu'un solide de notre système subisse un choc absorbant sur un élément statique. Cela signifie que le solide perdra toute sa vitesse. Il suffit d'ajouter une contrainte (non satisfaite) qui impose que ce solide ait une liaison rotule avec le sol. Le solide et le sol ont donc des vitesses différentes. En appliquant la correction en vitesse sur ce nouveau système, des forces vont s'appliquer entre les joints et entre le sol et le solide qui subit le choc, de sorte que le solide perde toute sa vitesse par rapport au sol et que les autres contraintes restent satisfaites. Le choc a donc bien été calculé.

Si nous voulons un rebond conservatif, il faut alors appliquer deux fois les forces de cette correction. Il devient alors intéressant de vérifier que cette méthode fonctionne bien en regardant la variation d'énergie totale du système.

Si nous voulons un rebond partiel, il suffit d'appliquer un coefficient entre 1 et 2.

IV-E - Cinématique inverse

La cinématique inverse désigne l'ensemble des méthodes de calcul des positions et rotations d'un modèle articulaire afin qu'il satisfasse des contraintes en plus. Il faut par exemple que notre bonhomme place sa main à une certaine position et nous voulons une valeur des variables de position et de rotation du bonhomme qui vérifie cette propriété.

Pour cela, nous pouvons poser une nouvelle contrainte à notre bonhomme et utiliser la méthode de Newton Raphson décrite ci-dessus pour trouver une solution. Il faut alors itérer la méthode pour trouver le résultat — rappelons que ce n'est pas une méthode directe, mais à convergence quadratique. Cela fonctionne à condition que l'initialisation des variables de position soit assez proche du résultat, sinon, la méthode ne convergera pas. Il faut, dans ce cas-là, utiliser des méthodes plus complexes pour trouver le résultat. Habituellement, la méthode de Denavit-Hartenberg est utilisée pour les systèmes articulaires (avec des rotules ou pivots) et arborescents. La méthode que nous allons présenter n'est pas une approche habituelle, mais est générale et assez intuitive et n'introduit pas de nouveaux calculs.

Dans notre cas, nous allons simplement borner lecalculé par Newton-Raphson à chaque pas de temps. La borne se fera sur les rotations : si le maximum est supérieur à 2, par exemple, on divise letrouvé par. Ceci est réalisé comme suit :

Alors, si la solution est éloignée (c'est-à-dire), l'itération utilisée sera du type descente de gradient, c'est-à-dire, on fixe la longueur du pas (ici égale à 2) et l'on suit simplement le gradient de la fonction afin de minimiser l'erreur sur la contrainte. Donner un pas élevé permettra que les variables sautent les creux de la fonction à satisfaire, de telle sorte que les variables ne restent pas dans une configuration où la méthode de Raphson ne converge pas. Lorsque l'on sera assez proche de la solution, alors l'itération utilisée sera du type Newton-Raphson, avec une convergence quadratique.

Dans notre cas, nous avons utilisé cette méthode pour que le bonhomme attrape son pied. Initialement, il se tenait debout. Il y a des longueurs de segments différents et la jambe doit se retourner pour que le bonhomme attrape son pied : la méthode de Newton ne va donc pas fonctionner directement, les deux poses sont trop différentes l'une de l'autre. L'erreur calculée ici est l'erreur quadratique en contrainte. Le système s'approche peu à peu de la solution, puis, lorsqu'il est assez proche, on observe une convergence quadratique. Finalement, le résultat a été trouvé après 26 itérations.
On peut alors exécuter une simulation dynamique…

Résultat de l'exécution d'une simulation

V - Résolution du système linéaire

Pour calculer l'accélération du système nous devons résoudre l'équation linéaire en. Il s'agit d'inverser une matrice creuse définie positive.

V-A - Propriété du problème

Solides liés de manière arborescente.

Dans cet exemple, plusieurs solides sont liés par des contraintes de façon arborescente (les contraintes ne forment pas de boucle). On peut faire correspondre chaque force à une contrainte ; ici, les forces peuvent être considérées comme les multiplicateurs de Lagrange. Pour mieux se représenter le système, admettons que les liaisons soient des liaisons rotules : ainsi, les contraintes entre deux solides sont des forces, des vecteurs de trois variables, et se représentent par des contraintes vectorielles qui forment trois équations scalaires.

Les équations que nous avons proviennent de contraintes de position entre deux solides. Les forces appliquées aux solides ont été remplacées à la place des inconnues de déplacement à l'intérieur des équations, de sorte que l'équation finale pour une contrainte fait intervenir toutes les forces qui s'appliquent sur les deux solides contraints, comme le montrent les flèches vertes. Cela ne fait pas beaucoup de termes par équation, c'est pour cela que nous pouvons parler de matrice creuse.

V-B - Algorithme des corps articulés

L'algorithme des corps articulé, plus connu sous le nom de « articulated body algorithm » ou la méthode de Featherstone est une méthode efficace pour résoudre ce problème lorsque les contraintes sont arborescentes. Sa complexité est linéaire en fonction du nombre d'inconnues, alors qu'une inversion classique à l'aide d'un pivot de Gauss a une complexité cubique en fonction du nombre d'inconnues.

Si vous n'êtes pas familier du pivot de Gauss, je vous invite à lire cet article sur les systèmes linéaires.

L'algorithme de Featherstone est une méthode directe pour inverser un système linéaire de forme arborescente. C'est en fait un pivot de Gauss par bloc (de 6 par 6) dans lequel on sait quels termes il faut éliminer, il se compose d'une descente et d'une remontée. Les pivots sont choisis sur la diagonale. La descente du pivot de Gauss sur une telle matrice lorsque les inconnues sont rangées dans le bon ordre (les pères après les fils) ne remplit pas la matrice. On sait exactement quels sont les termes modifiés. Nous n'avons donc pas besoin de considérer une matrice pleine. La remontée se fait aussi sans problème puisque l'on sait quels termes sont non nuls. La partie sur le pivot de Gauss creux permet de comprendre comment un tel algorithme peut fonctionner.

L'algorithme est accessible ici.

V-C - Élargissement de l'algorithme

Pour traiter le cas où il y a des boucles dans les contraintes (le système n'est plus arborescent), on peut utiliser d'autres algorithmes tels que la méthode de Gauss Seidel, très utilisée pour ce genre d'équations ou une combinaison de celle-ci et de l'algorithme des corps articulés. On peut aussi penser à la méthode des gradients conjugués, puisque notre matrice à inverser est symétrique définie positive.

Par exemple si on note le système

, avecle système arborescent, on peut inverser le système par blocs. Il suffit alors d'inverseret le complément de Schur de notre matrice :(qui est toujours inversible, carest définie positive).

On peut alors inverser le complément de Schur avec la méthode du gradient conjugué ou un pivot de Gauss si sa dimension est petite. On a alors une complexité maximum deavecla dimension deetle nombre de termes non nuls danset .

On peut aussi utiliser une variante de la méthode de Gauss Seidel, une méthode itérative qui profite de la propriété creuse des matrices, ce qui permet d'inverser des gros systèmes là où les méthodes directes prendraient trop de temps. Il y a aussi une version spécialement conçue pour les problèmes avec des contraintes d'inégalité (Projected Gauss Seidel et dérivées), très utilisée :

avec initialementetla partie triangulaire inférieure, la partie triangulaire supérieure etla diagonale de notre matrice initiale.

Il suffit de poser, et de suivre l'algorithme précédent pour résoudre notre système avec une vitesse de convergence plus élevée.

V-D - Pivot de Gauss creux

Habituellement, le pivot de Gauss n'est pas utilisé pour résoudre des systèmes linéaires creux, car, lors de la descente du pivot de Gauss, la matrice se remplit et ne devient plus du tout creuse, mais pas dans notre cas.

Lorsque l'on réalise la descente du pivot de Gauss, on élimine les termes sous la diagonale par des additions de lignes de la matrice. On choisit un pivot sur une ligne, puis on l'utilise pour annuler les termes qui sont sur la colonne du pivot. La complexité du pivot de Gauss vient principalement de la descente. Ainsi, pour que notre pivot de Gauss creux soit efficace, il faut que la matrice reste creuse lors de la descente du pivot de Gauss.

Éliminer une inconnue dans le système d'équations revient à enlever un terme aux équations et à le remplacer par d'autres. Pour que le système d'équations ne grandisse pas lorsque l'on élimine les inconnues, il faut que l'équation utilisée pour éliminer l'inconnue soit petite et que cette équation soit substituée dans un nombre minimum d'équations. Dans les deux exemples suivants, nous voyons que les matrices après descente du pivot de Gauss ne se sont pas remplies. Dans le cas d'une chaîne et d'une boucle, le pivot de Gauss creux est efficace.

Chaîne de solides liés.

Si nous n'avons qu'une chaîne de contraintes, la matrice à inverser sera tridiagonale par bloc de 3 par 3, car les équations ne feront intervenir que trois termes à la fois.

Descente du pivot de Gauss

Boucle de solides liés.

Lorsqu'il y a une boucle, la matrice n'est pas tridiagonale, mais le pivot de Gauss reste efficace. Ici, les équations ont seulement un terme en plus par rapport à la descente du pivot de Gauss tridiagonal.

Descente du pivot de Gauss.

De manière plus générale, la substitution d'une inconnue ne modifiera que les équations des contraintes voisines, car les équations ne font intervenir que des contraintes voisines. Ainsi, si nous substituons une zone d'inconnues, les équations modifiées seront toutes celles en bordure de cette zone et le nombre de termes dans ces équations sera équivalent au nombre d'équations de la bordure.

Cela explique pourquoi, dans notre cas, cet algorithme est efficace. Si nous avons substitué les équations de toute une branche de contrainte, il n'y aura qu'une seule équation en bordure, là où la branche a été coupée. Au contraire, dans un maillage d'éléments finis, prenons par exemple un cube, supposons que la moitié des inconnues du cube aient déjà été substituées. Alors, la surface en bordure de ces équations est beaucoup plus importante. Ainsi, nous aurons des complexités différentes :

- pour un corps articulé arborescent ;

- pour un maillage d'éléments finis 3D.

De plus, cela nous renseigne quant au choix de l'ordre des pivots.

Solides liés de manière arborescente.

Deux branches du système peuvent être substituées sans interférer entre elles si elles ne sont pas voisines. Ainsi les équations ne grossissent pas.
Par exemple, dans le système ci-contre, si nous substituons toutes les inconnues en dehors du cercle vert, les équations pour f2,f 3 et f4 sont devenues plus petites (deux termes au lieu de trois).
Au contraire, si nous substituons toutes les inconnues à l'intérieur du cercle vert, les équations pour f1, f5 et f9 sont devenues plus grosses, puisque ces équations font intervenir toutes celles en bordures, à savoir f1, f5 et f9.

Il faut donc substituer les inconnues en partant des extrémités des branches.

V-E - Implémentation du pivot de Gauss creux

Comme évoqué précédemment, les équations seront rangées sous forme de listes chaînées ordonnées de matrices trois par trois. À chaque terme (matrice trois par trois) sera associé le numéro de sa force ; les termes nuls ne font pas partie de la liste. Les équations sont rangées dans un tableau. Notre matrice creuse est donc un tableau de listes de termes.

Exemple du système à inverser pour cinq corps articulés bouclés.

Les termes sont chaînés par ordre croissant de la colonne du terme. L'utilisation de listes chaînées ordonnées nous permet d'additionner deux lignes facilement. Pour additionner deux lignes (deux blocs de trois lignes), il suffit de parcourir en même temps les chaînes des deux équations. On insère alors les éléments par ordre croissant dans une nouvelle liste.

Si le terme de l'équation 1 se place avant le terme de l'équation 2, on place le terme de l'équation 1 à la suite de la liste et on réitère avec le prochain terme de l'équation 1. On fait de même si c'est le terme de l'équation 2 qui se place avant. Si les termes des deux équations sont de la même colonne, on les additionne avant de les insérer. Ainsi, il suffit de parcourir une fois les deux équations dans l'ordre pour les additionner.

Il reste un problème embêtant : savoir quelle équation contient un terme qui doit être éliminé par le pivot. On pourrait parcourir l'ensemble des équations, mais la complexité serait mauvaise. La solution est de ne pas substituer directement l'inconnue du pivot dans toutes les équations, mais de le faire uniquement lorsque l'équation arrive : il suffit alors de parcourir l'équation ; les termes ayant le numéro de leur colonne, on trouve les équations correspondant aux pivots à éliminer. Puis on élimine les termes dans l'ordre.

Descente du pivot de Gauss conventionnel : on parcourt les équations de haut en bas et on élimine les termes sous le pivot.

Descente du pivot de Gauss creux : on parcourt les équations de haut en bas et on élimine les termes à gauche du pivot de l'équation en utilisant les pivots précédents.

De même pour la remontée du pivot de Gauss, on élimine les termes à droite du pivot en utilisant les pivots en dessous de l'équation.

VI - Dynamique inverse

Pour commander les actionneurs d'un robot, il peut être intéressant de savoir quelle force appliquer aux articulations pour réaliser directement le mouvement et avoir un asservissement très précis. La cinématique inverse permet de trouver un paramétrage des variables d'un robot pour qu'il effectue un mouvement, la dynamique inverse calcule les forces à appliquer à partir de ce paramétrage pour que le robot réalise réellement le mouvement.

Bras articulé.

Le problème que nous voulons résoudre est le suivant : connaissant un paramétrage des coordonnées dans le temps du système, nous voulons trouver quelles forces doivent appliquer les actionneurs pour réaliser ce déplacement.

En d'autres termes, connaissant les positions, vitesse et accélération,, entre, nous allons calculer les forces au niveau des joints qui engendrent ce déplacement.

Nous cherchons ici des forces qui travaillent, pas seulement des forces de contrainte. Nous ne pouvons donc pas considérer les forces inconnues sous forme d'un multiplicateur de Lagrange.

En pratique, ce genre de problème arrive sur des contraintes arborescentes. Ainsi, nous allons considérer que tous les joints ont trois inconnues de force et trois inconnues de couple. Dans ce cas-là, il y a alors autant de variables de déplacement que de contraintes inconnues. Autrement dit, le système est inversible.

Voici les équations obtenues.

La forceest la force à la racine de l'élément, elle est comptée positivement dans notre expression, alors que les actions des forces des fils sont soustraites. Cela signifie que nous avons choisi que les forces s'appliquent du père sur le fils.

En résolvant ce système, nous obtenons les forces de contraintes et les forces d'action en même temps, on récupère alors les forces d'action en projetant les forces obtenues. Nous pouvons alors commander nos actionneurs.

cinématique inverse

void Solveur::inversedyn(){
int k,i;
Matr9d *equa1; // équation courante
Matr9d *equa2;
Vect3d *secondterme; // second terme équation courante

// initialisations
constante.setend(listejoints.end()*2);
for(k=0;k<listejoints.end();k++){constante(k)=Vect3d();}
sys.setend(listejoints.end()*listejoints.end()*4);
for(k=0;k<sys.end();k++)sys(k)=Matr9d();

for(k=0;k<listejoints.end();k++){
Solide &partk = *(*listejoints(k)).part2;
equa1 = &sys(2*k*listejoints.end()*2);
equa2 = &sys((2*k+1)*listejoints.end()*2);
secondterme = &constante(k*2);

secondterme[0] += partk.masse*partk.acceleration;
secondterme[1] += partk.getinertie()*partk.a_rot + (partk.w^partk.getinertie()*partk.w);
for(i=0;i<partk.liens.end();i++){
Joint &link = *partk.liens(i);
if(link.is_pere(partk)){
equa1[link.iforce*2] += Matr9d(1);
equa2[link.iforce*2+1] += Matr9d(1);
equa2[link.iforce*2] += (-link.AG1()).v();
}else{
equa1[link.iforce*2] += Matr9d(-1);
equa2[link.iforce*2+1] += Matr9d(-1);
equa2[link.iforce*2] += link.AG2().v();
}
}
}
gauss(sys, constante);
}

VII - Structures de données

Le moteur physique a été développé en C++. À ce stade du développement, nous pouvons déjà diviser le moteur en plusieurs parties différentes :

- un module avec différentes classes permettant de faire des calculs avec des vecteurs en trois dimensions, avec des matrices 3 par 3, utiles pour les rotations, et les quaternions ;

- une classe solide, qui permet de stocker toutes les caractéristiques des solides dans une classe, comme l'orientation, la position, la liste des liaisons qui font intervenir le solide, sa masse, son inertie ;

- une classe liaison ou « joint », qui enregistre toutes les caractéristiques relatives aux contraintes ;

- une classe qui permet d'orchestrer le système de solides et de contraintes que nous appellerons « corps ».

Pour un moteur physique, il vaut mieux utiliser des vecteurs en double précision, voire plus : lorsque les calculs impliquent des inversions de matrices et beaucoup d'itérations, les erreurs d'arrondi deviennent importantes.

Une boucle dans la fonction principale permet d'appeler les fonctions d'intégration, de gestion des impulsions, de correction et d'affichage.

À ce niveau, on peut penser au parallélisme pour accélérer les calculs. La première transformation à effectuer pour résoudre les équations différentielles est la parallélisation spatiale : il peut être intéressant de séparer les ensembles de solides liés pour les traiter sur des cœurs différents. On peut aussi penser à effectuer le pivot de Gauss creux en partant de différentes branches qui peuvent être effectuées sur des processeurs différents tant que les équations n'interfèrent pas. D'autres approches plus récentes visent à paralléliser l'intégration de l'équation différentielle. RK4 appelle quatre fois de suite la fonction dérivée pour un pas de temps, mais il y a des variantes parallèles de RK4 qui permettent d'appeler quatre fois la fonction dérivée en parallèle, puis de trouver le pas suivant grâce à une combinaison de ces appels. Il y a aussi d'autres algorithmes comme Parareal plus généraux pour les méthodes d'intégrations parallèles.

VII-A - Articulation

La classe qui représentera l'articulation sera appelée Joint(qui vient de l'anglais). Celle-ci permet de définir chaque contrainte. L'approche que nous avons choisie est un joint défini par sa fonction contrainte. Nous allons donc utiliser une fonction virtuelle (ou pointeur de fonction) qui donnera cette fonction contrainte. Pour accélérer les calculs, cette fonction virtuelle calculera les gradients et les autres termes utilisés. Ces variables seront stockées dans la classe joint pour plus de simplicité.

Une autre approche moins flexible aurait été de considérer les liaisons comme une liaison encastrée à laquelle on supprime des degrés de liberté au choix. Il n'y a alors plus besoin de pointeur de fonction, il y a seulement les équations d'une liaison encastrée auxquelles on supprime des lignes et colonnes. Certains moteurs physiques pour jeux comme bullet proposent de telles liaisons. Cependant, elles ne permettent pas de décrire autant de liaisons qu'une fonction contrainte quelconque comme la nôtre. Cela reste un point intéressant au niveau de l'interface utilisateur, choisir la contrainte en choisissant les degrés de liberté de la liaison.

classe Joint

class Joint{
public:
Joint(Solide *p1,Solide *p2,Vect3d const&bar1A1,Vect3d const&A2bar2);
inline bool is_pere(Solide &s){return &s==part1;}

// précalcule les gradients, les erreurs de contrainte, et le second terme
virtual void updatevar();
// remplissage de l'équation associée à la contrainte
void fillequa(Matr9d *equa);

Vect3d A()const; // donne la position d'attache du joint, utile pour l'affichage
Vect3d AG2()const; // donne la position d'attache du joint par rapport au solide 2
Vect3d AG1()const; // donne la position d'attache du joint par rapport au solide 1
void corrigevar(); // corrige la contrainte pour un système arborescent

//variables de stockage
Vect3d errc; // valeur de la fonction contrainte (doit valoir 0)
Vect3d errv; // dérivée de la fonction contrainte (doit valoir 0)
Vect3d terme2; // second terme: J'u' qui apparait dans l'équation
Matr9d Gq1,Gq2,Gp1,Gp2; // les gradients

//caractéristiques de la liaison (constantes)
Matr9d base0; // si la contrainte s'exprime dans une certaine base
Quaternion quat0;

int iforce; // numéro de la force associé au joint
int libertes; // degrés de liberté (entre 0 et 3)

Solide* part1; // pointeur vers le solide 1
Solide* part2; // pointeur vers le solide 2
Vect3d G1A1; // position de la liaison dans la base locale au solide 1
Vect3d G2A2; // position de la liaison dans la base locale au solide 2
};

VII-B - Solide

La classeSoliden'aura pas de méthodes importantes associées. Cette classe est surtout une structure pour stocker l'état du solide au cours du temps (position, orientation, vitesses) et stocker ses caractéristiques physiques, comme son moment d'inertie, sa masse et ses liaisons avec les autres solides.

Pour simplifier le code, des variables secondaires ont été ajoutées pour utiliser la méthode RK4 directement. Bien évidemment, il vaut mieux éviter de procéder de la sorte pour du code de production, afin de séparer clairement l'affichage des calculs.

solide

class Solide{
public:
Solide(double m,Vect3d const& dimensions);
inline Matr9d getinertie()const {return R*inertie0*transpo(R);}

inline Vect3d A()const{if(liens.end()>0){return (*liens(0)).A();} else return p;}

/** caractéristiques constantes **/
double masse;
Matr9d inertie0;
dynamictab<Joint*> liens;
int numero; // numéro du solide

/** caractéristiques variables **/
// termes pour stocker l'itération actuelle ou pour calculer les accélérations
Quaternion q;// passage de la base terrestre à la base locale
Matr9d R;// passage de la base terrestre à la base locale
Vect3d w;// dans la base terrestre, vitesse de rotation
Vect3d p;// dans la base terrestre, position du barycentre
Vect3d v;// dans la base terrestre, vitesse du barycentre

// termes pour stocker les accélérations
Vect3d acceleration;// dans la base terrestre, accélération du barycentre
Vect3d a_rot;// dans la base terrestre, accélération en rotation

// termes pour stocker l'itération précédente
Quaternion q0;
Vect3d w0;
Vect3d p0;
Vect3d v0;

// termes pour stocker la somme (k1 + 2*k2 + 2*k3 + k4)
Quaternion q1;
Vect3d w1;
Vect3d p1;
Vect3d v1;
};

VII-C - solveur

La classeSolveurest en fait un système de solides sous contrainte. On peut ajouter n'importe quelle contrainte ou solide au système.

Elle comprend un tableau de solides et de joints, nous y mettrons aussi les fonctions d'affichage graphique. Lorsque nous avons des contraintes arborescentes, nous avons tendance à penser que des fonctions récursives pourraient être plus efficaces, mais ce n'est pas le cas. En rangeant les éléments de façon ordonnée, par exemple les fils derrière les pères, certains algorithmes récursifs deviennent alors plus faciles à écrire et plus efficaces, car ils deviennent itératifs : par exemple, il suffit de parcourir le tableau des solides et de corriger la position du solide fils, pour chaque solide, pour corriger de façon directe tout le système comme évoqué dans IV.A.

Cette classe comprendra toutes les fonctions de résolution, comme les méthodes d'intégration, les corrections et l'impulsion.

Nous avons ici choisi d'utiliser un pivot de Gauss pour inverser les systèmes d'équations.

classe Solveur

// pivot de Gauss
void gauss(dynamictab<Matr9d> &s, dynamictab<Vect3d> &v);

class Solveur
{
public:
inline Solveur(){}
Solveur(double poids); // constructeur
~Solveur(); // destructeur

void accelere_rot(Vect3d const& dw); // permet de mettre en rotation tout le bonhomme

int add(Solide *s);// ajoute un solide
int link(unsigned int p1,unsigned int p2, Vect3d const&d1, Vect3d const&d2);// ajoute une liaison rotule

void solvesystemeforce(); // calcule les multiplicateurs de Lagrange et les accélérations (forces)
void incrementeforce(double dt); // Euler
void incrementeforce2(double dt); // RK2
void incrementeforce4(double dt); // RK4

void corrigevars();// corrige pour des liaisons rotule arborescente
double corrigepos();// correction en position (une itération de la méthode de Newton), retourne l'erreur
void corrigev();// corrige en vitesse
void impulsion(); // génère une impulsion si le solveur sort de la zone
void corrige();// corrige la position et la vitesse en effectuant plusieurs itérations (méthode de Newton)

double geterreur(); // retourne la norme cartésienne de la contrainte en position
void inversedyn(); // dynamique inverse
double getenergie(); // retourne l'énergie totale du système
void recentre(); // positionne le bonhomme au centre et supprime sa vitesse de groupe
void affiche(); // affiche les segments

// quantité totale du système
double massetot;
Vect3d barycentretot;
Vect3d vtot;

// liste des joints et solides
dynamictab<Solide*> listepartie;
dynamictab<Joint*> listejoints;

// système à inverser préaloué
dynamictab<Vect3d> constante;
dynamictab<Matr9d> sys;
};

VIII - Sources

Les sources sont accessibles ici. Le programme dépend de SDL 1.2 et d'OpenGL.

IX - Remerciements

Merci à dourouc05 et à Littlewhite pour les relectures de cet article plutôt long.

Merci à ClaudeLeloup pour ses relectures orthographiques.