La présente invention concerne un procédé et un dispositif de simulation interactive du contact entre au moins un premier objet déformable et au moins un second objet avec un pas prédéterminé de temps d'échantillonnage d'un modèle simulé. On a déjà proposé d'effectuer une simulation de mesures d'interpénétration entre un objet rigide et un objet déformable à partir d'estimations volumiques ou de distance, notamment pour des applications de chirurgie virtuelle où un outil virtuel rigide de chirurgie coopère avec un organe virtuel déformable du corps humain. Toutefois, selon ces méthodes la relation entre la mesure d'interpénétration et les forces de contact de réaction n'ont pas de base physique et des forces artificielles peuvent s'appliquer à des nœuds des maillages des objets qui ne sont pas en contact, ce qui nuit à la fiabilité, ou les forces de contact ne répondent pas aux conditions du problème de Signorini. La présente invention vise à remédier aux inconvénients précités et à permettre de réaliser une simulation interactive en temps réel du contact entre des objets, dont au moins certains sont déformables, de façon simplifiée et économique tout en respectant les contraintes des lois physiques qui régissent les contacts, de telle manière que les contacts simulés entre des objets soient fiables et qu'ainsi la stabilité de la simulation soit garantie. Ces buts sont atteints grâce à un procédé de simulation interactive du contact entre au moins un premier objet déformable et au moins un second objet avec un pas prédéterminé de temps d'échantillonnage d'un modèle simulé,

caractérisé en ce que :

(a) on calcule au préalable les paramètres décrivant les caractéristiques physiques de chacun des objets, tels que la géométrie et la mécanique des matériaux de chacun des objets, et on stocke ces paramètres dans une mémoire,

(b) au début de chaque pas de temps d'échantillonnage du modèle simulé, on procède au niveau de chaque objet à une analyse en temps réel du comportement propre de l'objet pour prédire les positions, vitesses et accélérations de cet objet selon un mouvement libre qui ne tient pas compte d'éventuels contacts ultérieurs,

(c) à chaque pas de temps d'échantillonnage du modèle simulé, on analyse en temps réel, au niveau d'une scène globale comportant les objets susceptibles de venir en contact, des paires d'objets qui sont détectés en intersection, et on établit une liste de groupes de collisions qui contient une chaîne d'objets en collision et une description des collisions,

(d) à chaque pas de temps d'échantillonnage du modèle simulé, on rapatrie en temps réel, pour chaque groupe de collisions, des paramètres représentant les caractéristiques physiques des objets et la description des collisions, de manière à déterminer, pour chaque cas, la solution au problème de Signorini qui gouverne le contact entre deux objets dans le cas d'un glissement relatif pur,

(e) à la fin de chaque pas de temps d'échantillonnage du modèle simulé, on procède au niveau de chaque objet à une visualisation en temps réel du comportement propre de l'objet suite à la collision, et

(f) l'ensemble des traitements en temps réel s'effectue avec un pas de temps de calcul plus court que le pas de temps d'échantillonnage du modèle simulé, de manière à définir une simulation interactive où l'utilisateur peut intervenir directement en cours de simulation.

Lors de l'étape a) de calcul préalable des paramètres décrivant les caractéristiques physiques de chacun des objets, on utilise pour les paramètres décrivant la mécanique des matériaux une description des déformations de type éléments finis, avec le remplissage et l'inversion de matrices, la résolution de systèmes d'équation et le stockage des données en mémoire. Selon un mode particulier de réalisation, chaque objet est décrit dans une configuration au repos comme un ensemble de triangles reproduisant sa surface et un ensemble de tétraèdres décrivant l'intérieur de l'objet. Avantageusement, chaque triangle est décrit par trois points, placés dans un ordre qui permet de calculer des normales qui sont invariablement dirigées vers l'extérieur de l'objet. De préférence, les déformations des objets sont interpolées par la méthode des éléments finis en utilisant un maillage tétraédrique linéaire. A chaque pas de temps de calcul, on intègre lors de l'étape b) au niveau d'un objet, les forces explicites appliquées sur l'objet, qui sont déjà connues au départ du pas de calcul, de manière à définir le mouvement qu'elles créent sur l'objet, tandis que la valeur des forces de contact implicites, qui dépendent elles-mêmes du mouvement des objets dans le pas de temps de calcul, est déterminée lors de l'étape d) de recherche au niveau d'une scène globale, de la solution au problème de Signorini. Lors de l'étape c) d'analyse au niveau d'une scène globale, on détecte géométriquement les intersections existantes entre les objets de la scène pour extraire des couples d'éléments d'objets en intersection, une longueur et une direction d'interpénétration entre les deux éléments d'un couple d'éléments d'objets. Selon une variante de réalisation, lors de l'étape c) d'analyse au niveau d'une scène globale, pour extraire des couples d'éléments d'objets en intersection, une longueur et une direction d'interpénétration entre les deux éléments d'un couple d'éléments d'objet, on prend également en compte un mouvement intermédiaire des objets entre le pas de calcul précédent et le pas de calcul courant, pour calculer une direction d'interférence privilégiée entre ces objets. Lors de l'étape d) de recherche de la solution au problème de Signorini, on reconstruit les points extrêmes d'application de la force de contact entre deux objets soumis à une collision lorsque ces points extrêmes d'application n'ont pas été déterminés à l'étape précédente. Selon un mode de réalisation particulier, lors de l'étape d) lors de l'étape d) de recherche de la solution au problème de Signorini, dans le cas d'une intersection segment-segment de deux objets en triangle, les deux points choisis pour constituer les points extrêmes d'application de la force de contact entre les deux objets soumis à une collision sont situés à l'intersection de chacun des deux segments avec le plan formé par la face du triangle en intersection. Selon un autre mode de réalisation particulier, lors de l'étape d) de recherche de la solution au problème de Signorini, dans le cas d'une intersection point-face de deux objets en triangle, un premier point choisi pour constituer un point extrême d'application de la force de contact entre les deux objets soumis à une collision est le point de l'intersection tandis que le second point extrême d'application de la force de contact entre les deux objets soumis à une collision est la projection du premier point extrême sur la face du triangle en intersection. Selon un aspect particulier de la présente invention, on utilise les coordonnées barycentriques pour répartir les déplacements et les forces des points d'application de la force de contact entre les points extrêmes d'application de la force de contact en effectuant une interpolation linéaire pour une modélisation en éléments finis. De façon plus particulière, on peut calculer la distance δ d'interpénétration entre les deux points extrêmes d'application de la force de contact dans le cas d'un contact segment-segment entre un premier segment et un second segment d'un second triangle à partir de l'équation suivante :

