Real-time Computational Algorithms

Contact-Response Algorithms

Iterative Predictor Corrector (IPC) algorithm

When simulating the contact response between two bodies, such as a rigid surgical tool contacts the soft tissue surface, in addition to elastodynamic equations associated with the simulation of the tissue, a system of inequalities representing the kinematic contact constraints needs to be solved as well. This is known as a Mixed Linear Complementarity Problem (MLCP). A standard way of solving MLCP proceeds by converting the MLCP into a standard linear complementarity problem (LCP) and then solving the LCP using an interactive solver like projected Gauss-Seidel method. This approach works well for rigid body contact but is computationally costly for modeling elastic contact occurred in deformable objects. To address this, the Iterative Constraint Anticipation (ICA) technique [OtTa09] uses a staggered iteration approach in which the linear equalities and complementarity inequalities are solved in a successive manner. The cost of ICA is O(n+m) per iteration, where nis the degree of freedom and m is the number of contacts. However, the ICA is not efficient when friction is considered. We have developed an Iterative Predictor Corrector (IPC) algorithm [ArDe15] that avoids modeling of the friction cone but instead uses the knowledge of intermediate solution of the MLCP solver to determine the unknowns of a frictional contact.

Modified Iterative Constraint Anticipation (MICA) algorithm

Although Linear projection constraints (LPC) [BaWi98] are commonly used for collision response, they are prone to sticking artifacts. To enforce no slip conditions between surfaces in static friction, we therefore enforce the immobility for all the nodes in static friction state using projection operators which are then incorporated into the ICA method, resulting in a modified interactive constraint anticipation (MICA) algorithm. To accelerate convergence, we further extend the MLCP method to a multilevel framework - Modified Multilevel ICA (M-MgICA) [Arikatla14], which demonstrates faster convergence compared to a single level MICA.

Contact-Response Algorithms image

Left: snapshots of an elastic block resting on a planar surface and pushed from left to slide against friction. Right: Convergence comparison of single and 3-level MICA.


  • [OtTa09] M. A. Otaduy, R. Tamstorf, D. Steinemann, and M. Gross, “Implicit contact handling for deformable objects,” Comput. Graph. Forum, vol. 28, no. 2, pp. 559–568, 2009.
  • [ArDe15] V. S. Arikatla and S. De, “An iterative predictor-corrector approach for modeling static and kinetic friction in interactive simulations,” Graph. Models, vol. 82, pp. 29–42, 2015.
  • [BaWi98] D. Baraff and A. Witkin, “Large steps in cloth simulation,” in Proceedings of the 25th annual conference on Computer graphics and interactive techniques  - SIGGRAPH ’98, 1998, pp. 43–54.
  • [Arikatla14] V. S. Arikatla, “Algorithms for Imposition of Constraints and Topology Changes in Interactive Simulations”, PhD. thesis, 2014.

Point Associated Finite Field

This is a meshfree computational environment based on the moving least squares approximation functions, compactly supported on spherical subdomains, used in a point collocation residual minimization technique. Advantages over traditional finite elements include not having to perform numerical integration, constant Jacobians which preclude distortions and smooth solutions.This is a meshfree computational environment based on the moving least squares approximation functions, compactly supported on spherical subdomains, used in a point collocation residual minimization technique. Advantages over traditional finite elements include not having to perform numerical integration, constant Jacobians which preclude distortions and smooth solutions.

Point Associated Finite Field image

Dynamic Point

Dynamic Point image

In real time virtual reality applications with haptic feedback, update rates have to be maintained close to 1000 Hz. In such applications the interactions are modeled as a point contact. In applications such as virtual surgery, long slender tools are used to interact with the organs and are better modeled as a line rather than a point interaction. Dynamic Point is a novel method that can compute efficiently the collision of line tools with rigid and deformable meshes. This method has constant computational complexity for arbitrary mesh sizes.

Constrained Particle Dynamics

