FR2837307A1 - Dynamic simulator for virtual reality objects comprises a computer modeling device that resolves kinematic equations relating to a number of objects using a first unresolved equation and a second constraint representing equation - Google Patents

Abstract

Computer modeling device for modeling the kinematic behavior of a number of objects comprises a calculator for implementation of a set of algorithms at each time step to resolve a contact between at least two objects. The calculation module resolves a first expression corresponding to the behavior of objects without constraints. The first equation is resolved using a second equation representing at least one constraint condition for the unknown quantity in the first equation. The invention also relates to a corresponding computer modeling method.

Description

<Desc / Clms Page number 1>

The invention relates to the field of virtual reality.

Rigid virtual object simulators are used, which can interact with each other or with their contact environment.

A contact between virtual objects is simulated by defining their dynamic behavior through data relating to the physical properties of these objects: friction, non-interpenetration constraints, for example.

The simulators known to date treat the contacts between virtual objects by the resolution of contact models. The unknowns of these models represent forces of stress, for a permanent contact, or pulses, for an instantaneous contact. Simulation methods of this kind are developed in the books - "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, or the corresponding electronic reference www.graphics.cornell.edu/pubs/1992/Bar92a.html.

- "Collision / Contact Models for the Dynamic Simulation of Complex Environments", D. Ruspini, 0. Khatib, IEEE / RSJ International Conference on Intelligence Robots and Systems: IROS'97, September 1997, Grenoble, France.

These methods have certain disadvantages, in particular algorithmic disadvantages, as will be seen below.

The present invention improves the situation.

A computer modeling device for working on input data, comprising on the one hand environment data, on the other hand object data representing a plurality of defined objects

<Desc / Clms Page number 2>

 by kinematic magnitude, the device comprising a calculation module, able to implement a set of algorithms at each time step, from previous values of the object data and environment data, the game of algorithms corresponding to solving a set of equations for a given contact between at least two objects.

According to a main characteristic of the invention, the calculation module is able, at each time step, and for at least two objects in contact, to implement the resolution of a first equation representing the minimization of a Euclidean norm. applied to a first predetermined expression having an unknown dependent on a dynamic magnitude of each object, this first predetermined expression corresponding to a behavior of the objects out of constraints, and said resolution of the first equation being under stress of a second equation representing at least one constraint condition on the unknown of the first equation and in that the set of algorithms corresponding to the resolution of the set of equations comprises an algorithm for projecting a point on a polyhedral cone.

The invention also relates to a computer modeling method, based on input data, comprising on the one hand environment data, on the other hand object data representing a plurality of objects defined by kinematic quantities. , the method implementing iteratively a set of algorithms corresponding to the resolution of a set of equations for a given contact between at least two objects, and comprising the following steps: a. to develop a first equation representing the minimization of a Euclidean norm applied to a first predetermined expression having an unknown dependent on a dynamic quantity of each object, this first predetermined expression corresponding to a behavior of the non-constrained objects, and a second equation representing less

<Desc / Clms Page number 3>

 a constraint condition on the unknown of the first equation. b. implementing said set of algorithms, comprising an algorithm for projecting a point on a polyhedral cone, for solving said first equation, said resolution of the first equation being under stress of said second equation The invention concerns also the software product, defining the elements used in the computer modeling device as defined above.

The invention also relates to the software product, comprising the functions used in the computer modeling method as defined above.

The invention further relates to a modeling software, comprising the software product as defined above.

Other characteristics and advantages of the invention will become apparent on examining the detailed description below, as well as the appended drawings in which: FIG. 1 diagrammatically illustrates a simulator according to the prior art, FIG. 2 schematically illustrates a simulator according to the invention, - Figure 3 is the representation of two objects in a continuous contact, - Figure 4 is the representation of two objects in an instantaneous contact, - Figure 5 is a flow chart showing the steps of the simulation. according to the invention.

<Desc / Clms Page number 4>

In addition, the detailed description is accompanied by annexes, in which: - Appendix A presents an example of the code used according to the invention, - Appendix B shows the continuation of the example of the code, - the appendix. 1 details the general matrices used, - Appendix 2 details the stress formulas according to the invention, - Appendix 3 details the modeling formulas according to the invention.

The drawings contain, for the most part, elements of a certain character. They can therefore not only serve to better understand the description, but also contribute to the definition of the invention, if any.

Most of the simulators known to date treat the contacts between virtual objects by the formulation of a classical contact model of the "Linear Complementarity Problem" type, called the LCP problem, then by the resolution of this contact model using appropriate algorithms quite complex. The unknowns of these models represent forces of stress, for a permanent contact, or pulses, for an instantaneous contact. In addition, the number of objects in the simulation and the number of constraints applied to these objects are coupled.

This formulation has drawbacks, in particular at the level of algorithmic complexity. The speed of resolution of this formulation as well as the accuracy of the solutions found are affected.

The formulation of the classical contact model is taken up to extract a particular algorithmic expression, solved by appropriate simple algorithms. This algorithmic expression has important advantages, especially for real-time simulation, such as the improvement of

<Desc / Clms Page number 5>

 calculation time of the solutions and the precision of the solutions.

In the following description, the term "kinematics" will be used to refer to movement and the term "dynamic" to refer to movement and the causes of movement (eg forces).

In the following description, the vectors are formally represented in capital letters and in bold. The vectors, matrices and operations used are shown in Appendix 1.

Point J1 is considered the center of gravity of object i. Point H is a point of contact between two objects i and j, and is called Hi and Hj when it is respectively attached to object i and object j.

Thus are defined the unknown vector X, the acceleration vector A and the velocity vector V in appendix 1-I. In each vector X, A and V representing N objects in contact, each object i is represented by a first component, respectively Xi (Ji), ai (Ji) and vi (Ji), translation vectors, and by a second component , respectively yi, ai and wi, rotation vectors. Each component is a vector of dimension equal to 3. If the simulation relates to N objects in contact, the vectors X, A and V are of dimension 6N.

Advantageously, the vectors used in the description are similarly formed: an acceleration vector without constraint Au a constrained acceleration vector Ac a vector designating an unknown stress acceleration vector Xa a constrained speed vector Vu a constrained speed vector Vc a vector designating an unknown constraint speed vector Xc

<Desc / Clms Page number 6>

 - a velocity vector before collision V- - a velocity vector after collision V + - a vector designating a velocity vector after unknown collision Xv - an unknown generalized vector Xg - a generalized unconstrained vector Xk - solution vectors each of an equation Xg * , Xa *, Xv *, Xc *.

Xk vectors; Xg and Xg * may each designate an acceleration vector or a suitable velocity vector. In Appendix 1-II, a generalized mass matrix M for N objects in contact is defined by a mass mi of the object i, an identity matrix Id of dimension 3 × 3 and the inertia matrix Ii of the object i, of 3x3 dimension. Thus, M is a matrix of dimension 6Nx6N, formed of N sub-matrices Mi.

The operations on the matrices are recalled in Appendix 1III.

Figure 1 shows a simulator of virtual objects according to the prior art. The simulator 1 is provided with a database 10, which can comprise: data on the objects, in particular position, speed and mass of the virtual objects in contact, and data on their environment, in particular external forces applied to the objects. virtual contacts.

To simulate the movement of objects, the simulator works at each time step to determine certain data on the objects, for example their position. Continuous contact is detected depending in particular on the positions and speeds of the objects.

The simulator uses a module 12 for continuous contact modeling according to stress forces. In the case of continuous contact between virtual objects, the module 12 models the

<Desc / Clms Page number 7>

 movement of virtual objects and uses, as unknowns, forces of constraint. Indeed, the movement of virtual objects in continuous contact is modified by forces of constraint.

A collision detector 14 recognizes the state of collision between virtual objects at a given instant. The collision detector may use a first so-called continuous collision detection technique or a second so-called discrete technique. In the event of a collision detected between virtual objects, the collision detector 14 triggers the instantaneous contact modeling module 13 according to pulses. This module 13 models the movement of the virtual objects by integrating the data on the objects and the data on their environment just before the collision, and using, as unknowns, pulses.

A computer 16 is provided with specific calculation algorithms for each modeling module 12 and 13. These algorithms calculate the unknown stress forces in the case of continuous contact modeling and the unknown pulses in the case of instantaneous contact modeling. .

In Figure 2 is presented a simulator 2 according to the invention.

