Patent Publication Number: US-10776532-B2

Title: Modified effective mass for parallel rigid body simulation

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application is related to concurrently filed patent application Ser. No. 13/774,906, entitled “PARALLEL LINEAR COMPLEMENTARITY SOLVER FOR RIGID BODY DYNAMICS,” having a filing date of Feb. 22, 2013, which is herein incorporated by reference in its entirety. 
     BACKGROUND 
     Rigid body dynamics is widely used in applications ranging from movies to engineering to video games. Piles of objects are particularly common, because ultimately, gravity pulls all rigid bodies to the ground. Some of the most visually interesting simulations involve destruction, such as projectile impacts and explosions, and these can generate large piles of debris. In mechanical engineering some of the most computationally challenging problems involve simulating interaction with large resting systems of soil particles or rocks. Piles require stable simulation of static friction, dynamic friction and resting contact, which present many challenges. 
     In large-scale or real-time simulations, the time available for computation can be small compared to the number of rigid body contacts. In these cases iterative methods are used to solve rigid body dynamics. However, iteration is often terminated well before convergence to maximize the rigid body count. This results in residual energy, which can cause artifacts when its distribution over the contacts varies too much from frame to frame. The effect of this can be seen when bodies coming to rest oscillate (jitter) or move around each other non-physically (swim). 
     A commonly used iterative algorithm for simulating interactions is serial projected Gauss-Seidel (PGS), which solves contacts in sequence. When the iteration is terminated, the last contact solved has no error and the other contacts have non-zero error, resulting in an uneven distribution of the residual energy. On single threaded implementations the distribution of this error can be made consistent from frame to frame by, for example, by ensuring that collisions are detected in the same order each frame and that addition and deletion of bodies do not disrupt the order. 
     Unfortunately, things are not as straightforward for parallel implementations. PGS has limited parallelism, as updates to bodies having multiple contacts must be serialized to ensure they are not lost and the connectivity between bodies can be complex, especially in piles. It is very challenging to avoid jittering and swimming with parallel PGS, and due to serialization its performance does not scale well as more threads are added. 
     SUMMARY 
     A computer implemented method and system for a parallel iterative rigid body solver that simulates complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts. The systems and methods deal with resolution of contacts after contacts have been detected. The parallel iterative solver and computer-implemented method for implementation distributes the residual among the contacts consistently from frame to frame eliminating swimming and visible jittering. This improves simulation quality in use cases where the iteration is stopped well before convergence without introducing swimming or visible jittering. It utilizes one thread per contact without requiring serialization, resulting in a low cost per iteration and providing excellent parallel scaling. In one embodiment, the method includes splitting the mass of each body between its contacts before solving the contacts independently. In another embodiment, the method splits the bodies so that each sub-body receives a block of contacts with the corresponding portion of the mass in order to accelerate convergence. The system solver and method provides for correct momentum propagation, static friction, dynamic friction and stable resting contact. 
     In one embodiment, a computer implemented method for solving linear complementarity problems for rigid body simulation is disclosed. The method includes providing, receiving and/or determining one or more contact constraints affecting an original object having an original mass. The method includes splitting the original object by a total number of the contact constraints into a plurality of sub-bodies. The method includes assigning a contact constraint to a corresponding sub-body. The method further includes solving contact constraints in isolation. The method also includes adding fixed joints between the sub-bodies, and enforcing them exactly, such that positions and/or velocities of each sub-body are identical. 
     In another embodiment, another method for solving linear complementarity problems for rigid body simulation or solving contact constraints for systems of rigid bodies is disclosed. The computer implemented method includes determining one or more patches associated with an original object having an original mass. The patch comprises a pair of colliding objects including the original object and the patch defines a plurality of contact constraints affecting the pair of colliding objects. The method includes splitting the original object by a total number of the patches into a plurality of sub-bodies. Each sub-body is assigned to a patch and comprises a sub-body mass defined by the original mass divided by the total number of the patches. The method also includes solving patches in parallel, such that for each patch, in sequential fashion contact constraints are solved approximately. The method also includes adding fixed joints between the sub-bodies, and enforcing them exactly, such that positions and/or velocities of each sub-body are identical. 
     In still another embodiment, a linear complementarity solver for rigid body simulation is disclosed. The solver includes a patch determining module for determining one or more patches associated with an original object having an original mass, wherein a patch comprises a pair of colliding objects including the original object and the patch defines a plurality of contact constraints affecting the pair of colliding objects. The solver includes an object splitter for splitting the original object by a total number of the patches into a plurality of sub-bodies. Each sub-body is assigned to a patch and comprises a sub-body mass defined by the original mass divided by the total number of patches. The solver also includes a contact constraint solver for solving patches in parallel, such that for each patch, in sequential fashion contact constraints are solved approximately. The solver also includes a fixed joint solver for enforcing fixed joint constraints exactly, such that positions and/or velocities of each sub-body are identical. 
     In another embodiment, a computer implemented method for solving linear complementarity problems for rigid body simulation is disclosed. The iterative method includes: determining a plurality of modified effective masses for a plurality of contacts between a plurality of bodies, wherein each modified effective mass term is based on a corresponding number of contacts. The method includes determining a plurality of relative velocities based on a plurality of body velocities determined from a last iteration or an initial state. The method further includes determining a plurality of impulse corrections based on the plurality of modified effective masses and the plurality of relative velocities in order to negate influences of the plurality of relative velocities. The method also includes determining a plurality of updated impulses based on the plurality of impulse corrections. The method further includes applying the plurality of updated impulses to the plurality of bodies based on a plurality of original masses of the plurality of bodies and the plurality of body velocities determined from the last iteration or initial state to determine a plurality of updated velocities of the plurality of bodies. 
     In still another embodiment, a computer implemented method for solving linear complementarity problems for rigid body simulation is disclosed. The iterative method includes in a projected iteration for determining a plurality of modified effective masses for a plurality of contacts between a plurality of bodies, dividing each body mass by a corresponding number of contacts affecting said body mass. 
     In still another embodiment, a computer system is disclosed that includes a processor and memory coupled to the processor having stored therein instructions that, if executed by the computer system, cause the computer system to execute a method. The method includes determining a plurality of modified effective masses for a plurality of contacts between a plurality of bodies, wherein each modified effective mass term is based on a corresponding number of contacts. The method includes determining a plurality of relative velocities based on a plurality of body velocities determined from a last iteration or an initial state. The method further includes determining a plurality of impulse corrections based on the plurality of modified effective masses and the plurality of relative velocities in order to negate influences of the plurality of relative velocities. The method also includes determining a plurality of updated impulses based on the plurality of impulse corrections. The method further includes applying the plurality of updated impulses to the plurality of bodies based on a plurality of original masses of the plurality of bodies and the plurality of body velocities determined from the last iteration or initial state to determine a plurality of updated velocities of the plurality of bodies. 
     These and other objects and advantages of the various embodiments of the present disclosure will be recognized by those of ordinary skill in the art after reading the following detailed description of the embodiments that are illustrated in the various drawing figures. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The accompanying drawings, which are incorporated in and form a part of this specification and in which like numerals depict like elements, illustrate embodiments of the present disclosure and, together with the description, serve to explain the principles of the disclosure. 
         FIG. 1  depicts a block diagram of an exemplary computer system suitable for implementing the present methods, in accordance with one embodiment of the present disclosure. 
         FIG. 2  is a block diagram of a exemplary system for performing rigid body simulation, in accordance with one embodiment of the present disclosure. 
         FIG. 3  is a block diagram of a linear complementarity, parallel iterative rigid body solver that implements a scalar method for simulating complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts, in accordance with one embodiment of the present disclosure. 
         FIG. 4  is a flow diagram illustrating a scalar method for solving linear complementarity iterative rigid body dynamics, in accordance with one embodiment of the present disclosure. 
         FIG. 5  is a diagram illustrating coloring for a parallel projected Gauss-Seidel solution for rigid body dynamics. 
         FIG. 6  is a diagram illustrating the splitting of the mass of an original object into one or more sub-bodies, in accordance with one embodiment of the present disclosure. 
         FIG. 7  is a block diagram of a linear complementarity, parallel iterative rigid body solver that implements a block method for simulating complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts, in accordance with one embodiment of the present disclosure. 
         FIG. 8  is a flow diagram illustrating a block method for solving linear complementarity iterative rigid body dynamics, in accordance with one embodiment of the present disclosure. 
         FIG. 9  is a block diagram of a linear complementarity, parallel iterative rigid body solver that implements a scalar method for simulating complex systems coming to rest using a modified effective mass that divides each body mass term in the effective mass by the number of contacts acting on the body, in accordance with one embodiment of the present disclosure. 
         FIG. 10  is a flow diagram illustrating a scalar method for solving linear complementarity iterative rigid body dynamics using a modified effective mass that divides each body mass term in the effective mass by the number of contacts acting on the body, in accordance with one embodiment of the present disclosure. 
         FIG. 11  is a diagram showing a system of three bodies with four contacts to illustrate partitioning of contacts by body pairs, in accordance with one embodiment of the present disclosure. 
     
    
    
     DETAILED DESCRIPTION 
     Reference will now be made in detail to the various embodiments of the present disclosure, examples of which are illustrated in the accompanying drawings. While described in conjunction with these embodiments, it will be understood that they are not intended to limit the disclosure to these embodiments. On the contrary, the disclosure is intended to cover alternatives, modifications and equivalents, which may be included within the spirit and scope of the disclosure as defined by the appended claims. Furthermore, in the following detailed description of the present disclosure, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure. However, it will be understood that the present disclosure may be practiced without these specific details. In other instances, well-known methods, procedures, components, and circuits have not been described in detail so as not to unnecessarily obscure aspects of the present disclosure. 
     Some portions of the detailed descriptions that follow are presented in terms of procedures, logic blocks, processing, and other symbolic representations of operations on data bits within a computer memory. These descriptions and representations are the means used by those skilled in the data processing arts to most effectively convey the substance of their work to others skilled in the art. In the present application, a procedure, logic block, process, or the like, is conceived to be a self-consistent sequence of steps or instructions leading to a desired result. The steps are those utilizing physical manipulations of physical quantities. Usually, although not necessarily, these quantities take the form of electrical or magnetic signals capable of being stored, transferred, combined, compared, and otherwise manipulated in a computer system. It has proven convenient at times, principally for reasons of common usage, to refer to these signals as transactions, bits, values, elements, symbols, characters, samples, pixels, or the like. 
     It should be borne in mind, however, that all of these and similar terms are to be associated with the appropriate physical quantities and are merely convenient labels applied to these quantities. Unless specifically stated otherwise as apparent from the following discussions, it is appreciated that throughout the present disclosure, discussions utilizing terms such as “splitting,” “solving,” “enforcing,” “assigning,” “determining,” or the like, refer to actions and processes (e.g., flowcharts  400  and  800  of  FIGS. 4 and 8 , respectively) of a computer system or similar electronic computing device or processor (e.g., system  100 ,  200 ,  300 , and  700  of  FIGS. 1, 3 and 7  respectively). The computer system or similar electronic computing device manipulates and transforms data represented as physical (electronic) quantities within the computer system memories, registers or other such information storage, transmission or display devices. 
       FIGS. 4 and 8  are flowcharts of examples of computer-implemented methods for processing data according to embodiments of the present invention. Although specific steps are disclosed in the flowcharts, such steps are exemplary. That is, embodiments of the present invention are well-suited to performing various other steps or variations of the steps recited in the flowcharts. 
     Embodiments of the present invention described herein are discussed within the context of hardware-based components configured for monitoring and executing instructions. That is, embodiments of the present invention are implemented within hardware devices of a micro-architecture, and are configured for monitoring for critical stall conditions and performing appropriate clock-gating for purposes of power management. 
     Other embodiments described herein may be discussed in the general context of computer-executable instructions residing on some form of computer-readable storage medium, such as program modules, executed by one or more computers or other devices. By way of example, and not limitation, computer-readable storage media may comprise non-transitory computer storage media and communication media. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform particular tasks or implement particular abstract data types. The functionality of the program modules may be combined or distributed as desired in various embodiments. 
     Computer storage media includes volatile and nonvolatile, removable and non-removable media implemented in any method or technology for storage of information such as computer-readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, random access memory (RAM), read only memory (ROM), electrically erasable programmable ROM (EEPROM), flash memory or other memory technology, compact disk ROM (CD-ROM), digital versatile disks (DVDs) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to store the desired information and that can accessed to retrieve that information. 
     Communication media can embody computer-executable instructions, data structures, and program modules, and includes any information delivery media. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, radio frequency (RF), infrared and other wireless media. Combinations of any of the above can also be included within the scope of computer-readable media. 
       FIG. 1  is a block diagram of an example of a computing system  100  capable of implementing embodiments of the present disclosure. Computing system  10  broadly represents any single or multi-processor computing device or system capable of executing computer-readable instructions. Examples of computing system  100  include, without limitation, workstations, laptops, client-side terminals, servers, distributed computing systems, handheld devices, or any other computing system or device. In its most basic configuration, computing system  100  may include at least one processor  110  and a system memory  140 . 
     Both the central processing unit (CPU)  110  and the graphics processing unit (GPU)  120  are coupled to memory  140 . System memory  140  generally represents any type or form of volatile or non-volatile storage device or medium capable of storing data and/or other computer-readable instructions. Examples of system memory  140  include, without limitation, RAM, ROM, flash memory, or any other suitable memory device. In the example of  FIG. 1 , memory  140  is a shared memory, whereby the memory stores instructions and data for both the CPU  110  and the GPU  120 . Alternatively, there may be separate memories dedicated to the CPU  110  and the GPU  120 , respectively. The memory can include a frame buffer for storing pixel data drives a display screen  130 . 
     The system  100  includes a user interface  160  that, in one implementation, includes an on-screen cursor control device. The user interface may include a keyboard, a mouse, and/or a touch screen device (a touchpad). 
     CPU  110  and/or GPU  120  generally represent any type or form of processing unit capable of processing data or interpreting and executing instructions. In certain embodiments, processors  110  and/or  120  may receive instructions from a software application or hardware module. These instructions may cause processors  110  and/or  120  to perform the functions of one or more of the example embodiments described and/or illustrated herein. For example, processors  110  and/or  120  may perform and/or be a means for performing, either alone or in combination with other elements, one or more of the monitoring, determining, gating, and detecting, or the like described herein. Processors  110  and/or  120  may also perform and/or be a means for performing any other steps, methods, or processes described and/or illustrated herein. 
     In some embodiments, the computer-readable medium containing a computer program may be loaded into computing system  100 . All or a portion of the computer program stored on the computer-readable medium may then be stored in system memory  140  and/or various portions of storage devices. When executed by processors  110  and/or  120 , a computer program loaded into computing system  100  may cause processor  110  and/or  120  to perform and/or be a means for performing the functions of the example embodiments described and/or illustrated herein. Additionally or alternatively, the example embodiments described and/or illustrated herein may be implemented in firmware and/or hardware. 
       FIG. 2  is a block diagram of a rigid body physics system or engine  200  for performing rigid body simulation, in accordance with one embodiment of the present disclosure. In particular, system  200  calculates the motion of a set of bodies such that when the bodies are displayed on a screen, they appear to obey the laws of physics and don&#39;t pass through each other. 
     As shown in  FIG. 2 , at each displayed frame, a physics system  200  executes a sequence of stages (e.g., a pipeline) to calculate new the positions and orientations of the bodies. For instance, the pipeline of system  200  includes a collision detector  210  that takes as inputs shapes and poses of objects, and outputs contacts that describe the interactions of contacting pairs of objects. Each contact can be represented as a constraint that constrains the relative motion of pairs of contacting bodies. A constraint generator  220  takes each contact and converts it to a constraint representation that can be solved by a solver. That is, generator  220  outputs constraints. Block  230  takes the initial velocities of the objects and applies external forces or impulses affecting the objects. A vector of velocities is output from block  230 . Solver  250  takes the set of constraints that the velocities, positions, and/or orientations of bodies must satisfy and calculates a set of corrections such that the bodies will satisfy the constraints once the corrections are applied to the bodies. In one embodiment, solver  250  calculates position corrections and velocity corrections. In another embodiment, solver  250  calculates only velocity corrections. In one embodiment, velocities are corrected with impulses. In another embodiment, velocities are corrected with forces. 
     In embodiments of the present invention, solver  250  solves a linear complementarity problem. The solver  250  is configured to process in parallel operations needed to simulate a very large number of rigid bodies at real-time frame rates (e.g., 60 Hz). As such, the desired frame rate is achieved using the parallel linear complementarity solver  250 , which is running on a parallel processor or processors (e.g., GPU, multi-core CPU, multi-threaded processor, etc.). 
     Scalar Method and System for a Rigid Body Solver 
       FIG. 3  is a block diagram of a linear complementarity, parallel iterative rigid body solver  300  for simulating complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts, in accordance with one embodiment of the present disclosure. In one embodiment, solver  300  performs the functions of solver  250  of  FIG. 2  within the rigid body physics system or engine  200 . 
     Solver  300  includes an object mass splitter  310 , for splitting an original object into a plurality of sub-bodies. The bodies are not split spatially, as performed by traditional finite element methods, but rather the mass is split so that all the sub-bodies have the same spatial extents, center of mass and initial velocity. In one embodiment, the splitter  310  also splits the original mass of the original object, such that each sub-body has the same mass distribution. More particularly, solver  300  is configured to split the original object by a total number of contact constraints that are affecting the original object. 
     Solver  300  also includes a contact constraint assigner  315 . Each contact constraint is assigned to a corresponding sub-body. In that manner, by splitting the original object into one or more sub-bodies, multiple computing threads are used to compute solve contact constraint in parallel. That is, contact constraint solver  320  is configured to solve contact constraints approximately for each sub-body. Further, each contact constraint is solved with respect to its corresponding sub-body in isolation. As such, contact constraint assigner  315  is configured to use a relative velocity at the contact point as defined by the contact constraint and calculate an impulse applied at that point to eliminate the relative velocity. 
     Solver  300  also includes a fixed joint solver  325  that enforces fixed joint constraints exactly, such that positions and orientations of each sub-body are identical. With multiple sub-bodies, applying an impulse at one contact point for one constraint can affect the velocity at many other contact points. That is, after the sub-bodies have been computed, the impulse results are combined to find the force on the original body. The combination of sub-bodies is modeled as a fixed joint that forces all the sub-bodies of the original object to have the same position, orientation, and velocity. That is, in one embodiment, by splitting the original object into a plurality of sub-bodies, the fixed joint is enforced or solved at each iteration by averaging the sub-body velocities, positions and/or orientations and applying the average velocity to each of the sub-bodies. The average velocity is used as the initial velocity of the sub-bodies in the next iteration to solve contact constraints and fixed joint constraints 
       FIG. 4  is a flow diagram  400  illustrating a computer implemented scalar method for solving linear complementarity iterative rigid body dynamics, in accordance with one embodiment of the present disclosure. In another embodiment, flow diagram  400  is implemented within a computer system including a processor and memory coupled to the processor and having stored therein instructions that, if executed by the computer system causes the system to execute a scalar method for solving linear complementarity iterative rigid body dynamics. In still another embodiment, instructions for performing a method are stored on a non-transitory computer-readable storage medium having computer-executable instructions for causing a computer system to perform a scalar method for solving linear complementarity iterative rigid body dynamics. In one embodiment, flow diagram  400  is implemented by solver  300  of  FIG. 3 , and/or solver  250  of  FIG. 2 . More particularly, flow diagram  400  is implemented to take a set of constraints that the velocities, positions, and/or orientations of objects must satisfy and calculate a set of forces and/or impulses such that the objects will satisfy the constraints once the forces and/or impulses are applied to the objects. 
     The method of flow diagram  400  is implemented after a set of constraints are defined, wherein contact constraints affect one or more objects or bodies. For instance, constraint generator  220  of  FIG. 2  is able to determine the constraints that affect one or more objects or bodies. Once the constraints are defined, the scalar method of flow diagram  400  is implemented to determine the positions, orientations, and/or velocities of the one or more objects or bodies. 
     At  410 , splitting of an original object is performed, wherein the original object is associated with an original mass. The original object is split into a plurality of sub-bodies by splitting the original object by a total number of contact constraints affecting the original object. Additionally, the sub-bodies all have the same spatial extent, center of mass, and initial velocity. 
     In one embodiment, each of the sub-bodies have the same mass distribution. That is, each of the sub-bodies have an equally divided mass that is defined by dividing the original mass by the total number of contact constraints. In another embodiment, the masses of each sub-body may be assigned to varying values, such that the masses of two sub-bodies may be different. 
     At  420 , assigning each contact constraint to a corresponding sub-body is performed. By splitting into sub-bodies, different computing threads are used to compute different contact constraints in parallel. As such, parallel computations may be performed for solving contact constraints and solving fixed joint constraints, as provided in  430  and  440  below. 
     At  430 , solving contact constraints is performed. In addition, the solving of contact constraints is performed in isolation for each sub-body. The affect one sub-body has on another sub-body is considered when solving for fixed joint constraints at  440 . 
     In one embodiment, the velocity Signorini condition is satisfied when solving the contact constraints. That is, a contact constraint is solved between a pair of sub-bodies, such that impulses are applied in pairs to each sub-body, that have equal magnitudes but opposite directions in order to satisfy Newton&#39;s third law of motion. In that manner, the magnitude of the pair of impulses is solved. Applying the pair of impulses changes the relative velocity of the pair of sub-bodies at the contact point. The corresponding contact constraint imposes requirements on both the magnitude and the resulting relative velocity, as follows: the resulting relative velocity must be zero or positive (separating); and the magnitude of the impulse must be zero or positive (repulsive), wherein the velocity must be zero when the impulse is repulsive, and the impulse must be zero when the relative velocity is separating. 
     At  440 , fixed joint constraints are solved exactly. More particularly, positions and/or velocities of each sub-body are solved such that they are substantially identical. For instance, in one embodiment, the enforcing of fixed joints is performed by averaging a plurality of sub-body velocities corresponding to the plurality of sub-bodies. As such, an average velocity is defined. Further, the average velocity is applied to each of the plurality of sub-bodies. 
     In one embodiment, the operations in  430  and  440  are iteratively performed. That is, for each iteration, given a starting average velocity from the previous iteration, a new average velocity and new position and/or new orientations are determined for a corresponding sub-body. These new values again correspond to an average velocity and position and/or orientation for the original object that can be used finally or for the next iteration. 
     In another embodiment, a frictional force is applied to determine the sub-body velocity at each iteration. Embodiments of the present disclosure support any means for applying friction. Embodiments of the present disclosure apply a coulomb friction model, though other friction models are supported. For instance, friction may be considered by applying a friction model, such that the more force an object is exerting on a surface, the harder it needs to be pushed to move it across the surface. As such, heavier objects are harder to push than lighter objects. In addition, the transition from sticking (e.g., when the object does not move) to sliding (e.g., when the object first moves) is considered. 
     The solution for the scalar method for performing linear complementarity in parallel for an iterative rigid body solver is described below. As previously described, the scalar method simulates complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts. Portions described below are equally adaptable for use in the block method for performing linear complementarity in parallel for an iterative rigid body solver that is further described below in relation to  FIGS. 7 and 8 , in other embodiments. 
     Rigid Body Contact Dynamics 
     The rigid body contact model used is summarized. Consider a scene having n rigid bodies with positions xϵ   6n , external forces f e ϵ   6n  and masses Mϵ   6n×6n . Collision detection identifies m contacts between the rigid bodies, represented by constraints Φ(x)≥0 with Jacobian ∂Φ∂x=Jϵ   m×6n . For simplicity of description constraint stabilization is omitted, so the unknowns are the forces λϵ   m  necessary to satisfy the time derivative of the constraints. Friction is treated separately, in one embodiment. 
     Contacts must satisfy the velocity Signorini condition: i.e., forces must not be attractive (e.g., λ≥0), velocities must move the system out of penetration 
               (           ∂   Φ       ∂   x       ⁢   v     ≥   0     )     ,         
