Physics-Based Animation: Affine Body Dynamics with Penalty Spring Contact

In this assignment you will learn to implement affine body dynamics with penalty springs for simulating stiff objects in contact. Rendered assignment output played back at high-speed

WARNING: Do not create public repos or forks of this assignment or your solution. Do not post code to your answers online or in the class discussion board. Doing so will result in a 20% deduction from your final grade.

Checking out the code and setting up the python environment

These instructions use Conda for virtual environment. If you do not have it installed, follow the instructions at the preceeding link for your operating system

Checkout the code git clone git@github.com:dilevin/pba-assignment-abd.git {ROOT_DIR}, where {ROOT_DIR}* is a directory you specify for the source code.

Next create a virtual environment and install relevant python dependencies install.

cd {ROOT_DIR}
conda create -n csc417  python=3.12 -c conda-forge
conda activate csc417
pip install -e . 

Optionally, if you have an NVIDIA GPU you might need to install CUDA if you want to use the GPU settings

conda install cuda -c nvidia/label/cuda-12.1.0

Assignment code templates are stored in the {ROOT_DIR}/assginment directory.

WINDOWS NOTE: If you want to run the assignments using your GPU you may have to force install torch with CUDA support using

pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu121

Installation without conda

If you are having too many problems with Conda or prefer to use another package manager, we recommend using UV. If you do not have it installed, follow the instructions at the preceeding link for your operating system

Next, create a virtual environment and install relevant python dependencies:

cd {ROOT_DIR}
uv venv
uv pip install -e . 

If you opted to use UV, you can run the examples using:

uv run python main.py <arguments-for-tests>

Tools You Will Use

1. NVIDIA Warp -- python library for kernel programming
2. PyTorch -- python library for array management, deep learning etc ...
3. SymPy -- python symbolic math package

Running the Assignment Code

cd {ROOT_DIR}
python main.py --scene=tests/{SCENE_PYTHON_FILE}.py

By default the assignment code runs on the cpu, you can run it using your GPU via:

python main.py --scene=tests/{SCENE_PYTHON_FILE}.py --device=cuda

Finally, the code runs, headless and can write results to a USD file which can be viewed in Blender:

python main.py --scene=tests/{SCENE_PYTHON_FILE}.py --usd_output={FULL_PATH_AND_NAME}.usd --num_steps={Number of steps to run}

Occasionaly on windows the assignment will fail to run the first time, but subsequent attempts should work fine

Assignment Structure and Instructions

1. You are responsible for implementing all functions found in the assignments subdirectory.
2. The tests subdirectory contains the scenes, specified as python files, we will validate your code against.
3. This Google Drive link contains output from the solution code that you can use to validate your code. The output consists of USD (Universal Scene Description) files which contain simulated results. These can be played back in any USD viewer. I use Blender. You can output your own simulations as USD files, load both files in blender and examine the simulations side-by-side.

Background and Resources

Affine body dynamics was introduced to graphics in this paper by Lei Lan. The online physics-based animation text book by Li and colleagues as a good overview of the technique. Penalty springs for contact have a long history in simulation. In this assignment we are using the formulation described in in this paper

This assignment draws on previous lectures on Newton's Method, Affine Body Dynamics and Contact Handling.

Affine-Body Kinematics

Affine Bodies are a substitute for Rigid Body mechanics, useful when simulating stiff, effectively rigid objects. The key change is to replace a rigid kinematic map, represented by a rotation and a translation with an affine transformation, represented by a linear transformation and a translation:

Flattening transformation $Q(t)$ to a $12\times 1$vector$\mathbf{q}(t)$ allows us to rewrite the kinematic equation above as

where $J(\mathbf{X})$ is a $3\times 12$ matrix called the kinematic Jacobian. That matrix has the following form (which maintains equivalence with the matrix-valued map):

Kinetic Energy

The kinetic energy of an affine body can be defined using the generalized velocity $\dot{\mathbf{q}}$ as

The mass matrix $M$ consists of volume integrals of second-order monomials. This integral can be converted to a surface integral via the divergence theorem for evaluation using only an input triangle mesh.