According to the invention, modeling for any type of contact between virtual objects is carried out according to a principle, called "GAUSS principle" or "principle of the least GAUSS constraints". This principle states that for rigid objects on which geometric constraints are exerted, a particular stress acceleration Xa *, over all the possible accelerations, minimizes the Ga function of Annex 3-I-al). This principle will be developed during the rest of the description and defined in a generalized way to an unknown X * that can be an acceleration or a speed in Annex 3-III-a1).

<Desc / Clms Page number 8>

Each modeling is based on constraints applied to the movement of virtual objects. Modeling is considered "consistent" if object positions and velocities are consistent with these constraints.

The simulator 2 comprises a module 20. This module 20 provides a database on the objects and data on their environment. Each simulated virtual object is assumed to be rigid, that is, to be non-deformable. Thus, the data on the objects of the module 20 define for each virtual object i a mass mi, a position Pi, a speed Vi and geometric constraints that must satisfy each object i. These geometric constraints, also called non-penetration constraints developed in Appendix 2-1, can be bilateral, for example for permanent connections, such as joints or joints, or unilateral, for example for a continuous contact between two objects. Thus, these object data define the state of the system at every moment.

As for the environmental data of the module 20, they correspond to the external actions applied to each object, such as gravity, wind, actions defined on the virtual objects by the user of the simulator.

Modeling modules are defined in FIG. 2: a modeling module 211 for modeling the continuous contact; a module 212 for modeling the instantaneous contact.

The modeling module 210 groups the modules 211 and 212 to define modelizations according to the GAUSS principle.

In known manner, a collision detector 24 triggers the module 212 which then defines the modeling of the instantaneous contact. This collision detector 24 advantageously uses a continuous collision detection technique that makes it possible to interpolate the successive positions of the objects and to calculate the moment at which a collision occurs between

<Desc / Clms Page number 9>

 some objects. This technique of continuous collision detection is presented in the following document: "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.

A calculator 26 makes it possible to search for the minimum value of an unknown vector X in a model proposed by the modeling module 210 according to the GAUSS principle.

In a first environment, said world of the first order, objects are not subject to acceleration or inertia, unlike a second environment, said world of the second order.

The data of the module 20 defines the order of the world in which the simulation is performed. For the modelizations of contact, the constraints defined by the module 210 are detailed in appendix 2.

In the case of continuous contact and a world of the second order, the module 211 defines geometric constraints, called non-penetration constraints, applied to accelerations. A non-penetration stress on the accelerations of two objects i and j in continuous contact at the point H in Appendix 2-I-a1) is obtained by derivation of a non-penetration stress on the position of the object. In the example, this formula is expressed in the center of gravity of the objects i and j, with as unknown acceleration possible constraints of rotation and translation.

A system of m unilateral constraints applied to the N objects in contact makes it possible to define the possible stress accelerations Xa by the formula in appendix 2-I-a2). It is defined as a matrix of geometric constraints on accelerations, of dimension mx6N. Ba is defined as a

<Desc / Clms Page number 10>

 vector of conditions applied to geometric constraints on accelerations.

In the proposed example, a bilateral constraint is expressed as two unilateral constraints.

In the case of a continuous contact and a world of the first order, the module 211 defines the geometric constraints applied to the speeds of the objects.

A non-penetration constraint for a continuous contact between two objects i and j is expressed by the formula in Appendix 2-I-bl). This formula is advantageously expressed at the center of gravity of the objects i and j, with as unknown speeds possible constraints of rotation and translation.

For m points of contact, a system of m unilateral constraints applied to the N objects in contact makes it possible to define the possible constraint speeds Xv by the formula in appendix 2I-b2). Cv is defined as a matrix of geometric constraints on velocities, of dimension mx6N.

In the case of an instantaneous contact between virtual objects, the constraints defined by the module 212 apply to the speeds. A point of contact H between two objects i and j is defined as Hi on object i and Hj on object j. At a given moment and for a point of contact H between two objects i and j, the hypothesis of Poisson defines a constraint on the velocities of the objects as indicated in the formula of appendix 2-II-1). In this formula, e represents the coefficient of restitution which is chosen according to the materials of the virtual objects in contact. This formula can be expressed in the center of gravity of objects, with as unknown speeds of rotation and translation.

For m points of contact, a system of m constraints makes it possible to define the possible collision velocities Xc by the formula II-2). This is defined as a matrix of contrasts

<Desc / Clms Page number 11>

 on velocities during collision (s), dimension mx6N.

Bc is defined as a vector of the conditions applied to the constraints on the velocities during collision (s).

The module 210 determines a modeling according to the principle of the least GAUSS constraints. Thus, this principle states a minimization of the Ga scalar function based on the unstressed acceleration Au described in Appendix 3-I-a1).

According to the invention, the module 210 is able to define a system of equations, described in Appendix 3-I-a2), to minimize the scalar function Ga as a function of the possible accelerations determined according to constraints on accelerations. The calculator 26 solves this system of equations.

Advantageously, in the case of a world of the first order, the module 211 defines a system of equations given in Appendix 3-II-b2) whose unconstrained acceleration is replaced by the unconstrained speed Vu. Function Ga is replaced by function Gv defined in Annex 3-II-bl). The calculator 26 looks, among the possible speeds, for the constrained speed closest to the speed without constraint, while minimizing the function Gv.

In the case of an instantaneous contact, the module 212 defines a system of equations given in Appendix 3-II-1). For N objects in instantaneous contact, Au = O translates the hypothesis of a negligible intensity of the external actions with respect to the intensity of the actions of the stresses during a collision.

Acceleration is considered constant during the collision. Also, the function Ga is replaced by the function Gc based on the speed before collision V-. The computer 26 looks, among the possible speeds after collision, the speed after collision that minimizes the function Gc.

According to one embodiment, the module 210 is arranged to determine a generalized modeling according to the principle of the least GAUSS constraints. The module 210 is able to define a system of generalized equations described in the appendix

<Desc / Clms Page number 12>

 3-III-a2), to minimize the scalar function Gg in appendix 3III-a1) generalized to a vector without constraint Xk, as a function of constraints generalized to the stress matrix Ck and to the vector of the conditions applied to the constraints Bk.

Thus, the generalized equation system comprises a predetermined primary equation (3-III-a2, Eql) representing the minimization of a primary expression, in particular the scalar function Gg, having for unknown a kinematic magnitude (Xg, 3-III- a2) of each object, and a secondary equation (3-IIIa2, Eq2) able to apply to the kinematic magnitude (Xg, 3III-a2) of each object specific constraints (Ck, 3III-a2, Eq2) respecting specified conditions (Bk, 3-IIIa2, Eq2).

This system of equations defines any contact between virtual objects so that the movement of these objects is calculated by the calculator 26. More precisely, a system of generalized equations states that among the possible unknowns Xg, defined by the equations of m constraints , the unknown Xg * minimizes the generalized Gg function. Assimilated at a kinetic distance, the square root of the Gg function is also minimized by the constrained unknown.

According to a particularly interesting aspect, the module 210 defines a simplified generalized modeling according to the principle of the least constraints of Gauss. From the generalized model of appendix 3-III-a2), the system of equations is reduced to the formula of equations Eql and Eq2 of appendix 3-III-bl) by making a change of variables. As presented in this appendix, the matrix M is positive definite and is decomposed according to the matrices D and F, thus making it possible to obtain the new equations formulated according to the variable X1. C1 is a matrix of dimension mx6N. By making a basic change to this system, the simplified system of generalized equations can be written as shown in Appendix 3-III-b2). X2 and E are defined as vectors

<Desc / Clms Page number 13>

 of size 6N + 1, C2 is a matrix of dimension mx (6N + 1).

On the other hand, the solution X2 * is assumed in a particular form as a function of a vector Q and a component h as defined in Annex 3-III-b2). The vector Q is defined as a vector of dimension 6N. The shape of the vector X1 * is deduced.

According to the invention, the Lagrange method makes it possible, from the simplified generalized equation system 3-III-b2), to obtain the Lagrange 3-III-b3 equation system.

Thus, the invention allows a simple formulation of the problem of contact between objects. This system of equation according to Lagrange facilitates calculations of the calculator by proposing the minimization of a Euclidean norm with as unknown the variable of Lagrange # satisfying a condition of positivity. The Lagrange variable represents a dynamic variable (for example a contact force or a pulse).