and a force should be applied at a contact only if that contact is not separating, which is expressed as λ≥0⊥Jv≥0. Putting this all together, the continuous model is described in the following equations.
 
 Mx″=J   T   λ+f   e   (1)
 
 x′=v   (2)
 
λ≥0 ⊥Jv≥ 0.  (3)
 
Time Stepping
 
     The model is discretized using a semi-implicit stepping scheme with time step h, and constraint impulse z=hλ, as described in the following equations:
 
 M ( v   new   −v   old )= J   T   z+hf   e   (4)
 
 x   new   −x=hv   new   (5)
 
 z≥ 0 ⊥Jv   new ≥0.  (6)
 
     The Signorini condition causes the discretized model for rigid body dynamics to be a linear complementarity (LCP) formulation rather than a linear system. As such, let x:=LCP(A, b) be defined in the following equation:
 
find  xϵ     m  such that, for  all i= 1  . . . m,  
 
 x   i ≥0 and ( Ax+b ) i ≥0 and
 
 x   i =0 or ( Ax+b ) i =0.  (7)
 
     As such, the discretized model is solved in the following equations:
 
 v   0   =v   old   +hM   −1   f   e   (8)
 
 q=Jv   0   (9)
 
 N=JM   −1   J   T   (10)
 
 z=LCP ( N,q )  (11)
 
     Two iterative methods can be used to solve the rigid body LCP: projected Gauss-Seidel (PGS) and projected Jacobi. Both these methods are described below. 
     Projected Guass-Seidel 
     The projected Guass-Seidel model solves each constraint individually in sequence. Let L, U, D be the lower triangle, upper triangle and diagonal of N respectively. Abstractly PGS generates the following iteration from an initial guess, z°, in the following equations: 
     
       
         
           
             
               
                 
                   
                     z 
                     
                       r 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       ( 
                       
                         
                           z 
                           r 
                         
                         - 
                         
                           
                             D 
                             
                               - 
                               1 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               q 
                               + 
                               
                                 Lz 
                                 
                                   r 
                                   + 
                                   1 
                                 
                               
                               + 
                               
                                 Uz 
                                 r 
                               
                             
                             ) 
                           
                         
                       
                       ) 
                     
                     + 
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       ( 
                       
                         x 
                         + 
                       
                       ) 
                     
                     i 
                   
                   ⁢ 
                   
                     = 
                     def 
                   
                   ⁢ 
                   
                     
                       max 
                       ⁡ 
                       
                         ( 
                         
                           0 
                           , 
                           
                             x 
                             i 
                           
                         
                         ) 
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     Forming N explicitly would take O(m 3 ) space and O(m 2 ) time, which is often infeasible as m can be very large (≥10000). The formation of N may be avoided by explicitly updating the body velocities at each iteration, maintaining the invariant v= M   −1 J T z∔v old , and calculating Jv when elements of Nz are required, which is shown in Algorithm 1, below: 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 1 PGS iteration optimized for N = JM −1  J T   
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                 1: for all contacts, i do 
               
               
                   
                 2:  z i =max(0, z i  − M ii   −1  (Jv) i   
               
               
                   
                 3:  v=M −1  J T z 
               
               
                   
                 4: end for 
               
               
                   
               
            
           
         
       
     
     Each contact i constrains two bodies, let b i1  be the index of the first body, and b i2  be the index of the second. Each row of J corresponds to a contact and the (six element wide) column blocks of J correspond to the bodies constrained by each row, so J is sparse with two groups of six contiguous elements per row. By taking advantage of this sparsity, the iteration cost becomes O(m) space and O(m) time, which is described in Algorithm 2, below: 
                             Algorithm 2 PGSiter(z, M, J, v)                                            1: for all contacts, i do           2:  t = z i             3:  z i =max(0, z i  − (M ii ) −1  (Jv) i             4.  for all j ∈ {1,2} do           5:   β = b i,j             6:   v β =v β  + M β   −1  (J iβ ) T  (z i  − t)           7:  end for           8: end for                    
Parallel Projected Gauss-Seidel
 
     Each body is often in contact with more than one other body, so if all the contacts are processed in parallel, then each body velocity will be updated by many contacts simultaneously, which in turn causes velocity updates to be lost, thereby adversely affecting convergence. One solution is contact coloring, which partitions (or colors) the contacts so that each body is referenced at most once in each partition. The PGS solver can then process each color sequentially, processing the contacts within each color in parallel. 
       FIG. 5  is a diagram illustrating coloring for a parallel projected Gauss-Seidel solution for rigid body dynamics. As shown, five blocks are shown in various contact with each other and a surface  540 . For instance, Block-A is in contact with surface  540  at R-0 and Block-D at G-3; Block-B is in contact with surface  540  at R-1, Block-D at B-4, and Block-E at G-5; Block-C is in contact with surface  540  at R-2, and Block-E at B-6; Block-D is in contact with Block-A at G-3, and Block-B at B-4; and Block-E is in contact with Block-B at G-5, and Block-C at B-6. Three colors are used to partition out the contacts so that each body is referenced at most in one partition, as follows: Red as “R,” Green as “G,” and Blue as “B.” As such, the constraints are processed by partition, such as, processing R-0, R-1, and R-2 first, then processing G-3 and G-5, and then processing B-4 and B-6 last. 
     The number of colors required depends on Δ, the maximum number of contacts per body over all the bodies. The minimum number of colors required is either Δ or Δ+1. A parallel algorithm is used to find the optimal coloring, which gives at most 2Δ−1 colors. In piles containing objects of very different sizes, or long thin objects like planks of wood, Δ can be quite high. In the example with planks of wood, there may be up to 10 contacts per body and the greedy coloring algorithm may produce 16 colors. Generally, the first few colors contain a large proportion of the contacts, but the following colors contain fewer and fewer contacts. As the colors must be processed sequentially, the time taken by each iteration is multiplied by the number of colors, reducing the speedup of parallel PGS versus serial PGS. As such, on architectures (like graphic processing units [GPUs]) with hundreds of cores, most of them are idle for most of the iteration, especially as back end colors are processed. In embodiments of the present disclosure, processors are used to process contact constraints for one or more sub-bodies, and as such this unused processing power is harnessed to reduce swimming and jitter. 
     Parallel Projected Jacobi 
     In the parallel projected Jacobi method, each thread takes a constraint and calculates the impulse required to satisfy the constraint completely. Formally, the iteration is calculated in the following equation:
 
 z   r+1 =( z   r   −D   −1 ( Nz   r   +q )) + ,  (14)
 
Mass Splitting
 
     In embodiments of the present disclosure, each body is split into sub-bodies so that each sub-body has one contact. In one implementation, the mass is split evenly between the sub-bodies. In another implementation, the mass is split unevenly. As such, each sub-body is solved independently, and then joined back together after having processed all the contacts and their contact constraints. 
     A common way in non-rigid simulation is to split bodies spatially, as is done when simulating deformables using a finite element method. However, this would cost time and memory and would lead to a complicated method. On the other hand, in embodiments of the present disclosure, a body is split into sub-bodies that occupy the same region of space as the original body. Each sub-body has the same mass distribution as the original body, in one embodiment. Hence, each sub-body moment of inertia is just a scalar fraction of the original moment of inertia. The body is split so that each sub-body has just one contact. Fixed joints are added to force all the sub-bodies of each body to have the same position and/or velocity. 
     Splitting the Bodies 
     The split system is constructed from the original system.  FIG. 6  is a diagram  600  illustrating the splitting of the mass of an original object into one or more sub-bodies, in accordance with one embodiment of the present disclosure. As such, original object  610  is split into sub-bodies  610 A and  610 B. A fixed joint constraint joins the sub-bodies back together. Also shown in  FIG. 6  is the original object  610  in contact with a non-moving block  640  and a moveable object, such as, a concrete block  630 . For instance, contact patch  650  represents the contact formed between original object  610  and non-moving block  640 . Contact patch  660  represents the contact formed between original object  610  and moving block  630 . Further, sub-bodies are formed and assigned to a contact or contact constraint. For instance, contact patch  650  represents the contact formed between sub-body  610 A and non-moving block  640 . Also, contact patch  660  represents the contact formed between sub-body  610 B and moving block  660 . 
     Continuing, let n β  be the number of constraints affecting original body β. We split body β into n β  sub-bodies, so the split system will have sub-bodies as described in Eqn. 15. 
     
       
         
           
             
               
                 
                   
                     n 
                     s 
                   
                   = 
                   
                     
                       ∑ 
                       
                         β 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                       n 
                       β 
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     The mass and inertia of body β is split evenly between its n B  sub-bodies, in one embodiment. As such, the inverse mass matrix of the sub-bodies of body β is described in Equation 16, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       W 
                       β 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 
                                   n 
                                   β 
                                 
                                 ⁢ 
                                 
                                   M 
                                   β 
                                   
                                     - 
                                     1 
                                   
                                 
                               
                             
                             
                               
                                   
                               
                             
                             
                               0 
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               0 
                             
                             
                               
                                   
                               
                             
                             
                               
                                 
                                   n 
                                   β 
                                 
                                 ⁢ 
                                 
                                   M 
                                   β 
                                   
                                     - 
                                     1 
                                   
                                 
                               
                             
                           
                         
                         ] 
                       
                       ∈ 
                       
                         
                           6 
                           ⁢ 
                           
                             n 
                             β 
                           
                           × 
                           6 
                           ⁢ 
                           
                             n 
                             β 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     The inverse mass matrix of the whole split system is described in Equation 17, as follows: 
     
       
         
           
             
               
                 
                   W 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               W 
                               1 
                             
                           
                           
                             
                                 
                             
                           
                           
                             0 
                           
                         
                         
                           
                             
                                 
                             
                           
                           
                             ⋱ 
                           
                           
                             
                                 
                             
                           
                         
                         
                           
                             0 
                           
                           
                             
                                 
                             
                           
                           
                             
                               W 
                               n 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     Each sub-body is just a portion of the mass of the original body, so the initial velocity of each sub-body of body β is set to the initial velocity of the original body β, as described in Equation 18, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       β 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 v 
                                 β 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 v 
                                 β 
                               
                             
                           
                         
                         ] 
                       
                       ∈ 
                       
                         
                           6 
                           ⁢ 
                           
                             n 
                             i 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
     
     The velocity vector for the whole split system is described in Equation 19, as follows: 
     
       
         
           
             
               
                 
                   
                     v 
                     S 
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               V 
                               1 
                             
                           
                         
                         
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               V 
                               n 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
           
         
       
     
     The unconstrained velocity, v S0  is constructed from v S  in the same way. 
     Joining the Bodies Back Together 
     A particular type of bilateral constraint is a fixed joint, a joint which forces the position and orientation of two bodies to be the same. Fixed joints are not commonly used, as it would be wasteful to calculate impulses to join two bodies when they could be replaced by a single body. However, during the derivation of the method, impulses will disappear in the final algorithm. Fixed joints are represented by the constraint function, as described in Equation 20, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       Φ 
                       FJ 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             
                               x 
                               1 
                             
                           
                         
                         
                           
                             
                               x 
                               2 
                             
                           
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         x 
                         1 
                       
                       - 
                       
                         x 
                         2 
                       
                     
                     = 
                     0 
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
           
         
       
     
     The Jacobian of the function is described in Equation 21, as follows: 
     
       
         
           
             
               
                 
                   
                     J 
                     FJ 
                   
                   = 
                   
                     
                       
                         ∂ 
                         
                           Φ 
                           FJ 
                         
                       
                       
                         ∂ 
                         x 
                       
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               I 
                             
                             
                               
                                 - 
                                 I 
                               
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     For each body β, fixed joints n β −1 are introduced to join the sub-bodies together. This gives a total number of fixed joints, m f , as expressed in Equation 22, as follows: 
     
       
         
           
             
               
                 
                   
                     m 
                     f 
                   
                   = 
                   
                     
                       ∑ 
                       
                         β 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             n 
                             β 
                           
                           - 
                           1 
                         
                         ) 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
     The Jacobian of the fixed joints of body β is described in Equation 23, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       F 
                       β 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                             
                               0 
                             
                             
                               … 
                             
                             
                               0 
                             
                           
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               0 
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                             
                               
                                   
                               
                             
                             
                               ⋮ 
                             
                           
                           
                             
                               ⋮ 
                             
                             
                               ⋮ 
                             
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               0 
                             
                           
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               0 
                             
                             
                               … 
                             
                             
                               0 
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                           
                         
                         ] 
                       
                       ∈ 
                       
                         
                           6 
                           ⁢ 
                           
                             n 
                             β 
                           
                           × 
                           6 
                           ⁢ 
                           
                             ( 
                             
                               
                                 n 
                                 β 
                               
                               - 
                               1 
                             
                             ) 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
           
         
       
     
     Also, the Jacobian of all the fixed joints in the split system is described in Equation 24, as follows: 
                   F   =       [           F   1                     0                       ⋱                       0                     F   n           ]     .             (   24   )               