δ = [a, b, cj [a 1 451 - l - β] V1 W. V. (1)

où : α et 1-α sont les coordonnées barycentriques sur le premier segment,

β et 1-β sont les coordonnées barycentriques sur le deuxième segment,

ai bj q sont les coordonnées de la direction ni d'interpénétration,

Wi et W2 sont les coordonnées du premier segment,

Vi et V2 sont les coordonnées du second segment.

On peut également calculer la distance δ d'interpénétration entre les deux points extrêmes d'application de la force de contact dans le cas d'un contact point-plan entre un point d'un second triangle et un plan d'un premier triangle à partir de l'équation suivante : "W1 " δ = [a, b, cj [a β y] W2 - V1 (2) W3

où : α, β et γ sont les coordonnées barycentriques sur le premier triangle,

ai bj q sont les coordonnées de la direction nj d'interpénétration,

Wi, W2, W3 sont les coordonnées du premier triangle,

Vi représente les coordonnées du point de contact constitué par un sommet du second triangle. Après avoir déterminé les points d'application des forces de contact entre deux objets en collision, lors de l'étape d) on procède au transfert des caractéristiques mécaniques des objets dans l'espace défini des contacts dans lequel on traite l'ensemble d'un groupe de m contacts avec n objets où m et n sont des entiers. De façon plus particulière, lors de l'étape d) on considère la masse et l'inertie d'un objet de façon globale, en son centre de gravité et on établit une relation instantanée entre les forces de contact fc dans la direction du contact, les accélérations δ"c dues aux contraintes dans la même direction et les accélérations libres δ"ιlbre dans la même direction connues lors de l'étape c) au niveau d'une scène globale, selon l'équation suivante :

<?'„ = J. MΓ'J/ £.+ <?'_ libre (3)

où : Jc est une matrice jacobienne m*6n qui transfère le mouvement instantané linéaire et angulaire dans l'espace des contacts,

3CT est la matrice transposée de Jc,