This Lagrange equation system comprises a first equation (3-III-b3, Eql) representing the minimization of a Euclidean norm applied to a first predetermined expression having an unknown depending on a dynamic quantity of each object (in particular kinematic) , this first predetermined expression corresponding to a behavior of the non-constrained objects, and said resolution of the first equation being under stress of a second equation (3-III-b3, Eq2) representing at least one constraint condition on the unknown from the first equation.

The calculator 26 solves the Lagrange system which represents one or the other of the modelizations of the modules 211 and 212. The problem of the resolution of this system of equations according to Lagrange consists in a problem of geometric projection of the vector E on a polyhedral cone of size 6N + 1 in the example. According to the invention, the calculator uses a simple algorithm for solving this system of equations. This algorithm is a projection algorithm of a point

<Desc / Clms Page number 14>

 on a polyhedral cone. Various algorithms known from the prior art can be advantageously used and are described in the following works: - K. G. Murty and Y. Fathi, "A Critical Index Algorithm for Nearest Point Problems on Simplicial Cones", 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.

In one embodiment of the invention, the algorithm for projecting a point on a polyhedral cone used is the Wilhelmsen algorithm, for example. An algorithm for projecting a point on a polyhedral cone advantageously makes it possible to find a solution to the modelizations considered to be consistent. In appendix 3-III-b4) is presented the formulation of the solution of vectors X2 * and Xg *.

As an alternative embodiment, from the generalized modeling of Appendix 3-III-a2), the system of equations is reduced directly to the system of equations Eql and Eq2 of Annex 3-III-b5). For this, a Lagrange function L (Xg, #) of appendix 3-III-b6) associated with the modelization of appendix 3-III-a2) is used and the variable Xg is eliminated after having posed equal to zero the derivative of the function L according to the variable Xg. In Appendix 3-III-b5), the matrices M, D and F are identical to the matrices defined in Appendix 3-111bl). Cl is the same matrix of dimension mx6N defined as being a function of matrices Ck, D and F. A vector K is chosen such that it corresponds to the equation defined in appendix 3-III-b5).

A vector S of dimension 6N is defined according to an equation of the same appendix.

<Desc / Clms Page number 15>

The new system of equations Eql, Eq2 can be written as indicated in this appendix 3-III-b5) corresponding to a Lagrange system. This system of equation having the same form as the system of equation 3-III-b3), the problem of its resolution consists of a problem of geometric projection of the vector S on a polyhedral cone of dimension 6N. The resolution is thus carried out using the same type of algorithm as before, in particular an algorithm for projecting a point on a polyhedral cone, for example the Wilhelmsen algorithm.

Part of the code according to the invention is presented in Annexes A and B.

In Annex A, the constraint management is presented according to three functions: - ConstrainGeneralVelocity () which allows the computation of the constraint velocities, - ConstrainGeneralAcceleration () which allows the computation of the constraint accelerations, - Resolvelmpact () which allows the computation of the rebound velocities during a collision.

In Appendix B, Wilhelmsen's algorithm is presented according to three other functions: - cFilter () which allows the transformation of the problem into a problem of the Wilhelmsen form, and calls the function cProjection (), - cProjection () which defines the Wilhelmsen algorithm, - cLinearCholesky () which solves a matrix equation AX = b, which is called by the Wilhelmsen algorithm.

According to the invention, the module 210 defines a group of contacts representing N moving objects forming a chain of contacts between the first object and the nth object. Thus groups of contacts are dynamically independent of one another and can be subject to independent simulations.

Figures 3 to 4 show virtual objects in contact to illustrate the contact modelings according to the invention.

<Desc / Clms Page number 16>

FIG. 3 represents a first object in contact with a second object Q, this second object Q being stationary. Thanks to the data on the environment of the objects, for example the gravity, the module of continuous contact modeling in relation with the computer determines the acceleration without constraint Au. This unrestrained acceleration Au represents the acceleration that the object would have if it were not subject to non-penetration constraints. In the example of FIG. 3, the unrestrained acceleration Au of the virtual object R corresponds to the severity of this same object. According to the non-penetration constraint indicated in appendix 2-I-a1) and applied to the virtual object R, the possible accelerations of this same object are oriented according to the slope LA of the virtual object Q. Among these possible accelerations, the The acceleration acceleration Ac of the virtual object R corresponds to the nearest possible acceleration of the acceleration without constraint Au.

The acceleration acceleration Ac, as defined by the Gauss principle of the least constraints in appendix 3-I-a2), is oriented in FIG.

Figure 4 provides an example of a collision between two virtual objects R and Q. In the example, the object Q is immobile.

Thus, the object R is animated with a speed V- being the speed before collision. At the instant Tcol the virtual object R comes into contact with the virtual object Q at a point H. Hypotheses are made concerning the instantaneous contact. Thus, a first hypothesis considers the duration of the collision as infinitesimal. In other words, the position of the object R remains unchanged, as well as the accelerations of stress.

Another hypothesis considers that the data on the environment of the objects or on the external actions, can be neglected during the collision compared to the intensity of the actions of constraint.

As indicated in Appendix 3-II-2), GAUSS's principle of least constraints is arranged so that the speed after collision is the unknown. Thus, the speed after collision V +, among the possible speeds after collision, minimizes the

<Desc / Clms Page number 17>

 Gc function, developed in Annex 3-II-1). These possible velocities after collision are determined by the collision response constraint according to the Poisson hypothesis detailed in Appendix 2-II. Thus, in a group of contacts, the possible speeds after collision are the speeds closest to the speeds before collision.

In FIG. 4, the set of possible speeds after collision is represented by the half-space situated above the line LC. The speed after collision V + is, according to GAUSS's principle of the least constraints, as close as possible to the speed before collision V-.

In a world of the first order and in the case of a continuous contact, GAUSS's principle of least constraints is to define the object constraint speed of a group of contacts. This contact group is subjected to a non-penetration constraint applied to the speeds of the objects and indicated in Annex 2-I-bl). This linear function makes it possible to determine the possible speeds of the objects of the group of contacts.

Finally, in the function Ga indicated in the appendix 3-I-al), the acceleration without constraint is replaced by the unconstrained speed Vu defining a new function Gv in the appendix 3-I-bl). Among the possible speeds, the constrained speed of the group of contacts, is chosen to minimize the function Gv. In other words, at any time, the contact group stress rate is chosen to correspond to the closest possible speed to the unconstrained speed.

According to the invention, GAUSS's principle of least constraints allows the computer to find the movement of the virtual objects in contact. Advantageously, the invention is capable of defining, in all cases of contact, a set of solutions for the accelerations and the possible speeds determined by the constraints applied to the objects.

<Desc / Clms Page number 18>

 virtual. In this set of solutions, a single solution can be determined.

According to another advantage of the invention, each contact model can be simplified in two steps. These two stages are a change of variable and a homogenization of the considered model. Thus, the contact models are reduced to a single model described in Appendix 3-111. Advantageously, the unknown of this unique model is found after geometric projection in a space of determined dimension, 6N or 6N + 1 where N is the number of objects of a group of contacts.

At the computer level, the problem solving by geometric projection is performed by an algorithm for projecting a point on a polyhedral cone known from the state of the art. Advantageously, this simple algorithm can be Wilhelmsen's algorithm.

FIG. 5 proposes a method for simulating rigid virtual objects according to the invention.

At a time T, the object data defines the N objects of the contact group. These data are in particular the position P, the velocity V, the generalized acceleration A in the case of a world of the second order (also called mode A), the position P, the velocity V in the case of a world of the first order (also called V mode). Thus, we define the virtual objects in their environment.

If constraints are applicable on an object or a group of contacts in step 102, a modeling of the continuous contact is carried out according to the invention. In order to be able to process this modeling numerically, a method of solving differential equations is applied to the modeling considered in step 104. Advantageously, this method of solving differential equations can be the Runge-Kutta method. Over a time interval dT, this resolution method applies the algorithm of

<Desc / Clms Page number 19>

 Wilhelmsen. The application of step 104 makes it possible to determine in step 108: the constrained acceleration Ac sought if the simulation is carried out in a world of the second order; the constrained speed Vc sought if the simulation is carried out in a world of the first order.

Advantageously, the time interval can be variable.

In step 102, if no constraint is applicable on an object or a group of contacts, in step 110 the generalized acceleration equals the unconstrained acceleration Au if the simulation is performed in a world of the second order and the generalized velocity equals the unconstrained velocity Seen if the simulation takes place in a world of the first order.

