Project for Dr. Chenfanfu Jiang's CIS 563 course: Making a PBD simulation from scratch
Switch branches/tags
Nothing to show
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Failed to load latest commit information.

Position Based Dynamics

View Demo Video HERE

What's Implemented:


  • bending constraint
  • stiffness [stretch and compress] constraint
  • collision constraints
  • implicit solver
  • kinematic collisions (ground)
  • kinematic collisions (moving sphere)
  • restricting positions

Additional features

  • compression vs stretching constraint
  • friction
  • velocity affecting of the mesh on collision with other objects

Project Structure:

I formulated the project to be simple and easily readable. With Eigen as the backbone for all data-structures and mathematics, I also used std::vectors for listings such as a collection of all constraints and when doing the original file i/o so that I would not need to dynamically change the size of my eigen matrices while filling in the initial matrices. BaseInfo was used as the backend for all main includes, creation of structs and data types for use in the simulation, and just main simulation info such as k-values, dimensions, and the number of simulation frames to run. The pbd file by itself is just the algorithm, allowing me to comment and uncomment specific lines such as friction and velocity damping and what file outputs I’ll be working with for that simulation run instead of working with commenting and uncommenting larger blocks of code.

I also made use of namespaces for repeatedly used calculations and all classes used for the simulation, which include TriangleMesh, Face, and Particles, as a way to more easily organize code sections. Additionally, I chose to store all vertex information, such as positions, velocities, weights, etc using Eigen Matrices in Particles removing the need to loop through each of the indices to update their velocities for force updates and other common math updates used in the algorithm. Instead, I could just multiply the matrices appropriately, making the code that much more readable. As it’s name states, TriangleMesh holds the mesh data, which includes the original face information (and its triangulations) for input and output of obj files while also holding a Particles object which is used for all data manipulation.

For the Constraints, I created a class hierarchy so that I could reference any lower constraint such as a FaceBendingConstraint or a StretchConstraint by just the name of their parent class ‘Constraint,’ allowing them all to be stored in one std::vector in the overarching Constraints class. Additionally, for each constraint, any reference to a specific particle is done so by index so that after any position/velocity/etc attribute gets updated, the particular constraint can still reference the newly updated value [also allows me to not need as many pointers and helps with memory management].


Current limitations of my design are that I have not implemented self-collision or a more robust kinematic collision with larger objects instead of just signed-distance functions. However, one of the main issues of the pbd version I implemented (the one from 2006) is that the collision constraints used for calculations depend on the mesh size and could allow for, depending on the velocity of the moving object, an object to stretch and “pop through” the mesh. This phenomenon occurs when an object with large velocity moves and pushes the mesh to the point where (if allowed and given the proper constraints) the mesh stretches s.t. the whole is as big as the sphere moving through it. Since the sphere is no longer in the way of any of the mesh’s vertices, no collisions are detected and the ball continues forward. One way to fix this is to actually implement breakage when an edge stretches past a particular threshold; another is to alter the collision detection so that it includes either the midpoints of the edges of the triangle or the closest point along the edges to the collision object (to do the closest points it would be the same as doing the midpoints except the delta p for each vertex would change to become a linear interpolation weight instead of just used as an average estimate).

Additionally, I considered this a limitation at first; however, ended up figuring out how to perform partially the same effect without tetrahedralization. I originally thought that the only way pbd could be extended fully to an obj mesh is that the mesh would need to be tetrahedralized and pbd performed as a “springed” system inside of it. This is largely the case for simulating stronger metallic objects that should remain entirely rigid; however, if the user is simulating more ballooning objects that can bend given the proper forces [as noted by the last clip on my vimeo video], the mesh can bounce and act as if it was tetrahedralized. As a note for the mesh in the video, its stiffness [stretch and compress] and bending constraints were both set to 1 and I had added an extremely large force on top of the existing gravitational force pushing it downwards. As the mesh falls and bounces appropriately against the ground, the force still pushes, bending one of the cow’s legs out of shape, which ultimately gets rapidly popped back into position, and the cow bounces to rest onto its side. Compared to the two previous clips, the two main reasons the cow tries to maintain its shape in this is that its stiffness and bending constraints are set to 1 and this clip includes a “velocity bounce” for when the simulated mesh hits a collision object (bouncing in the direction of the collision normal).

Lastly, I am not currently updating my constraints in a random order, which means I am always updating the stretch constraints all at the same time then the bending constraints etc. This can lead to a skewing of the mesh as it is solved. Simply selecting the constraints randomly when updating will even out this skewing appropriately.

Not really a limitation but more something I found was that the number of constraint update iterations needed in the simulation for just stretch and bend constraints could be between 1-3; however, when including collision constraints, 5+ were required to remove artifacts created during the collision updates. This could partly be because I’m not currently updating my constraints randomly; however, it’s mainly because that’s a feature of how the pbd algorithm is implemented and is an effect of the finickyness of the k-values and the constraints chosen.

Future Extensions:

I’d ultimately like to implement complex collision constraints (ie between two meshes) and to finish my implementation of the balloon constraint and tearing constraints.