Constrained Particle Dynamics (CPD) is a method for real-time simulation of a system of particles and continua, based on explicit and implicit time integration of the governing set of dynamical equations.

Most of the popular methods for the simulation of dynamic systems in computer graphics are force based. Accelerations are computed from forces and then velocities and positions are integrated. The advantage of CPD is that it  works directly on positions. The instability introduced by explicit time integration is avoid. In addition, all constraints can be handled in an uniform way, including the collision constraints. The objective is to establish a set of techniques that are as simple, general, robust and efficient as the position-based dynamics, but are physically accurate.

Constrained Cutting Methods
Two frames of real-time cloth simulation. FPS: ~60; Mesh: 2500 vertices; number of Constraints: 15000(7600 distance + 7400 dihedral).


Constrained Spline Dynamics

Constrained spline dynamics utilizes a smooth discretization based on b-splines for real-time simulation of rod-like objects under dynamic loading. As against discrete methods, the formulation enables accurate simulation of rods with much fewer discretization points, while exhibiting better accuracy in capturing the strain fields such as the curvature and bi-normal vector fields even in scenarios of complex bent configurations. The discretization is formulated using the nodal-quantities as degrees of freedom, thus having a simple and robust collision treatment.

In-plane and out-of-plane deformed configurations under bending imposed by constraints simulated with 12 sub-divisions.


Physics-driven neural networks-based simulation system (PhyNNeSS)

PhyNNeSS image

The primary goal of PhyNNess is to predict mechanical properties of soft tissues, based on neural networks. The basic idea of this research is to implement elastography using deep neural networks. Mechanical properties in this research, for example, can be Young's modulus and shear modulus of soft tissues. Because tumors are generally harder than normal tissues, elasticity map of tissues delivers useful information for distinguishing normal tissues and tumors. To train neural networks, ultrasound images (B-mode images) is used because how the tissue is deformed can be observed from ultrasound images of tissues. The displacement field of tissues then serve as inputs to the neural network and output of the neural network is mechanical properties of interest.

Voxel-based cutting methods

An Interactive voxel-based cutting algorithm for deformable objects

Cutting simulation of deformable objects has been an active research topic in physics-based simulation for decades. Efficient and robust incorporating cuts into deformable objects is a challenging task. 

To allow for real-time simulation of deformable objects with changing topology, introduced by cuts for instance, we proposed an efficient voxel-based cutting technique for deformable objects. Besides cutting, this structure is also exploited to simulate real-time tissue interactions like pressing and dragging scenarios.


Multigrid methods

Algorithms for cutting of deformable objects simulated using multilevel methods

Multilevel cutting is complicated by the need to reflect fine grid topology changes at the coarser grid level which requires modification of the coarse grid operators. We consider a two-level formulation and address the real time issue by developing two separate strategies for updating the map relating the coarse and fine grid operators for structured as well as unstructured meshes.  Considerable speedups of this multigrid approach compared to a preconditioned conjugate gradient (PCG) solver are observed.

Multigrid algorithms for multiphysics problem

A taxonomy of virtual surgery has been divided into five different generations, namely systems representing accurately the geometry of the organs at a macroscopic level, physical dynamics of the tissue at a continuum level, physiology, microscopic phenomenon and biochemical phenomenon. The higher the generation for virtual surgery, the more physical fidely, yet the more complex of the problem. The system matrix of this type of coupled multiphysics problem is generally a big sparse matrix with a special block structure representing each of the underlying phyiscs and frequently has a very large condition number. Current best practices for the solution of approach for this type of multiphysics problem are almost exclusively based on multigrid methods and preconditioned Krylov subspace iterative methods. Multigrid methods have optimal arithmetic complexity, ideally with O(N) floating point operations and therefore applicable to multiphysics problem with large number of degrees of freedom.

Left: Monolithic blocked system matrix consisting of three physics; Middle: three level multigrid hierarchy; Right: linear scaling of algebraic and geometric multigrid methods