If no collision is detected in step 114, the new speed V1 and position P1 are calculated for a time T equal to T + dT. If, on the other hand, a collision is detected at time Tcol equal to T + a.dT with a less than 1, at step 116: the new speeds V- and position Pi are calculated before the collision if the simulation is carried out in a world of the second order; the new position Pi is calculated before the collision if the simulation is carried out in a world of the first order.

In step 117, the process returns to step 100 if the simulation is in a first-order world. In the case of a simulation in a world of the second order, the instantaneous contact is modeled according to the invention in step 118. In step 120, the Wilhelmsen algorithm is applied to the modeling of the instantaneous contact. This algorithm makes it possible to find a unique solution for the speed after collision V + in step 124. Returning to step 100, the movement of the virtual objects is calculated no time after no time.

In an exemplary embodiment of the invention, a continuous collision detection technique is applied to the objects.

<Desc / Clms Page number 20>

This technique assigns rigid fixed displacements to the objects to interpolate the successive positions of these objects and thus to detect any collisions and moments of impact between these objects.

Thus, the simulation according to the invention provides, among other things, an improvement as regards the numerous disadvantages of the prior simulation methods.

The invention is applicable in many industrial fields and allows, for example, to present virtual products in dynamic simulation, to simulate in real time a virtual assembly of parts on assembly / disassembly lines, to simulate operating tests. The invention can also be applied in the field of virtual games.

Of course, the invention is not limited to the embodiment described above by way of example, it extends to other variants.

Advantageously, other algorithms can solve the projection problem concerning the modeling of the continuous contact.

<Desc / Clms Page number 21>

Appendix A Al-ConstrainGeneralVelocityO void cScene :: ConstrainGeneralVelocity (pcContactGroup cg, double * intmd, double * conmd) {
This function is spread over the speed space to control the intentional movement of the objects so that they respect the constraints on their movement. It is called more often, during the same time step, by the function of cRGVelocity. NB: function does not update the list of contacts, nor their
It's their number. nobjectsis the total number of objects in the group of contacts (mobilesandimmobiles). nmobile is the number of moving objects in the group of contacts. ncontacts is the number of contacts. Contact is // an array of pointers to classrooms. Contact, tactless contacts. invidestabletotalfullnumber // nobjects, which matches the object's ID to the mobile object list (when the object is immobile, its place is -1). mobp is an integer array of nmobile size, which gives the // inverse match: the ID number of the mobile object whose motion is and filter. intm is a table // doubletetaille 6 * nmobile representing the intentional movement of mobile objects (enradians). conm is // a double table of size 6 * nmobile representing the final filter movement of moving objects (in // radians).

 NB: in this function, we also protect the safety distance inti; // contact group variables int ncontacts = cg-> GetNContactsO; pcContact * contact = new pcContact [ncontacts]; for (i = O; i <ncontacts; i ++) contact [i] = cg-> ContactPointer (i); // A-1-a First step: we transform the contacts into constraints on the initialization movement of the array, we suppose that it can not be more than 50 // NB: we suppose that the list of contacts (contact) is ajour (Update function was called) intnee = 0; // number of constraints intnvf = 0; number of constraints vf intnfv = 0; number of constraintsfv int spacedim = 6 * cg-> GetNObjects (); double * cmat = new double [(spacedim + l) * ncontacts]; double * intmh = new double [spacedim + 1]; // intentional movementgeneralizedhomogenized double * conmh = new double [spacedim + l]; // homogenized constraintgeneralized motion for (i = O; i <spacedim; i ++) intmh [i] = intmd [i]; intmh [spacedim] = 1.0;

<Desc / Clms Page number 22>

switch (contaet[i]->GetTypeO) { case VF~CONSTRAINT: nvf++; break ; case EE~CONSTRAINT: nee++ ; break ; case FV~CONSTRAINT: nfv++ ; break ; default : break ; } // il y a 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<spacedim+ l;j++) cmat[idim+j]=0; // coefficients relatifs aupremierobjet cmat[offl]=-n[0]; cmat[offl+l]=-n[l]; cmat[offl+2]=-n[2]; cmat[offl+3]=-cprod[0]; cmat[offl+4]=-cprod[l]; cmat[offl+5]=-cprod[2]; for (i = O; i <ncontacts; i ++) {// we add a constraintcontact // count constraints

switch (contaet [i] -> GetTypeO) {case VF ~ CONSTRAINT: nvf ++; break; case EE ~ CONSTRAINT: nee ++; break; case FV ~ CONSTRAINT: nfv ++; break; default: break; } // there are two objects involved in a contact // Let two mobile objects, one movable object and one hard object // have a constraint on the relative movement

cMatrix4gml (contact [i] ->GetO10->GlobalMatrix0); cVertexpivotgtobl = gml * Contact [i] -> Get01 () ->GetPivotO; cVertex gcl (contact [i] ->GetLocationn-pivotglobi); cVertex n (contact [i] ->GetNormalo); calculation of the coefficients of the constraint // n: normal vector at the plane of stress in the global coordinate system // NB: n is unitary assumption // coefficients of the constraint cVertex cprod (gcl ^ n); int idim = i * (spacedim + l); // position of constraint in the table

int indldim = cg-> InvId (contact [i] ->Get0l-> GetIdn) * 6; int offl = idim + indldim; for (int j = O; j <spacedim + l; j ++) cmat [idim + j] = 0; // coefficients for the first object cmat [offl] = - n [0]; cmat [offl + l] = - n [l]; cmat [offl + 2] = - n [2]; cmat [offl + 3] = - PROCOD [0]; cmat [offl + 4] = - PROCOD [l]; cmat [offl + 5] = - PROCOD [2];

<Desc / Clms Page number 23>

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<spacedim;j++) cmat[idim+j]=-cmat[idim+j]; //A-1-b- Deuxième étape : on contraint le mouvement // transformation du mouvement intentionneldes objets (passage des degrés aux radians) cFilter(spacedim+l,ncontacts,cmat,intmh,conmh); // appela Wilhelmsen for (i=O;i<spacedim;i++) conmd[i]=conmh[i]; delete [] cmat; delete [] intmh; delete [] conmh; delete [] contact; }

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 // coefficients related to the second object

