Projet image 2013 : Animation temps-réel de liquides (Cascade)

De Ensiwiki
Aller à : navigation, rechercher
Project schedule.png
Titre du projet Projet Image 2013 : Animation temps-réel de liquides (Cascade)
Cadre Projets de spécialité
Page principale Projets de spécialité image

Encadrants Marie-Paule Cani, Pierre-Luc Manteaux


L'équipe

ExempleSPH.png

Anne-Hermine Allain (MMIS)

Jean-Philippe Bauchet (MMIS)

Loïc Ciccone (MMIS)

Héloïse Guillemaud (MMIS)

Contexte des travaux et problématique

(La version originale du sujet est disponible ici.)

La problématique de ce projet de spécialité s'inscrit dans le cadre de la production de scènes virtuelles animées. Des jeux vidéo à la chimie en passant par la médecine et la thermodynamique, ingénieurs et chercheurs recourent de plus en plus régulièrement à la simulation de systèmes complexes comportant plusieurs centaines, voire plusieurs milliers de particules. On impose ici la contrainte temps-réel, c'est-à-dire l'absence de calculs en amont face à l'introduction d'une perturbation ; les calculs se font en même temps que l'affichage. Dans notre cas, il peut s'agir d'un robinet que l'on ouvre, d'une barrière qui se lève,... entraînant la propagation d'un liquide dans un milieu.

De par sa thématique, ce projet prend donc la forme d'un travail de recherche. Quelles sont les limites de la contrainte temps-réel ? Est-il possible de simuler en temps-réel de l'eau tombant dans un bassin ?

Objectifs du projet

Nos objectifs pour ce projet couvraient trois aspects : modélisation physique, rendu graphique, et optimisations.

Assimilant le fluide à un système de particules, dotées de caractéristiques propres telles que masse volumique ou densité, il nous fallait en premier lieu implémenter un modèle physique décrivant le fluide et ses caractéristiques globales, mais aussi les forces extérieures auxquelles sont soumises les particules qui le composent, pour en déduire les équations du mouvement qui leurs sont associées. Nous disposions pour cela de deux modèles :

  • Le premier, "Smoothed Particle Hydrodynamics" (SPH), consiste à considérer les particules comme une discrétisation du champ de pression et de masse volumique dans le fluide. À partir de ces points, il est possible de définir toutes les grandeurs physiques du système à n'importe quel endroit par interpolation. On peut ainsi en déduire les expressions des forces de pression, de viscosité, de tension de surface auxquelles est soumise une particule donnée, par somme des contributions des autres particules, celle-ci étant d'autant plus importante que l'autre particule est proche de la particule considérée, et nulle si elle est trop éloignée. Les équations du mouvement s'obtiennent alors par intégration par rapport au temps, selon la deuxième loi de Newton.
  • Le second, "Adaptively Restrained Particle Simulation" (ARPS), utilise des principes relatifs à la mécanique hamiltonienne et reformule les équations décrivant le mouvement des particules dans le fluide. La différence majeure par rapport au SPH est que les particules possèdent des états (elles peuvent être actives ou inactives). Si le moment d'une particule devient trop faible, elle peut être mise en sommeil et ne plus intervenir dans le calcul des autres équations du mouvement, améliorant ainsi la vitesse d'exécution du programme.

Nous reviendrons plus en détail sur ces modèles ultérieurement.

Un autre aspect du projet concernait le rendu graphique de la simulation. En effet, notre fluide peut être représenté de différentes manières :

  • Une première méthode basique consiste à dessiner des billes centrées en chacune des particules du système ;
  • On peut aussi utiliser le concept de surfaces implicites, vu en cours de graphique 3D. On dispose d'une fonction d'évaluation à support fini et centré en l'origine, considérant la distance d'un point de l'espace par rapport à une particule, et d'une isovaleur. Il faut alors représenter la surface séparant les points de l'espace où cette fonction d'évaluation est inférieure à l'isovaleur, de ceux où elle est supérieure : c'est la définition d'une surface implicite, que permettent de représenter certains algorithmes comme l'algorithme des "Marching Cubes" [ LC87 ]. Ici, il serait question de donner un aspect plus réaliste à notre fluide, en dessinant des gouttes d'eau, le film à la surface du fluide, et les rides qui le parcourent ;
  • D'autres techniques peuvent être utilisées comme le point splatting : une image ou une texture est projetée sur les objets de telle sorte que la forme de l'image varie doucement de son centre vers ses bords.