Assigning the Contacts to the Split Bodies
 
     Contacts are assigned to sub-bodies on a one-to-one basis. That is, contacts in the set of n β  contacts acting on body β are assigned, such that each contact is assigned to a corresponding sub-body. The original Jacobian is partitioned into m×n blocks, such that J αβ ϵ^   I×6  represents the effect of contact α on body β, as described in Equation 25, as follows: 
     
       
         
           
             
               
                 
                   J 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               J 
                               11 
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               
                                 1 
                                 ⁢ 
                                 n 
                               
                             
                           
                         
                         
                           
                             ⋮ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               J 
                               
                                 m 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               mn 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   25 
                   ) 
                 
               
             
           
         
       
     
     Each J β  has n β  nonzero rows corresponding to the contacts that affect body β. Let r j  be the index of the jth nonzero row. In one implementation, the jth constraint of body β is assigned to the jth sub-body of body β, as expressed in Equation 26, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         C 
                         αβ 
                       
                       ] 
                     
                     
                       1 
                       , 
                       j 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             J 
                             αβ 
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               α 
                             
                             = 
                             
                               r 
                               j 
                             
                           
                         
                       
                       
                         
                           0 
                         
                         
                           
                             otherwise 
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   26 
                   ) 
                 
               
             
           
         
       
     
     The expression, [M αβ ]ij is used to represent element i, j of block α, β of any block matrix M. Later, the inverse map of Equation 26 is needed. That is, given a split contact Jacobian, C, the corresponding unsplit contact Jacobian, J needs to be found. This is accomplished using a couple more definitions. For instance, let s α,1  be the rank of contact α in its first body&#39;s constraint list, and s a,2  be the rank in its second. As such, the inverse map is now expressed in Equation 27, as follows: 
     
       
         
           
             
               
                 
                   
                     J 
                     αβ 
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               [ 
                               
                                 C 
                                 αβ 
                               
                               ] 
                             
                             
                               1 
                               , 
                               
                                 s 
                                 
                                   α 
                                   , 
                                   1 
                                 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   1 
                                 
                               
                             
                             = 
                             β 
                           
                         
                       
                       
                         
                           
                             
                               [ 
                               
                                 C 
                                 αβ 
                               
                               ] 
                             
                             
                               1 
                               , 
                               
                                 s 
                                 
                                   α 
                                   , 
                                   2 
                                 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   2 
                                 
                               
                             
                             = 
                             
                               β 
                               . 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   27 
                   ) 
                 
               
             
           
         
       
     
     Putting all the bodies back together, the contact Jacobian for the whole split system is expressed in Equation 28, as follows: 
     
       
         
           
             
               
                 
                   C 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               C 
                               11 
                             
                           
                           
                             … 
                           
                           
                             
                               C 
                               
                                 1 
                                 ⁢ 
                                 n 
                               
                             
                           
                         
                         
                           
                             ⋮ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               C 
                               
                                 m 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               C 
                               mn 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   28 
                   ) 
                 
               
             
           
         
       
     
     The impulse required to enforce the contact constraints on the split bodies (i.e., ignoring the fixed joints) is now expressed in Equation 29, as follows:
 
 LCP ( CWC   T   ,q ).  (29)
 
A Mixed LCP to Approximately Solve Contacts and Exactly Solve Fixed Joints
 
     A system is presented that contains contact constraints and bilateral joint constraints, or fixed joint constraints. In particular, the system is represented with a mixed LCP (MLCP) formulation, as described in Equation 30, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         
                           
                             
                               z 
                               C 
                             
                           
                         
                         
                           
                             
                               z 
                               F 
                             
                           
                         
                       
                       ] 
                     
                     = 
                     
                       MLCP 
                       ⁡ 
                       
                         ( 
                         
                           
                             [ 
                             
                               
                                 
                                   
                                     N 
                                     CC 
                                   
                                 
                                 
                                   
                                     N 
                                     CF 
                                   
                                 
                               
                               
                                 
                                   
                                     N 
                                     FC 
                                   
                                 
                                 
                                   
                                     N 
                                     CC 
                                   
                                 
                               
                             
                             ] 
                           
                           , 
                           
                             [ 
                             
                               
                                 
                                   
                                     q 
                                     C 
                                   
                                 
                               
                               
                                 
                                   
                                     q 
                                     F 
                                   
                                 
                               
                             
                             ] 
                           
                         
                         ) 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   30 
                   ) 
                 
               
             
           
         
       
     
     More particularly, the system containing the contact constraints and bilateral joint constraints is described in the following equations, as follows:
 
( N   CC   z   C   +N   CF   z   F   +q   C )≥0 ⊥z   C ≥0  (31)
 
( N   FC   z   C   +N   FF   z   F   +q   F )=0  (32)
 
     The impulse required to simultaneously satisfy the contacts on the split bodies and the joints that hold them together is expressed in Equation 33, as follows: 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             z 
                             C 
                           
                         
                       
                       
                         
                           
                             z 
                             F 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     MLCP 
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             [ 
                             
                               
                                 
                                   C 
                                 
                               
                               
                                 
                                   F 
                                 
                               
                             
                             ] 
                           
                           ⁢ 
                           
                             
                               W 
                               ⁡ 
                               
                                 [ 
                                 
                                   
                                     
                                       C 
                                     
                                   
                                   
                                     
                                       F 
                                     
                                   
                                 
                                 ] 
                               
                             
                             T 
                           
                         
                         , 
                         
                           [ 
                           
                             
                               
                                 q 
                               
                             
                             
                               
                                 0 
                               
                             
                           
                           ] 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   33 
                   ) 
                 
               
             
           
         
       
     
     This is useful because z C  is a solution of the original (i.e., unsplit) system. That is, Equation 34 shows the Theorem to be proved, as follows: 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             z 
                             C 
                           
                         
                       
                       
                         
                           
                             z 
                             F 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       
                         MLCP 
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 [ 
                                 
                                   
                                     
                                       C 
                                     
                                   
                                   
                                     
                                       F 
                                     
                                   
                                 
                                 ] 
                               
                               ⁢ 
                               
                                 
                                   W 
                                   ⁡ 
                                   
                                     [ 
                                     
                                       
                                         
                                           C 
                                         
                                       
                                       
                                         
                                           F 
                                         
                                       
                                     
                                     ] 
                                   
                                 
                                 T 
                               
                             
                             , 
                             
                               [ 
                               
                                 
                                   
                                     q 
                                   
                                 
                                 
                                   
                                     0 
                                   
                                 
                               
                               ] 
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         
 
                       
                       ⟹ 
                       
                         
 
                       
                       ⁢ 
                       
                         z 
                         C 
                       
                     
                     = 
                     
                       
                         LCP 
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 JM 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 J 
                                 T 
                               
                             
                             , 
                             q 
                           
                           ) 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   34 
                   ) 
                 
               
             
           
         
       
     
     More particularly, suppose that 
             [           z   C               z   F           ]         