int ind2dim = cg-> InvId (contact [i] ->Get02->GetId); if (ind2dim> -l) {// the second object ismobile also ind2dim * = 6; int off2 = idim + ind2dim; cmat [off2] = n [O]; cmat [off2 + l] = n [t]; 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: for face / point contact, consider the relativeoppose motion: if (#ntact [i] -> GetTypeO == FV ~ CONSTRAINT) {for (j = O; j <spacedim; j ++) cmat [idim + j] = - cmat [idim + j]; // A-1-b- Second step: we constrain the motion // transformation of the intentional movement of objects (transition from degrees to radians) cFilter (spacedim + l, ncontacts, cmat, intmh, conmh); // called Wilhelmsen for (i = O; i <spacedim; i ++) conmd [i] = conmh [i]; delete [] cmat; delete [] intmh; delete [] conmh; delete [] contact; }

A-2- ConstrainGeneralAccelerationp void cScene :: ConstrainGeneralAcceleration (pcContactGroup cg, double * intmd, double * conmd) {
This function is placed in the accelerator space to control the intentional movement of the objects so that they respect the constraints on their movement. It is called several times, during the same
It is not time, by the decRGAcceleration function. NB: this function did not update the listescontacts, // nor their number. nobjects is the total number of objects in the contact group (mobile and immobile). he

<Desc / Clms Page number 24>

int n#ntacts=cg->GetNContactsO; pcContact *contact=new pcContact[ncontacts]; for (i=O;i<ncontacts;i++) contact[i]=cg->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; nombre de 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<spacedim;i++) intm [i]=intmd[i]; intm[spacedim]=1.0; // attention a cette lignequand ongérera la distancenon euclidienne cnumber=-l; // lapremière contrainte a le numéro 0 nmobile is the number of moving objects in the contact group. ncontacts is the number of contacts. He is an array of pointers to classes, contact, and contact. invidestun array of integers // detaillenobjects, which matches the ID of the object to themobilehostlist // (when the object is still, saplace is -1). mobp is an integer array of nmobile size, which gives the // inverse match: the ID number of the mobile object whose motion has and filters. vdint is a size 6 * nmobile double table representing the intentional movement of mobile objects (in radians). vdfiltest // a double table of size 6 * nmobile representing the final filter movement of moving objects (in Il radians). Lafriction is drained by adding constraints to the set // unilateral stresses. Therefore, even if the number of contact points is known, there is not yet a number of constraints. NB: in this function, we also overwrite the security distance; // variables related to the decontact group

int n # ntacts = cg->GetNContactsO; pcContact * contact = new pcContact [ncontacts]; for (i = O; i <ncontacts; i ++) contact [i] = cg-> ContactPointer (i); // A- 2-a .. First step: we transform the contacts into constraints on the // initialization motion of the array, we suppose that it can not get more than 50 // NB: we suppose that the list of contacts (contact) is ajour ( the function Update) was called intnee = 0; // number of constraints ee int nvf = 0; number of constraints vf int nfv = 0; // number of constraintsfv int nconstraints; // number of constraints at the end int cnumber; // number of the unilateral constraint under construction // count of the final number of constraints nconstraints = ncontacts; creation ofprimary variables (before Euclidean distance) ll these variables are homogenized int spacedim = 6 * cg-> GetNObjects (); // number of degrees of freedom in the contactgroup double * cmat = new double [(spacedim + l) * nconstraints]; double * intm = new double [spacedim + l]; double * conm = new double [spacedim + l]; for (i = O; i <spacedim; i ++) intm [i] = intmd [i]; INTM [spacedim] = 1.0; pay attention to this line when the Euclidean distancenon will be increased cnumber = -l; // the first constraint has the number 0

<Desc / Clms Page number 25>

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 for (i = O; i <ncontacts; i ++) {// we examine the contacts one by one if (contact [i] -> Error ()) {cout "The contact" i "is an incorrect edge / edge contact (edges parallel). "endl;break; enumeration of the constraints switch (contact [i] -> GetTypeO) {case VF ~ CONSTRAINT: nvf ++; break; case EE ~ CONSTRAINT: nee ++; break; caseFV ~ CONSTRAINT: nfv ++; break; default: break; } // A-2-b Creation of the non-interpenetration constraint cnumber ++; // there are two objects involved in a contact
Let two mobile objects, one mobile object, and one object of the design lion have a constraint on the relative motion 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; // compute coefficients of the constraint // n: normal vector at the plane of constraint in the globalphere // NB: n are unitary // NB: the coefficients of beforedxdy dz wxwy wz (which here represent the accelerations) are

<Desc / Clms Page number 26>

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<spacedim+l;j++) cmat[idim+j]=0; // coefficients relatifs aupremier objet cmat[offl]=-n[0]; cmat[offl+l]=-n[l]; cmat[offl+2]=-n[2]; cmat[offl+3]=-cprod[0]; cmat[offl+4]=-cprod[l]; cmat[offl+5]=-cprod[2]; // coefficients relatifs audeuxième objet

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. // the same as for speeds. Only the term constantchange. gm1 = (O1-> GlobalMatrix ()); pivotglobl = (gml * ol -> GetPivot ());

gel = (contact [i] -> GetLocationO-pivotglob 1); n = (contact [i] ->GetNormaIO); // coefficients of the constraint cVertex cprod (gcl ^ n); int idim = cnumber * (spacedim + l); // position of the constraint in the table

int indldim = cg-> InvId (ol-> GetldO) * 6; int offl = idim + indldim; for (int j = O; j <spacedim + l; j ++) cmat [idim + j] = 0; // coefficients related to the first object cmat [offl] = - n [0]; cmat [offl + l] = - n [l]; cmat [offl + 2] = - n [2]; cmat [offl + 3] = - PROCOD [0]; cmat [offl + 4] = - PROCOD [l]; cmat [offl + 5] = - PROCOD [2]; // relative coefficients for the second object

int ind2dim = cg-> Invld (o2->Getldo); if (ind2dim> -l) {// the second object ismobile too ind2dim * = 6; int off2 = idim + ind2dim; cmat [off2] = n [0]; cmat [off2 + l] = n [t]; 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]; // calculation of the constant term / / NB: this term contains centripetal acceleration (w ^ (w ^ GC)) # n, for both objects.

 // It also holds 2 * (va-vb) lndot, where np is the derivative with respect to the time of the normal vector.

 Finally, it contains a term function of the distance to the safety distance (DISPLACEMENT + contact [i] -> Penetration ()),

<Desc / Clms Page number 27>

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<spacedim;j++) cmat[idim+j]=-cmat[idim+j]; } // A-2-c. Deuxieme étape : on contraint le mouvement

cFilter(spacedim+l,nconstraints,cmat,intm,conm); // appel a Wilhelmsen for (i=O;i<spacedim;i++) conmd [i]=conm[i]; delete [] cmat; delete [] intm; delete [] conm; delete [] vt; delete [] dfriction; delete [] striction; delete [] contact ; It acts as a spring in the methods of vanity. cVertex wl, w2; // centripetal acceleration of the object 1 wl = (ol->GetRotVelocityO); double dis = - ((w1 ^ (w1 ^ gcl)) n); if (ind2dim> -l) {// the second object ismobile // acceleration centripetede the object2

w2 = (o2-> GetRotV elocityO); say + = ((w2 ^ (w2 ^ gc2)) n); // displacement of the normal vector

cVertex ndot (contact [i] ->GetNormalDoto); // derived from the time of the normal vector vcl = (ol-> GetTransVelocityO + (wl "gel; // velocity of the contact point in object 1 vc2 = (o2-> GetTransVelocityO + (w2Agc2)); // speed of point of contact in object 2 dis = 2. 0 * ((vc1-vc2) #ndot);} maintaining the distance of safety

double delta = DISPLACEMENT + contact (i) ->Penetration0; dis + = (KDISPLACEMENT * delta); cmat [idim + spacedim] = dis; // NB: for contacts of typeface / point, we consider the relative motion opposes:

if (contact (i) -> GetType == FV CONSTRAINT) {for (j = O; j <spacedim; j ++) cmat [idim + j] = - cmat [idim + j];} // A-2-c Second step: we constrain the movement

cFilter (spacedim + l, nconstraints, cmat, INTM, conm); // call to Wilhelmsen for (i = O; i <spacedim; i ++) conmd [i] = conm [i]; delete [] cmat; delete [] intm; delete [] conm; delete [] vt; delete [] dfriction; delete [] narrowing; delete [] contact;

<Desc / Clms Page number 28>

} A-3. ResolveImpact () void cScene :: ResolveImpact (pcContactGroup cg) {
This function is located in the velocity space to calculate the motion of objects colliding.

 Implicitly, it is therefore a simultaneous calculation of pulses. Explicitly, one directly calculates // the movement of the objects. This function is similar to the other functions of the constraintsNB //: this function has not updated the list of contacts, nileurnumber. It is assumed, therefore, that the colliding function of contacts can safely be called to find out whether the objects collide at that point of contact or not.

 It is the total number of objects in the group of contacts (mobilesandimmobiles). nmobile is the number of moving objects in the group of contacts. ncontacts is the number of contacts.

cVertex tspeed( mobp[i]->GetTrans V elocityO); cVertex rspeed(mobp[i]->GetRotVelocity); int i6=i*6; intm [i6]=tspeed[0]; There is a pointer to the contactclasses, to the contact.mustestabletheintenancefaces, which matches the object's ID to the movableobjectlist Il (when the estimemobile object, its place is -1). mobp is an array of pcObjectdetaillenmobile. mobp [i] gives // a pointer to the object whose place is j in the current list of moving objects. NB: useful to manage // also the safety distance in this function? int int nmobile = cg-> GetNObjects (); pcObject * mobp = new pcObject [nmobile]; // mobile objectpointer for (i = O; i <nmobile; i ++) mobp [i] = objects [cg-> Id (i)]; int ncontacts = cg-> GetNContacts (); pcContact * contact = new pcContact [ncontacts]; for (i = O; i <ncontacts; i ++) contact (i] = cg-> ContactPointer (i); // A-3-a First step: we transform the contacts into constraints on the motion // NB: we assumes that the contact list is up to date (the Update function has been called) intnee = 0; // number of constraints int nvf = 0; // number of constraints vf int nfv = 0; number of constraintsfv int spacedim = 6 * nmobile double * cmat = new double [(spacedim + l) * ncontacts]; double * intm = new double [spacedim + l]; double * conm = new double [spacedim + l]; // we store velocities of objects in a common vector // make it big array of matches mobp for (i = O; i <nmobile; i ++) {

cVertex tspeed (mobp [i] -> GetTrans V elocityO); cVertex rspeed (mobp [i] ->GetRotVelocity); int i6 = i * 6; intm [i6] = tspeed [0];

<Desc / Clms Page number 29>

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; 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; // compute the array of constraints for (i = O; i <ncontacts; i ++) {// add a constraint by contact // enumeration of the constraints switch (contact [i] -> GetType ()) {case VF ~ CONSTRAINT: nvf ++ ; break; case EE ~ CONSTRAINT: nee ++; break; case FV ~ CONSTRAINT: nfv ++; break; default: break;
There are two objects involved in a contact // Two moving objects, one moving object and one decor object // one has a constraint on the relative movement

cMatrix4gml (contact [i] ->GetO10->GlobalMatrix0); cVertexpivotgIobl = gml * Contact [i] -> Get01 () ->GetPivotO; cVertex gcl (contact (i) ->GetLocationQ-pivotglobl); cVertex n (contact [i] ->GetNormalo); // compute coefficients of the constraint // n: normal vector at the constraint plane in the global coordinate system // NB: nestsuppose unit coefficients of the constraint cVertex cprod (gc1 ^ n); int idim = i * (spacedim + l); // position of the constraint in the table

intindldim = CG-> Invid (contact [i] ->GetO10-> GetId0) * 6; int offl = idim + indldim;

<Desc / Clms Page number 30>

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<spacedim;j++) cmat[idim+j]=-cmat[idim+j]; } for (int j = O; j <spacedim + l; j ++) cmat [idim + j] = O; // coefficients related to the first object cmat [offl] = - n [0]; cmat [offl + l] = - n [l]; cmat [offl + 2] = - n [2]; cmat [offl + 3] = - PROCOD [0]; cmat [offl + 4] = - PROCOD [l]; cmat [offl + 5] = - PROCOD [2]; // coefficients relative to the second object

