Average Vector Field Integration for St. Venant-Kirchhoff Deformable Models
In this post I will explain the main idea behind the paper Average Vector Field Integration for St. Venant-Kirchhoff Deformable Models. I will try to give some context before explaining the main idea of the method in case you’re not too familiar with physics-based simulation. At the end of the post I will include a more technical summary of the method in case you’re already familiar with physics-based simulation and just want to implement the method without going through all the derivations presented in the paper.
Solving the equations of motion
The conventional story in physics is that differential equations describe the evolution of mechanical systems over time. However, these equations can only be solved analytically for very simple systems. In practice, solving the equations of motion for systems with complex potential energies and geometries requires numerical integration methods; in other words, algorithms to run the simulation rather than compact mathematical formulas that describe the entire motion. The fact that position-dependent forces act on a physical system during a given time interval is often expressed with an integral that looks like the following:
where f is a force function and x(t) is the time-dependent trajectory of the system. Intuitively, you can think of a spring connecting two particles. If the spring is deformed, the particles will be under the effect of restoring forces (f) trying to bring the spring back to its rest shape. This integral appears in the equations of motion of the system because a force acting on a particle implies changes in its velocity, which in turn implies changes in its position. To run a simulation, we need to replace this integral with actual numbers and there are different ways to do so. One option is to evaluate the force at the current frame, assume that it remains constant for the entire frame duration (h = t1 - t0) and use that to update the velocity of the system. Explicit methods such as forward Euler use this approach, but they are unstable in the sense that they can “explode” very easily unless we use very small time steps.
The backward Euler method approximates the force integral with an evaluation of the force at the next frame rather than the current frame.
At first this might look like it doesn’t make any sense, since the state of the system at the next frame (x1) is exactly the thing we’re trying to compute, but if you do some algebra to manipulate the original equations of motion, replacing the force integral with this rule, you end up with a different system of equations. Solving this is more involved than just applying explicit operations like in forward Euler (this is why backward Euler is considered an implicit method), especially when the forces are non-linear functions, but it can be done. The system of non-linear equations is often reformulated as a numerical optimization problem which can be solved with gradient-based algorithms.
Backward Euler is very stable in practice even with very large time steps. It doesn’t explode, but it introduces artificial damping as shown in the following animation. There is an initial potential energy due to elastic deformation in the ears of the bunny, but after running backward Euler, the total energy dissipates and the motion stops. Sometimes this is a desired effect, but it’s important to note that this is an artifact of the integration method over which we have little control.
Simulating the bunny with other implicit methods such as implicit midpoint and Newmark-beta avoids energy dissipation, but suffers from a different problem, the total energy can unboundedly increase, which manifests visually as “explosions”. This problem is even worse in explicit methods such as forward Euler, which requires very small time steps for stability.
Main idea of the method
As we have seen, the choice of different quadrature rules to evaluate the force integral can affect the simulation results significantly. The Average Vector Field (AVF) method can guarantee energy conservation (the paper has a proof, adapted from previous work, in the appendix) at the cost of evaluating exactly a particular approximation of the force integral, which corresponds to replacing the unknown trajectory with a straight line:
At this point we have just replaced one integral with another. We still need to actually compute it somehow. In the paper we explain in more detail how we can introduce a quadrature rule and transform the resulting non-linear equations into a minimization problem that can be solved with gradient-based optimization, which allows us to exactly conserve the total energy of the bunny and other simulation scenarios, avoiding both artificial damping and explosions.
Although the method is very general in principle and the paper provides some additional notes to make it work for arbitrary polynomial energies, we focus on the St. Venant-Kirchhoff (StVK) material, a quartic energy model that can produce interesting animations using finite element discretizations. The bunny was simulated using a tetrahedral mesh with the StVK material. We also adapted this material model to mass-spring systems to simulate cloth. Damping models can also be added to our method. Since damping is not a side effect of the integrator to achieve stability and we don’t have to worry about adding artificial damping to avoid explosions, we have more freedom to experiment with damping parameters to achieve different effects.
Method
I will omit some derivations explained in the paper in order to provide a summary of the method containing only a few equations, relevant for implementation purposes.
At the beginning of each frame, we know the positions of the mesh nodes x0, as well as their velocities v0. We want to compute x1 and v1, which correspond to the updated positions and velocities of the mesh nodes after h units of time.
Computing x1 requires solving the following optimization problem. Finding any local minimum of g(x) is enough to guarantee energy conservation. If you’re familiar with the variational formulation of backward Euler, popular in computer animation, the following formulas will probably look familiar to you. The formulas are similar to the backward Euler formulas, except for some additional terms and different multiplicative constants.
M is the mass and E is the potential energy function (a polynomial of degree four or less), including potentials due to elasticity and gravity.
After computing x1, we can compute v1 using the following rule:
Citation
@article {
rojas2018-avf-stvk,
title = {Average Vector Field Integration for {St. Venant-Kirchhoff} Deformable Models},
author = {Junior Rojas and Tiantian Liu and Ladislav Kavan},
journal = {IEEE Transactions on Visualization and Computer Graphics},
year = {2018}
}