is a solution of the split system, resulting in Equation 35, as follows:
 
 C   v   S   +q≥ 0 ⊥z   C ≥0  (35)
 
     In Equation 35, the term v S  is the velocity as expressed in Equation 36, as follows:
 
 v   S   =W ( C   T   z   C   +F   T   z   F ).  (36)
 
     Let e(v S ) be the vector consisting of the velocity of the first sub-body for each body. It is shown that z C  is a solution of the unsplit system. This is accomplished in two stages. First, it is shown that if impulse z C  is applied to the unsplit system, then the velocity of the bodies is e(v S ). Then, it is shown that the velocities e(v S ) satisfy the constraints of the unsplit system. Then, applying impulse z C  to the unsplit system will give the following velocity for each body β, in Equation 37, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 M 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 J 
                                 T 
                               
                               ⁢ 
                               
                                 z 
                                 C 
                               
                             
                             ) 
                           
                           β 
                         
                         = 
                           
                         ⁢ 
                         
                           
                             W 
                             β 
                           
                           ⁢ 
                           
                             n 
                             β 
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 k 
                                 = 
                                 1 
                               
                               
                                 n 
                                 b 
                               
                             
                             ⁢ 
                             
                               
                                 ( 
                                 
                                   
                                     C 
                                     β 
                                     T 
                                   
                                   ⁢ 
                                   
                                     z 
                                     C 
                                   
                                 
                                 ) 
                               
                               k 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               W 
                               β 
                             
                             ⁡ 
                             
                               ( 
                               
                                 I 
                                 - 
                                 
                                   
                                     
                                       
                                         F 
                                         β 
                                         T 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             F 
                                             β 
                                           
                                           ⁢ 
                                           
                                             W 
                                             β 
                                           
                                           ⁢ 
                                           
                                             F 
                                             β 
                                             T 
                                           
                                         
                                         ) 
                                       
                                     
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁢ 
                                   
                                     F 
                                     β 
                                   
                                   ⁢ 
                                   
                                     W 
                                     β 
                                   
                                 
                               
                               ) 
                             
                           
                           ⁢ 
                           
                             C 
                             β 
                             T 
                           
                           ⁢ 
                           
                             z 
                             C 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               [ 
                               
                                 
                                   W 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         
                                           C 
                                           T 
                                         
                                         ⁢ 
                                         
                                           z 
                                           C 
                                         
                                       
                                       + 
                                       
                                         
                                           F 
                                           T 
                                         
                                         ⁢ 
                                         
                                           z 
                                           F 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                                 β 
                               
                               ] 
                             
                             1 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   37 
                   ) 
                 
               
             
           
         
       
     
     This results in Equation 38, as follows:
 
 M   −1   J   T   z   C   =e ( v   S ).  (38)
 
     Next, the velocities of the split system is shown to satisfy the constraints of the unsplit system in Equation 39, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               Cv 
                               S 
                             
                             ) 
                           
                           α 
                         
                         = 
                           
                         ⁢ 
                         
                           
                             ∑ 
                             
                               β 
                               = 
                               1 
                             
                             n 
                           
                           ⁢ 
                           
                             
                               C 
                               αβ 
                             
                             ⁢ 
                             
                               V 
                               β 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             ∑ 
                             
                               β 
                               = 
                               1 
                             
                             n 
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 j 
                                 = 
                                 1 
                               
                               
                                 n 
                                 β 
                               
                             
                             ⁢ 
                             
                               
                                 
                                   [ 
                                   
                                     C 
                                     αβ 
                                   
                                   ] 
                                 
                                 
                                   1 
                                   , 
                                   j 
                                 
                               
                               ⁢ 
                               
                                 v 
                                 β 
                               
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               ( 
                               
                                 Je 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     v 
                                     S 
                                   
                                   ) 
                                 
                               
                               ) 
                             
                             α 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   39 
                   ) 
                 
               
             
           
         
       
     
     This results in the theorem in Equation 40, as follows:
 
 Cv   S   +q≥ 0 ⊥z   C ≥0
 
   JM   −1   J   T   z   C   +q≥ 0 ⊥z   C ≥0.  (40)
 
     Continuing, an MLCP can be solved by interleaving PGS iterations for the contact constraints, and with a direct (i.e., exact) solver for the fixed joints, as shown in Algorithm 3, as follows: 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 3 MLCP solver with iterative contacts  
               
               
                 and exact joints 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                 1: z C =0 
               
               
                   
                 2: z F =0  
               
               
                   
                 3: for all iterations do 
               
               
                   
                 4:  for all contacts do 
               
               
                   
                 5:   z C  = max (0, z C  − D CC   −1 (q C  + N CC z C  + N CF z F )) 
               
               
                   
                 6:  end for 
               
               
                   
                 7:  z F  = linsolve(N FF , q F  + N FC Z C ) 
               
               
                   
                 8: end for 
               
               
                   
               
            
           
         
       
     
     Proof of convergence for the iterative MLCP with exact joints for Algorithm 3 is described. In particular, an LCP iteration is applied, wherein λ=1, and ω=1 in Equation 41, as follows: 
     
       
         
           
             
               
                 
                   
                     E 
                     = 
                     
                       [ 
                       
                         
                           
                             
                               D 
                               CC 
                               
                                 - 
                                 1 
                               
                             
                           
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                           
                             
                               N 
                               FF 
                               
                                 - 
                                 1 
                               
                             
                           
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     K 
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 L 
                                 CC 
                               
                             
                             
                               0 
                             
                           
                           
                             
                               
                                 N 
                                 FC 
                               
                             
                             
                               0 
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   41 
                   ) 
                 
               
             
           
         
       
     
     A new projection operator P is defined that applies (+) only to z C . This gives the iteration, as described in Equations 42 and 43, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       z 
                       C 
                       
                         r 
                         + 
                         1 
                       
                     
                     = 
                     
                       
                         ( 
                         
                           
                             z 
                             C 
                             r 
                           
                           - 
                           
                             
                               D 
                               CC 
                               
                                 - 
                                 1 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   q 
                                   C 
                                 
                                 + 
                                 
                                   
                                     L 
                                     CC 
                                   
                                   ⁢ 
                                   
                                     z 
                                     C 
                                     
                                       r 
                                       + 
                                       1 
                                     
                                   
                                 
                                 + 
                                 
                                   
                                     
                                       [ 
                                       
                                         
                                           
                                             
                                               U 
                                               CC 
                                             
                                           
                                           
                                             
                                               N 
                                               CF 
                                             
                                           
                                         
                                       
                                       ] 
                                     
                                     ⁡ 
                                     
                                       [ 
                                       
                                         
                                           
                                             
                                               z 
                                               C 
                                             
                                           
                                         
                                         
                                           
                                             
                                               z 
                                               F 
                                             
                                           
                                         
                                       
                                       ] 
                                     
                                   
                                   r 
                                 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                       + 
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   42 
                   ) 
                 
               
             
             
               
                 
                   
                     z 
                     F 
                     
                       r 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       z 
                       F 
                       r 
                     
                     - 
                     
                       
                         
                           N 
                           FF 
                           
                             - 
                             1 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               q 
                               F 
                             
                             + 
                             
                               
                                 [ 
                                 
                                   
                                     
                                       
                                         N 
                                         FC 
                                       
                                     
                                     
                                       
                                         N 
                                         FF 
                                       
                                     
                                   
                                 
                                 ] 
                               
                               ⁡ 
                               
                                 [ 
                                 
                                   
                                     
                                       
                                         z 
                                         C 
                                         
                                           r 
                                           + 
                                           1 
                                         
                                       
                                     
                                   
                                   
                                     
                                       
                                         z 
                                         F 
                                         r 
                                       
                                     
                                   
                                 
                                 ] 
                               
                             
                           
                           ) 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   43 
                   ) 
                 
               
             
           
         
       
     
     This is the iteration that is evaluated by Algorithm 3. To prove convergence, for all y≥0, the following must be shown, as described in Equation 44, as follows:
 
( z   P   −z ) T   E   −1 ( z   P   −y )≥0.  (44)
 
     This is proved in Equation 45, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         ( 
                         
                           
                             z 
                             P 
                           
                           - 
                           z 
                         
                         ) 
                       
                       T 
                     
                     ⁢ 
                     
                       
                         E 
                         
                           - 
                           1 
                         
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             z 
                             P 
                           
                           - 
                           y 
                         
                         ) 
                       
                     
                   
                   = 
                   
                     
                       
                         
                           
                             [ 
                             
                               
                                 
                                   
                                     
                                       z 
                                       C 
                                       + 
                                     
                                     - 
                                     
                                       z 
                                       C 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     
                                       z 
                                       F 
                                     
                                     - 
                                     
                                       z 
                                       F 
                                     
                                   
                                 
                               
                             
                             ] 
                           
                           T 
                         
                         ⁡ 
                         
                           [ 
                           
                             
                               
                                 
                                   D 
                                   CC 
                                 
                               
                               
                                 0 
                               
                             
                             
                               
                                 0 
                               
                               
                                 
                                   N 
                                   FF 
                                 
                               
                             
                           
                           ] 
                         
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             
                               
                                 
                                   z 
                                   C 
                                   + 
                                 
                                 - 
                                 
                                   y 
                                   C 
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   z 
                                   F 
                                 
                                 - 
                                 
                                   y 
                                   F 
                                 
                               
                             
                           
                         
                         ] 
                       
                     
                     = 
                     
                       
                         
                           
                             ( 
                             
                               
                                 z 
                                 C 
                                 + 
                               
                               - 
                               
                                 z 
                                 C 
                               
                             
                             ) 
                           
                           T 
                         
                         ⁢ 
                         
                           
                             D 
                             CC 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 z 
                                 C 
                                 + 
                               
                               - 
                               
                                 y 
                                 C 
                               
                             
                             ) 
                           
                         
                       
                       ≥ 
                         
                       ⁢ 
                       0 
                     
                   
                 
               
               
                 
                   ( 
                   45 
                   ) 
                 
               
             
           
         
       
     
     Continuing, considering the specific MLCP arising from contact Jacobian C and joint Jacobian F, as described in Equation 46, as follows: 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             z 
                             C 
                           
                         
                       
                       
                         
                           
                             z 
                             F 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       MLCP 
                       ⁡ 
                       
                         ( 
                         
                           
                             
                               [ 
                               
                                 
                                   
                                     C 
                                   
                                 
                                 
                                   
                                     F 
                                   
                                 
                               
                               ] 
                             
                             ⁢ 
                             
                               
                                 W 
                                 ⁡ 
                                 
                                   [ 
                                   
                                     
                                       
                                         C 
                                       
                                     
                                     
                                       
                                         F 
                                       
                                     
                                   
                                   ] 
                                 
                               
                               T 
                             
                           
                           , 
                           
                             [ 
                             
                               
                                 
                                   
                                     Cv 
                                     
                                       S 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       0 
                                     
                                   
                                 
                               
                               
                                 
                                   0 
                                 
                               
                             
                             ] 
                           
                         
                         ) 
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   46 
                   ) 
                 
               
             
           
         
       
     
     As with PGS, the MLCP is optimized in one embodiment by introducing an auxiliary variable, as described in Equation 47, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       v 
                       S 
                     
                     = 
                     
                       
                         
                           W 
                           ⁡ 
                           
                             [ 
                             
                               
                                 
                                   
                                     C 
                                     T 
                                   
                                 
                                 
                                   
                                     F 
                                     T 
                                   
                                 
                               
                             
                             ] 
                           
                         
                         ⁡ 
                         
                           [ 
                           
                             
                               
                                 
                                   z 
                                   C 
                                 
                               
                             
                             
                               
                                 
                                   z 
                                   F 
                                 
                               
                             
                           
                           ] 
                         
                       
                       + 
                       
                         v 
                         
                           S 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           0 
                         
                       
                     
                   
                   ⁢ 
                   
                       
                   
                 
               
               
                 
                   ( 
                   47 
                   ) 
                 
               
             
           
         
       
     
     This results in Algorithm 4, as described below: 
                             Algorithm 4 Velocity MLCP solver                                            1: for all iterations do           2:  v S  = PGSiter(z C , CWC T , C, v S )           3:  z F  = linsolve(FWF T , Fv S )           4:  v S  = v S  + WF T  z F              5: end for                    
Solving the Split System
 
     The split system includes contact impulses that satisfy the original system, but the system can be reduced in embodiments of the present disclosure. The system is reduced in size through transformation. Algorithm 4 is applied to the split system, which results in Algorithm 5, provided below. This allows for early stoppage of the iterative process to give an approximately solution to the contact constraints, while exactly solving the fixed joint constraints. 
     The first part applies a single iteration of PGS to all the contacts. Each sub-body has exactly one constraint, so the contacts can be safely computed concurrently on m threads. Equations 48 and 49 update the velocity of body j of contact α, as follows:
 
 i=s   α,j   (48)
 
[ v   β   S ] i   =[v   β   S ] i   +[W   β ] i   [C   α,β ] i   T ( z   α   C   −t ).  (49)
 
     The fixed joints could be solved with a direct linear solver, but this would take too much memory and time. However, in embodiments of the present disclosure, the fixed joints can be solved using a lower-cost analytical solution, as expressed in Equation 50, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 F 
                                 β 
                               
                               ⁢ 
                               
                                 W 
                                 β 
                               
                               ⁢ 
                               
                                 F 
                                 β 
                                 T 
                               
                             
                             ) 
                           
                           
                             - 
                             1 
                           
                         
                         = 
                           
                         ⁢ 
                         
                           
                             n 
                             β 
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               
                                 M 
                                 β 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                   + 
                                   
                                     S 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                             
                               - 
                               1 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             n 
                             β 
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               
                                 M 
                                 β 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                   - 
                                   
                                     
                                       n 
                                       β 
                                       
                                         - 
                                         1 
                                       
                                     
                                     ⁢ 
                                     
                                       S 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             n 
                                             β 
                                           
                                           - 
                                           1 
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   54 
                   ) 
                 
               
             
           
         
       
     
     Equation 50 shows that the fixed joints can be enforced by averaging the sub-body velocities of each body at each iteration. This is proven in the following discussion, wherein for notational convenience, let S(n,m) be a matrix of size (6n, 6m) whose elements are all one, and let I(n,m) be the identify matrix of the same size, where it is defined that S(n)=S(n,n) and I(n)=I(n,n). This results in Equation 51, as follows:
 
 F   β   =[S ( n   β −1,1)− I ( n   β −1)].  (51)
 
     The following proof makes use of the identity, as described in Equation 52, as follows:
 
 S ( n,m ) S ( m,p )= mS ( n,p ).  (52)
 
     What is proven is shown in Equation 53, as a Theorem, as follows:
 
( I ( n   β )− F   β   T ( F   β   W   β   F   β   T ) −1   F   β   W   β )= n   β   −1   S ( n   β ).  (53)
 
     Equation 54 starts the proof, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 F 
                                 β 
                               
                               ⁢ 
                               
                                 W 
                                 β 
                               
                               ⁢ 
                               
                                 F 
                                 β 
                                 T 
                               
                             
                             ) 
                           
                           
                             - 
                             1 
                           
                         
                         = 
                           
                         ⁢ 
                         
                           
                             n 
                             β 
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               
                                 M 
                                 β 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                   + 
                                   
                                     S 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                             
                               - 
                               1 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             n 
                             β 
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               
                                 M 
                                 β 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           n 
                                           β 
                                         
                                         - 
                                         1 
                                       
                                       ) 
                                     
                                   
                                   - 
                                   
                                     
                                       n 
                                       β 
                                       
                                         - 
                                         1 
                                       
                                     
                                     ⁢ 
                                     
                                       S 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             n 
                                             β 
                                           
                                           - 
                                           1 
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   54 
                   ) 
                 
               
             
           
         
       
     
     In addition, Equation 55 shows the following:
 
 F   β   T ( I ( n   β −1)− n   β   −1   S ( n   β −1)) F   β   =I ( n   β )− n   β   −1   S ( n   β ).  (55)
 
     As such, Equation 56 results, wherein the velocities are averaged, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       I 
                       ⁢ 
                       
                         ( 
                         
                           n 
                           β 
                         
                         ) 
                       
                     
                     - 
                     
                       
                         
                           
                             F 
                             β 
                             T 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 F 
                                 β 
                               
                               ⁢ 
                               
                                 W 
                                 β 
                               
                               ⁢ 
                               
                                 F 
                                 β 
                                 T 
                               
                             
                             ) 
                           
                         
                         
                           - 
                           1 
                         
                       
                       ⁢ 
                       
                         F 
                         β 
                       
                       ⁢ 
                       
                         W 
                         β 
                       
                     
                   
                   = 
                   
                     
                       
                         I 
                         ⁡ 
                         
                           ( 
                           
                             n 
                             β 
                           
                           ) 
                         
                       
                       - 
                       
                         
                           n 
                           β 
                           
                             - 
                             1 
                           
                         
                         ⁢ 
                         
                           M 
                           β 
                         
                         ⁢ 
                         
                           
                             F 
                             β 
                             T 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 I 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       n 
                                       β 
                                     
                                     - 
                                     1 
                                   
                                   ) 
                                 
                               
                               - 
                               
                                 
                                   n 
                                   β 
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   S 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         n 
                                         β 
                                       
                                       - 
                                       1 
                                     
                                     ) 
                                   
                                 
                               
                             
                             ) 
                           
                         
                         ⁢ 
                         
                           F 
                           β 
                         
                         ⁢ 
                         
                           W 
                           β 
                         
                       
                     
                     = 
                     
                       
                         n 
                         β 
                         
                           - 
                           1 
                         
                       
                       ⁢ 
                       
                         
                           S 
                           ⁡ 
                           
                             ( 
                             
                               n 
                               β 
                             
                             ) 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   56 
                   ) 
                 
               
             
           
         
       
     
     As a corollary, the fixed joints of body β is solved by setting the velocity of all its sub-bodies iϵ[1 −  . . . {dot over (n)} β ] as follows in Equation 57: 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         
                           ( 
                           
                             
                               v 
                               S 
                             
                             - 
                             
                               
                                 
                                   
                                     WF 
                                     T 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       FWF 
                                       T 
                                     
                                     ) 
                                   
                                 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 Fv 
                                 S 
                               
                             
                           
                           ) 
                         
                         β 
                       
                       ] 
                     
                     i 
                   
                   = 
                   
                     
                       n 
                       β 
                       
                         - 
                         1 
                       
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         
                           n 
                           β 
                         
                       
                       ⁢ 
                       
                         
                           
                             [ 
                             
                               V 
                               β 
                             
                             ] 
                           
                           k 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   57 
                   ) 
                 
               
             
           
         
       
     
     Continuing, fixed joints are enforced by averaging the sub-body velocities of each body at each iteration, as applied in Algorithm 5, as follows: 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 5 Parallel PGS contacts with averaging for joints 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                  1: p = diag(CWC T ) −1   
               
               
                   
                  2: for all iterations do 
               
               
                   
                  3:  {Apply one iteration of PGS to the sub-body contacts} 
               
               
                   
                  4:  for all contacts, α do {in parallel} 
               
               
                   
                  5:   t = z Cα   
               
               
                   
                  6:   z Cα  = max(0,z Cα  − P α  (Cv S ) α ) 
               
               
                   
                  7:   for all j ∈ {1,2} do 
               
               
                   
                  8:    β = b a,j   
               
               
                   
                  9:    i = s α,j   
               
               
                   
                 10:    [v β   S ] i  = [v β   S ] i  + [W β ] i [C α,β ] i   T  (z α   C  − t) 
               
               
                   
                 11:   end for  
               
               
                   
                 12:  end for 
               
               
                   
                 13:  {Solve the fixed joints by velocity averaging} 
               
               
                   
                 14:  for all bodies, β do 
               
               
                   
                 15:   a = n β   −1  Σ k=1   nβ  [V β ] k   
               
               
                   
                 16:   for all sub-bodies, i do {all sub-bodies get the same vel.} 
               
               
                   
                 17:    [v β   S ] i  = a 
               
               
                   
                 18:   end for 
               
               
                   
                 19:  end for 
               
               
                   
                 20: end for 
               
               
                   
               
            
           
         
       
     
     Also, the sub-body velocities of each body are equal at the end of each iteration, so in one embodiment, the sub-body velocity need not be explicitly stored for every sub-body. In fact, in one embodiment, the masses of the sub-bodies are scaled or calculated dynamically, as the corresponding body mass is scaled before the sub-body mass is used. Using these optimizations and making use of the definitions of C, V S  and W, Algorithm 6 is provided below. In one embodiment, the algorithm is parallelizable, wherein lines 24 to 32 can be implemented using a parallel segmented reduction. 
                             Algorithin 6 Final mass splitting solver, scalar version                                                 1:   {Assemble mass splitting preconditioner}            2:   for all patches, α do            3:    p α  = 0            4:    for all j ∈ {1, 2} do            5:     β = b α,j              6:     {divide mass by number of number of bodies it touches}            7:     
         m   split     =       M   β       n   β           
            8:     p α  = p α  + J αβ  * m split   −1  * J βα   T              9:    end for           10:    p α  = p α   −1             11:   end for           12:   {Solver}           13:   for all iterations do           14:    {Solve each sub-body contact}           15:    for all contacts, α do {in parallel}           16:     t = z Cα             17:     z Cα  = max (0, z Cα  − p α (Jv) α )           18:     for all j ∈ {1, 2} do           19:      β = b αj             20:      D αj  = v β  + n β M β   −1 [J αβ ] T (z α   C  − t)           21:     end for           22:    end for           23:    {Fixed joints, implement with parallel segmented reduction}           24:    for all bodies, β do           25:     a = 0           26:     for all j ∈ {1, 2} do           27:      for all i such that b ij  = β do           28:       a = a + D ij             29:      end for           30:     end for           31:     v β  = n β   −1 a           32:    end for           33:   end for                    