Illustration des méthodes de rendu issues de [MCG03]. À gauche: rendu à partir de billes. Au milieu: rendu à l'aide des "splatting images". À droite: rendu à l'aide des "Marching Cubes".


Enfin, un sujet de ce type mobilisant un certain nombre de ressources de l'ordinateur, et dans l'approche temps-réel qui est voulue ici, nous avons également cherché à optimiser le plus possible notre programme. Ce travail d'optimisation se décline en les quelques points suivants :

  • Utilisation d'une table de hachage pour la recherche des plus proches voisins, afin d'éviter des parcours superflus de la liste des particules qui conduiraient à l'ajout de termes nuls dans l'expression des forces ;
  • Choix des algorithmes, des structures de données et des structures de classe les plus performantes ;
  • Implémentation de l'ARPS ;
  • etc...


En guise de démonstration de notre programme, notre objectif final est l'intégration de tous ces éléments dans une scène complexe. C'est ainsi que notre groupe a choisi d'animer une chute d'eau dans un bassin, visible à la fin de cette page.

Planning et Organisation

La première étape a été de prendre connaissance de l'état actuel des recherches dans ce domaine. Tous les membres de l'équipe ont travaillé en même temps sur ce travail de documentation, de manière à tous avoir le même niveau de connaissances concernant ce projet. Ce travail a été assez long, et s'est finalement étalé sur une grande moitié du projet, car il y avait de nombreux articles, et nous avons dû faire la synthèse de tous ceux-ci, en prenant les éléments les plus intéressants de chaque. Ensuite, après avoir précisément défini nos objectifs, nous nous sommes réparti le travail en fonction des compétences et des goûts de chacun (une ou deux personnes sur chaque tâche en général). Cette parallélisation des tâches nous a, de plus, permis d'optimiser le temps de travail dont nous disposions. Cela dit, si l'un d'entre nous rencontrait des difficultés sur sa partie, les autres membres de l'équipe restaient toujours disponibles pour l'aider.

En ce qui concerne la communication, nous nous retrouvions tous les jours, du lundi au samedi, dans les locaux de l'Ensimag. Cela nous permettait de rester informés de l'avancement de nos coéquipiers et de pouvoir être présent en cas de blocage. Grâce à cette communication constante, nous n'avons rencontré aucun problème d'organisation, ou tension particulière au sein du groupe. Avec nos encadrants, nous communiquions régulièrement par mail et organisions une réunion tous les 3 à 4 jours environ. Nous étions donc bien suivis et redirigés vers la bonne voie quand cela s'avérait nécessaire. De plus, ils avaient constamment accès à notre code (sur dépôt Git) en cas de besoin.

Présentation des modèles existants

Modèles physiques

SPH

Noyau de Monaghan. La courbe représentée est celle ode la fonction W(r) dans le cas unidimensionnel. Il est clair que plus un point de l'espace est proche de la particule, plus l'importance accordée à la valeur de la grandeur A en ce point est grande. Si la particule est trop éloignée, l'influence est en revanche nulle.
Comme énoncé précédemment, le principe du modèle SPH est que les grandeurs scalaires initialement définies en des points discrets de l'espace, peuvent ensuite être évaluées n'importe où dans cet espace, par interpolation. Pour cela, nous considérons que chaque point de l'espace où les grandeurs sont définies influe sur tous les points compris dans son voisinage. Soit A une grandeur. Dans le cadre du SPH, nous avons :


A(r) = \sum_{j = 1}^{n} A_j \frac{m_j}{\rho_j} W(r - r_j)

n est le nombre de particules, m_j la masse de la j-ème particule, \rho_j sa masse, et A_j sa valeur pour la grandeur scalaire recherchée. W représente pour sa part un coefficient de pondération pour chaque terme de la somme, qui est d'autant plus élevé que la norme de r - r_j est faible. Il s'agit d'une fonction de r, appelée noyau de lissage. Compte tenu de la contrainte d'égale influence de part et d'autre de la particule, W doit être paire, de support fini et centrée en zéro.

Pour le choix de la fonction W, plusieurs expressions sont possibles : gaussienne, noyau de Monaghan [Mon05]...

Le calcul des forces de pression, de viscosité et de tension de surface s'appuie sur le même schéma : ce sont des sommes pondérées de contributions de toutes les particules suffisamment proches.




