During the article where we created a Hookean particle system, we used cumulative external and internal forces from which accelerations are calculated based on Newton’s second law to simulate physical effects such as stretching, shearing, and bending.
However, a quasi-new methodology exists that removes the layers of velocity and acceleration, which is called position-based dynamics.
This method will allow us to solve our system by working directly on the positions of the particles by satisfying constraint functions.
With this methodology, we can simulate fabric as well as fluids, deformations, etc!
To begin, I would like to highlight the advantages and disadvantages compared to simulations based on external and internal forces.
- Advantages :
- Quick, easy to set up and controllable.
- It avoids “overshooting” problems encountered when using explicit time integration.
- It is able to work with arbitrary constraints as long as the gradient of the constraint function can be determined.
- The stiffness of the model depends on the stiffness parameter as well as the magnitude of the time step and the number of iterations of the solver.
for all vertices i initialize xᵢ = xᵢ_₀, vᵢ = vᵢ_₀, wᵢ = 1/mᵢ while simulating for all particles 𝑖 (pre-solve) 𝐯ᵢ ← 𝐯ᵢ + ∆𝑡*wᵢfₑₓₜ(𝐱ᵢ) 𝐩ᵢ ← 𝐱ᵢ + ∆𝑡*𝐯ᵢ for solverIterationCount for all constraints 𝐶 (solve) solve(𝐶, ∆𝑡) (solve for 𝐱ᵢ) for all particles 𝑖 (post-solve) 𝐱ᵢ ← 𝐱ᵢ + ∆𝐱ᵢ 𝐯ᵢ ← (𝐱ᵢ − 𝐩ᵢ)/∆t
Let’s go through this pseudo-code together to better understand it. Suppose we are using a particle that is constrained to a wire, like what you might see on a ‘bead rollercoaster’.
To start the simulation, we need to initialize our attributes. Each particle has three attributes: its mass (m), its position (𝐱) and its speed (𝐯).
Let’s imagine from this scenario that there are two contributing velocities - the velocity of gravity and an initialized velocity in the direction of the wire - which gives a final velocity vector in the direction shown in the image.
After initializing them, we can start updating them in the following steps.
Before imposing the wire constraints, the particle moves in the direction of the velocity vector for a given time interval, Δt.
To do this, we take a Euler step. We have already learned about this numerical method months ago.
Brilliant 🎉! We managed to approximate the new position of the particle! However… it’s no longer attached to the wire! Let’s fix that.
Onward to the beauty of position-based dynamics. The iterative solver manipulates the position of the particle in order to satisfy the constraints.
There are two types of constraints:
- Equality (bilateral) - Equality constraints are constraints that must always be enforced. An example of this type of constraint is a distance constraint where, say, :
- Inequality (unilateral) - An inequality constraint uses the comparison operator or . There are therefore several solutions that can satisfy such a constraint. For example, here is this particle which is falling to the ground. It will continue to fall as long as :
In addition, an important point is that these constraints belong to a group of constraints considered nonlinear. Here is the difference between nonlinear and linear constraints:
Linear - Linear inequality constraints have the form . Linear equalities have the form .
Nonlinear - Nonlinear inequality constraints have the form , where . Similarly, nonlinear equality constraints also have the form .
The difference between linear and nonlinear systems is beyond the scope of this article, but we may rediscover them soon!
As I said before, the job of the solver is to move the particle in order to satisfy the constraints.
The problem that needs to be solved comprises of a set of M equations for the 3N unknown position components. However, the number of constraints is not necessarily equal to . In this case, the system is considered as asymmetrical.
Knowing these details, we can consider the system that we have to solve as an asymmetric and nonlinear system with equalities and inequalities. How scary! 👻
The process is iterative as we can see by seeing the loop using
solverIterationCount. Also, each constraint is solved individually using a nonlinear Gauss-Seidel solver.
We will see more in a following article. For now, my goal is to present all this knowledge in a high level way.
After satisfying the constraints and moving the particle back to the wire, we need to solve for the new velocity. Finally, the loop is repeated!
- Position based Dynamics - Matthias Müller, Bruno Heidelberger, Marcus Hennix, John Ratcliff
- A Survey on Position Based Dynamics, 2017 - Jan Bender, Matthias Müller and Miles Macklin
- What is “first order” and “second order” in time? - StackExchange
- OPF Equality and Inequality Constraints - PowerWorld
- Nonlinear Constraints - MathWorks
- Linear Constraints - MathWorks