Block Method and System for a Rigid Body Solver
 
     The scalar rigid body solver described in  FIGS. 3 and 4  so far meets the jitter and load balancing requirements. To make the solution provided by the solver descried in  FIGS. 3 and 4  converge more quickly, constraints are solved in blocks, in one embodiment as is described in  FIGS. 7 and 8 . In the scalar method described previously in relation to  FIGS. 3-6 , the bodies are split so that each sub-body had exactly one constraint. In the present embodiment, the collision detection system takes pairs of potentially colliding bodies and finds their intersections. For convex pairs, the intersection set has one connected component, for non-convex pairs the intersection can have many components. Each component of each intersection is referred to as a patch. In the block solver of the present embodiment, the bodies are split so that so that there is one sub-body per patch. As such, each sub-body now has multiple contacts, so contacts cannot be processed in parallel. Instead, each patch is processed in parallel, and the contents are solved within each path sequentially on its thread using PGS. 
       FIG. 7  is a block diagram of a linear complementarity, parallel iterative rigid body solver  700  that implements a block method for simulating complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts, in accordance with one embodiment of the present disclosure. In one embodiment, solver  700  performs the functions of solver  250  of  FIG. 2  within the rigid body physics system or engine  200 . 
     Solver  700  includes a patch determining module  710 . For a particular original object, a pair of potentially colliding bodies defines an intersection set of contacts. That is, a patch defines a plurality of contact constraints affecting the pair of colliding objects. One of the pair of colliding objects includes the original object, which is associated with an original mass. In one embodiment, the intersection set has one connected component (e.g., for convex pairs). In another embodiment, the intersection set has multiple components (e.g., for non-convex pairs). 
     Solver  700  includes an object and mass splitter  720 , for splitting an original object into a plurality of sub-bodies. An original object is split so that there is one sub-body per patch. As such, each sub-body has multiple contacts. Instead of processing the contacts in parallel, patches are processed in parallel. As such, object splitter  720  is configured to split the original object by a total number of patches. Also, each sub-body has a mass, which is defined by the original mass divided by the total number of patches. 
     Contact constraint solver  730  is configured to solve contact constraints approximately for each sub-body. Further, each sub-body has multiple contacts and/or multiple contact constraints. In one embodiment, each contact constraint within a patch is solved approximately and sequentially. By solving for the contact constraints within a patch, a sub-body position and/or velocity is determined. 
     Solver  700  also includes a fixed joint solver  740  that enforces fixed joint constraints exactly, such that positions of each sub-body are identical. The fixed joint solver  740  is necessary because applying an impulse at a contact of the original body results in velocities at other contacts of the same body. The goal is for each of the sub-bodies to behave exactly like the original body, such that when an impulse is applied at one sub-body, it is necessary to change the velocities of other sub-bodies of the same body. The requirement that all the sub-bodies move together is modeled as a fixed joint, and is enforced by the fixed joint solver  740 . In one embodiment, the fixed joint is enforced or solved at each iteration by averaging the sub-body velocities and applying the average velocity to each of the sub-bodies. The average velocity is used as the initial velocity of the sub-bodies in the next iteration to solve contact constraints and fixed joint constraints 
       FIG. 8  is a flow diagram  800  illustrating a computer implemented block method for solving linear complementarity iterative rigid body dynamics, in accordance with one embodiment of the present disclosure. In another embodiment, flow diagram  400  is implemented within a computer system including a processor and memory coupled to the processor and having stored therein instructions that, if executed by the computer system causes the system to execute a block method for solving linear complementarity iterative rigid body dynamics. In still another embodiment, instructions for performing a method are stored on a non-transitory computer-readable storage medium having computer-executable instructions for causing a computer system to perform a block method for solving linear complementarity iterative rigid body dynamics. In one embodiment, flow diagram  800  is implemented by solver  700  of  FIG. 7 , and/or solver  250  of  FIG. 2 . More particularly, flow diagram  800  is implemented to take a set of constraints that the positions and/or velocities of objects must satisfy and calculate a set of forces and/or impulses such that the objects will satisfy the constraints once the forces and/or impulses are applied to the objects. 
     The method of flow diagram  800  is implemented after a set of constraints are defined, wherein contact constraints affect one or more objects or bodies. For instance, constraint generator  220  of  FIG. 2  is able to determine the constraints that affect one or more objects or bodies. Once the constraints are defined, the block method of flow diagram  800  is implemented to determine the positions and/or velocities of the one or more objects or bodies. 
     At  810 , determining one or more patches associated with an original object is performed. The original object has an original mass. Also, each patch comprises a pair of colliding objects, one of which includes the original object. Further, each patch defines a plurality of contact constraints that affect the pair of colliding objects. As such, each patch is associated with multiple contacts. 
     At  820 , splitting the original object is performed. In particular, the original object is split into a plurality of sub-bodies, such that there is one sub-body per patch. As such, the number of patches defines the number of sub-bodies. In addition, the mass for each sub-body is defined by dividing the original mass for the original object by the total number of patches. 
     At  830 , solving patches for their corresponding contact constraints is performed. In particular, for each patch, contact constraints are solved sequentially. That is, for each contact constraint of a patch, an impulse is calculated to satisfy the contact constraint. This is accomplished by solving for a sub-body velocity vector of the sub-body due to the contact constraint. In one embodiment, the relative velocity vector is constrained to be resting or separating (non-negative). Also, the impulses are constrained to be non-attractive (non-negative). In addition, impulse should be zero at separating contacts. This is performed until a final velocity and position is determined for the sub-body associated with the patch. 
     At  840 , enforcing fixed joints exactly is performed, such that fixed joint constraints are solved exactly. More particularly, positions and/or velocities of each sub-body are solved such that they are substantially identical. For instance, in one embodiment, the enforcing of fixed joints is performed by averaging a plurality of sub-body velocity vectors corresponding to the plurality of sub-bodies. As such, an average velocity vector is defined. Further, the average velocity vector is assigned to each of the plurality of sub-bodies. 
     In one embodiment, the operations in  830  and  840  are iteratively performed. That is, for each iteration, given a starting average velocity vector from the previous iteration, a new average velocity vector and new position are determined for a corresponding sub-body. These new values again correspond to an average velocity vector and/or position for the original object that can be used finally or for the next iteration. 
     In another embodiment, a frictional force is applied to the sub-bodies at each iteration. Embodiments of the present disclosure support any means for applying friction to the block method for solving linear complementarity iterative rigid body dynamics. Embodiments of the present disclosure apply a coulomb friction model, though other friction models are supported. For instance, friction may be considered by applying a friction model, such that the more force an object is exerting on a surface, the harder it needs to be pushed to move it across the surface. As such, heavier objects are harder to push than lighter objects. In addition, the transition from sticking (e.g., when the object does not move) to sliding (e.g., when the object first moves) is considered. 
     The solution for the block method for performing linear complementarity in parallel for an iterative rigid body solver is described below. As previously described, the block method simulates complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts. The block method converges more quickly than the scalar method previously described in relation to  FIGS. 3-6 . 
     For definition, let b i,1   P  and b i,2   P  be the two bodies constrained by patch i. Also, let n β   P  be the number of patches affecting body β in the unsplit system, as described in Equation 58, as follows:
 
 n   β   P   =|{i:b   i,1   P   =βvb   i,2   P =β}|.  (58)
 
Splitting the Bodies in the Block Method
 
     The bodies are split into sub-bodies so that each sub-body has one patch. The total number of sub-bodies is described in Equation 59, as follows: 
     
       
         
           
             
               
                 
                   
                     n 
                     B 
                   
                   = 
                   
                     
                       ∑ 
                       
                         β 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                       
                         n 
                         β 
                         P 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   59 
                   ) 
                 
               
             
           
         
       
     
     Again, the mass of each sub-body is split evenly between its sub-bodies, in one embodiment, as described in Equation 60, as follows: 
                     W   β   P     =       [             n   β   P     ⁢     M   β     -   1                         0                       ⋱                       0                       n   β   P     ⁢     M   β     -   1               ]     ∈         6   ⁢     n   β   P     ×   6   ⁢     n   β   P         .               (   60   )               
Joining the Bodies Back Together
 
     The sub-bodies of each body are joined together with n β   P −1 fixed joints. The 
     Jacobian of the fixed joints of body β is described in Equation 61, as follows: 
                     F   β   P     =       [           I   6           -     I   6           0       …       0             I   6         0         -     I   6                       ⋮           ⋮       ⋮                   ⋱       0             I   6         0       …       0         -     I   6             ]     ∈         6   ⁢     n   β   P     ×   6   ⁢     (       n   β   P     -   1     )         .               (   61   )               