M est une matrice bloc diagonale correspondant à la masse et l'inertie de n objets du groupe de contacts. Lors de l'étape d) pour le transport des caractéristiques mécaniques locales, on établit une relation entre : • la différence de déplacement (Uk1) des points du maillage déformable représentant l'objet i à l'instant k, entre la déformation libre (U'uiure) et la déformation contrainte (U^c) soit Uk' = U1K1C - U'k.ϋbre

• Les positions relatives des objets, libres et contraintes, dans l'espace des contacts : δπbre et δc. δ = ∑W Nc1 Ukj + δϋbre (4)

où Nc est une matrice de passage de l'espace des déplacements du maillage vers l'espace des déplacements aux contacts. De même, on établit une relation entre les forces dans l'espace des contacts /c et les forces dans l'espace des forces de déformation Fk. Fk = (NcO1Zc (5)

De façon plus particulière, lors de l'étape d) on établit une relation linéaire instantanée de caractérisation des déformations ou déplacement de contact δc à partir des forces de contact fc et des déplacements libres δiibre dus aux mouvements libres intégrant uniquement les forces connues de façon explicite en début de pas de temps de calcul, selon l'équation suivante :

Sc = [∑i=1n N; A (IU (Nc!)τ ] fc+ δlibre (6)

où : N[. est une matrice de passage de l'espace des déplacements du maillage vers l'espace des déplacements de contacts,

(N n^\)T est la matrice transposée de N'c

A est une matrice permettant de définir la déformation de l'objet au niveau local, de telle sorte que si Uk représente le vecteur du déplacement dans le repère local de l'objet à l'instant courant et Uk-i représente le vecteur de déplacement dans le repère local de l'objet au pas de calcul précédent dont les valeurs instantanées sont connues au début du pas de calcul courant, on a :

UK = A (Uk-i) Fk + b (Uk-i) (7)

où Fk est un vecteur représentant les forces extérieures appliquées à l'objet exprimées dans le repère local, et b est un vecteur qui a une valeur dans l'espace des déplacements et qui dépend du modèle de déformation de l'objet. Dans un cas plus général, lors de l'étape d) on établit une relation instantanée de caractérisation des déformations ou déplacements de contact δc à partir des forces de contact fc et des déplacements libres δπbre dus aux mouvements libres intégrant uniquement les forces connues de façon explicite en début de pas de temps de calcul, selon l'équation suivante :

Sc = [θ dt2 Jc MT1JZ+V N; A (U,,) (Nc')τ ] /β+ δUhre (8)

où : Jc est une matrice jacobienne m*6n qui transfère le mouvement instantané linéaire et angulaire dans l'espace des contacts,

3CT est la matrice transposée de Jc'

M est une matrice bloc diagonale correspondant à la masse et l'inertie des n objets du groupe de contacts,

θ est une constante dépendant de la méthode d'intégration en temps,

N^ est une matrice de passage de l'espace des déplacements du maillage vers l'espace des déplacements de contacts,

(N^.)τ est la matrice transposée de N^.,

A est une matrice permettant de définir la déformation de l'objet au niveau local, de telle sorte que si Uk représente le vecteur du déplacement dans le repère local de l'objet à l'instant courant et Uk-i représente le vecteur de déplacement dans le repère local de l'objet au pas de calcul précédent dont les valeurs instantanées sont connues au début du pas de calcul courant, on a :

UK = A (Un) Fk + b (Un) (7)

où Fk est un vecteur représentant les forces extérieures appliquées à l'objet exprimées dans le repère local, et b est un vecteur qui a une valeur dans l'espace des déplacements et qui dépend du modèle de déformation de l'objet.

On peut ainsi remarquer que dans le cas d'un objet rigide indéformable, on a Uk = Uk-i qui traduit l'absence de modification en fonction du temps du vecteur Uk.

Avantageusement, le procédé selon l'invention comprend en outre une étape de couplage avec un module d'interface haptique pour produire un retour de sensation haptique sur un dispositif mécanique par lequel un opérateur manipule les objets dans une scène virtuelle.

L'invention concerne également un dispositif de simulation interactive du contact entre au moins un premier objet déformable et au moins un second objet avec un pas prédéterminé de temps d'échantillonnage d'un modèle simulé, caractérisé en ce qu'il comprend :

(a) un module de calcul préalable des paramètres décrivant les caractéristiques physiques de chacun des objets, tels que la géométrie et la mécanique des matériaux de chacun des objets,

(b) une mémoire de stockage des paramètres préalablement calculés dans le module de calcul,

(c) un module de couplage avec une interface utilisateur comprenant un dispositif mécanique tenu par un utilisateur lui permettant d'exercer virtuellement des forces sur lesdits objets dans une scène du modèle simulé,

(d) un écran de visualisation pour visualiser lesdits objets représentés sous forme de maillages,

(e) une unité centrale de traitement associée à des moyens d'entrée, comprenant au moins

el) un module d'analyse d'objet pour analyser en temps réel au niveau de chaque objet le comportement propre de l'objet pour prédire les positions, vitesses et accélérations de cet objet selon un mouvement libre qui ne tient pas compte d'éventuels contacts ultérieurs,

e2) un modèle d'analyse d'une scène globale comportant les objets susceptibles de venir en contact, pour analyser en temps réel des paires d'objets qui sont détectés en interaction et établir une liste de groupes de collisions qui contient une chaîne d'objets en collision et une description des collisions,

e3) un module de rapatriement en temps réel, pour chaque groupe de collisions, des paramètres représentant les caractéristiques physiques des objets et la description des collisions pour déterminer, pour chaque cas, la solution au problème de Signorini qui gouverne le contact entre deux objets dans le cas d'un glissement relatif pur,

e4) un module de traitement de chaque objet pour permettre en temps réel au niveau de chaque objet une visualisation en temps réel du comportement propre de l'objet suite à une collision, et

e5) des moyens de détermination d'un pas de calcul plus court que le pas de temps d'échantillonnage du modèle simulé de manière à définir une simulation interactive.