int ind2dim = cg-> InvId (contact [i] ->Get02n->GetId); if (ind2dim> -1) {// the second mobile object is ind2dim * = 6; int off2 = idim + ind2dim; cmat [off2] = n [0]; cmat [off2 + l] = n [t]; 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]; } // calculate the constant coefficient double e = 0.34; // coefficient of restitution (good value: 0.38)

double nspeed = -e * (contact [i] ->RelativeVelocityoln); // normal speed // relative to the contact point if (nspeed> DISPLACEMENT) cmat [idim + spacedim] = nspeed; // elsecmat [idim + spacedim] = O; else

ecmat [idim + spacedim] = 0.1 * KDISPLACEMENT * (DISPLACEMENT + contact [i] ->PenetrationO); // NB: for face / point contact, we consider the relative motion opposite: if (contact [i] -> GetTypeO == FV ~ CONSTRAINT) {for (j = O; j <spacedim; j ++) cmat [idim + j] = - cmat [idim + j]; }

<Desc / Clms Page number 31>

 // A-3-b. Second step: constrain the movement cFilter (spacedim + l, ncontacts, cmat, intm, conm); // called Wilhelmsen // modification of the object velocities for (i = O; i <nmobile; i ++) {inti6 = i * 6; cVertex tspeed (conm [i6], conm [i6 + 1], conm [i6 + 2]); cVertex rspeed (conm [i6 + 3], conm [i6 + 4], conm [i6 + 5]); MOBP [i] -> SetTransVelocity (Tspeed); MOBP [i] -> SetRotVelocity (rspeed); } delete [] cmat; delete [] intm; delete [] conm; }

<Desc / Clms Page number 32>

if ((scal)&&(scal≥-P~EPSILON)&&(scal≤P~EPSILON)) { // cout "risque de bugpour scal dans cProjection !!!!!!!!!!!!!" #endl; // scal=0; } Appendix B Wilhelmsen algorithm Bl.cProjectionO void cProjection (int n, int m, double c [], double v [], double s []) {
This function uses Wilhelmsen's algorithm to find the most prochedupointpoint in the positive polyhedralone vector by m vectors. These vectors are stored in the matrix matrix n * m. All vectors are in a space of size n. NB: the jccoordordonnee of the vectorciesten c [i * n + j] Thus, when one reads the coordinates in the order, one does only one multiplication. We suppose that the // vectors ci are all of norm 1. we start by setting the norm of all the constrained vectors to CNORME, in order to avoid the errors of rounding. NB: it would be necessary to take account of the fact that one projects the point // v, which is very particular (from almost all) int i, j; int mp; // cardinal of the set F int mm = m * m; // initialize variables for (i = O; i <m; i ++) {innervtag [i] = 0; } for (i = 0; i <(n * m); i ++) c [i] * = CNORME; // we change the norm of constrained vectors // Bla. step 0 // we calculate all the innerv of a hit and take the largest int found = 0; int index = 0; double scalmax = 0; for (i = O; i <m; i ++) {// calculation of the scalar product v.ci int in = i * n; double scat = 0; for (j = O; j <n; j ++) scal + = v [j] * c [in + j];

if ((scal) && (scal≥-P ~ EPSILON) && (scal≤P ~ EPSILON)) {// cost "risk of bug for scal in cProjection !!!!!!!!!!!!!"#endl; // scal = 0; }

<Desc / Clms Page number 33>

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<n;i++) s[i]=O; // on nettoie la memoireet on arrete return ; //a partird'ici, l'originen'estpas lepoint leplusproche de v on peut passer a l'etape 1. L'indice du vecteur choisiest indice. innerv [i] = scal; // store the result

if ((scal> P EPSILOI &&(scal> scalmax)) {found = l; index = i; scalmax = scal; // NB: compared to the old ones, the innerv coefficients are multiplied by CNORME if (! found) {// the origin is the closest point of v for (i = O; i <n; i ++) s [i] = O; // clean the memory and stop return; // from here, the origin is not the closest point of v on can go to step 1. The index of the chosen vector is index.

// initialize the double rhomin variables; double comprhomin; for (i = O; i <mm; i ++) innertag [i] = 0; for (i = O; i <m; i ++) e [i] = O; e [index] = l; // the set e1 is the singleton c (index) double temp = innerv [subscript]; lambda [index] = Temp; int indicen = index * n; for (j = O; j <n; j ++) p [j] = temp * c [indicen + j], // calculation of the current project of {// start of the outer loop // B-1-b. step 1 // calculating the difference vector, and its norm, to see sivest in the subconect for (i = O; i <n; i ++) d [i] = CNORME2 * v [i] -p [i]; // v is not in the current sub-cone. We are looking for a point on the wrong side of the hyperplane. it is not possible to test the points which are in the set of active constraints, since they are in the hyperplane. This allows moreover to avoid the errors of rounding // In addition, one calculates them all, to take the maximum

<Desc / Clms Page number 34>

if ((scal>P EPSILON)&&(scal>scalmax)) { found=1; indice=i ; scalmax=scal; } } } if (!found) { // le projeté courant est le pointsolution for (i=O;i<n;i++) s[i]=INVCNORME2*p[i]; break ; } else // a partir d'ici, le projeté courant n'est pas le point solution lion peut passer a l'étape2. L'indice du vecteurchoisi est indice. index = 0; found = 0; double scalmax = 0; for (i = O; i <m; i ++) {if (! (e [i])) {// if the vectori is not part of the active set
It calculates the scalar product d.ci int in = i * n; double scal = 0; for (j = O; j <n; j ++) scal + = d [j] * c [in + j];

if ((scal> P EPSILON) &&(scal> scalmax)) {found = 1; index = i; scalmax = scal; }}} if (! found) {// the current project is the pointsolution for (i = O; i <n; i ++) s [i] = INVCNORME2 * p [i]; break; } else // from here, the current project is not the solution point lion can go to step2. The index of the vector chosen is index.

// we initialize the variables // creation of the setF 1 [0] = index; // put the selected vector in the list lambda [index] = 0; mp = 1; for (i = O; i <m; i ++) {if (e [i]) {
1 [mp] = i; mp ++; }}
It creates barycentric coordinates of R // at the beginning R = pj // NB: all the lambda are multiplied by CNORME for (i = O; i <mp; i ++) lambdar [i] = lambda [l [i]]; // the current project is for the moment unchanged

<Desc / Clms Page number 35>