Assigning the Contacts to the Split Bodies
 
     Let the total number of patches be p and the number of contacts in patch i be m i . The rows of J are partitioned into blocks as well as the columns, so that J iβ   ^P ϵ   m     i     ×6  represents the effect of patch i on body β as shown in Equation 62, as provided below: 
     
       
         
           
             
               
                 
                   J 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               J 
                               11 
                               P 
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               
                                 1 
                                 ⁢ 
                                 n 
                               
                               P 
                             
                           
                         
                         
                           
                             ⋮ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               J 
                               
                                 p 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                               P 
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               pn 
                               P 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   62 
                   ) 
                 
               
             
           
         
       
     
     In addition, each patch of a number of patches n β   P  that is acting on body β is assigned to its own sub-body, as described in Equation 63, as follows: 
                       [     C   αβ   P     ]     ij     =     {             [     J   αβ   P     ]     ij             if   ⁢           ⁢     r   j       =   α             0         otherwise   .                     (   63   )               
Solving the Block Split System
 
     A block-sparse MLCP is defined in Equation 64, as follows: 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             z 
                             C 
                             P 
                           
                         
                       
                       
                         
                           
                             z 
                             F 
                             P 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       MLCP 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           ( 
                           
                             
                               [ 
                               
                                 
                                   
                                     
                                       N 
                                       CC 
                                       P 
                                     
                                   
                                   
                                     
                                       N 
                                       CF 
                                       P 
                                     
                                   
                                 
                                 
                                   
                                     
                                       N 
                                       FC 
                                       P 
                                     
                                   
                                   
                                     
                                       N 
                                       FF 
                                       P 
                                     
                                   
                                 
                               
                               ] 
                             
                             , 
                             
                               [ 
                               
                                 
                                   
                                     q 
                                   
                                 
                                 
                                   
                                     0 
                                   
                                 
                               
                               ] 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
 
                         
                         [ 
                         
                           
                             
                               
                                 N 
                                 CC 
                                 P 
                               
                             
                             
                               
                                 N 
                                 CF 
                                 P 
                               
                             
                           
                           
                             
                               
                                 N 
                                 FC 
                                 P 
                               
                             
                             
                               
                                 N 
                                 FF 
                                 P 
                               
                             
                           
                         
                         ] 
                       
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 C 
                                 P 
                               
                             
                           
                           
                             
                               
                                 F 
                                 P 
                               
                             
                           
                         
                         ] 
                       
                       ⁢ 
                       
                         
                           
                             
                               W 
                               P 
                             
                             ⁡ 
                             
                               [ 
                               
                                 
                                   
                                     
                                       C 
                                       P 
                                     
                                   
                                 
                                 
                                   
                                     
                                       F 
                                       P 
                                     
                                   
                                 
                               
                               ] 
                             
                           
                           T 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   64 
                   ) 
                 
               
             
           
         
       
     
     In Equation 64, the term N CC   P  is a block diagonal with blocks defined in Equation 65, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           B 
                           αα 
                         
                         = 
                           
                         ⁢ 
                         
                           
                             ∑ 
                             
                               k 
                               = 
                               1 
                             
                             n 
                           
                           ⁢ 
                           
                             
                               C 
                               
                                 α 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 k 
                               
                               P 
                             
                             ⁢ 
                             
                               
                                 
                                   W 
                                   k 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     C 
                                     
                                       α 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       k 
                                     
                                     P 
                                   
                                   ) 
                                 
                               
                               T 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             ∑ 
                             
                               k 
                               = 
                               1 
                             
                             n 
                           
                           ⁢ 
                           
                             
                               n 
                               k 
                               P 
                             
                             ⁢ 
                             
                               J 
                               
                                 α 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 k 
                               
                               P 
                             
                             ⁢ 
                             
                               
                                 
                                   M 
                                   k 
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     J 
                                     
                                       α 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       k 
                                     
                                     P 
                                   
                                   ) 
                                 
                               
                               T 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               n 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   1 
                                 
                               
                             
                             ⁢ 
                             
                               J 
                               
                                 α 
                                 , 
                                 
                                   b 
                                   
                                     α 
                                     , 
                                     1 
                                   
                                 
                               
                               P 
                             
                             ⁢ 
                             
                               
                                 
                                   M 
                                   
                                     b 
                                     
                                       α 
                                       , 
                                       1 
                                     
                                   
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     J 
                                     
                                       α 
                                       , 
                                       
                                         b 
                                         
                                           α 
                                           , 
                                           1 
                                         
                                       
                                     
                                     P 
                                   
                                   ) 
                                 
                               
                               T 
                             
                           
                           + 
                           
                             
                               n 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   2 
                                 
                               
                             
                             ⁢ 
                             
                               J 
                               
                                 α 
                                 , 
                                 
                                   b 
                                   
                                     α 
                                     , 
                                     2 
                                   
                                 
                               
                               P 
                             
                             ⁢ 
                             
                               
                                 
                                   
                                     M 
                                     
                                       b 
                                       
                                         α 
                                         , 
                                         2 
                                       
                                     
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       J 
                                       
                                         α 
                                         , 
                                         
                                           b 
                                           
                                             α 
                                             , 
                                             2 
                                           
                                         
                                       
                                       P 
                                     
                                     ) 
                                   
                                 
                                 T 
                               
                               . 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   65 
                   ) 
                 
               
             
           
         
       
     
     Let L αα  and D αα  be the lower triangle and diagonal of B αα . Applying the LCP iteration, with λ=1 and ω=1 results in Equations 66-68, as follows: 
     
       
         
           
             
               
                 
                   
                     E 
                     = 
                     
                       [ 
                       
                         
                           
                             
                               D 
                               11 
                             
                           
                           
                             
                                 
                             
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                         
                           
                             
                                 
                             
                           
                           
                             ⋱ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             0 
                           
                           
                             
                                 
                             
                           
                           
                             
                               D 
                               pp 
                             
                           
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                           
                             … 
                           
                           
                             0 
                           
                           
                             
                               N 
                               FF 
                               P 
                             
                           
                         
                       
                       ] 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   and 
                 
               
               
                 
                   ( 
                   66 
                   ) 
                 
               
             
             
               
                 
                   
                     K 
                     = 
                     
                       [ 
                       
                         
                           
                             
                               L 
                               11 
                             
                           
                           
                             
                                 
                             
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                         
                           
                             
                                 
                             
                           
                           
                             ⋱ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             0 
                           
                           
                             
                                 
                             
                           
                           
                             
                               L 
                               PP 
                             
                           
                           
                             0 
                           
                         
                         
                           
                             
                               
                                 ( 
                                 
                                   N 
                                   PC 
                                   P 
                                 
                                 ) 
                               
                               1 
                             
                           
                           
                             … 
                           
                           
                             
                               
                                 ( 
                                 
                                   N 
                                   FC 
                                   P 
                                 
                                 ) 
                               
                               p 
                             
                           
                           
                             0 
                           
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   And 
                 
               
               
                 
                   ( 
                   67 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         
                           z 
                           
                             C 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             1 
                           
                           
                             r 
                             + 
                             1 
                           
                         
                       
                     
                     
                       
                         ⋮ 
                       
                     
                     
                       
                         
                           z 
                           Cp 
                           
                             r 
                             + 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           z 
                           F 
                           
                             r 
                             + 
                             1 
                           
                         
                       
                     
                   
                   = 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 z 
                                 1 
                                 r 
                               
                               - 
                               
                                 
                                   D 
                                   11 
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       q 
                                       1 
                                     
                                     + 
                                     
                                       
                                         L 
                                         11 
                                       
                                       ⁢ 
                                       
                                         z 
                                         1 
                                         
                                           r 
                                           + 
                                           1 
                                         
                                       
                                     
                                     + 
                                     
                                       
                                         U 
                                         11 
                                       
                                       ⁢ 
                                       
                                         z 
                                         1 
                                         r 
                                       
                                     
                                     + 
                                     
                                       
                                         
                                           ( 
                                           
                                             N 
                                             CF 
                                             P 
                                           
                                           ) 
                                         
                                         1 
                                       
                                       ⁢ 
                                       
                                         z 
                                         F 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                             ) 
                           
                           + 
                         
                       
                     
                     
                       
                         ⋮ 
                       
                     
                     
                       
                         
                           
                             
                               ( 
                               
                                 
                                   z 
                                   p 
                                   r 
                                 
                                 - 
                                 
                                   
                                     D 
                                     pp 
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         q 
                                         p 
                                       
                                       + 
                                       
                                         
                                           LL 
                                           pp 
                                         
                                         ⁢ 
                                         
                                           z 
                                           p 
                                           
                                             r 
                                             + 
                                             1 
                                           
                                         
                                       
                                       + 
                                       
                                         
                                           U 
                                           pp 
                                         
                                         ⁢ 
                                         
                                           z 
                                           p 
                                           r 
                                         
                                       
                                       + 
                                       
                                         
                                           
                                             ( 
                                             
                                               N 
                                               CF 
                                               P 
                                             
                                             ) 
                                           
                                           p 
                                         
                                         ⁢ 
                                         
                                           z 
                                           F 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ) 
                             
                             + 
                           
                           - 
                         
                       
                     
                     
                       
                         
                           
                             
                               ( 
                               
                                 N 
                                 FF 
                                 P 
                               
                               ) 
                             
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 α 
                                 = 
                                 1 
                               
                               P 
                             
                             ⁢ 
                             
                               
                                 ( 
                                 
                                   N 
                                   FC 
                                   P 
                                 
                                 ) 
                               
                               ⁢ 
                               
                                 
                                   α 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       z 
                                       C 
                                     
                                     ) 
                                   
                                 
                                 α 
                                 
                                   r 
                                   + 
                                   1 
                                 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   68 
                   ) 
                 
               
             
           
         
       
     
     Algorithm 7 provides an implementation of this iteration. In one embodiment, the blocks of z C  are computed in parallel. 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 7 Block MLCP solver with iterative contacts  
               
               
                 and exact joints 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                  1: z C =0 
               
               
                   
                  2: z F =0 
               
               
                   
                  3: for all iterations do 
               
               
                   
                  4:  for all blocks, α do {in parallel} 
               
               
                   
                  5:   for all contacts in block, i do 
               
               
                   
                  6:    t = [z Cα ] i  − [D αα ] ii   −1  ([q α ] i  + [B αα z Cα ] i  + [(N CF   P z F ) α ] i ) 
               
               
                   
                  7:    [z Cα ] i  = max(0,t) 
               
               
                   
                  8:   end for 
               
               
                   
                  9:  end for 
               
               
                   
                 10:  z F  = linsolve(N FF ,q F  + N FC z C ) 
               
               
                   
                 11: end for 
               
               
                   
               
            
           
         
       
     
     Lines 4-7 of Algorithm 7 are equivalent to running PGS on each block independently. As an improvement, Algorithm 8 is a velocity version of Algorithm 7, as described below: 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 8 Velocity MLCP block solver 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                 1: for all iterations do 
               
               
                   
                 2:  for all blocks, α do {in parallel} 
               
               
                   
                 3:   v S  = PGSiter(z α ,D αα , [C α1  · · · C αn ], v S ) 
               
               
                   
                 4:  end for  
               
               
                   
                 5:  z F  = linsolve(N FF ,N FC z C )  
               
               
                   
                 6:  v S  = v S  + WF T z F    
               
               
                   
                 7: end for 
               
               
                   
               
            
           
         
       
     
     The fixed joints are solved by velocity averaging, as previously described in relation to the scalar method for solving rigid body dynamics. Algorithm 9 provides the velocity averaging for the fixed joints below. 
     
       
         
           
               
             
               
                   
               
               
                 Algorithm 9 Block solver with velocity averaging for joints 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
            
               
                   
                  1: for all iterations do 
               
               
                   
                  2:  for all patches, α do {in parallel) 
               
               
                   
                  3:    for all contacts in patch, i do {sequentially on thread α} 
               
               
                   
                  4:     t = [z Cα ] i   
               
               
                   
                  5:    [z Cα ] i  = max (0, [z Cα ] i  − [D αα ] ii   −1 [C α1  · · · C αn ]v i   S   
               
               
                   
                  6:     for all j ∈ {1,2} do 
               
               
                   
                  7:     β = b i,j   
               
               
                   
                  8:     k = s α,j   
               
               
                   
                  9:     [v β   S ] k  = [v β   S ] k  + [W β ] k [C α,β ] k   T (z cα  − T) 
               
               
                   
                 10:    end for 
               
               
                   
                 11:    end for 
               
               
                   
                 12:   end for 
               
               
                   
                 13:   {Solve the fixed joints by velocity averaging} 
               
               
                   
                 14:   for all bodies, β do 
               
               
                   
                 15:    a = n β   −1 Σ k=1   nβ [V β ] k   
               
               
                   
                 15:   for all sub-bodies, i do {all sub-bodies get the same vel.} 
               
               
                   
                 17:     [v β   S ] i  = a 
               
               
                   
                 18:   end for  
               
               
                   
                 19:   end for 
               
               
                   
                 20: end for 
               
               
                   
               
            
           
         
       
     
     For the definitions of C P  and W P , Equation 69 is provided below:
 
diag( B   αα )=diag 1 ([ J   α1   P    . . . J   αn   P   ]M   −1   [J   α1   P    . . . J   αn   P ] T ).  (69)
 
     Finally, the sub-body velocities are eliminated in the same way as previously described in the scalar method for solving rigid body dynamics, as described in Algorithm 10, as follows: 
                             Algorithm 10 Final mass splitting solver, block version                                                 1:   {Assemble mass splitting preconditioner}            2:   for all patches, α do            3:    for all contacts in patch, i do            4:     [p α ] i  = 0            5:     for all j ∈ {1, 2} do            6:      β = b α,j   P              7:      {divide body&#39;s mass by its patch count}            8:      
         m   split     =       M   β       n   β   P           
            9:      [p α ] i  = [p α ] i  + [J αβ ] i  * m split   −1  * [J βα ] i   T             10:     end for           11:     p α  = p α   −1             12:    end for           13:   end for           14:   {Solver}           15:   for all iterations do           16:    for all patches, α do {in parallel}           17:     β 1 , β 2  = b α1   P , b α2   P             18:     D α1 , D α2  = v β     1   , v β     2               19:     for all contacts in patch, i do {sequentially on thread α}           20:      t = [z Cα ] i             21:      [z Cα ] i  = max (0, [z Cα ] i  − [p α ] i (J α , β     1   D α1  + J α , β     2   D α2 ))           22:      for all j ∈ {1, 2} do           23:       k = s α,j             24:       D αj  = D αj  + M βj   −1 J α,βj   T (z Cα  − t)           25:      end for           26:     end for           27:    end for           28:    {Fixed joints, implement with parallel segmented reduction}           29:    a = 0           30:    for all j ∈ {1, 2} do           31:     for all α such that b αj  = β do           32:      a = a + D αj             33:     end for           34:    end for           35:    v β  = n β   −1 a           36:   end for                    
Modified Effective Mass for Parallel Rigid Body Simulation
 
     Other embodiments of the present invention disclose a parallel iterative rigid body solver that avoids common artifacts at low iteration counts by implementing a projected Jacobi-based method that has guaranteed convergence like PGS, but freedom from jitter and parallel scaling like Jacobi. In large or real-time simulations, iteration is often terminated before convergence to maximize scene size. If the distribution of the resulting residual energy varies too much from frame to frame, then bodies close to rest can visibly jitter. Projected Gauss-Seidel (PGS) distributes the residual according to the order in which contacts are processed, and preserving the order in parallel implementations is very challenging. In contrast, Jacobi-based methods provide order independence, but have slower convergence. Embodiments of the present invention accelerate projected Jacobi by dividing each body mass term in the effective mass by the number of contacts acting on the body, while using the full body mass to apply impulses. 
     Additionally, acceleration of the method is achieved by solving contacts in blocks, providing wallclock performance competitive with PGS while avoiding visible artifacts. In particular, convergence of the method is extended by solving blocks of contacts with PGS, and using projected Jacobi to combine the blocks. Again, the effective mass of the Jacobi iteration is modified, and more particularly dividing the mass terms of the effective mass by the number of contact blocks per body, and again using the full mass to apply impulses. The improved convergence combined with its lower cost per iteration, the method&#39;s wall-clock performance is competitive with parallel PGS. 
       FIG. 9  is a block diagram of a linear complementarity, parallel iterative rigid body solver  900  for simulating complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts, in accordance with one embodiment of the present disclosure. In one embodiment, solver  900  performs the functions of solver  250  of  FIG. 2  within the rigid body physics system or engine  200 . 
     As shown rigid body solver  900  includes an effective mass module  910  that determines a plurality of effective masses for a plurality of contacts between a plurality of bodies, wherein the body mass term in the effective mass is based on and/or related to the number of contacts acting on the body. Also, solver  900  includes a relative velocity module  920  for determining the relative velocities of the contacts based on body velocities of a last iteration or an initial state. Solver  900  also includes an impulse correction solver  930  that is configured for determining for each contact an impulse correction that is based on the determined relative velocity. Further, solver  900  includes an impulse correction applicator  940  for applying the impulse determined from a corresponding impulse correction. The impulse constraining module  950  applies contact constraints to the plurality of impulses to allow for separation of relative velocities and to modify negative relative velocities such that penetration of bodies is disallowed. Solver  900  includes a constrained impulse applicator that applies the constrained impulses to determine updated velocities for the plurality of bodies. 
       FIG. 10  is a flow diagram  1000  illustrating a scalar method for solving linear complementarity iterative rigid body dynamics wherein a body mass term in the effective mass is based on the number of contacts acting on the body, in accordance with one embodiment of the present disclosure. In another embodiment, flow diagram  1000  is implemented within a computer system including a processor and memory coupled to the processor and having stored therein instructions that, if executed by the computer system causes the system to execute a scalar method for solving linear complementarity iterative rigid body dynamics wherein a body mass term in the effective mass is based on the number of contacts acting on the body. In still another embodiment, instructions for performing a method are stored on a non-transitory computer-readable storage medium having computer-executable instructions for causing a computer system to perform a scalar method for solving linear complementarity iterative rigid body dynamics wherein a body mass term in the effective mass is based on the number of contacts acting on the body. In one embodiment, flow diagram  1000  is implemented by solver  900  of  FIG. 9 , and/or solver  250  of  FIG. 2 . 
     The method of flow diagram  1000  is implemented after a set of constraints are defined, wherein contact constraints affect one or more objects or bodies. For instance, constraint generator  220  of  FIG. 2  is able to determine the constraints that affect one or more objects or bodies. Once the constraints are defined, the scalar method of flow diagram  1000  is implemented to determine the positions, orientations, and/or velocities of the one or more objects or bodies. 
     The method outlined in  FIG. 10  is performed iteratively in order to converge values for velocities and/or positions of a plurality of objects. Specifically, at  1010 , the method includes determining a plurality of modified effective masses for a plurality of contacts between a plurality of bodies. The distinctive feature of the present embodiment provides for determining each modified effective mass based on a number of contacts affecting the corresponding body. In particular, for each body mass, the modified effective mass is determined by dividing the corresponding body mass by a corresponding number of contacts affecting the corresponding body mass. For example, the effective mass module  910  of  FIG. 9  is configured to determine the plurality of modified effective masses. 
     The method includes at  1020 , determining a plurality of relative velocities based on a plurality of body velocities that is determined from a last iteration and/or an initial state. That is, relative velocities for the plurality of contacts associated with pairs of bodies are determined based on the velocities of the bodies as determined in the last iterative step. In another implementation, the relative velocities are determined based on an initial state of the body velocities, such as when performing the first iterative step. That is, an initial velocity is assigned to each of the plurality of bodies. 
     In one embodiment, the plurality of relative velocities are adjusted by an additional velocity component. For instance, the component provides for restitution or elasticity of the bodies, such as, through the application of restitution parameters, in one embodiment. In another embodiment, the component provides for position correction of bodies. 
     The method includes at  1030 , determining a plurality of impulse corrections based on the plurality of modified effective masses and the plurality of relative velocities in order to negate influences of the plurality of relative velocities. That is a plurality of relative impulses are determined to negate the relative velocities of the contacts, in order to bring the bodies to rest. The relative impulse for each contact point is determined based on the negative relationship between the effective mass and corresponding relative velocity. 
     The method includes at  1040 , determining a plurality of updated impulses based on the plurality of impulse corrections. Specifically, for each contact, the impulse correction is subtracted from the impulse determined from the last iteration, and/or its initial state, to determine a corresponding updated impulse for the current iteration. 
     The method includes at  1050 , applying the plurality of updated impulses to the plurality of bodies based on a plurality of original masses of the plurality of bodies and the plurality of body velocities determined from the last iteration and/or initial state to determine a plurality of updated velocities of the plurality of bodies. Of particular note, the full body mass is used to apply impulses. 
     In another embodiment, the method includes constraining the plurality of updated impulses to determine a plurality of constrained, updated impulses. For example, a contact constraint provides for separation of relative velocities, and disallows the penetration of bodies associated with negative relative velocities. For instance, the contact constraint modifies negative relative velocities for corresponding contacts by assigning each corresponding negative relative impulse to a value of zero. As such, the plurality of constrained, updated impulses is applied to the plurality of bodies based on the original masses of the plurality of bodies, and the plurality of velocities of the bodies relating to the a last iteration. In this manner, a plurality of updated velocities for the plurality of bodies is determined that considers contact constraints. 
     In one embodiment, the method outlined in  FIG. 10  is performed in order to determine positions of the plurality of bodies. For instance, instead of determining relative velocities, the plurality of impulses is applied to the plurality of bodies to determine a plurality of new positions of the plurality of bodies based on the plurality of original masses of the plurality of bodies and the plurality of velocities related to a last iteration or initial state. 
     In one embodiment, the operations performed in  1020 ,  1030 ,  1040 , and  1050  are performed iteratively to converge the velocities and positions of the plurality of bodies. That is, for each iteration, given a velocity from the previous iteration, an updated velocity and updated position and/or updated orientations are determined for a corresponding body. These updated values again correspond to a velocity and position and/or orientation for the object that can be used finally or for the next iteration. In particular, at decision step  1060 , the method determines whether the velocities and/or positions converge sufficiently. For instance, threshold parameters related to velocities and/or distance may be defined to determine whether convergence is met. In another implementation, the number of iterations determines whether or not convergence is satisfied. If not, the process iteratively proceeds back to  1020 . Otherwise, the process ends. 
     In another embodiment, a frictional force is applied to determine the velocities at each iteration. Embodiments of the present disclosure support any means for applying friction. Embodiments of the present disclosure apply a coulomb friction model, though other friction models are supported. For instance, friction may be considered by applying a friction model, such that the more force an object is exerting on a surface, the harder it needs to be pushed to move it across the surface. As such, heavier objects are harder to push than lighter objects. In addition, the transition from sticking (e.g., when the object does not move) to sliding (e.g., when the object first moves) is considered. 
     In still another embodiment, the process outlined in flow diagram  1000  is preformed in parallel, such that operations outlined in  1010 ,  1020 ,  1030 ,  1040 , and  1050  to include determining a plurality of modified effective masses, determining a plurality of relative velocities, determining a plurality of impulse corrections, determining a plurality of updated impulses, and applying the plurality of updated impulses are performed in parallel. 
     The solution for the scalar method for performing linear complementarity in parallel for an iterative rigid body solver, wherein a body mass term in the effective mass is based on the number of contacts acting on the body, is described below. As previously described, the scalar method simulates complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts. Portions described below are equally adaptable for use in the block method for performing linear complementarity in parallel for an iterative rigid body solver, in other embodiments. 
     As previously described, a rigid body contact model is described in a scene having n rigid bodies with positions xϵ   6n , external forces f e ϵ   6n  and masses Mϵ   6n×6n . Collision detection identifies m contacts between the rigid bodies, represented by constraints Φ(x)≥0 with Jacobian ∂Φ/∂x=Jϵ   m×6n . For simplicity of description constraint stabilization is omitted, so the unknowns are the forces λϵ   m  necessary to satisfy the time derivative of the constraints. Friction is treated separately, in one embodiment. Also, contacts must satisfy the Signorini condition: i.e., forces must not be attractive (e.g., λ≥0), velocities must move the system out of penetration 
               (           ∂   Φ       ∂   x       ⁢   v     ≥   0     )     ,         
and a force should be applied at a contact only if that contact is not separating, which is expressed as λ≥0⊥Jv≥0. The continuous model was described previously in Equations 1-3.
 
     Also previously described in relation to Equations 4-11, the model is discretized for purposes of time stepping. Two iterative methods are used in the current model to solve the symmetric LCP: projected Guass-Seidel (PGS) and projected Jacobi. Both methods are described below. 
     Projected Jacobi 
     Projected Jacobi solves each constraint in isolation and then merges the results, allowing efficient parallel implementation. Abstractly, it generates the following iteration, as described in Equation 70, as follows:
 
 z   r+1 =( z   r   −D   −1 ( Nz   r   +q )) +   (70)
 
     In Equation 70, D=diag(N) and (x + ) i =max(0, x i ). In the rigid body model, D i   −1  can be interpreted as being the effective mass along the contact normal, that is, the mapping between constraint space velocity and constraint space impulse at contact i, as defined in Equation 71, as follows:
 
 E   ii   =D   ii   −1 =(Σ j=1   n   J   ij   M   j   −1   J   ij   T ) −1   (71)
 
     Rather than implementing Equation 71 directly, applications typically update velocities at each iteration in order to reduce storage and computation costs, wherein velocities are determined in Equation 72, as follows:
 
 z   r+1 =( z   r   −E ( q+Jv   r )) + 
 
 v   r+1   =v   r   +M   −1   J   T ( z   r+1   −z   r ).  (72)
 
     Convergence is guaranteed if 2D−N is positive definite, but many rigid body systems don&#39;t converge with Jacobi. 
     Modified Jacobi Methods 
     There are a number of ways to modify Jacobi to guarantee convergence. One possibility is to add a line search step, as provided in Equation 73, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       z 
                       
                         r 
                         + 
                         0.5 
                       
                     
                     = 
                     
                       
                         ( 
                         
                           
                             z 
                             r 
                           
                           - 
                           
                             
                               D 
                               
                                 - 
                                 1 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   Nz 
                                   r 
                                 
                                 + 
                                 q 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                       + 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       d 
                       r 
                     
                     = 
                     
                       
                         z 
                         
                           r 
                           + 
                           0.5 
                         
                       
                       - 
                       
                         z 
                         r 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     α 
                     = 
                     
                       
                         
                           arg 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           min 
                         
                         
                           
                             0 
                             ≤ 
                             α 
                           
                           
                             0 
                             ≤ 
                             
                               z 
                               + 
                               
                                 α 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   d 
                                   r 
                                 
                               
                             
                           
                         
                       
                       ⁢ 
                       
                         f 
                         ⁡ 
                         
                           ( 
                           
                             z 
                             + 
                             
                               α 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 d 
                                 r 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       f 
                       ⁡ 
                       
                         ( 
                         z 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           1 
                           2 
                         
                         ⁢ 
                         
                           z 
                           T 
                         
                         ⁢ 
                         Nz 
                       
                       + 
                       
                         
                           z 
                           T 
                         
                         ⁢ 
                         q 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       z 
                       
                         r 
                         + 
                         1 
                       
                     
                     = 
                     
                       
                         z 
                         r 
                       
                       + 
                       
                         α 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           
                             d 
                             r 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   73 
                   ) 
                 
               
             
           
         
       
     
     Another method is to scale D −1  using ρ(N), wherein the spectral radius of N is defined in Equation 74, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       z 
                       
                         r 
                         + 
                         1 
                       
                     
                     = 
                     
                       
                         ( 
                         
                           
                             z 
                             r 
                           
                           - 
                           
                             ω 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 D 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     Nz 
                                     r 
                                   
                                   + 
                                   q 
                                 
                                 ) 
                               
                             
                           
                         
                         ) 
                       
                       + 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     ω 
                     &lt; 
                     
                       2 
                       
                         ρ 
                         ⁡ 
                         
                           ( 
                           N 
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   74 
                   ) 
                 
               
             
           
         
       
     
     Another solution proposes a Jacobi modification for both cloth simulation that can be adapted to the rigid body iteration of Equation 72. Let n b     β    be the number of bodies contacting body β and divide the impulse applied to body β by n n     β   , as follows in Equation 75:
 
 z   r+1 =( z   r   −E ( q+Jv   r )) + 
 
 v   r+1   =v   r   −M   −1 diag( n   1   −1   , . . . ,n   n   −1 ) J   T ( z   r+1   −z   r ).  (75)
 
     However, this solution does not specify the convergence conditions for the method and has not been proved. 
     Projected Gauss-Seidel 
     PGS solves each constraint individually in sequence. Let L be the (strict) lower triangle of N. Abstractly, PGS generates the following iteration, as defined in Equation 76, as follows:
 
 z   r+1 =( z   r   −D   −1 ( q+Lz   r+1 +( N−L ) z   r )) +   (76)
 
     Applications generally use an equivalent iteration that updates the velocity after processing each constraint, as defined in Equation 77:
 
 z   i   r+1 =( z   i   r   −E   ii ( q   i   +J   i   v   rm+i−1 )) + 
 
 v   rm+i   =v   rm+i−1   +M   −1   J   i   T ( z   i   r+1   −z   i   r ).  (77)
 
     The size of z is m (the number of contacts), and the indexing of v is arranged so that the velocity used to calculate the first element of each iteration (i=1) is the velocity calculated after the last element (i=m) of the previous iteration. 
     Parallel Projected Gauss-Seidel 
     Ideally, processing of each contact is conducted in parallel on a separate thread. However, each body can be in contact with more than one other body, so if PGS were parallelized naively then multiple threads could try to update the same body simultaneously, causing velocity updates to be lost and convergence to be affected. As previously described in relation to  FIG. 5 , a common solution is contact coloring, wherein partitions (or colors) are assigned to the contacts so that each body is referenced at most once in each partition. The PGS solver can then process each color sequentially, processing the contacts within each color in parallel. 
     The number of colors required depends on Δ, the maximum number of contacts per body over all the bodies. From Vizing&#39;s theorem, the minimum number of colors required is either Δ, or Δ+1. Finding the optimal coloring is NP-hard so a parallel greedy algorithm is used in one embodiment, which gives at most 2Δ-1 colors. In piles of objects, Δ can be quite high and generally, the first few colors contain a large proportion of the contacts, but the following colors contain fewer and fewer contacts. As the colors must be processed sequentially, the time taken by each iteration is multiplied by the number of colors, reducing the speedup of parallel PGS versus serial PGS. Also, on architectures (such as GPUs) with hundreds of cores, most of them are idle for most of the iteration. Embodiments of the present invention harness this unused processing power to reduce jitter. 
     Mass Splitting Solver 
     Embodiments of the present invention distinctly take the Jacobi algorithm and divide each body mass term in the effective mass (as defined in Equation 71) by the number of contacts affecting it, but use the full mass to apply the impulses. Let n i  be the number of contacts involving body i, then the mass splitting algorithm is provided in Equation 78, as follows:
 
 z   r+1 =( z   r   −E   S ( q+Jv   r )) + 
 
 v   r+1   =v   r   +M   −1   J   J ( z   r+1   −z   r )
 
 E   ii   S =(Σ j=1   n   n   j   J   ij   M   j   −1   J   ij   T ) −1   (78)
 
     This is different from the method outlined in Equation 75 where the impulse application is scaled, but the effective mass uses the full mass. 
     Block Mass Splitting Solver 
     A standard method for parallelizing solvers is to partition a system into blocks, solving the blocks on separate threads using one method and then combining the results in an outer iteration, possibly using a different method. Embodiments of the present invention improve upon the convergence of the basic method by solving contacts in blocks, using PGS to solve the contacts within the blocks, and projected Jacobi to combine the blocks. Solving all the contacts in a block on a single thread allows for the use of serial PGS, which converges faster than Jacobi. The collision detection system typically generates more than one (point) contact to represent the area of contact between each body pair, so a natural choice is to partition the contacts by body pair, as is shown in  FIG. 11  by way of an example. In  FIG. 11 , the blocks shown in the top left corner show a system  1100  of three bodies with four contacts. The contacts are grouped into two blocks of two. In the split system, body c is split into two sub-bodies, on contacting body a and one contacting body b. In the top right, a contact Jacobian  1120  is shown for the original system. On the bottom, a contact Jacobian  1130  is shown for the split system. Again the effective mass is modified to accelerate convergence, but now the mass of each body β is divided by n b     β    the number of body pairs it belongs to, as shown in Equation 79: 
     
       
         
           
             
               
                 
                   
                     z 
                     α 
                     
                       r 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 z 
                                 α 
                                 r 
                               
                               - 
                               
                                 
                                   
                                     E 
                                     ^ 
                                   
                                   α 
                                 
                                 ( 
                                 
                                   
                                     q 
                                     α 
                                   
                                   + 
                                   
                                     
                                       L 
                                       α 
                                     
                                     ⁢ 
                                     
                                       z 
                                       α 
                                       
                                         r 
                                         + 
                                         1 
                                       
                                     
                                   
                                   + 
                                   
                                     
                                       ( 
                                       
                                         
                                           N 
                                           α 
                                         
                                         - 
                                         
                                           L 
                                           α 
                                         
                                       
                                       ) 
                                     
                                     ⁢ 
                                     
                                       z 
                                       α 
                                       r 
                                     
                                   
                                   + 
                                   
                                     
                                       ∑ 
                                       
                                         γ 
                                         ≠ 
                                         α 
                                       
                                       
                                           
                                       
                                     
                                     ⁢ 
                                     
                                       
                                         N 
                                         αγ 
                                       
                                       ⁢ 
                                       
                                         z 
                                         γ 
                                         r 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                             ) 
                           
                           + 
                         
                         ⁢ 
                         
                           
 
                         
                         [ 
                         
                           
                             E 
                             ^ 
                           
                           α 
                         
                         ] 
                       
                       ii 
                     
                     = 
                     
                       
                         
                           ( 
                           
                             
                               ∑ 
                               
                                 β 
                                 = 
                                 1 
                               
                               n 
                             
                             ⁢ 
                             
                               
                                 
                                   
                                     n 
                                     
                                       b 
                                       β 
                                     
                                   
                                   ⁡ 
                                   
                                     [ 
                                     
                                       J 
                                       αβ 
                                     
                                     ] 
                                   
                                 
                                 i 
                               
                               ⁢ 
                               
                                 
                                   
                                     M 
                                     β 
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁡ 
                                   
                                     [ 
                                     
                                       J 
                                       αβ 
                                     
                                     ] 
                                   
                                 
                                 i 
                                 T 
                               
                             
                           
                           ) 
                         
                         
                           - 
                           1 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   79 
                   ) 
                 
               
             
           
         
       
     
     Here, α and γ are used to index contact blocks and β is used to index bodies. [M αβ ] i,j  represents element i, j of block α, β of matrix M. In particular, [J αβ ] i  represents the effect of contact i from contact block α on body β. Derivation of this solver is provided below, and Algorithm 11 provides an exemplary pseudocode, in embodiments of the present invention. 
     Derivation and Convergence 
     Unlike the scalar algorithm, for the block solver, the original system is transformed into a system of split masses and fixed joints. As described below, solving this system with PGS interleaved with an exact fixed joint solver converges and is operationally equivalent to running the block solver on the original system. Given that this interleaved solver provably converges, then the block algorithm must also converge. As such, solutions of the split system are also solutions of the original system, as will be further described below. This also demonstrates why scaling the terms of the effective mass, but not the applied impulses is proper, as previously determined. 
     The split system is constructed by splitting each original body into sub-bodies, one for each contacting body. Each sub-body has the same spatial extent as the original body, but the mass of the original body is divided equally between the sub-bodies. Fixed joints are added to ensure that the position and orientation of all the sub-bodies are equal. See  FIG. 6  for an example. 
     The key insight is that fixed joints have a closed form solution: the final velocity of every sub-body should be the average of all the sub-body velocities. 
     Splitting the Bodies 
     Below is described in detail how to construct the split system from the original system. Let n b     β    be the number of bodies contacting original body β. The body β is split into n b     β    sun-bodies, giving a total of sub-bodies described in Equation 80, as follows: 
     
       
         
           
             
               
                 
                   
                     n 
                     s 
                   
                   = 
                   
                     
                       ∑ 
                       
                         β 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                       
                         n 
                         
                           b 
                           β 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   80 
                   ) 
                 
               
             
           
         
       
     
     The mass and inertia of body β is split evenly between its sub-bodies, so the inverse mass matrix of the sub-bodies of body β is described in Equation 81, as follows:
 
 W   β   n   b     β   diag( M   β   −1   , . . . ,M   β   −1 )ϵ   6nb     β     6nb     β     (81)
 
and the inverse mass matrix of the whole split system is described in Equation 82, as follows:
 
 W =diag( W   1   , . . . ,W   n ).  (82)
 
     Each sub-body is just a portion of the mass of the original body, so the initial velocity of each sub-body V β  is set to the initial velocity of the original body β, giving whole system velocity v s , as is shown in Equation 83: 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       β 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 v 
                                 β 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 v 
                                 β 
                               
                             
                           
                         
                         ] 
                       
                       ∈ 
                       
                         
                           6 
                           ⁢ 
                           
                             n 
                             i 
                           
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       v 
                       s 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 V 
                                 1 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 V 
                                 n 
                               
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   83 
                   ) 
                 
               
             
           
         
       
     
     The unconstrained velocity, v s0  is constructed from v s , in the same way. 
     Assigning the Contacts to the Split Bodies 
     Let the total number of contact blocks (body pairs) be p and the number of contacts in block i be m i . The original Jacobian is partitioned into p×n blocks, such that J αβ ϵ   m     f     ×6  represents the effect of the contacts in contact block α on body β, as is shown in Equation 84: 
     
       
         
           
             
               
                 
                   J 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               J 
                               11 
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               
                                 1 
                                 ⁢ 
                                 n 
                               
                             
                           
                         
                         
                           
                             ⋮ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               J 
                               
                                 p 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               J 
                               pn 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   84 
                   ) 
                 
               
             
           
         
       
     
     In each column block β of J, there are n b     β    non-zero blocks, corresponding to the contact blocks that affect body β. Let r j,β  be the index of the jth non-zero block in column block β of J. The contact Jacobian is built for the split system, C, in such a way that each sub-body has exactly one contact block, as is shown in Equation 85: 
     
       
         
           
             
               
                 
                   C 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               C 
                               11 
                             
                           
                           
                             … 
                           
                           
                             
                               C 
                               
                                 1 
                                 ⁢ 
                                 n 
                               
                             
                           
                         
                         
                           
                             ⋮ 
                           
                           
                             
                                 
                             
                           
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               C 
                               
                                 p 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               C 
                               pn 
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   85 
                   ) 
                 
               
             
           
         
       
     
     Each block C αβ  has n b     β    sub-blocks, one for each sub-body of β, and [C αβ ] 1,j  represents the effect of contact block α on sub-body j of original body β. The jth non-zero contact block of body β in J is assigned to sub-body j of body β in C, as described in Equation 86: 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         C 
                         αβ 
                       
                       ] 
                     
                     
                       1 
                       , 
                       j 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             J 
                             αβ 
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               α 
                             
                             = 
                             
                               r 
                               
                                 j 
                                 , 
                                 β 
                               
                             
                           
                         
                       
                       
                         
                           0 
                         
                         
                           
                             otherwise 
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   86 
                   ) 
                 
               
             
           
         
       
     
     Let b α,1  and b α,2  be the bodies constrained by contact block α and s α,1  be the rank of contact block α in its first body&#39;s constraint list, and s α,2  be the rank in its second, so that Equation 87 defines: 
     
       
         
           
             
               
                 
                   
                     J 
                     αβ 
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               [ 
                               
                                 C 
                                 αβ 
                               
                               ] 
                             
                             
                               1 
                               , 
                               
                                 s 
                                 
                                   α 
                                   , 
                                   1 
                                 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   1 
                                 
                               
                             
                             = 
                             β 
                           
                         
                       
                       
                         
                           
                             
                               [ 
                               
                                 C 
                                 αβ 
                               
                               ] 
                             
                             
                               1 
                               , 
                               
                                 s 
                                 
                                   α 
                                   , 
                                   2 
                                 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 b 
                                 
                                   α 
                                   , 
                                   2 
                                 
                               
                             
                             = 
                             
                               β 
                               . 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   87 
                   ) 
                 
               
             
           
         
       
     
     For an example of how to construct C from J, see  FIG. 11 . 
     Joining the Bodies Back Together 
     A fixed joint forces the position and orientation of two bodies to be the same, as shown in Equation 88: 
     
       
         
           
             
               
                 
                   
                     
                       
                         Φ 
                         FJ 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             
                               
                                 x 
                                 1 
                               
                             
                           
                           
                             
                               
                                 x 
                                 2 
                               
                             
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           x 
                           1 
                         
                         - 
                         
                           x 
                           2 
                         
                       
                       = 
                       0 
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   88 
                   ) 
                 
               
             
           
         
       
     
     and has Jacobian in Equation 89: 
     
       
         
           
             
               
                 
                   
                     J 
                     FJ 
                   
                   = 
                   
                     
                       
                         ∂ 
                         
                           Φ 
                           FJ 
                         
                       
                       
                         ∂ 
                         x 
                       
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               I 
                             
                             
                               
                                 - 
                                 I 
                               
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   89 
                   ) 
                 
               
             
           
         
       
     
     For each body β a number of fixed joints n b     β   −1 are introduced to join the sub-bodies together, giving a total of 
     
       
         
           
             
               
                 
                   
                     m 
                     f 
                   
                   = 
                   
                     
                       ∑ 
                       
                         β 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             n 
                             
                               b 
                               β 
                             
                           
                           - 
                           1 
                         
                         ) 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   90 
                   ) 
                 
               
             
           
         
       
     
     The Jacobian of the fixed joints of body β is defined in Equation 91, as follows: 
     
       
         
           
             
               
                 
                   
                     
                       F 
                       β 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                             
                               0 
                             
                             
                               … 
                             
                             
                               0 
                             
                           
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               0 
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                             
                               
                                   
                               
                             
                             
                               ⋮ 
                             
                           
                           
                             
                               ⋮ 
                             
                             
                               ⋮ 
                             
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               0 
                             
                           
                           
                             
                               
                                 I 
                                 6 
                               
                             
                             
                               0 
                             
                             
                               … 
                             
                             
                               0 
                             
                             
                               
                                 - 
                                 
                                   I 
                                   6 
                                 
                               
                             
                           
                         
                         ] 
                       
                       ∈ 
                       
                         ℝ 
                         
                           6 
                           ⁢ 
                           
                             n 
                             
                               b 
                               β 
                             
                           
                           × 
                           6 
                           ⁢ 
                           
                             ( 
                             
                               
                                 n 
                                 
                                   b 
                                   β 
                                 
                               
                               - 
                               1 
                             
                             ) 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   91 
                   ) 
                 
               
             
           
         
       
     
     The Jacobian of all fixed joints in the split system is defined in Equation 92, as follows:
 
 F =diag( F   1   , . . . ,F   m     f   ).  (92)
 
