FR 2837307 A1 20030919 FR 0203253 A 20020315 Simulateur dynamique d'objets virtuels L'invention se rapporte au domaine de la réalité virtuelle. On y utilise des simulateurs d'objets virtuels rigides, qui peuvent interagir entre eux ou avec leur environnement par contact. Un contact entre objets virtuels est simulé en définissant leur comportement dynamique à travers des données relatives aux propriétés physiques de ces objets : frottements, contraintes de non-interpénétration, par exemple. Les simulateurs connus à ce jour traitent les contacts entre objets virtuels par la résolution de modèles de contact. Les inconnues de ces modèles représentent des forces de contrainte, pour un contact permanent, ou des impulsions, pour un contact instantané. Des méthodes de simulation de ce genre sont développées dans les ouvrages - "Fast contact force computation for non penetrating rigid bodies", D. BARAFF, in SIGGRAPH 94 Conference Proceedings, Annual Conference Series, pp 23-34, ACM SIGGRAPH, Addison Wesley, 1994, ou à la référence électronique correspondante www.graphics.cornell.edu/pubs/1992/Bar92a.htlm. - "Collision/Contact Models for the Dynamic Simulation of Complex Environments", D. Ruspini, 0. Khatib, IEEE/RSJ International Conference on Intelligence Robots and Systems s:IROS'97, September 1997, Grenoble, France. Ces méthodes présentent certains inconvénients, en particulier des inconvénients algorithmiques, comme on le verra plus loin. La présente invention vient améliorer la situation. L'invention concerne un dispositif informatique de modélisation destiné à travailler sur des données d'entrée, comprenant d'une part des données d'environnement, d'autre part des données d'objet représentant une pluralité d'objets définis par des grandeur cinématiques, le dispositif comprenant un module de calcul, apte à mettre en oeuvre un jeu d'algorithmes à chaque pas de temps, à partir de valeurs précédentes des données d'objet et des données d'environnement, le jeu d'algorithmes correspondant à la résolution d'un jeu d'équations pour un contact donné entre au moins deux objets. Selon une caractéristique principale de l'invention, le module de calcul est apte, à chaque pas de temps, et pour au moins deux objets en contact, à mettre en oeuvre la résolution d'une première équation représentant la minimisation d'une norme euclidienne appliquée à une première expression prédéterminée ayant une inconnue dépendante d'une grandeur dynamique de chaque objet, cette première expression prédéterminée correspondant à un comportement des objets hors contraintes, et ladite résolution de la première équation s'effectuant sous contrainte d'une deuxième équation représentant au moins une condition de contrainte sur l'inconnue de la première équation et en ce que le jeu d'algorithmes correspondant à la résolution du jeu d'équations comprend un algorithme de projection d'un point sur un cône polyédrique. L'invention concerne également un procédé de modélisation informatique, à partir de données d'entrée, comprenant d'une part des données d'environnement, d'autre part des données d'objet représentant une pluralité d'objets définis par des grandeurs cinématiques, le procédé mettant en oeuvre itérativement un jeu d'algorithmes correspondant à la résolution d'un jeu d'équations pour un contact donné entre au moins deux objets, et comprenant les étapes suivantes: a. élaborer une première équation représentant la minimisation d'une norme euclidienne appliquée à une première expression prédéterminée ayant une inconnue dépendante d'une grandeur dynamique de chaque objet, cette première expression prédéterminée correspondant à un comportement des objets hors contraintes, et une deuxième équation représentant au moins une condition de contrainte sur l'inconnue de la première équation. b. mettre en oeuvre ledit jeu d'algorithmes, comprenant un algorithme de projection d'un point sur un cône polyédrique, pour la résolution de ladite première équation, ladite résolution de la première équation s'effectuant sous contrainte de ladite deuxième équation L'invention concerne également le produit logiciel, définissant les éléments utilisés dans le dispositif informatique de modélisation tel que défini ci-dessus. L'invention concerne par ailleurs le produit logiciel, comprenant les fonctions utilisés dans le procédé de modélisation informatique tel que défini ci-dessus. L'invention concerne en outre un logiciel de modélisation, comprenant le produit logiciel tel que défini ci-dessus. D'autres caractéristiques et avantages de l'invention apparaîtront à l'examen de la description détaillée ci-après, ainsi que des dessins annexés sur lesquels: - la figure 1 illustre schématiquement un simulateur selon la technique antérieure, - la figure 2 illustre schématiquement un simulateur selon l'invention, - la figure 3 est la représentation de deux objets en un contact continu, - la figure 4 est la représentation de deux objets en un contact instantané, - la figure 5 est un ordinogramme représentant les étapes de la simulation selon l'invention. De plus, la description détaillée est accompagnée d'annexes, dans lesquelles : - l'annexe A présente un exemple du code utilisé selon l'invention, - l'annexe B présente la suite de l'exemple du code, - l'annexe 1 détaille les matrices générales utilisées, - l'annexe 2 détaille les formules de contraintes selon l'invention, - l'annexe 3 détaille les formules de modélisation selon l'invention. Les dessins contiennent, pour l'essentiel, des éléments de caractère certain. Ils pourront donc non seulement servir à mieux faire comprendre la description, mais aussi contribuer à la définition de l'invention, le cas échéant. La plupart des simulateurs connus à ce jour traitent les contacts entre objets virtuels par la formulation d'un modèle de contact classique de type "Problème à Complémentarité Linéaire" (traduction de l'anglais Linéarité Complementarity Problem), appelé encore problème LCP, puis par la résolution de ce modèle de contact en utilisant des algorithmes appropriés assez complexes. Les inconnues de ces modèles représentent des forces de contrainte, pour un contact permanent, ou des impulsions, pour un contact instantané. De plus, le nombre des objets de la simulation et le nombre des contraintes appliquées à ces objets sont couplés. Cette formulation présentent des inconvénients, en particulier au niveau de la complexité algorithmique. La vitesse de résolution de cette formulation ainsi que la précision des solutions trouvées en sont affectées. La formulation du modèle de contact classique est reprise pour en extraire une expression algorithmique particulière, résolue par des algorithmes simples adéquats. Cette expression algorithmique présente des avantages importants, notamment pour la simulation en temps réel, comme l'amélioration du temps de calcul des solutions et de la précision des solutions. Dans la description qui suit, on emploiera le terme "cinématique" pour désigner ce qui relève du mouvement et le terme "dynamique" pour désigner ce qui relève du mouvement et des causes de ce mouvement (par exemple les forces). Dans la description qui suit, les vecteurs sont formellement représentés en majuscule et en gras. Les vecteurs, matrices et opérations utilisées sont indiqués en annexe 1. Le point J1 est considéré comme le centre de gravité de l'objet i. Le point H est un point de contact entre deux objets i et j, et est appelé Hi et Hj lorsqu'il est respectivement rattaché à l'objet i et à l'objet j. Ainsi sont définis le vecteur inconnu X, le vecteur accélération A et le vecteur vitesse V en annexe 1-I. Dans chaque vecteur X, A et V représentant N objets en contact, chaque objet i est représenté par une première composante, respectivement Xi(Ji), ai(Ji) et vi(Ji), des vecteurs de translation, et par une deuxième composante, respectivement yi, ai et wi, des vecteurs de rotation. Chaque composante est un vecteur de dimension égale à 3. Si la simulation porte sur N objets en contact, les vecteurs X, A et V sont de dimension 6N. Avantageusement, les vecteurs utilisés dans la description sont formés de manière semblable: - un vecteur accélération sans contrainte Au - un vecteur accélération contrainte Ac - un vecteur désignant un vecteur accélération contrainte inconnu Xa - un vecteur vitesse sans contrainte Vu - un vecteur vitesse contrainte Vc - un vecteur désignant un vecteur vitesse contrainte inconnu Xc - un vecteur vitesse avant collision V- - un vecteur vitesse après collision V+ - un vecteur désignant un vecteur vitesse après collision inconnu Xv - un vecteur généralisé inconnu Xg - un vecteur généralisé sans contrainte Xk - des vecteurs solutions chacun d'une équation Xg*, Xa*, Xv*, Xc*. Les vecteurs Xk; Xg et Xg* peuvent désigner chacun un vecteur accélération ou un vecteur vitesse adéquat. En annexe 1-II, une matrice M de masse généralisée pour N objets en contact est définie par une masse mi de l'objet i, une matrice identité Id de dimension 3x3 et la matrice d'inertie Ii de l'objet i, de dimension 3x3. Ainsi, M est une matrice de dimension 6Nx6N, formée de N sous-matrices Mi. Les opérations sur les matrices sont rappelées en annexe 1III. La figure 1 représente un simulateur d'objets virtuels selon la technique antérieure. Le simulateur 1 est muni d'une base de données 10, qui peut comprendre: - des données sur les objets, notamment position, vitesse et masse des objets virtuels en contact, et - des données sur leur environnement, notamment forces extérieures appliquées aux objets virtuels en contact. Pour simuler le mouvement des objets, le simulateur travaille à chaque pas de temps pour déterminer certaines données sur les objets, par exemple leur position. Un contact continu est détecté en fonction notamment des positions et des vitesses des objets. Le simulateur utilise un module 12 de modélisation de contact continu selon des forces de contrainte. En cas de contact continu entre objets virtuels, le module 12 modélise le mouvement des objets virtuels et utilise, comme inconnues, des forces de contrainte. En effet, le mouvement des objets virtuels en contact continu est modifié par des forces de contrainte. Un détecteur de collision 14 reconnaît l'état de collision entre objets virtuels à un instant donné. Le détecteur de collision peut utiliser une première technique de détection de collision dite continue ou une deuxième technique dite discrète. En cas de collision détectée entre objets virtuels, le détecteur de collision 14 déclenche le module 13 de modélisation de contact instantané selon des impulsions. Ce module 13 modélise le mouvement des objets virtuels en intégrant les données sur les objets et les données sur leur environnement juste avant la collision, et en utilisant, comme inconnues, des impulsions. Un calculateur 16 est muni d'algorithmes de calcul spécifiques pour chaque module de modélisation 12 et 13. Ces algorithmes calculent les forces de contrainte inconnues dans le cas de la modélisation de contact continu et les impulsions inconnues dans le cas de la modélisation de contact instantané. Sur la figure 2 est présenté un simulateur 2 selon l'invention. Selon l'invention, une modélisation pour tout type de contact entre objets virtuels est effectuée selon un principe, dit "principe de GAUSS" ou "principe des moindres contraintes de GAUSS". Ce principe énonce que pour des objets rigides sur lesquels s'exercent des contraintes géométriques, une accélération contrainte particulière Xa*, sur l'ensemble des accélérations possibles, minimise la fonction Ga de l'annexe 3-I-al). Ce principe sera développé au cours de la suite de la description et défini de manière généralisée à une inconnue X* pouvant être une accélération ou une vitesse dans l'annexe 3-III-al). Chaque modélisation est basée sur des contraintes appliquées au mouvement des objets virtuels. La modélisation est considérée comme "consistante" si les positions et les vitesses des objets sont compatibles avec ces contraintes. Le simulateur 2 comprend un module 20. Ce module 20 met à disposition une base de données sur les objets et de données sur leur environnement. On suppose chaque objet virtuel i simulé comme rigide, c'est à dire comme non déformable. Ainsi, les données sur les objets du module 20 définissent pour chaque objet virtuel i une masse mi, une position Pi, une vitesse Vi et des contraintes géométriques que doit satisfaire chaque objet i. Ces contraintes géométriques, également appelées contraintes de non-pénétration développées en annexe 2-1, peuvent être bilatérales, par exemple pour des liaisons permanentes, comme des joints ou des articulations, ou unilatérales, par exemple pour un contact continu entre deux objets. Ainsi, ces données sur les objets définissent l'état du système à chaque instant. Quant aux données sur l'environnement du module 20, elles correspondent aux actions extérieures appliquées sur chaque objet, comme par exemple la gravité, le vent, des actions définies sur les objets virtuels par l'utilisateur du simulateur. Des modules de modélisation sont définis sur la figure 2 : - un module 211 de modélisation du contact continu, - un module 212 de modélisation du contact instantané. Le module 210 de modélisation regroupe les modules 211 et 212 pour définir des modélisations selon le principe de GAUSS. De façon connue, un détecteur de collision 24 déclenche le module 212 qui définit alors la modélisation du contact instantané. Ce détecteur de collision 24 utilise avantageusement une technique de détection de collision continue qui permet d'interpoler les positions successives des objets et de calculer le moment auquel s'effectue une collision entre certains objets. Cette technique de détection de collision continue est présentée dans le document suivant : " An Algebraic Solution to the problem of Collision Detection for Rigid Polyhedral Objects", S.Redon, A. Khedddar and S. Coquillart, Proceding of International Conference of Robotics and Automation, pp 3733-3738, April 2000. Un calculateur 26 permet de chercher la valeur minimale d'un vecteur inconnu X dans une modélisation proposée par le module 210 de modélisation selon le principe de GAUSS. Dans un premier environnement, dit monde du premier ordre, les objets ne sont soumis ni à une accélération ni à une inertie, contrairement à un deuxième environnement, dit monde du second ordre. Les données du module 20 définissent l'ordre du monde dans lequel la simulation est effectuée. Pour les modélisations de contact, les contraintes définies par le module 210 sont détaillées en annexe 2. Dans le cas d'un contact continu et d'un monde du second ordre, le module 211 définit des contraintes géométriques, dites contraintes de non pénétration, appliquées aux accélérations. Une contrainte de non pénétration sur les accélérations de deux objets i et j en contact continu au point H en annexe 2-I-al) est obtenue par dérivation d'une contrainte de non pénétration sur la position de l'objet. Dans l'exemple, cette formule s'exprime au centre de gravité des objets i et j, avec comme inconnues les accélérations contraintes possibles de rotation et de translation. Un système de m contraintes unilatérales appliquées aux N objets en contact permet de définir les accélérations contraintes possibles Xa par la formule en annexe 2-I-a2). Ca est définie comme une matrice des contraintes géométriques sur les accélérations, de dimension mx6N. Ba est défini comme un vecteur des conditions appliquées aux contraintes géométriques sur les accélérations. Dans l'exemple proposé, une contrainte bilatérale est exprimée sous forme de deux contraintes unilatérales. Dans le cas d'un contact continu et d'un monde du premier ordre, le module 211 définit les contraintes géométriques appliquées aux vitesses des objets. Une contrainte de non pénétration pour un contact continu entre deux objets i et j s'exprime par la formule en annexe 2-I-bl). Cette formule s'exprime avantageusement au centre de gravité des objets i et j, avec comme inconnues les vitesses contraintes possibles de rotation et de translation. Pour m points de contact, un système de m contraintes unilatérales appliquées aux N objets en contact permet de définir les vitesses contraintes possibles Xv par la formule en annexe 2I-b2). Cv est définie comme une matrice des contraintes géométriques sur les vitesses, de dimension mx6N. Dans le cas d'un contact instantané entre objets virtuels, les contraintes définies par le module 212 s'appliquent aux vitesses. Un point de contact H entre deux objets i et j est défini comme Hi sur l'objet i et Hj sur l'objet j . A un instant donné et pour un point de contact H entre deux objets i et j, l'hypothèse de Poisson définit une contrainte sur les vitesses des objets comme indiquée dans la formule de l'annexe 2-II-1). Dans cette formule, e représente le coefficient de restitution qui est choisi en fonction des matériaux des objets virtuels en contact. Cette formule peut s'exprimer au centre de gravité des objets, avec comme inconnues les vitesses de rotation et de translation. Pour m points de contact, un système de m contraintes permet de définir les vitesses après collision possibles Xcpar la formule II-2). Cc est définie comme une matrice des contrain- tes sur les vitesses lors de collision(s), de dimension mx6N. Bc est défini comme un vecteur des conditions appliquées aux contraintes sur les vitesses lors de collision(s). Le module 210 détermine une modélisation selon le principe des moindres contraintes de GAUSS. Ainsi, ce principe énonce une minimisation de la fonction scalaire Ga basée sur l'accélération non contrainte Au décrite en annexe 3-I-al). Selon l'invention, le module 210 est apte à définir un système d'équations, décrit en annexe 3-I-a2), pour minimiser la fonction scalaire Ga en fonction des accélérations possibles déterminées suivant des contraintes sur les accélérations. Le calculateur 26 résout ce système d'équations. Avantageusement, dans le cas d'un monde du premier ordre, le module 211 définit un système d'équations donné dans l'annexe 3-II-b2) dont l'accélération sans contrainte est remplacée par la vitesse sans contrainte Vu. La fonction Ga est remplacée par la fonction Gv définie dans l'annexe 3-II-bl) . Le calculateur 26 cherche, parmi les vitesses possibles, la vitesse contrainte la plus proche de la vitesse sans contrainte, en minimisant la fonction Gv. Dans le cas d'un contact instantané, le module 212 définit un système d'équations donné dans l'annexe 3-II-1). Pour N objets en contact instantané, Au=O traduit l'hypothèse d'une intensité négligeable des actions extérieures par rapport à l'intensité des actions des contraintes durant une collision. L'accélération est considérée comme constante durant la collision. Aussi, la fonction Ga est remplacée par la fonction Gc basée sur la vitesse avant collision V-. Le calculateur 26 cherche, parmi les vitesses possibles après collision, la vitesse après collision qui minimise la fonction Gc. Selon un mode de réalisation, le module 210 est agencé pour déterminer une modélisation généralisée selon le principe des moindres contraintes de GAUSS. Le module 210 est apte à définir un système d'équations généralisées décrit en annexe 3-III-a2), pour minimiser la fonction scalaire Gg en annexe 3III-al) généralisée à un vecteur sans contrainte Xk, en fonction de contraintes généralisées à la matrice des contraintes Ck et au vecteur des conditions appliquées aux contraintes Bk. Ainsi, le système d'équation généralisé comprend une équation primaire prédéterminée (3-III-a2, Eql) représentant la minimisation d'une expression primaire, notamment la fonction scalaire Gg, ayant pour inconnue une grandeur cinématique (Xg, 3-III-a2) de chaque objet, et une équation secondaire (3-IIIa2, Eq2) apte à appliquer à la grandeur cinématique (Xg, 3III-a2) de chaque objet des contraintes spécifiques (Ck, 3III-a2, Eq2) respectant des conditions déterminées (Bk, 3-IIIa2, Eq2). Ce système d'équations définit un contact quelconque entre objets virtuels afin que le mouvement de ces objets soit calculé par le calculateur 26. Plus précisément, un système d'équations généralisées énonce que parmi les inconnues Xg possibles, définies par les équations de m contraintes, l'inconnue Xg* minimise la fonction Gg généralisée. Assimilée à une distance cinétique, la racine carré de la fonction Gg est également minimisée par l'inconnue contrainte. Selon un aspect particulièrement intéressant, le module 210 définit une modélisation généralisée simplifiée selon le principe des moindres contraintes de Gauss. A partir de la modélisation généralisée de l'annexe 3-III-a2), le système d'équations se réduit à la formule des équations Eql et Eq2 de l'annexe 3-III-bl) en effectuant un changement de variables. Comme présenté dans cette annexe, la matrice M est définie positive et se décompose selon les matrices D et F, permettant ainsi d'obtenir les nouvelles équations formulées selon la variable X1. C1 est une matrice de dimension mx6N. En procédant à un changement de base sur ce système, le système d'équations généralisées simplifié peut s'écrire comme indiqué en annexe 3-III-b2). X2 et E sont définis comme des vecteurs de dimension 6N+1, C2 est une matrice de dimension mx(6N+1). D'autre part, la solution X2* est supposée sous une forme particulière en fonction d'un vecteur Q et d'une composante h telle que définie en annexe 3-III-b2). Le vecteur Q est défini comme un vecteur de dimension 6N. La forme du vecteur X1* en est déduite. Selon l'invention, la méthode de Lagrange permet, à partir du système d'équation généralisé simplifié 3-III-b2), d'obtenir le système d'équation selon Lagrange 3-III-b3). Ainsi, l'invention permet une formulation simple du problème de contact entre objets. Ce système d'équation selon Lagrange facilite les calculs du calculateur en proposant la minimisation d'une norme euclidienne avec comme inconnue la variable de Lagrange # vérifiant une condition de positivité. La variable de Lagrange représente une variable dynamique (par exemple une force de contact ou une impulsion). Ce système d'équation selon Lagrange comprend une première équation (3-III-b3, Eql) représentant la minimisation d'une norme euclidienne appliquée à une première expression prédéterminée ayant une inconnue dépendante d'une grandeur dynamique de chaque objet (notamment cinématique), cette première expression prédéterminée correspondant à un comportement des objets hors contraintes, et ladite résolution de la première équation s'effectuant sous contrainte d'une deuxième équation (3-III-b3, Eq2) représentant au moins une condition de contrainte sur l'inconnue de la première équation. Le calculateur 26 résout le système selon Lagrange qui représente l'une ou l'autre des modélisations des modules 211 et 212. Le problème de la résolution de ce système d'équations selon Lagrange consiste en un problème de projection géométrique du vecteur E sur un cône polyédrique de dimension 6N+1 dans l'exemple. Selon l'invention, le calculateur utilise un algorithme simple pour la résolution de ce système d'équations. Cet algorithme est un algorithme de projection d'un point sur un cône polyédrique. Divers algorithme connus de l'art antérieur peuvent être avantageusement utilisés et sont décrits dans les ouvrages suivants : - K. G Murty and Y. Fathi, "A critical Index Algorithm for Nearest Point Problems on Simplicial Cônes", Mathematical Programming, 23 (1982) 206-215; - P. Wolfe, "Algorithm for a Least Distance Programming Problem", Mathematical Programming Study 1, (1974) 190-205; - P.Wolfe, "Finding the Nearest Point in a Polytope" Mathematical Programming 11 (1976) 128-149; - D. R Whilhelmsenn "A Nearest Point Algorithm for Convex Polyhedral Cones and Applications to Positive Linear Approximations", Mathematics of Computation, 30 (1976) 48-57. Dans un mode de réalisation de l'invention, l'algorithme de projection d'un point sur un cône polyédrique utilisé est l'algorithme de Wilhelmsen par exemple. Un algorithme de projection d'un point sur un cône polyédrique permet avantageusement de trouver une solution aux modélisations considérées comme consistantes. En annexe 3-III-b4) est présentée la formulation de la solution des vecteurs X2* et Xg*. En variante de réalisation, à partir de la modélisation généralisée de l'annexe 3-III-a2), le système d'équations est réduit directement au système d'équations Eql et Eq2 de l'annexe 3-III-b5). Pour cela, une fonction de Lagrange L(Xg,#) de l'annexe 3-III-b6) associée à la modélisation de l'annexe 3-III-a2) est utilisée et la variable Xg est éliminée après avoir posé égale à zéro la dérivée de la fonction L selon la variable Xg. Dans l'annexe 3-III-b5), les matrices M, D et F sont identiques aux matrices définies en annexe 3-111bl). Cl est la même matrice de dimension mx6N définie comme étant fonction des matrices Ck, D et F. On choisit un vecteur K tel qu'il réponde à l'équation définie en annexe 3-III-b5). Un vecteur S de dimension 6N est défini selon une équation de la même annexe. Le nouveau système d'équations Eql, Eq2 peut s'écrire comme indiqué dans cette annexe 3-III-b5) correspondant à un système de Lagrange. Ce système d'équation ayant la même forme que le système d'équation 3-III-b3), le problème de sa résolution consiste en un problème de projection géométrique du vecteur S sur un cône polyédrique de dimension 6N. La résolution s'effectue ainsi en utilisant le même type d'algorithme que précédemment, en particulier un algorithme de projection d'un point sur un cône polyédrique, par exemple l'algorithme de Wilhelmsen. Une partie du code selon l'invention est présenté dans les annexes A et B. Dans l'annexe A, la gestion des contraintes est présentée suivant trois fonctions : - ConstrainGeneralVelocity() quipermet le calcul des vitesses contraintes, - ConstrainGeneralAcceleration() quipermet le calcul des accélérations contraintes, -Resolvelmpact() qui permet le calcul des vitesses de rebond, lors d'une collision. Dans l'annexe B, l'algorithme de Wilhelmsen est présentée suivant trois autres fonctions : - cFilter() qui permet la transformation du problème en un problème de la forme de Wilhelmsen, et appelle la fonction cProjection(), - cProjection() qui définit l'algorithme de Wilhelmsen, - cLinearCholesky() qui résout une équation matricielle A.X=b, et qui est appelée par l'algorithme de Wilhelmsen. Selon l'invention, le module 210 définit un groupe de contacts représentant N objets mobiles formant une chaîne de contacts entre le premier objet et le énième objet. Ainsi, des groupes de contacts sont dynamiquement indépendants les uns des autres et peuvent faire l'objet de simulations indépendantes. Les figures 3 à 4 représentent des objets virtuels en contact afin d'illustrer les modélisations de contact selon l'invention. La figure 3 représente un premier objet en contact avec un deuxième objet Q, ce deuxième objet Q étant immobile. Grâce aux données sur l'environnement des objets, par exemple la gravité, le module de modélisation de contact continu en relation avec le calculateur détermine l'accélération sans contrainte Au . Cette accélération sans contrainte Au représente l'accélération qu'aurait l'objet s'il n'était pas sujet aux contraintes de non-pénétration. Dans l'exemple de la figure 3, l'accélération sans contrainte Au de l'objet virtuel R correspond à la gravité de ce même objet. Selon la contrainte de non pénétration indiquée en annexe 2-I-al) et appliquée à l'objet virtuel R, les accélérations possibles de ce même objet sont orientées suivant la pente LA de l'objet virtuel Q. Parmi ces accélérations possibles, l'accélération contrainte Ac de l'objet virtuel R correspond à l'accélération possible la plus proche de l'accélération sans contrainte Au. L'accélération contrainte Ac, telle que définie par le principe de GAUSS des moindres contraintes en annexe 3-I-a2), est orientée sur la figure 3. La figure 4 propose un exemple de collision entre deux objets virtuels R et Q. Dans l'exemple, l'objet Q est immobile. Ainsi, l'objet R est animé d'une vitesse V- étant la vitesse avant collision. A l'instant Tcol l'objet virtuel R entre en contact avec l'objet virtuel Q en un point H. Des hypothèses sont faites concernant le contact instantané. Ainsi, une première hypothèse considère la durée de la collision comme infinitésimale. En d'autres termes, la position de l'objet R reste inchangée, ainsi que les accélérations de contrainte. Une autre hypothèse considère que les données sur l'environnement des objets ou sur les actions extérieures, peuvent être négligées durant la collision par rapport à l'intensité des actions de contrainte. Comme indiqué en annexe 3-II-2), le principe des moindres contraintes de GAUSS est aménagé de sorte que la vitesse après collision soit l'inconnue. Ainsi, la vitesse après collision V+, parmi les vitesses possibles après collision, minimise la fonction Gc, développée en annexe 3-II-1). Ces vitesses possibles après collision sont déterminées par la contrainte de réponse de collision selon l'hypothèse de Poisson détaillée en annexe 2-II. Ainsi, dans un groupe de contacts, les vitesses possibles après collision sont les vitesses les plus proches des vitesses avant collision. Sur la figure 4, l'ensemble des vitesses possibles après collision est représenté par le demi-espace situé au-dessus de la ligne LC. La vitesse après collision V+ est, selon le principe des moindres contraintes de GAUSS, la plus proche possible de la vitesse avant collision V-. Dans un monde du premier ordre et dans le cas d'un contact continu, le principe des moindres contraintes de GAUSS s'attache à définir la vitesse de contrainte d'objets d'un groupe de contacts. Ce groupe de contact est soumis à une contrainte de non-pénétration appliquée aux vitesses des objets et indiquée en annexe 2-I-bl). Cette fonction linéaire permet de déterminer les vitesses possibles des objets du groupe de contacts. Enfin, dans la fonction Ga indiquée dans l'annexe 3-I-al), l'accélération sans contrainte est remplacée par la vitesse sans contrainte Vu définissant une nouvelle fonction Gv dans l'annexe 3-I-bl). Parmi les vitesses possibles, la vitesse contrainte du groupe de contacts, est choisie pour minimiser la fonction Gv. En d'autres termes, à tout moment, la vitesse contrainte du groupe de contact est choisie pour correspondre à la vitesse possible la plus proche de la vitesse sans contrainte. Selon l'invention, le principe des moindres contraintes de GAUSS permet au calculateur de trouver le mouvement des objets virtuels en contact. Avantageusement, l'invention est apte à définir, dans tous les cas de contact, un ensemble de solutions pour les accélérations et les vitesses possibles déterminées grâce aux contraintes appliquées aux objets virtuels. Dans cet ensemble de solutions, une solution unique peut être déterminée. Selon un autre avantage de l'invention, chaque modèle de contact peut être simplifié en deux étapes. Ces deux étapes sont un changement de variable et une homogénéisation du modèle considéré. Ainsi, les modèles de contact sont réduits à un modèle unique décrit en annexe 3-111. Avantageusement, l'inconnue de ce modèle unique est trouvée après projection géométrique dans un espace de dimension déterminée, 6N ou 6N+1 où N est le nombre d'objets d'un groupe de contacts. Au niveau du calculateur, la résolution du problème par projection géométrique est effectuée par un algorithme de projection d'un point sur un cône polyédrique connu de l'état de la technique. Avantageusement, cet algorithme simple peut être l'algorithme de Wilhelmsen. La figure 5 propose une méthode de simulation d'objets virtuels rigides selon l'invention. A un instant T, les données sur les objets définissent les N objets du groupe de contact. Ces données sont notamment la position P, la vitesse V, l'accélération généralisée A dans le cas d'un monde du second ordre (appelé aussi mode A), la position P, la vitesse V dans le cas d'un monde du premier ordre (appelé aussi mode V). Ainsi, on définit les objets virtuels dans leur environnement. Si des contraintes sont applicables sur un objet ou un groupe de contacts à l'étape 102, une modélisation du contact continu est effectuée selon l'invention. Afin de pouvoir traiter numériquement cette modélisation, une méthode de résolution d'équations différentielles est appliquée à la modélisation considérée à l'étape 104. Avantageusement, cette méthode de résolution d'équations différentielles peut être la méthode de Runge-Kutta. Sur un intervalle de temps dT, cette méthode de résolution applique quatre fois de suite l'algorithme de Wilhelmsen. L'application de l'étape 104 permet de déterminer à l'étape 108: - l'accélération contrainte Ac cherchée si la simulation s'effectue dans un monde du second ordre; - la vitesse contrainte Vc cherchée si la simulation s'effectue dans un monde du premier ordre. Avantageusement, l'intervalle de temps peut être variable. A l'étape 102, si aucune contrainte n'est applicable sur un objet ou sur un groupe de contacts, à l'étape 110 l'accélération généralisée égale l'accélération non contrainte Au si la simulation s'effectue dans un monde du second ordre et la vitesse généralisée égale la vitesse non contrainte Vu si la simulation s'effectue dans un monde du premier ordre. Si aucune collision n'est détectée à l'étape 114, les nouvelles vitesse V1 et position P1 sont calculées pour un temps T égal à T+dT. Si au contraire une collision est détectée au temps Tcol égal à T+a.dT avec a inférieur à 1, à l'étape 116: - les nouvelles vitesse V- et position Pi sont calculées avant la collision si la simulation s'effectue dans un monde du second ordre; - le nouvelle position Pi est calculée avant la collision si la simulation s'effectue dans un monde du premier ordre. A l'étape 117, le procédé retourne à l'étape 100 si la simulation s'effectue dans un monde du premier ordre. Dans le cas d'une simulation dans un monde du second ordre, le contact instantané est modélisé selon l'invention à l'étape 118. A l'étape 120, l'algorithme de Wilhelmsen est appliqué à la modélisation du contact instantané. Cet algorithme permet de trouver une solution unique pour la vitesse après collision V+ à l'étape 124. En retournant à l'étape 100, le mouvement des objets virtuels est calculé pas de temps après pas de temps. Dans un exemple de réalisation de l'invention, une technique de détection de collision continue est appliquée aux objets. Cette technique attribue aux objets des déplacements rigides fixés pour interpoler les positions successives de ces objets et ainsi détecter les collisions éventuelles et les instants d'impact entre ces objets. Ainsi, la simulation selon l'invention apporte entre autre une amélioration quant aux nombreux désavantages des méthodes de simulation antérieures. L'invention est applicable dans de nombreux domaines industriels et permet, par exemple, de présenter des produits virtuels en simulation dynamique, de simuler en temps réel un assemblage virtuel de pièces sur des chaînes de montage/ démontage, de simuler des tests de fonctionnement. L'invention peut également s'appliquer dans le domaine des jeux virtuels. Bien entendu, l'invention n'est pas limitée à la forme de réalisation décrite précédemment à titre d'exemple, elle s'étend à d'autres variantes. Avantageusement, d'autres algorithmes peuvent résoudre le problème de projection concernant la modélisation du contact continu. Annexe A Gestion descontraintes A-l- ConstrainGeneralVelocityO void cScene::ConstrainGeneralVelocity(pcContactGroup cg, double *intmd, double *conmd) { Cette fonction seplacedans l'espace des vitesses pourcontraindre lemouvement intentionneldesobjetsafin qu'ils respectent lescontraintessurleurmouvement. Elle estappeléeplusieursfois, au coursd'un même pas // de temps, par la fonction de cRGVelocity. NB : fonction ne met pas a jour la liste des contacts, ni leur Il leurnombre. nobjectsest lenombre totald'objetsdans legroupe decontacts (mobilesetimmobiles). nmobile est le nombre d'objets mobiles dans legroupe de contacts. ncontacts est lenombre de contacts. Contact est // un tableau depointeursversdesclassescContact, de taillencontacts. invidestuntableau d'entiersdetaille // nobjects, quidonne la correspondance entre l'ID de l'objet etsaplace dans la liste d'objetsmobiles (quand l'objet est immobile, sa place est -1). mobp est un tableau d'entiers de taille nmobile, qui donne la //correspondance inverse : lenuméro ID de l'objet mobiledontlemouvement aet filtre. intm estun tableaude //doubledetaille 6*nmobile représentantle mouvement intentionneldesobjetsmobiles (enradians). conm est //un tableau de double de taille 6*nmobile représentant le mouvement filtre final des objets mobiles (en //radians). Il NB :dans cettefonction, ongère aussi la distance de sécurité inti ; //variables relativesau groupe de contact int ncontacts=cg->GetNContactsO; pcContact *contact=new pcContact[ncontacts]; for (i=O;i ContactPointer(i); //A-1-a Première étape : on transforme les contacts en contraintes sur le mouvement initialisation du tableau, on suppose qu'ilne peut pasy en avoir plus de50 // NB: on suppose que la liste de contacts (contact) est ajour (on afait appela lafonction Update) int nee=0; // nombrede contraintes ee intnvf=0; nombre de contraintes vf intnfv=0; nombre de contraintesfv int spacedim=6*cg->GetNObjects(); double *cmat=new double[(spacedim+l)*ncontacts]; double *intmh=new double[spacedim+1]; //mouvement intentionnelgénéraliséhomogénéisé double *conmh=new double[spacedim+l]; //mouvement contraintgénéralisehomogénéisé for (i=O;i for (i=O;i switch (contaet[i]->GetTypeO) { case VF~CONSTRAINT: nvf++; break ; case EE~CONSTRAINT: nee++ ; break ; case FV~CONSTRAINT: nfv++ ; break ; default : break ; } // il ya deux objets impliques dans un contact // Soitdeux objetsmobiles, soit un objet mobile et un objetdudecor // on a une contrainte sur lemouvement relatif cMatrix4gml(contact[i]->GetO10->GlobalMatrix0); cVertexpivotgtobl=gml*contact[i]->Get01()->GetPivotO; cVertex gcl(contact[i]->GetLocationn-pivotglobi); cVertex n(contact[i]->GetNormalo); calcul descoefficients de la contrainte // n: vecteurnormal auplan de contrainte dans le repèreglobal // NB: n est suppose unitaire // coefficientsde la contrainte cVertex cprod(gcl^n); int idim=i*(spacedim+l); // position de la contraintedans le tableau int indldim=cg->InvId(contact[i]->Get0l->GetIdn)*6; int offl=idim+indldim; for (int j=O;j // coefficients relatifs au deuxième objet int ind2dim=cg->InvId(contact[i]->Get02->GetId); if(ind2dim>-l) { // ledeuxième objet estmobile aussi ind2dim*=6 ; int off2=idim+ind2dim; cmat[off2]=n[O]; cmat[off2+l]=n[l]; cmat[off2+2]=n[2]; cMatrix4 gm2(contact[i]->Get020->GlobalMatrixO); cVertex pivotglob2=gm2*contact[i]->Get02.->GetPivot; cVertex gc2( contact[i] ->GetLocationO-pivotglob2); cVertex cprod2(gc2^n); cmat[off2+3]=cprod2[0]; cmat[off2+4]=cprod2[l]; cmat[off2+5]=cprod2[2]; } cmat[idim+spacedim ]=DISPLACEMENT +contact[ i]-> PenetrationO; // cmat[idim+spacedim]=2.0*DISPLACEMENT; // NB:pour les contactde type face/point, on considère lemouvement relatifoppose : if (#ntact[i]->GetTypeO==FV~CONSTRAINT) { for (j=O;j A-2- ConstrainGeneralAccelerationp void cScene::ConstrainGeneralAcceleration(pcContactGroup cg, double *intmd, double *conmd) { Il Cette fonctionse placedansl'espacedesaccélérations pourcontraindrelemouvement intentionneldesobjets ll afin qu'ils respectentlescontraintessurleurmouvement. Elleestappeléeplusieurs fois, au coursd'un même Il pasdetemps,par la fonction decRGAcceleration. NB : cette fonction nemet pasa jour la listedescontacts, // ni leur nombre. nobjects est le nombre total d'objets dans le groupede contacts (mobiles et immobiles). Il nmobile est le nombre d'objets mobiles dans le groupe de contacts. ncontacts est le nombre de contacts. Il contact estun tableau de pointeursversdes classescContact, de taillencontacts. invidestun tableau d'entiers //detaillenobjects, quidonne la correspondance entre l'ID de l'objet etsaplace dans la listed'objetsmobiles //(quand l'objet est immobile, saplace est -1). mobp est un tableau d'entiers de taille nmobile, qui donne la //correspondance inverse : lenuméro ID de l'objetmobile dont lemouvement a et filtre. vdintestun tableau de double de taille 6*nmobile représentantlemouvement intentionneldesobjetsmobiles (en radians). vdfiltest // un tableau de double de taille 6*nmobile représentant le mouvement filtrefinal des objets mobiles (en Il radians). Lafriction estgerée en rajoutant des contraintes a l'ensemble //decontraintesunilatérales. Donc même sionconnaîtaudépartlenombre de pointsdecontacts, onneconnaît Il pas encore lenombre de contraintes // NB:dans cettefonction, ongère aussi la distancede sécurité inti ; //variables relativesau groupe decontact int n#ntacts=cg->GetNContactsO; pcContact *contact=new pcContact[ncontacts]; for (i=O;i ContactPointer(i); // A- 2-a.. Premiere étape : on transforme lescontacts en contraintes sur le mouvement // initialisation du tableau, on supposequ'ilnepeutpasy en avoirplus de 50 // NB:on suppose que la listede contacts (contact) est ajour (on afait appela lafonction Update) int nee=0; // nombre de contraintes ee int nvf=0; nombrede contraintes vf int nfv=0; // nombrede contraintesfv int nconstraints; //nombrede contraintes aufinal int cnumber; //numéro de la contrainte unilatérale en coursde construction //denombremement du nombre finalde contraintes nconstraints=ncontacts ; création desvariablesprimaires (avantpassage a la distance euclidienne) ll ces variablessonthomogénéisées int spacedim=6*cg->GetNObjects(); //nombre de degrésde libertédans legroupe de contact double *cmat=new double[(spacedim+l)*nconstraints]; double *intm=new double[spacedim+l]; double *conm=new double[spacedim+l]; for (i=O;i for (i=O;i Error()){ cout "Le contact " i " estun contact arête/arête incorrect(arêtes parallèles). " endl; break ; dénombrement des contraintes switch (contact[i]->GetTypeO) { case VF~CONSTRAINT: nvf++; break ; caseEE~CONSTRAINT: nee++ ; break ; caseFV~CONSTRAINT: nfv++ ; break ; default: break ; } // A-2-b Création de la contrainte de non-interpénétration cnumber++ ; // ilya deux objets impliquesdans un contact Soit deux objetsmobiles, soit un objetmobile et un objetdu decor lion a une contraintesur le mouvement relatif pcObject ol,o2; cMatrix4 gml,gm2; cVertex pivotglobl,pivotglob2; cVertex gcl,gc2; cVertex vcl,vc2; cVertex n ; cVertex cprod2 ; ol=contact[i]->GetO10; o2=contact[i]->Get02; // calculdes coefficientsde la contrainte // n: vecteurnormal auplan de contraintedans lerepèreglobal // NB: n estsuppose unitaire // NB: les coefficientsdevantdxdy dz wxwy wz (qui représentent ici les accelerations) sont // les mêmes quepour les vitesses. Seul le terme constantchange. gm1=(o1->GlobalMatrix()); pivotglobl=(gml *ol ->GetPivot()); gel =( contact[i]->GetLocationO-pivotglob 1); n=(contact[i]->GetNormaIO); // coefficients de la contrainte cVertex cprod(gcl^n); int idim=cnumber*(spacedim+l);// position de la contrainte dans le tableau int indldim=cg->InvId(ol->GetldO)*6; int offl=idim+indldim; for (int j=O;j int ind2dim=cg->Invld(o2->Getldo); if(ind2dim>-l) { // le deuxième objetestmobile aussi ind2dim*=6 ; int off2=idim+ind2dim; cmat[off2]=n[0]; cmat[off2+l]=n[l]; cmat[off2+2]=n[2]; gm2=(o2->GlobalMatrix()); pivotglob2=(gm2 *o2->GetPivotO); gc2=( contact[ i]->GetLocationO-pivotglob2); cprod2=(gc2"n); cmat[off2+3]=cprod2[0]; cmat[off2+4]=cprod2[l]; cmat[off2+5]=cprod2[2]; // calculdu terme constant // NB: ce terme contientl'accélération centripete (w^(w^GC))#n, pour les deuxobjets. // Ilcontient aussi2*(va-vb)lndot, ou np est la derivee par rapport au tempsdu vecteurnormal. Il Enfin, il contient un terme fonction de l'ecart a la distance de securite (DISPLACEMENT+contact[i]->Penetration()), //quiagit comme un ressort dans les methodes depenalite. cVertex wl,w2; //acceleration centripete de l'objet 1 wl=(ol->GetRotVelocityO); double dis=-((w1^(w1^gcl))n); if(ind2dim>-l) { // le deuxieme objetestmobile //acceleration centripetede l'objet2 w2=( o2->GetRotV elocityO); dis+= ((w2^(w2^gc2))n); //deplacementdu vecteurnormal cVertex ndot(contact[i]->GetNormalDoto); // dérivée par rapport au temps du vecteur normal vcl=(ol->GetTransVelocityO+(wl "gel ; // vitesse du point de contact dans l'objet 1 vc2=(o2->GetTransVelocityO+(w2Agc2)); // vitesse du point de contact dans l'objet 2 dis-=2. 0*((vc1-vc2)#ndot); } maintien de la distance desecurite double delta=DISPLACEMENT+contact(i]->Penetration0; dis+=(KDISPLACEMENT* delta); cmat[idim+spacedim]=dis; // NB:pour les contactde typeface/point, on considere le mouvement relatifoppose : if (contact(i]->GetType==FV CONSTRAINT) { for (j=O;j cFilter(spacedim+l,nconstraints,cmat,intm,conm); // appel a Wilhelmsen for (i=O;i } A-3. ResolveImpact() void cScene::ResolveImpact(pcContactGroup cg) { Cette fonction se placedans l'espace des vitesses pour calculer le mouvement d'objets entrant en collision. Implicitement, c'est donc un calcul simultan# d'impulsions. Explicitement, on calcule directement // le mouvement des objets. Cette fonction esttressemblable auxautres fonctionsdegestions decontraintesNB //: cette fonctionnemet pasa jourla listedescontacts, nileurleurnombre. On suppose par consequent que // l'on peutfaire appel en toute securite a lafonction Colliding des contacts, pour savoir si les objets entrent // en collision en cepointde contactou non. Il nobjectsestlenombre totald'objetsdans legroupe decontacts(mobilesetimmobiles). nmobile estlenombre d'objets mobiles dans legroupe de contacts. ncontacts est le nombre decontacts. Il contactestuntableaude pointeursversdesclassescContact, de taillencontacts.mlistestuntableaud'entiers // de taillenobjects, quidonne la correspondance entre l'ID de l'objet etsaplace dans la listed'objetsmobiles Il(quandl'objetestimmobile, sa placeest-1). mobp estun tableau de pcObjectdetaillenmobile. mobp[i] donne // un pointeur vers l'objet dont la place est j dans la liste courante des objets mobiles. NB : utile de gerer // aussi la distance de securite dans cette fonction ? int int nmobile=cg->GetNObjects(); pcObject *mobp=new pcObject[nmobile]; //mobile objectpointer for (i=O;i Id(i)]; int ncontacts=cg->GetNContacts(); pcContact *contact=new pcContact[ncontacts]; for (i=O;i ContactPointer(i); // A-3-a. Premiere étape : on tranforme les contacts en contraintes sur le mouvement // NB:on suppose que la liste de contactsest a jour (on a fait appela lafonction Update) int nee=0; // nombre de contraintesee int nvf=0; // nombrede contraintes vf int nfv=0; nombre de contraintesfv int spacedim=6*nmobile; double *cmat=new double[(spacedim+l)*ncontacts]; double *intm=new double[spacedim+l]; double *conm=new double[spacedim+l]; // on stocke les vitesses desobjetsdans un vecteurcommun // on lefaitgraceau tableau des correspondances mobp for (i=O;i cVertex tspeed( mobp[i]->GetTrans V elocityO); cVertex rspeed(mobp[i]->GetRotVelocity); int i6=i*6; intm [i6]=tspeed[0]; intm[i6+l]=tspeed[l]; intm [i6+2]=tspeed[2]; intm[i6+3]=rspeed[0]; intm [i6+4]=rspeed[1]; intm[i6+5]=rspeed[2]; intm[spacedim]=1.0; //calculdu tableau des contraintes for (i=O;i GetType()) { caseVF~CONSTRAINT: nvf++ ; break ; caseEE~CONSTRAINT: nee++ ; break ; caseFV~CONSTRAINT: nfv++ ; break ; default : break ; Il ilya deux objets impliques dans un contact // Soitdeux objets mobiles, soit un objetmobile etun objetdu decor // on a une contrainte sur lemouvement relatif cMatrix4gml(contact[i]->GetO10->GlobalMatrix0); cVertexpivotgIobl=gml*contact[i]->Get01()->GetPivotO; cVertex gcl(contact(i]->GetLocationQ-pivotglobl); cVertex n(contact[i]->GetNormalo); // calculdes coefficientsde la contrainte // n: vecteurnormal auplan de contrainte dans le repere global // NB :nestsuppose unitaire coefficients dela contrainte cVertex cprod(gc1^n); int idim=i*(spacedim+l); // positionde la contrainte dans le tableau intindldim=cg->InvId(contact[i]->GetO10->GetId0)*6; int offl=idim+indldim; for (intj=O;j int ind2dim=cg->InvId(contact[i]->Get02n->GetId); if(ind2dim>-1) { // ledeuxieme objetest mobile aussi ind2dim*=6 ; int off2=idim+ind2dim; cmat[off2]=n[0]; cmat[off2+l]=n[l]; cmat[off2+2]=n[2]; cMatrix4gm2(contact[i]->GetO20->GlobalMatrix0); cvertex pivotglob2=gm2*contact[i]->GetO2o->GetPivoto; cVertex gc2( contact[ i]->GetLocationO-pivotglob2); cVertex cprod2(gc2^n); cmat[off2+3]=cprod2[0]; cmat[off2+4]=cprod2[l]; cmat[off2+5]=cprod2[2]; } // calculdu coefficient constant double e=0.34; //coefficient de restitution (bonne valeur : 0.38) double nspeed=-e*(contact[i]->RelativeVelocityoln); // vitesse normale //relative aupointdecontact if(nspeed>DISPLACEMENT) cmat[idim+spacedim]=nspeed; // elsecmat[idim+spacedim]=O; else ecmat[idim+spacedim ]=0.1 *KDISPLACEMENT* (DISPLACEMENT +contact[i]-> PenetrationO); // NB:pour les contactde type face/point, on considere le mouvement relatifoppose : if(contact[i]->GetTypeO==FV~CONSTRAINT) { for (j=O;j //A-3-b. Deuxieme étape : on contraint lemouvement cFilter(spacedim+l,ncontacts,cmat,intm,conm); //appela Wilhelmsen //modificationdes vitessesdesobjets for (i=O;i SetTransVelocity(tspeed); mobp[i]->SetRotVelocity(rspeed); } delete [] cmat; delete [] intm ; delete [] conm; } Annexe B Algorithme de Wilhelmsen B-l.cProjectionO void cProjection(int n, int m, double c[], double v[], double s[]) { Ilcette fonction utilisel'algorithmede Wilhelmsen pourtrouverle pointsle plus prochedupointvdans lecone 1/ polyédriquepositifengendre par m vecteurs. Ces vecteurssontstockes dans la matrice cdetaille n*m Tous les vecteurssontdans unespace dedimension n. NB :la jemecoordonnee du vecteurciesten c[i*n+j]Ainsi, quand on lit les coordonnees dans l'ordre, on nefait qu'une multiplication. On suppose deplus que les// vecteurs ci sont tous de norme 1. on commence par mettre a CNORME la norme de tous les vecteurs // contraintes, poureviter les erreursd'arrondi. NB :il faudrait tenircompte dufaitque l'onprojette lepoint // v, qui est tresparticulier (des0partoutpresque) int i,j; int mp; // cardinaldel'ensembleF int mm=m*m; //initialisationdes variables for (i=O;i if ((scal)&&(scal≥-P~EPSILON)&&(scal≤P~EPSILON)) { // cout "risque de bugpour scal dans cProjection !!!!!!!!!!!!!" #endl; // scal=0; } innerv[i]=scal; // on stocke le résultat if ((scal>P EPSILOI&&(scal>scalmax)) { found=l; indice=i ; scalmax=scal; // NB:par rapport aux anciens, lescoefficients innerv sontmultipliespar CNORME if(!found){ // l'origine estlepoint leplusproche de v for (i=O;i // on initialise les variables double rhomin; double comprhomin; for (i=O;i indice=0; found=0; double scalmax=0; for (i=O;i Il calculduproduit scalaire d.ci int in=i*n; double scal=0; for (j=O;j if ((scal>P EPSILON)&&(scal>scalmax)) { found=1; indice=i ; scalmax=scal; } } } if(!found) { // le projeté courant est le pointsolution for (i=O;i // on initialise les variables // création de l'ensembleF 1[0]=indice; // on met le vecteurchoisi dans la liste lambda[indice]=0; mp=1 ; for (i=O;i 1 [mp]=i; mp++ ; } } Il création des coordonnées barycentriques de R // au débutR=pj // NB: tous les lambda sontmultiplies par CNORME for (i=O;i do { //débutde la boucle interne if(mp==0){ for (i=O;i n) coût "pb : sous-ensemble de contraintes non indépendant (mp>n)" endl; for (i=0;i if (!(innertag[1[j]*m+I[i]j)) int in=l [i]*n; int jn=l [j]*n; double scat=0; for (intk=O;k innertag[1(j]*m+1[i]]=1; I/ on note le fait que le calcul a deja ete fait } submat(j *mp+i]=inner[1[j] * m+1[i]); submat[i*mp+j]=submat[j *mp+i]; } // NB:ne pas oublier les coefficients de la diagonale for (i=O;i cLinearCholesky(mp,submat,subvec,proj); //il faut tenir compte de la multiplicationpar CNORME for (i=0;i for (i=0;i break ; } //le projetéorthogonal n'est pas dans le nouveau sous-cone lion cherche alors un nouveau point R sur la frontière du cône //etsurlesegment quirelie le point R ancien au projeté orthogonal //determination des rho[i] for (i=O;i if ((diff)&&(diff>-P EPSILON)&&(diff≤P EPSILON)) { } if(diff>0) rho[i]=lambdar[i]/diff; else rho [i]=1.0; } //determination du minimum des rho[i] rhomin=rho[0]; for (i=l;i if gam[i])&&(gam[i]>-P ~EPSILON)&&(gam[i]≤P - EPSILON {} if(gam[i]>P~EPSILON) { Il pouréviter les erreurs d'arrondi ltemp [mptemp]=l[i]; lambdartemp[mptemp]=gam[i]; mptemp++; } } mp=mptemp; for (i=O;i for (i=O;i } Ilpuison resout Ux=y for (i=(n-1);i≥O;i--) sum=x[i]; for (k=(i+l);k Il NB:lajeme coordonnee du vecteurci est enc[i*n+j] //Ainsi, quandon lit les coordonnees dans l'ordre, on nefaitqu'une multiplication int i,j; //Avant de resoudre le probleme dual, il faut faire un changementde variables // NB:pour le moment, on suppose que la matrice des masses est l'identite double *c2=new double[n*m]; double *en=new double[n]; // vecteur (0,0,0,...,0,1) de taille n double *x2=new double[n]; // solution duprobleme dual for (i=0;i for (i=O;i // calcul dex x (x=(1lh).P.(llsqrt(D))A+Xu) double invh=1.0/x2[n-l]; for (i=0;i Annexe 1 Formules générales : Annexe 2 Formules de contraintes : I. dan al) (a(Ht)-a,(H,))+2x(v,(Hr)-v,(H,))dt0 a2) {CaxXa#Ba} bl) (v1(H1)-v1(H1)).n#0 b2) {CvxXv#0} II. 1) (v; (H,)-v 1 (H, .n ex (v- 1 (H,) - v- 1 (H, .n 2) {CcxXc#Bc} Annexe 3 Formules de modélisation : Revendications 1. Dispositif informatique de modélisation destiné à travailler sur des données d'entrée (20), comprenant d'une part des données d'environnement, d'autre part des données d'objet représentant une pluralité d'objets définis par des grandeur cinématiques, le dispositif comprenant un module de calcul (26), apte à mettre en oeuvre un jeu d'algorithmes à chaque pas de temps, à partir de valeurs précédentes des données d'objet et des données d'environnement, le jeu d'algorithmes correspondant à la résolution d'un jeu d'équations pour un contact donné entre au moins deux objets, caractérisé en ce que le module de calcul est apte, à chaque pas de temps, et pour au moins deux objets (R, Q) en contact, à mettre en oeuvre la résolution d'une première équation (3III-b3, Eql) représentant la minimisation d'une norme euclidienne appliquée à une première expression prédéterminée ayant une inconnue dépendante d'une grandeur dynamique de chaque objet, cette première expression prédéterminée correspondant à un comportement des objets hors contraintes, et ladite résolution de la première équation s'effectuant sous contrainte d'une deuxième équation (3-III-b3, Eq2) représentant au moins une condition de contrainte sur l'inconnue de la première équation et en ce que le jeu d'algorithmes correspondant à la résolution du jeu d'équations comprend un algorithme de projection d'un point sur un cône polyédrique. 2. Dispositif selon la revendication 1, caractérisé en ce que ladite première équation (3-III-b3, Eql) représentant la minimisation d'une norme euclidienne est tirée d'une équation primaire prédéterminée (3-III-a2, Eql) représentant la minimisation d'une expression primaire ayant pour inconnue une grandeur cinématique (Xg, 3-III-a2) de chaque objet. 3. Dispositif selon l'une des revendications 1 et 2, caractérisé en ce que ladite deuxième équation (3-III-b3, Eq2) est tirée d'une équation secondaire (3-III-a2, Eq2) apte à appliquer à la grandeur cinématique (Xg, 3-III-a2) de chaque objet des contraintes spécifiques (Ck, 3-III-a2, Eq2) respectant des conditions déterminées (Bk, 3-III-a2, Eq2). 4. Dispositif selon la revendication 2, caractérisé en ce que ladite équation primaire prédéterminée (3-III-a2, Eql) est apte à astreindre la grandeur cinématique (Xg, 3-III-a2) de chaque objet à se rapprocher le plus possible de la grandeur cinématique particulière (Xk, 3-III-a2, Eql) qu'aurait chaque objet s'il n'était pas soumis à des contraintes. 5. Dispositif selon l'une des revendications 3 et 4, caractérisé en ce que ladite grandeur cinématique comprend l' accélération contrainte cherchée (Xa, 3-I-a2), ladite grandeur cinématique particulière comprend une accélération sans contrainte (Au, 3-I-a2, Eql), lesdites contraintes spécifiques comprennent des premières contraintes géométriques (Ca, 3-Ia2, Eq2). 6. Dispositif selon l'une des revendications 3 et 4, caractérisé en ce que ladite grandeur cinématique comprend la vitesse contrainte cherchée (Vv, 3-I-b2), ladite grandeur cinématique particulière comprend la vitesse sans contrainte (Vu, 3-I-b2, Eql), lesdites contraintes spécifiques comprennent des deuxièmes contraintes géométriques (Cv, 3-I-b2, Eq2). 7. Dispositif selon l'une des revendications 3 et 4, caractérisé en ce que ladite grandeur cinématique comprend la vitesse après collision cherchée (Xe, 3-II-2) , ladite grandeur cinématique particulière comprend la vitesse avant collision (V, 3-II-2, Eql), lesdites contraintes spécifiques comprennent les contraintes sur la vitesse après collision cherchée (Cv, 3-II-2, Eq2). 8. Dispositif selon l'une des revendications 1 à 6, caractérisé en ce que ledit jeu d'algorithmes correspondant à la résolution d'un jeu d'équations, constitué de la première équation et la deuxième équation, comprend un algorithme de résolution d'équations différentielles propre à être appliqué en correspondance avec l'algorithme de projection d'un point sur un cône polyédrique (104). 9. Dispositif selon la revendication 8, caractérisé en ce que ledit algorithme de résolution d'équations différentielles est l'algorithme de Runge-Kutta. 10. Dispositif selon l'une des revendications 1 à 8, caractérisé en ce que ledit algorithme de projection d'un point sur un cône polyédrique est l'algorithme de Wilhelmsen. 11. Dispositif selon l'une des revendications 1 à 10, caractérisé en ce que ladite première équation représente une équation matricielle décrivant une pluralité d'équations linéaires. 12. Procédé de modélisation informatique, à partir de données d'entrée (20), comprenant d'une part des données d'environnement, d'autre part des données d'objet représentant une pluralité d'objets définis par des grandeurs cinématiques, le procédé mettant en oeuvre itérativement un jeu d'algorithmes correspondant à la résolution d'un jeu d'équations pour un contact donné entre au moins deux objets, et comprenant les étapes suivantes: a. élaborer une première équation (3-III-b3, Eql) représentant la minimisation d'une norme euclidienne appliquée à une première expression prédéterminée ayant une inconnue dépendante d'une grandeur dynamique de chaque objet, cette première expression prédéterminée correspondant à un comportement des objets hors contraintes, et une deuxième équation (3-III-b3, Eq2) représentant au moins une condition de contrainte sur l'inconnue de la première équation. b. mettre en oeuvre ledit jeu d'algorithmes, comprenant un algorithme de projection d'un point sur un cône polyédrique, pour la résolution de ladite première équation (3-III-b3, Eql), ladite résolution de la première équation s'effectuant sous contrainte de ladite deuxième équation (3-III-b3, Eq2). 13. Procédé de modélisation informatique selon la revendication 12, caractérisé en ce que l'étape a. comprend le fait de tirer ladite première équation (3-III-b3, Eql) d'une équation primaire prédéterminée (3-III-a2, Eql) représentant la minimisation d'une expression primaire ayant pour inconnue une grandeur cinématique (Xg, 3-III-a2) de chaque objet. 14. Procédé de modélisation informatique selon l'une des revendications 12 et 13, caractérisé l'étape a. comprend le fait de tirer ladite deuxième équation (3-III-b3, Eq2) d'une équation secondaire (3-III-a2, Eq2) apte à appliquer à la grandeur cinématique (Xg, 3-III-a2) de chaque objet des contraintes spécifiques (Ck, 3-III-a2, Eq2) respectant des conditions déterminées (Bk, 3-III-a2, Eq2). 15. Procédé de modélisation informatique selon la revendication 13, caractérisé en ce que l'équation primaire prédéterminée (3-III-a2, Eql) de l'étape a. comprend d'astreindre la grandeur cinématique (Xg, 3-III-a2) de chaque objet à se rapprocher le plus possible de la grandeur cinématique particulière (Xk, 3-III-a2, Eql) qu'aurait chaque objet s'il n'était pas soumis à des contraintes. 16. Procédé de modélisation informatique selon l'une des revendications 14t 15caractérisé en ce que ladite grandeur cinématique de l'étape a. comprend l'accélération contrainte cherchée (Xa, 3-I-a2), ladite grandeur cinématique particulière comprend une accélération sans contrainte (Au, 3-I-a2, Eql), lesdites contraintes spécifiques comprennent des premières contraintes géométriques (Ca, 3-I-a2, Eq2). 17. Procédé de modélisation informatique selon les revendications 14 et 15, caractérisé en ce que ladite grandeur cinématique de l'étape a. comprend la vitesse contrainte cherchée (Vv, 3-I-b2), ladite grandeur cinématique particulière comprend la vitesse sans contrainte (Vu, 3-I-b2, Eql), lesdites contraintes spécifiques comprennent des deuxièmes contraintes géométriques (Cv, 3-I-b2, Eq2). 18. Procédé de modélisation informatique selon les revendica- tions 14 et 15, caractérisé en ce que ladite grandeur cinéma- tique de l'étape a. comprend la vitesse après collision cherchée (Xc, 3-11-2), ladite grandeur cinématique particu- lière comprend la vitesse avant collision (V-, 3-11-2, Eql), lesdites contraintes spécifiques comprennent les contraintes sur la vitesse après collision cherchée (Cv' 3-11-2, Eq2). 19. Procédé de modélisation informatique selon l'une des revendications 12 à 17, caractérisé en ce que l'étape b. comprend de mettre en oeuvre un algorithme de résolution d'équations différentielles propre à être appliqué en corres- pondante avec l'algorithme de projection d'un point sur un cône polyédrique (104). 20. Procédé de modélisation informatique selon la revendica- tion 19, caractérisé en ce que ledit algorithme de résolution d'équations différentielles est l'algorithme de Runge-Kutta. 21. Procédé de modélisation informatique selon l'une des revendications 12 à 20, caractérisé en ce que ledit algorithme de projection d'un point sur un cône polyédrique est l'algorithme de Wilhelmsen. 22. Produit logiciel, comprenant le module de calcul du dispositif informatique de modélisation revendiqué dans l'une des revendications 1 à 11. 23. Produit logiciel, comprenant les fonctions utilisées dans le procédé de modélisation informatique revendiqué dans l'une des revendications 12 à 21. 24. Un logiciel de modélisation, comprenant le produit logiciel revendiqué dans l'une des revendications 22 et 23.