if (!(innertag[1[j]*m+I[i]j)) int in=l [i]*n; int jn=l [j]*n; double scat=0; for (int k=O;k<n;k++) scal+=c[in+k]*c[jn+k]; // NB:lesnouveauxinnersontmultiplies par CNORME2 par rapportauxanciensdemême pourlescoefficientssubmat inner[l[j]*m+l[i]]=scal; // onstocke lerésultat

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<mp;i++) submat[i*mp+i]=CNORME2; //calculdu sous-vecteur. for (i=O;i<mp;i++) subvec[i]=innerv[l[i]]; //calculduprojetéorthogonal do {// start of the inner loop if (mp == 0) {for (i = O; i <m; i ++) {e [i] = 0; lambda [i] = 0; } break; } // B-1c, step 2 2 // we calculate the orthogonal projection of the new set of constraints // we want to get the beta / compute coefficients of the sub-matrix. This sub-matrix is symmetrical. There are only 1 on the diagonal. Some elements have already been calculated and stored in the inner array. if (mp> n) cost "pb: non-independent subset of constraints (mp>n)"endl; for (i = 0; i <(mp-1); i ++) {for (j = (i + 1); j <mp; j ++) {

if (! (innertag [1 [j] * m + I [i] j)) int in = l [i] * n; int jn = l [j] * n; double scat = 0; for (int k = 0; k <n; k ++) scal + = c [in + k] * c [jn + k]; // NB: thenewinfinalsareincreased by CNORME2 in comparison with the former for the coefficients in the inner subgroup [l [j] * m + l [i]] = scal; // onstocke the result

innertag [1 (j] * m + 1 [i]] = 1; I / we note the fact that the calculation has already been done}

submat (j * mp + i) = inner [1 [j] * m + 1 [i]); submat [i * mp + j] = submat [j * mp + i]; } // NB: do not forget the diagonal coefficients for (i = O; i <mp; i ++) submat [i * mp + i] = CNORME2; // calculation of the sub-vector. for (i = O; i <mp; i ++) subvec [i] = innerv [l [i]]; // calculduprojetéorthogonal

<Desc / Clms Page number 36>

for (i=0;i<n;i++) p[i]+=Iambda[t[j]]*c[jn+i]; } // on sortde la boucle interne // ieon retourne au début de la boucle externe cLinearCholesky (mp, SUBMAT, subvec, proj); // consider multiplication by CNORME for (i = 0; i <mp; i ++) proj [i] * = CNORME2; // look now where is the projected int neg = 0; i = 0; do if (proj [i] <0) neg = 1; i ++; } while ((! neg) && (i <mp)); if (! neg) // the orthogonal projection is in the new sub-cone // we retain only the related constraints // we add the barycentric coordinates for (i = 0; i <m; i ++) {e [i] = 0 ; lambda [i] = 0; } for (i = 0; i <mp; i ++) {if (proj [i] == 0) {e [l [i]] = 0; } else {e [I [i]] = l; lambda [l [i]] = proj [i]; }} // we calculate the coordinates of the current new project for (i = 0; i <n; i ++) p [i] = O; for (j = 0; j <mp; j ++) {intjn = 1 [j] * n;

for (i = 0; i <n; i ++) p [i] + = Iambda [t [j]] * c [jn + i]; } // get out of the inner loop // ie return to the beginning of the outer loop

<Desc / Clms Page number 37>

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<mp;i++) { if(rho[i]<rhomin) rhomin=rho [i]; } // calculdes gam[i] comprhomin=l .0-rhomin; for (i=O;i<mp;i++) gam[i]=comprhomin*lambdar[i]+rhomin*proj[i]; lionpeut maintenant lancer l'étape3pour réduire mp // B-l-d. étape 3 //mise ajourde la liste des indices (nouveauF) int mptemp=0; //nouveau cardinal for (i=O;i<mp;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<mp;i++) l[i]=ltemp[i]; break; } // theorthogonalproject is not in the new sub-cone, it then searches for a new point R on the cone boundary // andonlyesquirements the oldR point to the orthogonal project // rho determination [i] for (i = O; i <mp; i ++) {double diff = lambdar [i] -proj [i];

if ((diff) &&(diff> -P EPSILON) && (diff≤P EPSILON)) {} if (diff> 0) rho [i] = lambdar [i] / diff; else rho [i] = 1.0; } // determine the minimum of rho [i] rhomin = rho [0]; for (i = l; i <mp; i ++) {if (rho [i] <rhomin) rhomin = rho [i]; } // calculation of gam [i] comprhomin = l .0-rhomin; for (i = O; i <mp; i ++) gam [i] = comprhomin * lambdar [i] + rhomin * proj [i]; lioncan now start step3to reduce mp // Bld. step 3 // update the index list (newF) int mptemp = 0; // new cardinal for (i = O; i <mp; i ++) {

if gam [i]) && (gam [i]> - P ~ EPSILON) && (gam [i] ≤P - EPSILON {} if (gam [i]> P ~ EPSILON) {It to avoid rounding errors ltemp [mptemp] = l [i]; lambdartemp [mptemp] = gam [i]; mptemp ++;}} mp = mptemp; for (i = O; i <mp; i ++) l [i] = ltemp [i];

<Desc / Clms Page number 38>

 for (i = O; i <mp; i ++) lambdar [i] = lambdartemp [i]; // we can now go back to step2} while (1); // end of the inner loop while (1); // end of outer loop lion cleans memory and stops return; } B-2. cLinearCholeskyO void cLinearCholesky (int n, double to [], double b [], x []) {// this function solves the equation ax = b, where a is an invertible nxn matrix // In a first time, we make a Cholesky decomposition delamatrice in table a, we access aa (i, j) by doing a [j * n + i]; inti, j, k; double sum; double * p = new double [n]; for (i = O; i <n; i ++) {for (j = i; j <n; j ++) {sum = a [i * n + j]; for (k = i-1; k≥O; k--) sum - = (a [i * n + k] * a [j * n + k]); if (i == j) {if (sum≤0.0) {cost "sum =" sum endl; } p [i] = sqrt (sum); } else a [j * n + i] = sum / p [i]; }} // solving the systemLUx = b first solves Ly = b (and stores the result inx) for (i = O; i <n; i ++) x [i] = 0.0; for (i = O; i <n; i ++) {sum = b [i]; for (k = i-1; k≥0; k--) sum- = a [i * n + k] * x [k]; x [i] = sum / p [i];