The Split System
 
     The split contact constraints and the fixed joints are concatenated to form system Jacobian [C T  F T ] T , and Equations 93 and 94 are defined, as follows: 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             N 
                             cc 
                           
                         
                         
                           
                             N 
                             cf 
                           
                         
                       
                       
                         
                           
                             N 
                             fc 
                           
                         
                         
                           
                             N 
                             ff 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             C 
                           
                         
                         
                           
                             F 
                           
                         
                       
                       ] 
                     
                     ⁢ 
                     
                       
                         W 
                         ⁡ 
                         
                           [ 
                           
                             
                               
                                 C 
                               
                             
                             
                               
                                 F 
                               
                             
                           
                           ] 
                         
                       
                       T 
                     
                   
                 
               
               
                 
                   ( 
                   93 
                   ) 
                 
               
             
             
               
                 
                   
                     D 
                     cc 
                   
                   = 
                   
                     
                       diag 
                       ⁡ 
                       
                         ( 
                         
                           N 
                           cc 
                         
                         ) 
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   94 
                   ) 
                 
               
             
           
         
       
     
     Impulses [z c   T z f   T ] T  are found that satisfy the contacts ad fixed joints simultaneously by solving the following mixed linear complementarity problem (MLCP), as shown in Equation 95.
 
 q+N   cc   T   +N   cf   T   z   f ≥0 ⊥z   c ≥0
 
 N   fc   T   z   c   +N   ff   T   z   f =0.  (95)
 
     This MLCP is solved by interleaving PGS iterations for the contacts with direct (exact) solves for the joints, in one embodiment, as shown in Equation 96:
 
 z   c   r+1 =( z   c   r   −D   cc   −1 ( q+N   cc   z   c   +N   cf   z   f )) + 
 
 z   f   r+1 =linsolve( N   ff   ,N   fc   z   c ).  (96)
 
     The iteration is transformed to a velocity iteration, as follows in Equations 97-100:
 
 v   s   0,1   =v   s0   (97)
 
 z   c     i     r+1 =( z   c     i     r   −D   cc     β     −1 ( q   i   +C   i   v   s   r,i )) +   (98)
 
 v   s   r,i+1   =v   s   v,i   +WC   i   T ( z   c     i     r+1   −z   c     i     r )  (99)
 
 v   s   r+1,1   =v   s   r,m+1   +WF   T  linsolve( N   ff   ,Fv   s   r,m+1 ).  (100)
 
     Equations 98 and 99 apply a single PGS iteration to all the contacts at each iteration. Each sub-body has exactly one constraint, so no coloring is required and the contacts can be safely computed concurrently on m threads. 
     Closed Form Solution for Fixed Joints 
     Equation 100 is solved using Equation 101: 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         
                           ( 
                           
                             v 
                             s 
                           
                           ) 
                         
                         β 
                       
                       ] 
                     
                     i 
                     
                       
                         r 
                         + 
                         1 
                       
                       , 
                       1 
                     
                   
                   = 
                   
                     
                       average 
                       ( 
                       
                         
                           
                             [ 
                             
                               
                                 ( 
                                 
                                   v 
                                   s 
                                 
                                 ) 
                               
                               β 
                             
                             ] 
                           
                           1 
                           
                             r 
                             , 
                             
                               m 
                               + 
                               1 
                             
                           
                         
                         ⁢ 
                         
                           
                             … 
                             ⁢ 
                             
                                 
                             
                             [ 
                             
                               
                                 ( 
                                 
                                   v 
                                   s 
                                 
                                 ) 
                               
                               β 
                             
                             ] 
                           
                           
                             n 
                             
                               b 
                               β 
                             
                           
                           
                             r 
                             , 
                             
                               m 
                               + 
                               1 
                             
                           
                         
                       
                       ) 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   101 
                   ) 
                 
               
             
           
         
       
     
     In other words, the fixed joints can be enforced simply by averaging the sub-body velocities of each body at each iteration. 
     Algorithm 11 provides pseudocode for the block mass splitting solver, below: 
                             Algorithm 11 Block Mass Splitting Solver                                             1: for all contact blocks, α do            2:  for all contacts in block, i do            3:   [p α ] i  = 0            4:   for all β ∈ {b α,1 ,b α,2 } do            5:    [p α ] i  = [p α ] i  + [J αβ ] i (n b     β   M β   −1 )[J βα ] i   T              6:   end for            7:   p α  = p α   −1              8:  end for            9: end for           10: for all iterations do           11:  for all contact blocks, α do {in parallel}           12:   β 1 , β 2  = b α1 , b α2             13:   D α1 , D α2  = v β     1   , v β     2               14:   for all contacts in block, i do {sequentially on thread α}           15:    t = [z Cα ] i             16:    [z Cα ] i  = max (0, [z Cα ] i  − [p α ] i (J α,β1  D α1 + J α,β2  D α2 ))           17:    for all j ∈ {1,2} do            18:     D αj  = D αj  + M β     j     −1 J α,βj   T  (z Ca  − t)           19:    end for           20:    end for           21:   end for           22:  a = 0           23:  for all j ∈ 1,2} do            24:   for all α such that by b αj  = β do           25:     a = a +D αj             26:    end for           27:  end for           28: end for                    
