This paper presents the most important methods used in computer graphics to produce realistic animations and simulations: finite element/difference/volume methods, mass-spring systems, meshfree methods, coupled particle systems and reduced deformable models based on modal analysis. Also, it discusses the benefits and drawbacks of each method, reports the improvement that made in other papers.
Time Integration
Time integration schemes are evaluated by two main criteria, their stability and their accuracy. In the field of physically based animation in Computer Graphics, stability is often much more important than accuracy.
$$ \begin{align*} \mathbf{x}(t+\Delta t) &= \mathbf{x}(t) + \Delta t \mathbf{v}(t) \\ \mathbf{v}(t+\Delta t) &= \mathbf{v}(t) + \Delta t F(\mathbf{v}(t), \mathbf{x}(t), t) \end{align*} $$Explicit methods are easy to implement but they are only conditionally stable, i.e. stable only if ∆t is smaller than a stability threshold. This problem can be solved by using implicit integration
$$ \begin{align*} \mathbf{x}(t+\Delta t) &= \mathbf{x}(t) + \Delta t \mathbf{v}(t + \Delta t) \\ \mathbf{v}(t+\Delta t) &= \mathbf{v}(t) + \Delta t F(\mathbf{v}(t + \Delta t), \mathbf{x}(t + \Delta t), t) \end{align*} $$This method is stable for large time steps $\Delta t$ but with the price of having to solve an algebraic system of equation at each time step.
$$ \begin{align*} \mathbf{v}(t+\Delta t) &= \mathbf{v}(t) + \Delta t F(\mathbf{v}(t), \mathbf{x}(t), t) \\ \mathbf{x}(t+\Delta t) &= \mathbf{x}(t) + \Delta t \mathbf{v}(t + \Delta t) \end{align*} $$This method is still explicit but more stable than forward Euler integration, without any additional computational overhead.
Lagrangian Mesh Based Methods
Finite Element Method
$$ \begin{equation}\label{fem} \rho \ddot{\mathbf{x}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f} \end{equation} $$$$ \tilde{\mathbf{x}}(\mathbf{m}, t) = \sum_i \mathbf{x}_i(t) \mathbf{b}_i(\mathbf{m}) $$$$ E = \int_V \boldsymbol{\epsilon}(\mathbf{m}) \cdot \boldsymbol{\sigma}(\mathbf{m}) d \mathbf{m} $$$$ \begin{equation}\label{nodal} \mathbf{f}_e = \mathbf{K}_e \mathbf{u}_e \end{equation} $$$$ \rho \ddot{\mathbf{x}} = \mu \Delta \mathbf{u} + (\lambda +\mu) \nabla (\nabla \cdot \mathbf{u}) $$Where the material constants $\lambda$ and $\mu$ can be computed directly from Young’s modulus and Poisson’s ratio.
$$ \mathbf{f} = \mathbf{RK}(\mathbf{R}^T \mathbf{x} - \mathbf{x}_0) $$Finite Differences Method
The method of Finite Differences is easier to implement than the general Finite Element Method. However, it is difficult to approximate the boundary of an arbitrary object with a regular mesh. Also, local adaptations are only possible with irregular meshes.
Finite Volume Method
In the explicit Finite Element Method, the forces acting on the nodes of an element are computed as the derivatives of the deformation energy with respect to the nodal positions.
$$ \mathbf{f}_A = \int_A \boldsymbol{\sigma} d \mathbf{A} $$$$ \mathbf{f}_A = A \boldsymbol{\sigma} \mathbf{n} $$Boundary Element Method
Roughly speaking, the integral form of the equation of motion is transformed into a surface integral by applying the Green-Gauss theorem. The method achieves substantial speedup because the three dimensional problem is reduced to two dimensions. However, the approach only works for objects whose interior is composed of a homogenous material. Also, topological changes (e.g. fractures) are more difficult to handle than in the explicit FEM approach where only local changes to the stiffness matrix or the connectivity of the elements need to be done.
Mass-Spring Systems
Instead of beginning with a PDE then discretizing in space, we begin directly with a discrete model.
$$ \mathbf{M} \ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{v}) $$where $\mathbf{M}$ is a $3n \times 3n$ diagonal mass matrix. Thus, mass-spring systems only require the solution of a system of coupled ODEs. This is done via time integration.
$$ \mathbf{f}_i = k_s(|\mathbf{x}_{ij}| - l_{ij}) \frac{\mathbf{x}_{ij}}{|\mathbf{x}_{ij}|} $$$$ \mathbf{f}_i = k_d (\mathbf{v}_j - \mathbf{v}_i) $$$$ \mathbf{f}_i = k_d \left( \frac{\mathbf{v}_{ij}^T \mathbf{x}_{ij}}{\mathbf{x}_{ij}^T \mathbf{x}_{ij}} \right) \mathbf{x}_{ij} $$This projects the velocity difference onto the vector separating the masses, and only allows a force along that direction.
Advantages:
- Particle systems tend to be intuitive and simple to implement.
- Computationally efficient and can handle large deformations with ease.
Limitations:
- Unlike the Finite Element and finite difference methods, which are built upon elasticity theory, mass-spring systems are not necessarily accurate.
- Most such systems are not convergent and the behavior is dependent on the mesh resolution and topology.
- Spring constants $k_s$ are often chosen arbitrarily so that one can say little quantitatively about the material being modeled.
- There is often coupling between the various spring types.
Time Integration:
$$ \begin{bmatrix} \Delta \mathbf{x} \\ \Delta \mathbf{v} \end{bmatrix} = \Delta t \begin{bmatrix} \mathbf{v}_n + \Delta \mathbf{v} \\ M^{-1} \mathbf{f}(\mathbf{x}_n +\Delta \mathbf{x}, \mathbf{v}_n +\Delta \mathbf{v}) \end{bmatrix} $$$$ \mathbf{f}(\mathbf{x}_n +\Delta \mathbf{x}, \mathbf{v}_n +\Delta \mathbf{v}) = \mathbf{f}_n + \frac{\partial \mathbf{f}}{\partial \mathbf{x}} \Delta \mathbf{x} + \frac{\partial \mathbf{f}}{\partial \mathbf{v}} \Delta \mathbf{v} $$$$ A \Delta \mathbf{v} = \mathbf{b} $$$$ \begin{align*} A &\equiv (I - \Delta t M^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{v}} - \Delta t^2 M^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}}) \\ \mathbf{b} &\equiv \Delta t M^{-1} (\mathbf{f}_n + \Delta t \frac{\partial \mathbf{f}}{\partial \mathbf{x}} \mathbf{v}_n) \end{align*} $$Lagrangian Mesh Free Methods
Loosely Coupled Particle Systems
Particles which interact with each other depending on their spatial relationship are referred to as spatially coupled particle systems. The interactions between particles evolve dynamically over time, thus, complex geometry and topological changes can be easily modeled using this approach.
$$ \phi_i = \sum_{j\ne i} \phi_{ij} $$$$ \mathbf{f}_i = -\nabla_{\mathbf{x}_i} \phi_i = -\sum_{j\ne i} \nabla_{\mathbf{x}_i} \phi_{ij} $$$$ \phi_{LJ}(d) = \frac{B}{d^n} - \frac{A}{d^m} $$$$ \phi(d) = \frac{-e_0}{m-n} \left( m \left( \frac{d_0}{d} \right)^n - n \left( \frac{d_0}{d} \right)^m \right) $$where $d_0$ is the equilibrium separation distance, and $−e_0$ is the minimal potential (the magnitude is called the dissociation energy). Increasing the dissociation energy increases the stiffness of the model, while with the exponents n and m the width of the potential well can be varied. Therefore, large dissociation energy and high exponents yield rigid and brittle material, while low dissociation energy and small exponents result in soft and elastic behavior of the object. This allows the modeling of a wide variety of physical properties ranging from stiff to fluid-like behavior.
Smoothed Particle Hydrodynamics
SPH is a Lagrangian technique where discrete, smoothed particles are used to compute approximate values of needed physical quantities and their spatial derivatives. Forces can be easily derived directly from the state equations. Furthermore, as a particle-based Lagrangian approach it has the advantage that mass is trivially conserved and convection is dispensable. This reduces both the programming and computational complexity and is therefore suitable for interactive applications. A drawback of SPH is that it is not possible so far to exactly maintain the incompressibility of material.
Do not actually understand this part
$$ A(\mathbf{x}) = \sum_j m_j \frac{A_j}{\rho_j} W(\mathbf{r}, h) $$where $\mathbf{r} = \mathbf{x} - \mathbf{x}_j$ and $W(\mathbf{r}, h)$ is a smoothing kernel with the properties $\int W(\mathbf{r}, h) d\mathbf{r} = 1$ and $\lim_{h\to 0} W(\mathbf{x}, h) = \delta(\mathbf{x})$, where $h$ is the support radius. For efficiency reasons, $h$ is usually chosen such that the kernel $W(\mathbf{r}, h)$ has compact support. The choice of smoothing kernel is very important. Often spline Gaussian kernels are used.
Point-Based Animation
$$ U = \frac{1}{2} (\epsilon \cdot \sigma) $$$$ \mathbf{f}_i = -\nabla_{\mathbf{u}_i} U = -\frac{1}{2} \nabla_{\mathbf{u}_i} (\epsilon_i \cdot E \cdot \epsilon_i) = -\sigma \cdot \nabla_{\mathbf{u}_i} \epsilon_i $$where Green’s nonlinear strain tensor is used for $\epsilon_i$, the spatial derivative of the displacement field $\nabla\mathbf{u}_i$ at $\mathbf{m}_i$ is needed. To guarantee zero elastic forces for rigid body modes, the approximation of $\nabla\mathbf{u}_i$ from the displacement vectors $\mathbf{u}_j$ of the neighboring particles must be at least first order accurate. Therefore, the Moving Least Squares Method with a linear basis is used which yields first order accurate interpolation of point-sampled functions (contrary to the SPH approximation which is not even zeroth-order accurate).
$$ \tilde{u}_j = u_i + \nabla u|_{\mathbf{m}_i} \cdot \mathbf{m}_{ij} $$$$ e_i = \sum_j (\tilde{u}_j - u_j)^2 w_{ij} $$$$ \left( \sum_j \mathbf{m}_{ij} \mathbf{m}_{ij}^T w_{ij} \right) \nabla u|_{\mathbf{m}_i} = \sum_j (u_j - u_i) \mathbf{m}_{ij} w_{ij} $$$$ \nabla u|_{\mathbf{m}_i} = \mathbf{A}^{-1} \left( \sum_j (u_j - u_i) \mathbf{m}_{ij} w_{ij} \right) $$where the inverse of $\mathbf{A} = \sum_j \mathbf{m}_{ij} \mathbf{m}_{ij}^T w_{ij}$ needs to be computed only once per particle and can be used for computing the derivation of $v$ and $w$ as well.
Linear Model Analysis
$$ \begin{equation}\label{motion} \mathbf{M} \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} = \mathbf{f}_{ext} \end{equation} $$$$ \mathbf{M} \Phi \Lambda = \mathbf{K} \Phi $$$$ \mathbf{u}(t) = \Phi \mathbf{q}(t) $$$$ \Phi^T \mathbf{M} \Phi \ddot{\mathbf{q}} + \Phi^T \mathbf{C} \Phi \dot{\mathbf{q}} + \Phi^T \mathbf{K} \Phi \mathbf{q} = \Phi^T \mathbf{f}_{ext} $$$$ \ddot{q}_i + (\alpha + \beta \lambda_i) \dot{q}_i + \lambda_i q_i = g_i $$where $\lambda_i$ is the $i$-th diagonal element of $\Lambda$.
We can solve this equation analytically.
Eulerian and Semi-Lagrangian Methods
Lagrangian techniques are a convenient way to define an object as a connected mesh of points or simply a cloud of points. The Eulerian point of view, on the other hand, looks at a stationary set of points and calculates how the material properties stored at those stationary grid points change over time. One of the drawbacks of changing to an Eulerian perspective is that the boundary of the object is no longer explicitly defined.
In Computer Graphics, fluids are often defined in an Eulerian framework.
$$ \begin{align} \nabla \cdot \mathbf{u} &= 0 \label{mass} \\ \mathbf{u}_t &= -(\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla \cdot (\nu \nabla \mathbf{u}) - \nabla p + \mathbf{f} \nonumber \end{align} $$These two equations represent the conservation of mass and momentum for an incompressible fluid. In this Eulerian formulation the quantities of the fluid are stored in a grid of cells. The velocity is stored on the cell faces and the pressure is stored at the center of the cells. This is commonly called a staggered (or MAC) grid. Notice that the position, $\mathbf{x}$, of the fluid is not defined; in this Eulerian frame work, the grid positions remain fixed.
Solving these two equations is usually done in several steps (breaking up a PDE into simpler pieces and solving the pieces separately is known as operator splitting). The $\mathbf{f}$ term is simply scaled by the time step and added to the current velocity. The advection term, $-(\mathbf{u} \cdot \nabla) \mathbf{u}$ is then solved. One popular way to solve the advection term is with the semi-Lagrangian technique because it is stable for large time step sizes.
$$ \tilde{\mathbf{u}} = \mathbf{u} + \Delta t (- (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla \cdot (\nu \nabla \mathbf{u}) + \mathbf{f}) $$$$ \begin{equation}\label{unew} \mathbf{u}^\text{new} = \tilde{\mathbf{u}} - \Delta t \nabla p \end{equation} $$$$ \Delta t \nabla^2 p = \nabla \cdot \tilde{\mathbf{u}} $$with which we can solve for $p$. Then plug it back $\eqref{unew}$ to find our final velocity. The final velocity is incompressible (divergence free), and for a constant density fluid this ensures mass conservation. Compressible fluids can also conserve mass, but their density must change to do so.