Avantageusement, le dispositif comprend des moyens de production d'un retour de sensation haptique sur l'interface utilisateur. Selon une caractéristique avantageuse, le pas de calcul correspond à une fréquence égale ou supérieure à environ 500 Hz. D'autres caractéristiques et avantages ressortiront de la description suivante de modes particuliers de l'invention, donnés à titre d'exemples, en référence aux dessins annexés, sur lesquels : - la Figure 1 est un schéma montrant les différentes étapes d'un procédé de simulation interactive du contact entre objets, selon l'invention, - la Figure 2 est un schéma montrant différents niveaux de traitement de l'interaction entre objets au cours de différentes étapes du procédé de simulation de la Figure 1, - les Figures 3A à 3C représentent trois exemples d'interaction entre deux objets virtuels représentés par des triangles, - la Figure 4 représente une interaction entre deux objets en triangle dans le cas d'une intersection segment-segment, - la Figure 5 représente une interaction entre deux objets en triangle dans le cas d'une intersection point-face, - la Figure 6 illustre schématiquement le cas d'une collision entre deux objets pour laquelle la configuration peut être définie à partir de seuls critères géométriques, - la Figure 7 illustre schématiquement le cas d'une collision entre deux objets pour laquelle la configuration est définie en tenant compte d'un mouvement intermédiaire, - la Figure 8 est un schéma-bloc montrant les constituants de base d'un dispositif de simulation interactive du contact entre objets, selon l'invention, - la Figure 9 montre un exemple de contact entre un objet virtuel déformable et un autre objet virtuel, et - les Figures 1OA à 1OC montrent trois positions relatives différentes entre un objet virtuel déformable, en forme de pince, et un objet virtuel rigide lors d'un processus de mise en place de l'objet virtuel déformable en forme de pince sur l'objet virtuel rigide. La Figure 8 illustre de façon schématique un exemple de dispositif permettant de mettre en œuvre l'invention et de réaliser la simulation interactive en temps réel du contact entre objets tout en permettant notamment d'avoir un retour de sensation haptique. Une unité centrale de traitement 100, qui peut être constituée à partir d'un ordinateur de type classique, permet d'effectuer les différents calculs nécessaires pour réaliser une simulation. Un écran de visualisation 107 relié à l'ordinateur 100 par une interface graphique permet l'affichage d'objets représentés sous la forme de maillage comprenant des nœuds ou sommets reliant des segments ou arêtes. Des informations peuvent être fournies à l'ordinateur 100 à partir d'une interface utilisateur classique 103 pouvant comprendre un clavier et par exemple une souris et constituant des moyens d'entrée. Un dispositif mécanique spécifique 104 tenu par un utilisateur relié par un module de couplage 101 à l'ordinateur 100 peut en outre être prévu pour permettre à l'utilisateur d'exercer virtuellement des forces sur les objets dans une scène d'un modèle simulé. Un tel dispositif mécanique 104 et le module de couplage 101 constituent une interface haptique qui permet à l'utilisateur d'exercer une sollicitation sur les objets virtuels de la scène et de recevoir en retour une simulation haptique qui est une réponse fournie par la simulation du contact entre objets. L'ordinateur 100 comprend de façon classiquement au moins un processeur, une mémoire permanente de stockage de programmes et de données et une mémoire de travail coopérant avec le processeur. Des supports externes de mémoire (disquettes, CD-ROM,...) ou un modem de liaison à un réseau peuvent naturellement être utilisés pour charger dans l'ordinateur des programmes ou données permettant d'effectuer tout ou partie du traitement de simulation. Sur la Figure 8, on a simplement représenté à titre symbolique un exemple de mémoire de stockage 102 coopérant avec le module 100 et pouvant être constituée par l'un ou l'autre des types de mémoires sus-mentionnés. D'une manière générale, au début d'un processus de simulation les paramètres décrivant la géométrie et la mécanique des matériaux des objets à simuler sont calculés dans l'unité centrale 100 et stockés dans une zone mémoire de la mémoire 102. Pour caractériser les déformations mécaniques des objets, il est utilisé lors du traitement par l'unité centrale 100 une description des déformations de type éléments finis. Cela se traduit par le remplissage et l'inversion de matrices, la résolution de systèmes d'équations et le stockage des données dans la mémoire 102 associée à l'unité centrale 100. Les positions et formes courantes des objets sont évaluées en fonction des sollicitations exercées et des lois mécaniques qui régissent les objets de la scène. Selon l'invention, pour garantir une simulation stable, il est pris en compte dans le calcul des contacts simulés entre objets les lois physiques qui régissent le contact. Pour permettre une simulation en temps réel, c'est-à-dire avec un délai très court et borné séparant la sollicitation exercée par l'utilisateur via l'interface haptique 104, 101 et la réponse fournie par l'unité centrale de traitement 100 à cette interface haptique 104, 101, le dispositif de simulation met en œuvre trois modules principaux qui sont sollicités itérativement à chaque pas de temps d'échantillonnage du modèle simulé. Par ailleurs, l'ensemble des traitements en temps réel s'effectue avec un pas de temps de calcul plus court que le pas de temps d'échantillonnage du modèle simulé. Les trois modules principaux du dispositif de simulation mettant en œuvre les différentes étapes du procédé de simulation se présentent essentiellement de la façon suivante : un premier module dit "mécanique" se situant au niveau de chaque objet et décrivant son comportement propre, permet de faire évoluer la position et la forme de l'objet suivant les forces et les lieux des forces exercés. Ce module est appelé en début de pas de calcul pour prédire les positions, vitesses et accélérations des objets sans tenir compte du contact puis sera de nouveau mobilisé pour tenir compte des forces calculées dans un troisième module dit "de traitement du contact". Un second module dit de "détection de collision", se situant au niveau de la scène globale établit des paires d'objets qui sont détectés en intersection. Ce module, de façon optionnelle, peut créer des mouvements intermédiaires entre les pas de la simulation pour savoir quand et comment les objets sont entrés en intersection. Ce module est régi avant tout par des lois géométriques optimisées qui permettent d'accélérer le calcul afin d'obtenir une chaîne d'objets en collision et une description des collisions. Un groupe de collision est ainsi un ensemble d'objets reliés entre eux par au moins une collision. Un objet entre dans un groupe s'il est en collision avec au moins un des objets du groupe. Une collision est décrite obligatoirement par la paire d'objets en collision et par le lieu de la collision en utilisant soit les éléments géométriques de base (par exemple deux triangles ou deux surfaces) en intersection, soit par un segment reliant les deux points qui sont localement les plus interpénétrés. Un troisième module dit "traitement du contact" est appelé par le module de "détection de collision" et fait appel en retour au module "mécanique". Pour chaque groupe de collisions, le module de traitement du contact rapatrie les caractéristiques physiques des objets et la description des collisions. Le module est apte à déterminer, pour chaque cas, la solution au problème de Signorini qui gouverne le contact entre deux objets dans le cas du glissement pur. L'invention permet d'effectuer une simulation interactive. Une simulation est définie par le pas de temps d'échantillonnage du modèle simulé et par le pas de temps de calcul de ce modèle. Le procédé selon l'invention met en œuvre un pas de temps de calcul qui est toujours inférieur au pas de temps choisi pour l'échantillonnage, ce qui permet d'avoir une simulation interactive où l'utilisateur va pouvoir intervenir directement en cours de simulation. La Figure 1 résume les principales étapes du procédé selon l'invention qui met en œuvre une boucle de simulation utilisant les trois modules principaux précités installés dans l'ordinateur 100 de la Figure 8. La Figure 2 illustre les différents niveaux de traitement entre objets au cours des différentes étapes du procédé de simulation. Une première étape de traitement 130 utilise le module dit "mécanique" et se situe au niveau de chaque objet (niveau objet 3). Les informations sont fournies à travers un module de couplage 120 à partir de l'interface utilisateur 110 ou interface haptique qui détermine la position et la forme de chaque objet (informations 135 élaborées à l'étape 130). Dans cette première étape 130, une modélisation 13 prend en compte chaque objet ou outil 201, 202, 203 de façon individuelle sans tenir compte d'interactions ultérieures éventuelles et permet de faire évoluer la position et la forme de l'objet suivant les forces et les lieux des forces exercés à partir de l'interface utilisateur 110. Une deuxième étape du traitement 140 utilise le module dit "détection de collision" et se situe au niveau d'une scène globale (niveau scène 4). Les informations 135 élaborées à l'étape 130 sont utilisées à l'étape 140 pour établir des paires d'objets qui sont détectés en intersection. Au cours de cette étape 130, il est élaboré une liste de groupes de collisions (informations 145) contenant une chaîne d'objets en collision et une description des collisions. Dans cette étape 140, une modélisation 14 prend ainsi en compte une paire d'objets en intersection tels que les objets 201, 202 au niveau d'une scène globale. Une troisième étape de traitement 150 utilise le module dit "traitement du contact" et se situe au niveau d'une scène globale (niveau scène 5). Les informations 145 élaborées à l'étape 140 ainsi que les informations 135 élaborées à l'étape 130 sont utilisées pour déterminer pour chaque cas la solution au problème de Signorini qui gouverne le contact entre deux objets dans le cas de glissement pur (information 155). Dans cette étape 150, une modélisation 15 prend ainsi en compte l'interaction entre deux objets tels que les objets 201, 202 au niveau d'une scène globale, les caractéristiques physiques des objets et la description des collisions étant rapatriées pour chaque groupe de collisions par le module "traitement de contact". La troisième étape de traitement 150 fournit des informations 155 concernant des forces et des lieux qui sont transmises au premier module "mécanique" au cours d'une quatrième étape de traitement 160 qui se situe à nouveau au niveau des objets (niveau objet 6). Dans cette étape 160, le résultat du traitement de simulation en temps réel peut être simplement normalisé dans une étape 170 de visualisation ou peut être transmis en retour à travers le couplage 120 vers l'interface utilisateur 110 pour conférer à l'utilisateur un retour de sensation haptique. Dans cette étape finale, une modélisation 16 prend ainsi à nouveau en compte chaque objet ou outil 201, 202, 203 individuellement tout en ayant tenu compte des contacts précédemment simulés. Au niveau de l'objet, dans une forme de réalisation préférée, on décrit l'objet comme un ensemble de triangles reproduisant sa surface et un ensemble de tétraèdres pour décrire son intérieur, le tout dans une configuration de repos. Cette configuration correspond à la forme de l'objet lorsqu'aucune force ne lui est appliquée. Avantageusement, les triangles sont décrits par trois points, placés dans un ordre qui permet au calcul des normales d'être invariablement dirigées vers l'extérieur de l'objet. Les surfaces des objets sont fermées de façon à pouvoir distinguer un extérieur d'un intérieur. Les déformations des objets sont interpolées par la méthode des éléments finis en utilisant un maillage tétraédrique linéaire. Le dispositif permet de simuler différentes lois de comportement sous réserve que l'on puisse en extraire localement, approximativement et pour un pas de calcul, une relation linéaire entre les forces exercées et les déplacements autour d'une configuration locale. Si Uk représente le vecteur de déplacement dans le repère local d'un objet à l'instant courant t et si Uk-i représente le vecteur de déplacement dans le repère local de l'objet au pas de calcul précédent t-1, on a la relation suivante :

Uk= A(Un) Fk + D(Uh-I) (7)

où : A est une matrice permettant de définir la déformation de l'objet au niveau local, Fk est un vecteur représentant les forces extérieures appliquées à l'objet exprimées dans le repère local, b est un vecteur qui a une valeur dans l'espace des déplacements et qui dépend du modèle de déformation de l'objet, et Uk-i est un vecteur dont les valeurs instantanées sont connues au début du pas de calcul, de l'instant courant t.

Avantageusement, une distinction est faite entre le mouvement global de l'objet décrit par la relation fondamentale de la dynamique et la déformation de l'objet autour d'une configuration courante décrite par la loi de déformation. Les forces appliquées sur l'objet se distinguent par leur caractère explicite ou implicite. Une force explicite est d'ores et déjà connue au départ du pas de calcul, et il faut l'intégrer pour connaître le mouvement qu'elle crée sur l'objet. Les forces de contact, au contraire, sont implicites au sens où elles dépendent elles-mêmes du mouvement des objets dans le pas de temps. Pour un pas de temps, on intègre donc les forces explicites et on passe alors au niveau de la scène globale pour pouvoir trouver la valeur des forces implicites. Après l'intégration des forces explicites au niveau de chaque objet, on appelle position et forme "libre" la configuration des objets dans la scène. Cette configuration est obtenue sans l'intervention des forces de contact. Ainsi, un "mouvement libre" est un mouvement intégrant uniquement les forces connues de façon explicite en début de pas de temps de calcul. Par suite on considère ici que le mouvement libre est le mouvement obtenu sur un pas de temps quand les forces de contact ne sont pas intégrées. Le système proposé met alors en œuvre un processus de détection de collision qui permet de tester géométriquement les intersections existantes entre les objets 223, 230 de la scène et les directions privilégiées pour sortir les objets de cette collision (Figures 6 et 7). Si l'on considère un objet 221 qui, dans une position 223, vient en interaction avec un autre objet 230, la direction privilégiée peut être calculée uniquement sur des critères géométriques (cas de la Figure 6) ou peut dépendre de la configuration par laquelle les objets 223, 230 sont entrés en collision en tenant compte d'un mouvement intermédiaire 222 d'au moins l'un des objets entre le pas de calcul précédent et le pas de calcul courant (cas de la Figure 7). Dans tous les cas le processus de détection des collisions permet d'extraire des couples d'éléments d'objets en intersection, une longueur et une direction d'interférence entre ces deux éléments. Dans le cas privilégié d'une description de la surface des objets par des triangles, un élément est soit un point, soit un segment, soit la face d'un triangle. La détection de collision peut prendre en compte trois cas canoniques d'intersection entre deux objets : intersection point/triangle, intersection segment/segment, intersection triangle/point. Le procédé peut également fournir un ensemble d'éléments d'objet en proximité qui pourrait potentiellement entrer en collision suite à l'intégration des forces de contact. Une distance et une direction séparant ces éléments sont alors calculées. Grâce à la description de toutes les interférences et proximités entre les objets, l'ensemble des groupes de collisions de la scène peut être construit. Chaque groupe va alors passer dans le module de contact. Après la détection d'une collision, un traitement dans le module de contact permet de déterminer la configuration du premier contact entre deux pas de temps de calcul. Si l'on prend en compte une description de la surface des objets par des triangles lorsqu'une zone d'interpénétration entre deux objets a été définie à un instant T du pas d'échantillonnage de la simulation, il est extrait une liste de triangles composant la paire d'objets. Si l'on a un objet rigide et un objet déformable, les coordonnées du triangle représentant l'objet déformable sont traduites dans le cadre de référence de l'objet rigide a des instants distincts d'échantillonnage de simulation T et T-I. Pour toute paire possible triangle/triangle il est effectué une interpolation linéaire du déplacement de trois points Di, D2, D3 du triangle déformable entre les positions initiale et finale aux instants discrets d'échantillonnage T et T-I. On peut ensuite opérer trois types de tests différents représentés sur les Figures 3A à 3C sur la paire comprenant un triangle déformable 20 défini par des sommets Di, D2, D3 et présentant diverses positions successives 21, 22, 23 entre les instants T-I et T et un triangle rigide 30 défini par des sommets Ri, R2, R3. Le test 1 (Figure 3A) correspond au cas où un plan de collision est formé par le triangle rigide et va établir une contrainte sur le point concerné du triangle déformable. Le test 2 (Figure 3B) correspond au cas où un plan de collision est formé par un segment rigide et un segment déformable à l'instant de collision t, et va établir une contrainte sur deux points de l'objet déformable. Le test 3 (Figure 3C) correspond au cas où un plan de collision est formé par l'objet déformable à l'instant de collision et va établir une contrainte sur trois points du triangle déformable. L'invention peut s'appliquer aussi bien au cas de collisions entre un objet rigide et un objet déformable qu'au cas de collisions entre deux objets déformables. Le module de contact va être décrit ci-dessous de façon plus détaillée en référence au cas préféré d'une description des objets en triangle. Le module de contact est appelé autant de fois qu'il y a de groupes de collisions dans la scène à l'instant du calcul. Le module de détection de collisions a stocké dans un espace mémoire pour chaque collision : - la normale, - la paire d'objets et les éléments concernés par la collision, - (éventuellement) les points d'application de la force de contact. Si l'algorithme de détection de collision ne donne pas les points d'application de la force de contact (cas illustré sur la Figure 6 d'une détection sans mouvement intermédiaire), il faut les reconstruire et dans tous les cas, il faut interpoler ces points d'applications par rapport au modèle de déformation choisi. Dans le cas préféré d'une description des objets en triangle, on sépare le problème en deux cas : soit on a deux éléments segments, soit on a un élément nœud et un élément face. Dans le cas d'une intersection segment/segment entre un triangle 40 défini par des sommets Pi, P2, P3 et un triangle 50 défini par des sommets Qi, Q2, Q3 (Figure 4), les deux points choisis 41, 51 sont situés à l'intersection de chacun des deux segments Qi, Q2, resp. Pi, P2 avec le plan formé par la face de l'autre triangle en intersection. Le vecteur reliant les deux points trouvés 41, 51 est appelé δ. Dans le cas d'une intersection Point/Face entre un triangle 60 défini par des sommets Pi, P2, P3 et un triangle 70 défini par des sommets Qi, Q2, Q3 (Figure 5), le premier point choisi est le point de l'intersection et le second est la projection de ce point sur la face de l'intersection. Le vecteur reliant les deux points trouvés est appelé δ. On peut utiliser l'algorithme d'intersection avec mouvement intermédiaire proposé ( par Xavier Provot (Collision and Self-collision handling in cloth model dedicated to design garment. Graphics Interface 1997, 177-189) qui permet d'obtenir une configuration approchée entre les deux triangles au moment de la collision. Pour répartir les déplacements et les forces de ces points, dans le cas privilégié d'une interpolation linéaire pour la modélisation en éléments finis (triangles - tétraèdres), on utilise les coordonnées barycentriques. Pour pouvoir calculer correctement les déformations sur les objets, il faut garantir la non-interpénétration. L'interpolation des éléments est donc utilisée de façon à n'avoir qu'une inconnue de force et de distance par contact. Pour simplifier le problème, on se place dans le cas d'une interpolation linéaire pour les éléments finis. On a ainsi dans le cas d'un contact segment/segment (Figure 4) la relation suivante :

δ = [a, b; cj [a 1 — a] I]-^Al (D

où alpha et 1-alpha sont les coordonnées barycentriques sur le premier segment Qi, Q2 et Beta, 1-Beta sur le deuxième segment Pi, P2, ai, bj, Q sont les coordonnées de la normale nt du triangle 40, Wi, W2 sont les coordonnées du premier segment Qi, Q2, Vi, V2 sont les coordonnées du deuxième segment Pi, P2. Dans le cas d'un contact point/plan (Figure 5), ai, bi , q sont les coordonnées de la normale n-, du triangle 60, la distance d'interpénétration δ entre les triangles 60 et 70 s'écrit :

"W1 - δ = [a, b, C1 ] [a β 7] W2 - V1 (2) W3

où : alpha, beta, gamma sont les coordonnées barycentriques sur le triangle 60 (addition = 1). Wi, W2, W3 sont les coordonnées du premier triangle 60, et Vi sont les coordonnées du point de contact 71 constitué par un sommet Qi du second triangle 70. On utilise une interpolation identique pour la force de contact. Une fois trouvé le point d'application des forces de contact, les caractéristiques mécaniques des objets sont transférées dans l'espace défini des contacts. Pour la suite, on suppose que l'on traite un groupe de m contacts avec n objets. Pour le transport des caractéristiques mécaniques globales dans le cas où l'on considère la masse et l'inertie de l'objet de façon globale, en son centre de gravité, on peut utiliser une matrice jacobienne classique, définie dans les travaux de Ruspini (Diego Ruspini & Oussama "A Framework for Multi-Contact Multi-Body Dynamic Simulation and Haptic Display", Proceedings of the 2000 IEEE/RSJ International Conférence on Intelligent Robots and Systems"). Grâce à cette matrice, on a une relation instantanée entre les forces de contact fc dans la direction du contact, les accélérations δ"c (contraintes) et δ'Ηbres dans la même direction :

d\ = J M'J τf> <?', libre (3)

où Jc est une matrice jacobienne m*6n qui transfère le mouvement instantané linéaire et angulaire dans l'espace des contacts, M est une matrice bloc diagonale correspondant à la masse et l'inertie des n objets du groupe de contacts et Jcτ est la matrice transposée de Jc. Pour le transport des caractéristiques locales, on peut utiliser l'interpolation définie pour les triangles. On a ainsi une relation entre les déplacements des points du maillage déformable et les déplacements dans l'espace des contacts et la même interpolation peut être utilisée entre les forces de contact et les forces sur les nœuds. En reprenant la relation linéaire instantanée de caractérisation des déformations, on obtient : δc = [∑i=; N; A (Uj (N;)T ] /c+ δlibre (6) où Nc1 représente la matrice de passage de l'espace des déplacements du maillage vers l'espace des déplacements aux contacts et A est une matrice comme défini plus haut. La modélisation du contact est choisie de façon à respecter au mieux les lois physiques. Par souci de simplification on peut néanmoins de préférence considérer que les contacts ne changeront pas de direction durant la résolution du calcul même si en pratique ce n'est pas strictement le cas. Le premier postulat du problème de Signorini est qu'il n'y a pas d'interférence entre les objets s'ils sont solides (les matières ne se mélangent pas). Ainsi, on souhaite qu'après la résolution du problème le déplacement au contact soit positif ou nul :

<?≥ 0 (9)

Le second postulat est que l'on est dans le cas d'un contact sans frottement, donc la force de contact est dirigée selon la normale : /c≥ O (10)

Le troisième postulat est que la force de contact est non nulle (fc≠ 0) si et seulement si il y a réellement un contact (<5fc =0). Cela crée une relation complémentaire entre les deux vecteurs :

«c J- Zc (H)

Pour pouvoir résoudre le problème de Signorini, on utilise le transport des caractéristiques mécaniques. Les effets des caractéristiques locales et globales sont additionnés en intégrant les accélérations pour se ramener uniquement à une relation entre déplacement et force dans l'espace des contacts. Avantageusement, on utilise une méthode numérique qui a tendance à conserver l'énergie comme la méthode des trapèzes (aussi nommée Tustin), pour intégrer les accélérations et qui peut s'exprimer de la façon suivante si l'on considère une grandeur X aux instants t (Xt) et t+1 (Xt+O :

X1+1 = X1 + Vi dt (x; + X1+1- ) X1+Z= X1^ y2 (It (X/' + X1+1" ) (12)

En utilisant cette méthode numérique et les équations (4) et (5), on obtient la relation suivante :

δe ≈ 0/4 dt2 JC M-1 jcτ+∑i=1n N;C OU (N;)T ] /c+ \bre (13)

Le coefficient 1/4 peut être différent si l'on utilise une autre méthode d'intégration des accélérations. Si le modèle mécanique choisi ne comporte pas de caractéristiques globales et si la masse et l'inertie sont intégrées au niveau local dans le modèle de déformation, on utilise uniquement l'équation (5) qui comporte déjà implicitement une intégration numérique en temps.

δc = [∑i=ιn N; A (u,,) (N;)T ] fc+δlibre (6)

Les postulats définis dans le problème de Signorini et l'équation linéaire instantanée créant une relation linéaire entre forces de contact et déplacements dans l'espace des contacts permettent de formuler le contact sous forme de problème de complémentarité linéaire (LCP). Pour ce type de problèmes, il existe de nombreux algorithmes de résolution (voir par exemple Murty, K.G., Linear Complementarity, Linear and Nonlinear Programming. Internet Edition 1997) qui sont capables de résoudre le problème en un temps compatible avec les performances demandées par l'haptique. A titre d'exemple, on peut effectuer un calcul à une fréquence de l'ordre de 500 Hz à 1000 Hz pour un nombre raisonnable de contacts (par exemple 30 à 40 contacts) en utilisant l'algorithme du Pivot principal sur un ordinateur de type PC Pentium IV 2GHz. La Figure 3 illustre l'interaction entre un objet déformable 80, tel qu'une pince avec un autre objet 90. Dans ce cas, l'objet déformable 80 est virtuellement attaché à l'interface haptique (dès lors qu'il est tenu par l'utilisateur) dans une zone d'un nœud O définissant un repère OxoyoZo. Pour de tels points on peut appliquer les conditions de Dirichlet. A chaque pas de simulation, la configuration de mouvement libre donne un espace de contact. Une résolution LCP donne les forces de contact fi non nulles et un mouvement contraint en est déduit. Ces forces (illustrées par la normale Ui de coordonnées ai bj q sur la Figure 9) sont transportées au point O pour créer la force et le couple sur l'interface haptique. Les Figures 1OA à 1OC montrent l'exemple d'un objet déformable 80 constitué par un clip mis en place sur un tube 91. On voit différentes déformations du clip 80 dans des positions 81, 82, 83 différentes par rapport au tube 91.