Friction
 
     In one embodiment, a pyramid approximation to the coulomb friction cone is implemented. To add this friction approximation to the discretized system, boxed LCPs (BLCPs) are represented, which have both upper and lower limits. Let x:=BLCP(A, b, 1, h) be defined as follows in Equation 102:
 
find  xϵ     m  such that, for  all i=l . . . m,  
 
 x   i   l   i  and ( Ax+b ) i ≥0 or
 
 x   i   =h   i  and ( Ax+b ) i ≤0 or
 
 l   i   &lt;x   i   &lt;h   i  and ( Ax+b ) i =0.  (102)
 
     All of the LCP algorithms used in this Application can be turned into BLCP algorithms by changing max(0, x i ) to max(l i , min(x i , h)). Let D be the Jacobian matrix of the discretized friction pyramid and let N fric =DM −1 D T . Now the discretized model can be written with friction as a pair of coupled complementarity problems, as provided in Equations 103 and 104:
 
 z=LCP ( N   s   J ( v   0   +M   −1   D   T   z   fric ))  (103)
 
 z   fric   =BLCP ( N   fric   ,D ( v   0   +M   −1   J   T   z ),−diag(μ) z ,diag(μ) z ).  (104)
 
     Thus, according to embodiments of the present disclosure, systems and methods are described that provide for a parallel iterative rigid body solver and method for implementing the same that uses a modified effective mass based on the number of contacts affecting a corresponding body to simulate complex systems coming to rest while avoiding artifacts such as, jittering and swimming at very low iteration counts. The linear complementarity solvers and methods for implementing the same for formulating rigid body dynamics provide for correct momentum propagation, static friction, dynamic friction, and stable resting contact. 
     While the foregoing disclosure sets forth various embodiments using specific block diagrams, flowcharts, and examples, each block diagram component, flowchart step, operation, and/or component described and/or illustrated herein may be implemented, individually and/or collectively, using a wide range of hardware, software, or firmware (or any combination thereof) configurations. In addition, any disclosure of components contained within other components should be considered as examples because many other architectures can be implemented to achieve the same functionality. 
     The process parameters and sequence of steps described and/or illustrated herein are given by way of example only and can be varied as desired. For example, while the steps illustrated and/or described herein may be shown or discussed in a particular order, these steps do not necessarily need to be performed in the order illustrated or discussed. The various example methods described and/or illustrated herein may also omit one or more of the steps described or illustrated herein or include additional steps in addition to those disclosed. 
     While various embodiments have been described and/or illustrated herein in the context of fully functional computing systems, one or more of these example embodiments may be distributed as a program product in a variety of forms, regardless of the particular type of computer-readable media used to actually carry out the distribution. The embodiments disclosed herein may also be implemented using software modules that perform certain tasks. These software modules may include script, batch, or other executable files that may be stored on a computer-readable storage medium or in a computing system. These software modules may configure a computing system to perform one or more of the example embodiments disclosed herein. One or more of the software modules disclosed herein may be implemented in a cloud computing environment. Cloud computing environments may provide various services and applications via the Internet. These cloud-based services (e.g., software as a service, platform as a service, infrastructure as a service, etc.) may be accessible through a Web browser or other remote interface. Various functions described herein may be provided through a remote desktop environment or any other cloud-based computing environment. 
     The foregoing description, for purpose of explanation, has been described with reference to specific embodiments. However, the illustrative discussions above are not intended to be exhaustive or to limit the invention to the precise forms disclosed. Many modifications and variations are possible in view of the above teachings. The embodiments were chosen and described in order to best explain the principles of the invention and its practical applications, to thereby enable others skilled in the art to best utilize the invention and various embodiments with various modifications as may be suited to the particular use contemplated. 
     Embodiments according to the present disclosure are thus described. While the present disclosure has been described in particular embodiments, it should be appreciated that the disclosure should not be construed as limited by such embodiments, but rather construed according to the below claims.