Force de pression

Nous avons utilisé comme expression de la force de pression :


f_i^{pressure} = - \rho_i \sum_j m_j ( \frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} ) \nabla W (r_i - r_j)

avec m_i masse, \rho_i masse volumique, P_i pression de la particule i.

Cette expression, dont l'intensité croît avec les P_i, traduit le fait que :

  • si des particules sont bien disposées régulièrement dans l'espace, aucune force de pression n'intervient ;
  • si des particules sont proches l'une de l'autre, la force de pression tend à les faire se repousser mutuellement ;
  • enfin, si à l'inverse des particules sont éloignées l'une de l'autre tout en restant dans leurs voisinages respectifs, alors elles vont s'attirer mutuellement.

Force de viscosité

L'expression de la force de viscosité est non nulle si et seulement si on a \langle v_i - v_j, r_i - r_j\rangle < 0. Dans ce cas on a :


f_i^{viscosity} = \rho_i \sum_j m_j \frac{A}{\rho_i + \rho_j} \frac{\langle v_i - v_j, r_i - r_j\rangle}{|r_i - r_j|^2 + \epsilon} \nabla W (r_i - r_j)

avec A une constante dépendant du fluide.

Cette force est inversement proportionnelle au carré de la distance à la particule, et traduit la résistance des particules du fluide à se mettre en mouvement.

Tension de surface

Nous avons pour expression de la tension de surface :


f_i^{tension} = - \sigma \nabla^2 c_i \frac {n_i}{| n_i |}

avec \sigma une constante, n_i la normale au fluide à la particule i et c_i une variable dite color field. Calculé dans tout l'espace par interpolation, le color field vaut 1 en le centre d'une particule et 0 ailleurs, et représente les variations de couleur sur un fluide : ainsi la surface d'un liquide aura-t-elle une couleur plus claire que le fond.

La tension de surface est proportionnelle à la divergence du vecteur n_i qui mesure la courbure d'une surface. Cette force illustre la rigidité de la surface d'un fluide lorsque celui-ci est soumis à une perturbation extérieure (c'est cette force qui permet à certains insectes de marcher sur l'eau). À l'inverse des forces de pression et de viscosité, il s'agit donc d'une force externe aux particules.

ARPS

La méthode de l'ARPS (Adaptively Restrained Particle Simulation) est une méthode générique pour un système de particules. Comme mentionné précédemment, il se base sur la mécanique hamiltonienne et lagrangienne, qui utilise pour chaque objet un certain nombre de degrés de liberté qui permettent de définir entièrement le mouvement. Par exemple pour un pendule plan, de tige rigide, il s'agit de l'angle par rapport à la normale et de la vitesse angulaire. Dans notre cas de particules se déplaçant dans l'espace en 3 dimensions, ils sont au nombre de 6 : 3 pour la position (dans les 3 dimensions), et 3 pour le moment (qui est égal au produit de la vitesse et de la masse).

Le principe de l'ARPS est de désactiver ou réactiver un certain nombre de ces degrés de liberté en fonction du mouvement de la particule à chaque pas de temps. Ici, il s'agit de restreindre le déplacement de la particule en fonction de son moment.

Méthodes de rendu graphique

Un exemple de marching cube : les points en vert vérifient f(x) \leq \alpha, les autres sont en rouge, les triangles à construire sont en orange et leurs normales en rouge.
Plusieurs méthodes de rendu graphique étaient possibles pour ce sujet. Le rendu sous forme de billes ou de points, vu précédemment, est facile à mettre en place et très efficace du point de vue des calculs ce qui le rend bien adapté au temps réel. En revanche, c'est une méthode peu esthétique, c'est pourquoi nous nous sommes également intéressés à la représentation de surfaces définies par une équation implicite (appelées métaballs).