<Desc / Clms Page number 39>

}
Then, let Ux = y for (i = (n-1); i≥O; i--) sum = x [i]; for (k = (i + l); k <n; k ++) sum- = a [k * n + i] * x [k]; x [i] = sum / p [i]; } delete [] p; } B-3. cFilterO void cFilter (int n, int m, double c [], double v [], x []) {// this function returns the closest point of point v in the form set to the intersection of m Ildemi -spaces The normals with half-spaces are vectors, stored in the region. Toutspasse in // dimension n. Dual Formulation is used to apply the Wilhelmsen algorithm.

 NB: the coordinate junction of the vector is in c [i * n + j] // Thus, when we read the coordinates in order, we only produce a multiplication int i, j; // Before solving the dual problem, we must make a change of variables // NB: for the moment, we assume that the mass matrix is the double identity * c2 = new double [n * m]; double * in = new double [n]; // vector (0,0,0, ..., 0,1) of size n double * x2 = new double [n]; // solution of the dual problem for (i = 0; i <(n-1); i ++) in [i] = O; in [n-l] = 1000; // computation of cl (cl = CP (l / sqrt (D)) // ci is stored in c2 // cl is of size m * (nl) // NB: icicl = c for (i = O; m; i ++) {int in = i * n; for (j = 0; j <(n-1); j ++) c2 [in + j] = c [in + j]; // computation of b1 (bl = B -CXu) // we must store -b1 in the matrix c2, that is -B + CXu // or -B is the last column of c // and Xu is v

<Desc / Clms Page number 40>

// calcul dex x (x=(1lh).P.(llsqrt(D))A+Xu) double invh=1.0/x2[n-l]; for (i=0;i<(n-l);i++) x[i]=invh*x2[i]+v[i]; x[n-l]=1.0; delete [] c2; delete [] en; delete [] x2; for (i = O; i <m; i ++) {// compute bl component by component int in = i * n; int ind = in + (nl); c2 [ind] = c [ind]; for (j = 0; j <(n-1); j ++) c2 [ind] + = c [in + j] * v [j]; // normalize the rows of the matrix cmat2 for (i = O; i <m; i ++) {// normalize line by line int idim = i * n; double ncnorme = 0; for (j = O; j <n; j ++) ncnorm + = c2 [idim + j] * c2 [idim + j]; ncnorme = sqrt (ncnorme); for (j = O; j <n; j ++) c2 [idim + j] / = ncormal; } // solving the dual problem cProjection2 (n, m, c2, en, lag); for (j = O; j <n; j ++) lag [j] / = in [nl]; // change of variables to get the solution for (i = 0; i <(nl); i ++) x2 [i] = - lag [i]; // on ax2 * = (A h) x2 [n-1] = 1.0-lag [n-1];

// compute dex x (x = (1lh) .P. (llsqrt (D)) A + Xu) double invh = 1.0 / x2 [nl]; for (i = 0; i <(nl); i ++) x [i] = invh * x2 [i] + v [i]; x [nl] = 1.0; delete [] c2; delete [] in; delete [] x2;

<Desc / Clms Page number 41>

Annex 1 General formulas:

<Desc / Clms Page number 42>

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.

Annex 2 Constraint formulas:

I. dan al) (a (Ht) -a, (H)) + 2x (v, (Hr) -v, (H)) dt0 a2) {CaxXa # Ba} b1) (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}

<Desc / Clms Page number 43>

Appendix 3 Modeling formulas:

<Desc / Clms Page number 44>

<Desc / Clms Page number 45>

Claims (24)

1. A computer modeling device for working on input data (20), comprising on the one hand environment data, and on the other hand object data representing a plurality of objects defined by magnitudes. kinematics, the device comprising a calculation module (26), able to implement a set of algorithms at each time step, from previous values of the object data and environment data, the game of algorithms corresponding to the resolution of a set of equations for a given contact between at least two objects, characterized in that the calculation module is capable, at each time step, and for at least two objects (R, Q) in contact, to implement the resolution of a first equation (3III-b3, Eql) representing the minimization of a Euclidean norm applied to a first predetermined expression having an unknown depending on a dynamic quantity of each object, this first predetermined expression corresponding to a behavior of the objects without constraints, and said resolution of the first equation being effected under constraint of a second equation (3-III-b3, Eq2) representing at least one constraint condition on the unknown of the first equation and in that the set of algorithms corresponding to the resolution of the set of equations comprises an algorithm for projecting a point on a polyhedral cone. 2. Device according to claim 1, characterized in that said first equation (3-III-b3, Eql) representing the minimization of a Euclidean norm is derived from a predetermined primary equation (3-III-a2, Eql) representing the minimization of a primary expression having for unknown a kinematic magnitude (Xg, 3-III-a2) of each object. 3. Device according to one of claims 1 and 2, characterized in that said second equation (3-III-b3, Eq2) is derived from a secondary equation (3-III-a2, Eq2) adapted to <Desc / Clms Page number 47>  apply to the kinematic magnitude (Xg, 3-III-a2) of each object specific constraints (Ck, 3-III-a2, Eq2) respecting specific conditions (Bk, 3-III-a2, Eq2). 4. Device according to claim 2, characterized in that said predetermined primary equation (3-III-a2, Eql) is able to constrain the kinematic magnitude (Xg, 3-III-a2) of each object to get as close as possible the particular kinematic magnitude (Xk, 3-III-a2, Eql) that each object would have if it were not subject to constraints. 5. Device according to one of claims 3 and 4, characterized in that said kinematic magnitude comprises the desired strain acceleration (Xa, 3-I-a2), said particular kinematic magnitude comprises an acceleration without constraint (Au, 3- I-a2, Eql), said specific constraints comprise first geometric constraints (Ca, 3-Ia2, Eq2). 6. Device according to one of claims 3 and 4, characterized in that said kinematic magnitude comprises the desired strain rate (Vv, 3-I-b2), said particular kinematic magnitude comprises the speed without constraint (Vu, 3-I -b2, Eql), said specific constraints comprise second geometric constraints (Cv, 3-I-b2, Eq2). 7. Device according to one of claims 3 and 4, characterized in that said kinematic magnitude comprises the desired speed after collision (Xe, 3-II-2), said particular kinematic magnitude comprises the speed before collision (V, 3- II-2, Eql), said specific constraints include the constraints on the desired speed after collision (Cv, 3-II-2, Eq2). 8. Device according to one of claims 1 to 6, characterized in that said set of algorithms corresponding to the resolution of a set of equations, consisting of the first equation and the second equation, comprises an algorithm of <Desc / Clms Page number 48>  resolution of differential equations adapted to be applied in correspondence with the projection algorithm of a point on a polyhedral cone (104). 9. Device according to claim 8, characterized in that said algorithm for solving differential equations is the Runge-Kutta algorithm. 10. Device according to one of claims 1 to 8, characterized in that said projection algorithm of a point on a polyhedral cone is the Wilhelmsen algorithm. 11. Device according to one of claims 1 to 10, characterized in that said first equation represents a matrix equation describing a plurality of linear equations. 12. A method of computer modeling, from input data (20), comprising on the one hand environment data, on the other hand object data representing a plurality of objects defined by kinematic quantities, the method implementing iteratively a set of algorithms corresponding to the resolution of a set of equations for a given contact between at least two objects, and comprising the following steps: a. constructing a first equation (3-III-b3, Eql) representing the minimization of a Euclidean norm applied to a first predetermined expression having an unknown depending on a dynamic quantity of each object, this first predetermined expression corresponding to a behavior of the objects out of constraints, and a second equation (3-III-b3, Eq2) representing at least one constraint condition on the unknown of the first equation. b. implementing said set of algorithms, comprising an algorithm for projecting a point on a polyhedral cone, for solving said first equation (3-III-b3, <Desc / Clms Page number 49>  Eql), said resolution of the first equation being under stress of said second equation (3-III-b3, Eq2). 13. Computer modeling method according to claim 12, characterized in that step a. comprises deriving said first equation (3-III-b3, Eql) from a predetermined primary equation (3-III-a2, Eq1) representing the minimization of a primary expression having a kinematic magnitude as unknown (Xg, 3 -III-a2) of each object. 14. Computer modeling method according to one of claims 12 and 13, characterized in step a. comprises deriving said second equation (3-III-b3, Eq2) from a secondary equation (3-III-a2, Eq2) adapted to apply to the kinematic magnitude (Xg, 3-III-a2) of each object specific constraints (Ck, 3-III-a2, Eq2) respecting certain conditions (Bk, 3-III-a2, Eq2). 15. Computer modeling method according to claim 13, characterized in that the predetermined primary equation (3-III-a2, Eql) of step a. comprises of constraining the kinematic magnitude (Xg, 3-III-a2) of each object to be as close as possible to the particular kinematic magnitude (Xk, 3-III-a2, Eql) that would have each object if n was not subject to constraints. 16. Computer modeling method according to one of claims 14t 15caractérisé in that said kinematic magnitude of step a. comprises the sought-after constraint acceleration (Xa, 3-I-a2), said particular kinematic magnitude comprises unconstrained acceleration (Au, 3-I-a2, Eql), said specific constraints include first geometric constraints (Ca, 3- I-a2, Eq2). 17. Computer modeling method according to claims 14 and 15, characterized in that said kinematic magnitude of step a. includes the desired strain velocity (Vv, 3-I-b2), said particular kinematic magnitude <Desc / Clms Page number 50>  comprises the unconstrained speed (Vu, 3-I-b2, Eql), said specific constraints comprise second geometrical constraints (Cv, 3-I-b2, Eq2).  18. Computer modeling process according to claims 14 and 15, characterized in that said cinematic magnitude of step a. includes the desired post-collision velocity (Xc, 3-11-2), said particular kinematic magnitude includes the pre-collision velocity (V-, 3-11-2, Eql), said specific constraints include velocity constraints after collision sought (Cv '3-11-2, Eq2).  19. Computer modeling method according to one of claims 12 to 17, characterized in that step b. comprises implementing an algorithm for solving differential equations adapted to be applied in correspondence with the projection algorithm of a point on a polyhedral cone (104).  20. The computer modeling method according to claim 19, characterized in that said algorithm for solving differential equations is the Runge-Kutta algorithm.  21. Computer modeling method according to one of claims 12 to 20, characterized in that said algorithm for projecting a point on a polyhedral cone is the Wilhelmsen algorithm.  22. Software product, comprising the calculation module of the computer modeling device claimed in one of claims 1 to 11.  A software product, comprising the functions used in the computer modeling method claimed in one of claims 12 to 21.  24. A modeling software, comprising the software product claimed in one of claims 22 and 23.