The surface integral is evaluated as a sum over triangle integrals. Each triangle integral can be performed using (barycentric integration)[https://en.wikipedia.org/wiki/Barycentric_coordinate_system]

Potential Energy

Unlike a true rigid body representation, affine bodies can still deform without additional constraints. To maintain rigid deformations in an affine simulation, we introduce an isometric potential energy which penalizes non rigid deformations. The simplest form of such an energy is the following

,

where $F\in\mathcal{R}^{3\times 3}$ is the deformation gradient of the affine kinematic map, $\kappa$ is a scalar material stiffness (in this assignment set to $1e8$)

You should be able to convince yourself that the deformation gradient for an affine body is constant for the entire object and so is the potential energy leading to the simple integration formula

Contact Resolution using Penalty Springs

The second part of this assignment involves implementing contact resolution using penalty springs. Penalty springs are springs that push contacting objects apart if, and only if they are in contact. In this assignment the collision detection algorithm (the part of the simulation that checks if two objects are in contact) has been implemented for you. You are responsible for implementing the penalty spring energy, gradient and hessian for active contacts.

The penalty spring energy is "one-sided" in the sense the energy is $0$when objects are not interpenetrating and increases quadratically as intepenetration increases. In the assignment this energy is only processed for active contacts, triangle-vertex contacts at which intepenetration is already occuring. This means you can ignore the$max$ operation in your implementations.

Backward Euler for Implicit-Integration

You will implement Backward Euler time integration to time step the affine body system with contact. Backward Euler proceeds by minimizing the following energy using Newton's Method

Newton's method proceeds by iteratively evaluating the Hessian and Gradient of the above energy, computing the Newton search direction and then updating the current configuration variables for all objects in the simulation.

In order to choose an approriate step-size, $\alpha$ you will need to implement backtracking line search to satisfy the sufficient decrease condition. Due to the stiff (re quickly increasing) energies present in contacting simulation of this kind, your simulation is likely to fail without a proper line search. Practically it is important to limit the maximum numer of line search iterations (in this assignment we set the max to 5) in order to prevent prohibitvely long simulation times.

Adding fixed constraints to Newton's Method

Fixing degrees of freedom (DOFs) in Newton’s Method can be accomplished by projecting them out of the solve, i.e., by guaranteeing that the update vector $\mathbf{d}$ is always zero for the pinned or fixed entries in $\mathbf{q}$.

At first glance, it might seem reasonable to simply zero out the corresponding entries in $\mathbf{d}$ after computing it. However, this approach ignores the physical coupling between the fixed and free DOFs: the solver would not correctly account for how fixed constraints influence the rest of the system. Instead, we solve for a reduced Newton step as follows. Let $P$ be a projection (or selection) matrix such that

$$P\mathbf{d} = \tilde{\mathbf{d}},$$

where $\tilde{\mathbf{d}} \in \mathbb{R}^{n-f}$ contains only the free DOFs (with $n$ total and $f$ fixed). Because $P$ selects the non-fixed components, we can also write

$$\mathbf{d} = P^T \tilde{\mathbf{d}},$$

where $\mathbf{d}$ contains zeros in the fixed positions.

Substituting this form into the Newton update and re-deriving the reduced system gives

$$\tilde{H} = P H P^T, \quad \tilde{\mathbf{g}} = P \mathbf{g},$$

where $H$ and $\mathbf{g}$ are the Hessian and gradient of the unconstrained energy. Using these reduced quantities in the Newton iteration yields a solver that automatically respects the prescribed fixed constraints without modifying the physics of the underlying system.

Some Debugging Hints

1. Start by implementing the transform_mesh method. This will cause scenes to be rendered correctly at startup and allow meshes to move during simulation.
2. Next work on the cube_compress.py example which doesn't require any contact forces.
3. Next debug cube_floor.py which uses contact forces but only between too objects. Finally move on to more complicated examples.
4. The examples in this assignment will run more slowly then the previous assignment due to collision detection overhead. You will notice the simulation slow down once objects come into contact. This is normal.
5. If you are using Visual Studio Code or Cursor, use the interactive debugger and python debugging console.
6. There is a new "Step" button in the interface which advances the simulation by a single time step.

Admissable Code and Libraries

You are allowed to use SymPy for computing formulas for integrals, derivatives and gradients. You are allowed to use any functions in the warp and warp.sparse packages. You ARE NOT allowed to use code from other warp packages like warp.fem. You are not allowed to use any of warps specialized spatial data structures for storing meshes, volumes or doing spatial subdivision. You cannot use code from any other external simulation library.

Hand-In

We will collect and grade the assignment using MarkUs

Late Penalty

The late penalty is the same as for the course, specified on the main github page.