Pour mémoire, une surface implicite S est définie par S = \{ x \in \mathbb{R}^3 | f(x) = \alpha \}, où f est une fonction (dite d'évaluation) et \alpha un réel appelé isovaleur.

Une surface implicite peut être dessinée grâce à l'algorithme des "Marching Cubes". Le principe est de discrétiser l'espace selon un maillage supposé régulier. Pour chaque point de ce maillage, on calcule la valeur de la fonction d'influence (somme des valeurs des fonctions d'évaluation de toutes les autres particules, puis on colore le sommet en "vrai" ou "faux" en fonction du résultat de la comparaison à \alpha. Ceci fait, on sélectionne des cubes formés de 8 points voisins et on trace les triangles qui correspondent à la configuration à laquelle le cube est associé. Comme chaque cube comprend 8 sommets et que chaque sommet peut prendre 2 valeurs possibles, il existe 2^8 = 256 configurations possibles (cf. un exemple à droite).

Implémentation des modèles, optimisations et résultats

Contribution aux modèles physiques

SPH

Au cours de ce projet, nous nous sommes attachés à implémenter le modèle SPH en synthétisant les contributions issues de [BT07], [Kel06] et [MCG03].

Un premier travail d'optimisation s'est avéré nécessaire à cette étape du projet. En effet, la formule du paragraphe précédent permet de constater que l'interpolation des grandeurs physiques et le calcul des différentes forces appliquées aux particules est en O(n^2), puisque pour chaque particule, on s'intéresse à toutes les autres y compris dans le cas où leur influence est nulle. Aussi nous a-t-il paru plus pertinent d'utiliser la propriété selon laquelle W était de support fini, et donc de construire une table de hachage adaptée à la recherche des plus proches voisins dans l'espace. L'implémentation se base sur une discrétisation de l'espace en une grille en 3 dimensions, dont les cubes sont appelés "voxels". Pour insérer une particule dans la table de hachage en 1 dimension, nous calculons ses coordonnées dans la grille de voxels, que nous transformons en clé entière grâce à une fonction de hachage fournie [Kel06]. Le principe de cette fonction est de donner une clé unique par particule dans la table. Pour la recherche de voisins, on utilise alors la grille pour définir un cube autour de la particule : les voisins sont recherchés parmi les clés correspondant à tous les nœuds compris dans le cube.

L'ajout de cette fonctionnalité a fait diminuer la complexité de O(n^2) à O(mn) avec m le nombre maximal de voisins, et a permis la simulation de systèmes comportant un nombre plus élevé de particules.

Ci-dessous, quelques captures d'écrans issues d'une simulation d'un fluide incompressible avec 200 particules (un nombre peu élevé).

Quelques étapes illustrant la chute d'un fluide au fond d'une boîte. À gauche, la situation initiale avec un empilement régulier de quatre colonnes de 50 particules, de positions perturbées aléatoirement. Au milieu, la simulation en cours. À droite, l'état final atteint par le système : les particules sont disposées régulièrement au fond de la boîte. On constate que la propriété de non-compressibilité est conservée.

Lorsque l'on s'intéresse à des simulations comportant un nombre plus élevé de particules, par exemple la formation d'une vague dans un bassin, on peut cependant observer quelques effets indésirables. A savoir que les particules les plus proches des bords du bassin ont tendance à s'empiler le long de ceux-ci. Même si cette situation est le plus souvent temporaire, l'affaissement de ces empilements prend du temps. Nous avons résolu ce problème en modifiant la gestion des collisions des particules contre les bords du bassin, pour certains cas précis. En augmentant le "coefficient de rebond" contre la paroi, on observe en effet une disparition de cet empilement, pour un résultat visuellement plus réaliste.

Voici les résultats obtenus avant et après cette modification :

Simulations obtenues
À gauche : sans appliquer la modification évoquée. On remarque un empilement indésirable des particules sur les bords de la boîte. Lien vers la vidéo
À droite : en appliquant cette modification. On remarque que le problème disparait. Lien vers la vidéo

ARPS

Pour ce projet, nous avons dû adapter le modèle de l'ARPS à notre système de particules de fluide, couplé à la SPH. Cela peut être vu comme un nouvel algorithme de mise à jour des particules, pour le calcul des forces et des autres grandeurs.

Adaptations à notre projet

Le grand avantage de l'ARPS ici est que cela permet de diminuer le nombre de calculs des forces à chaque pas de temps. En effet, cette fois chaque particule possède un attribut correspondant à la somme totale des forces s'exerçant sur elle. Elle est mise à jour à chaque pas de temps, mais elle n'a besoin de s'effectuer que si les particules sont actives : si des particules sont immobiles, les forces d'interaction entre elles ne sont pas modifiées. À des fins de simplicité, nous avons utilisé une matrice pour récapituler les forces, qui est triangulaire puisque les forces sont symétriques.

Après mise à jour des forces, nous allons pour chaque particule mettre à jour sa vitesse (assimilée au moment dans notre cas, puisque toutes les particules ont la même masse). Il est à noter que cette mise à jour concerne bien cette fois toutes les particules : même si elles sont immobiles, nous laissons leur vitesse évoluer, ce qui permet de changer leur état actif ou inactif. C'est la caractérisation de la méthode d'ARPS, et sa grande force.

Il suffit ensuite de calculer un certain critère sur cette valeur de la vitesse, que nous faisons intervenir dans la mise à jour de la position : si le moment est trop faible, la particule ne bouge pas. Le changement d'état actif/inactif est déterminé par des seuils que l'on peut paramétrer, et qui peuvent modifier grandement le comportement global du système. Malheureusement, aucune technique n'a été trouvée aujourd'hui pour déterminer des "bonnes" valeurs pour ces seuils, et il faut les déterminer empiriquement en fonction du résultat que l'on veut obtenir.

En dernier lieu, nous avons pensé à utiliser la même technique pour la mise à jour de la masse volumique. En effet, la valeur de la masse volumique d'une particule dépend des contributions de ses voisines, uniquement en fonction de leur distance. Nous pouvons donc intégrer le calcul de cette grandeur dans la boucle de mise à jour des forces, qui ne concerne que les particules actives. Cela évite de parcourir toutes les particules, et pour chacune de recalculer leur voisinage.

Résumé de l'algorithme

Illustration de l'ARPS : les particules en bleu sont actives, celles en rouge sont inactives, celles en vert sont en transition entre les deux états.
Cette vidéo est une application de l'ARPS sur l'exemple de la vague.

Boucle d'animation : à chaque pas de temps, nous faisons

  1. Mise à jour des forces et des grandeurs des particules
    • Au premier pas de temps : pour toutes les particules,
    1. Calculer la masse volumique, et la pression
    2. Calculer les forces d'interaction avec les autres particules à cet instant initial
    3. Ajouter les forces constantes (gravité, qu'on ne recalculera plus)
    • Pour les autres pas de temps : uniquement pour les particules actives,
    1. En se basant sur les anciennes positions et vitesses,
      1. Ôter les anciennes contributions à la masse volumique
      2. Enlever les forces qui s'appliquaient à la particule au pas de temps précédent
    2. Mettre à jour la table des voisins avec les nouvelles positions
    3. En se basant sur les positions et vitesses du pas de temps présent,
      1. Ajouter les contributions des particules à la masse volumique, et mettre ainsi à jour la pression
      2. Ajouter les forces d'interaction s'exerçant au pas de temps présent
  2. Réinitialisation de la liste des particules actives
  3. Mise à jour de la position et de la vitesse : pour toutes les particules,
    1. Mettre à jour la vitesse :  \vec v (t+dt) = \vec v (t) + dt * \vec F
    2. Calculer le critère d'activité correspondant :  \rho
    3. Si  \rho < 1 , insérer la particule dans la liste des actives
    4. Mettre à jour la position ; si  \rho = 1 , elle n'est pas modifiée
    5. Si la particule est active, gérer les collisions :
      1. S'il y a collision, modifier la position et la vitesse en conséquence
      2. Calculer une nouvelle fois le critère  \rho
      3. Si  \rho = 1 , retirer la particule de la liste des actives

Difficultés rencontrées

Nous avons rencontré quelques difficultés pour adapter l'ARPS à notre sujet. En effet, un point crucial de cet algorithme pour obtenir un système stable, est que les anciennes forces enlevées, dans un premier temps, soient bien exactement égales à celles qui opéraient au pas de temps précédent.

  • Il a donc fallu bien s'assurer que toutes les forces d'interaction, calculées entre deux particules, étaient symétriques. Cela concerne toutes les forces, qu'elles soient externes ou internes, tant qu'elles s'appliquent entre deux particules. Ainsi par exemple, nous avons rendu la force de tension de surface symétrique entre deux particules, par remplacement des grandeurs de chaque particule par la moyenne des deux.
  • Une autre difficulté a donc été la gestion des collisions. Nos particules ne peuvent pas vraiment être considérées comme des entités au sens propre, il n'y a pas de collision entre elles : l'équilibre des positions est atteint grâce aux forces d'interaction qui agissent entre les particules. Les collisions concernent les obstacles extérieurs au fluide, comme les parois d'une boîte. Nous les gérons une fois que la collision est détectée, en repositionnant la particule au point de collision, avec ajout d'une vitesse "miroir" de la vitesse précédant la collision, diminuée par un coefficient de restitution (de l'énergie). Cette gestion se fait donc juste après la mise à jour des positions, qui n'est effective que pour les particules actives. Or, il arrive qu'après gestion des collisions, qui change la vitesse, la particule redevienne inactive, et ce dans le même pas de temps. Ces changements ont été assez compliqués à gérer, mais la mise en place de la matrice des forces a bien résolu le problème.

Résultats

Le principal gain de cet algorithme en terme de performances se fait dans la mise à jour des forces, qui ne s'applique qu'aux particules actives. Dans le cas de notre scène finale, nous avons placé des bassins, dans lesquels de nombreuses particules finiront par atteindre un état stable, et donc inactif.

Cependant après quelques essais, nous avons été déçus du résultat, car il y avait trop peu de particules inactives; la vitesse d'exécution pour la même scène était même plus longue avec l'ARPS que sans! Nous avons donc décidé de forcer un peu la main, pour rendre plus facilement les particules inactives. Le passage d'un état à un autre se décide à partir de la valeur du critère, qui vaut entre 0 et 1 compris : une particule est inactive si et seulement si le critère vaut 1, elle est complètement active s'il vaut 0, et en transition sinon. Pour que plus de particules deviennent inactives, nous avons donc triché un peu, et mis un seuil de 1 - \epsilon, au-dessus duquel nous rendons les particules inactives. Le résultat est concluant, puisque nous avons avec cette technique gagné du temps d'exécution par rapport au SPH simple : rassurant!

État final avec l'ARPS. À gauche sans l'amélioration, à droite avec notre amélioration

Comparaison des deux méthodes

L'implémentation de l'ARPS avait pour objectif d'augmenter la vitesse d'exécution du programme. Nous avons lancé plusieurs simulations sur la scène de cascade finale, avec différents nombres de particules, pour mesurer l'éventuel gain de temps de l'ARPS. Grâce à notre petite correction, nous avons obtenu des résultats convaincants.

Comparaison des temps de calcul
Configurations SPH ARPS Gain de temps
5.000 particules
10.000 itérations
42 min 32 min 25%
10.000 particules
10.000 itérations
A planté 4h07 /
2.000 particules
50.000 itérations
2h39 1h52 30%

Ces simulations ont été lancées avec une seuil inférieur d'activation de 1, et d'inactivation de 2. Nous remarquons même que pour notre simulation la plus importante, avec la SPH elle n'est pas allée au bout, alors qu'avec l'ARPS elle s'est terminée sans problème.

Rendu graphique

Exemple de blob obtenu grâce à l'algorithme des "marching cubes", assimilable à la fusion de deux gouttes d'eau.
Dessiner un fluide sous la forme de sphères ou de points centrés en les particules composant le système fut aisé. Toutefois, nous avons cherché à obtenir un résultat plus réaliste en dessinant non pas des particules, mais la surface implicite délimitant le fluide de l'environnement extérieur.

Dans notre cas, cette surface S est définie à partir des particules P_i du fluide, qui jouent ici le rôle d'un squelette de points. Pour tout point de l'espace, la fonction d'influence f_i associée à la i-ème particule du fluide est une fonction décroissante de la distance d(M, P_i) à la particule. En sommant les f_i en tout point de l'espace pour obtenir la fonction d'évaluation f, ce sont toutes les surfaces locales que l'on entend fusionner en une surface lisse qui délimiterait les contours du fluide.

L'algorithme des "marching cubes" présente l'avantage de dessiner des surfaces dans l'espace pour tout degré de résolution et de posséder plus de réalisme. Toutefois, en l'état, il possède un coût en temps non négligeable qui le rend inadapté à une application temps-réel.

Nous avons toutefois refondu le code de l'algorithme en vue d'une optimisation. En effet, dans notre cas, la fonction d'influence f associée à une particule est de support fini. A chaque fois qu'une particule se déplace, ce ne sont donc que certaines particules qui sont affectées par une augmentation ou un déficit d'influence. Ainsi, il est inutile de recalculer la valeur de cette fonction pour toute la grille. Au lieu de considérer tous les points du maillage, nous avons donc préféré considérer toutes les particules en mouvement, les cases sur lesquelles leur mouvement influe et seulement celles-ci. En d'autres termes, pour chaque point dans le voisinage de départ de la particule, on va soustraire la valeur de la fonction d'influence en ce point (puisque la particule s'éloigne de ces points) et pour chaque point dans le voisinage d'arrivée, on ajoute la valeur de la fonction d'influence (puisque la particule s'en rapproche).

Cette méthode donne certes de meilleures performances, mais elles restent insuffisantes pour une application en temps-réel dans le cas d'un système complexe. Autrement dit, dans le cas d'une cascade ou d'un gros bloc de particules, l'animation doit faire l'objet d'un pré-calcul.

Au sujet du choix de la fonction d'évaluation, plusieurs familles de fonctions ont été traitées : inverse de carré, noyau exponentiel ou puissance, courbe de Hermite,... Les rendus ne sont pas toujours optimaux. Dans le rendu final, nous avons sélectionné un noyau de Cauchy de la forme f(x) = \frac{a}{1 + b x^2}, avec a et b paramètres réels.

Un exemple de résultat lors de la formation d'une vague dans une boîte :

Surface implicite obtenue lors de la formation d'une vague dans un récipient (2500 particules).
Lien vers la vidéo

Cette visualisation est correcte et visuellement acceptable dans le sens de la hauteur. Le résultat n'est cependant pas très beau sur les côtés : on aimerait que visuellement, le fluide donne l'impression d'être clairement délimité par la boîte. C'est pourquoi nous avons également implémenté l'algorithme dit des "marching squares" pour représenter les bords du fluide. Son fonctionnement est très similaire à celui de l'algorithme des "marching cubes" : pour chaque carré délimité par quatre points adjacents du maillage, on regarde si la valeur de la fonction d'influence en ces points est inférieure ou supérieure à l'isovaleur qui a été utilisée pour représenter la surface du fluide. On colore chaque sommet par rapport au résultat de cette configuration, puis on trace le jeu de triangles qui correspond à la configuration du carré. Chaque sommet pouvant prendre deux valeurs (0 ou 1) et chaque carré possédant quatre sommets, il existe 2^4 = 16 configurations possibles.

Voici le résultat que l'on obtient avec l'intégration de cette fonctionnalité :

Surface implicite obtenue lors de la formation d'une vague dans un récipient en intégrant la modification avec les "marching squares" (2500 particules).
Lien vers la vidéo

Le résultat est bien plus réaliste ainsi !

Scène finale : écoulement d'eau dans une cascade

Pour notre scène de démonstration, nous voulions bien mettre en évidence les gains de l'ARPS. Nous avons donc choisi l'écoulement d'une cascade, avec déversements de bassins. En effet dans ces bassins, il est facile d'avoir un certain nombre de particules inactives, qui vont augmenter les performances du programme par rapport au SPH traditionnel.

Au début de la simulation, nous faisons tomber un grand nombre de particules dans le premier bloc. Afin de réguler le flux, nous avons mis en place un palier physique dans la boîte. Les particules en attente de tomber sont donc largement inactives là encore. Voici les résultats de nos simulations :

Problèmes rencontrés

La principale difficulté de ce projet résidait dans la recherche des bons paramètres pour que le fluide ait un comportement réaliste (masse volumique, viscosité, rigidité, etc.). En effet, dans les documents de recherche que nous avons étudiés, des valeurs théoriques sont mentionnées mais aucune information sur la façon dont elles ont été obtenues n'est donnée, si tant est que des valeurs soient données. Cela conduisait souvent, dans notre implémentation, à des résultats fantaisistes (explosion du fluide par exemple). Il a donc fallu que nous trouvions les bonnes valeurs nous-mêmes. Et si parfois cela est faisable grâce à quelques calculs, d'autres fois il a fallu "tâtonner" pour arriver au bon résultat ; un travail fastidieux et décourageant.

D'autre part, nous programmions en grande partie sur des PC de l'Ensimag et avons été victimes de nombreux manques de logiciels. Nous aurions, par exemple, souhaité faire du calcul parallèle sur GPU pour accélérer les calculs de façon importante ou encore importer des modèles de 3ds Max dans notre code pour obtenir une belle cascade dans notre scène finale, mais les librairies nécessaires n'étaient pas installées.

Enfin, le sujet étant très large et recouvrant différents domaines, il nous a fallu faire des choix quant aux objectifs que nous allions nous fixer. Il était parfois difficile d'estimer le temps qu'allaient nous prendre les différentes tâches ou encore d'abandonner certaines pistes à cause de la contrainte du temps-réel.

Conclusions

Ce sujet portant sur des domaines très variés de l'informatique et des mathématiques (modèles physiques, rendu, surfaces implicites, optimisations, etc.), chacun d'entre nous a su y trouver de quoi s'épanouir. Cela nous a permis de garder un rythme régulier tout au long de ce projet et de bien avancer ; nous possédions déjà un premier livrable théorique (SPH et rendu par metaballs fonctionnels - mais non optimisés) une semaine avant le rendu final.

La scène de la cascade présentée ci-dessus (accessible via ce lien) représente un bon condensé de ce que nous avons implémenté durant ce mois de projet: un modèle physique cohérent avec une bonne interaction entre les particules, une méthode de rendu par metaballs et des optimisations efficaces (cf. comparaison des temps de calcul avec et sans ARPS).

Bien que nous regrettions certains points présentés dans les problèmes rencontrés, notamment l'impossibilité d'utiliser certaines librairies et d'aller plus loin dans le sujet, nous sommes globalement très satisfaits par ce sujet. Outre les connaissances qu'il nous a apportées sur l'animation de fluides, l'un des grands domaines de recherche en graphique 3D, nous avons été séduits par sa diversité et son aspect concret (le but étant la création d'une cascade). De plus, contrairement à la plupart des projets de l'Ensimag où nous sommes souvent très guidés, avec des squelettes et des bouts de codes imposés, nous sommes enchantés d'avoir été libres sur ce travail et d'avoir su créer une scène complexe en partant de zéro.

Le sujet de ce projet étant très large, nous pouvons, pour le futur, imaginer de multiples autres choses à développer : application à d'autres fluides (plus visqueux par exemple, qui colle aux parois), ajouter de la transparence, trouver d'autres optimisations, etc. En particulier, nous avons commencé l'implémentation de l'ajout d'écume sur notre fluide, en changeant la couleur des particules en fonction de leur nombre de voisins et de leur vitesse. Vous pouvez voir cette vidéo ici ou ici.

Paramètres physiques utilisés

Se référer au document de [Kel06] pour les explications des paramètres.

Valeurs utilisées dans nos simulations
Paramètres Valeurs
Pas de temps 0.001 s
Accélération de la gravitation (0, 0, -9.82) m.s-2
Température 293.15 K
Pression atmosphérique 101 325 Pa
Densité au repos 160 kg.m-3
Masse 0.02 kg
Coefficient de viscosité 3.5 Pa.s
Tension de surface 0.0728 N.m-1
Seuil pour la tension de surface 7.065
Rigidité du gaz 3 J
Coefficient de restitution 0.3
Nombre de particules dans le noyau 20
Rayon du noyau 0.03 m
Célérité du son dans l'eau 88.5 m.s-1
Constante alpha (viscosité) 0.08

Références

[AR12] Svetlana Artemova and Stephane Redon. "Adaptively restrained particle simulations". Physical review letters, 109(19) :190201, 2012.

[AR12] (Supplément) Svetlana Artemova and Stephane Redon. "(Supporting information for) Adaptively restrained particle simulations". Physical review letters, 109(19) :190201, 2012.

[BT07] Markus Becker and Matthias Teschner. "Weakly compressible SPH for free surface flows". In Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer anima-tion, pages 209-217. Eurographics Association, 2007.

[Kel06] Micky Kelager. "Lagrangian fluid dynamics using Smoothed Particle Hydrodynamics". University of Copenhagen : Department of Computer Science, 2006.

[LC87] William E Lorensen and Harvey E Cline. Marching cubes: A high resolution 3D surface construction algorithm. In ACM Siggraph Computer Graphics, volume 21, pages 163-169. ACM, 1987.

[MCG03] Matthias Müller, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications". In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 154-159. Eurographics Association, 2003.

[Mon05] Joe J Monaghan. "Smoothed particle hydrodynamics". Reports on progress in physics, 68(8) :1703, 2005.

Documents additionnels

Fichier:Soutenance Animation Temps Reel Liquides.pdf

Fichier:Sources Animation Temps Reel Liquides.tar.gz