Skip to content

Rigid Body Simulation

Defining the Transformation

Reviewing cloth simulation and deformable objects, we simulate a deformation from reference XX to world space xx as

ΔxxXFΔX\Delta x\approx \underset{F}{\frac{\partial x}{\partial X}}\Delta X

so that we can have the strain

ΔxTΔxΔXTΔX=ΔXFTFΔXΔXTΔX=ΔXT(FTFI)ΔX\begin{align*} \Delta x^T\Delta x - \Delta X^T\Delta X &= \Delta X F^TF\Delta X - \Delta X^T\Delta X\\ &= \Delta X^T(F^TF-I)\Delta X \end{align*}

where FTFIF^TF-I is the Green Lagrange Strain.

For a rigid object, the object cannot be deformed so that FTFI=0FTF=IFO(3)F^TF-I = 0\Rightarrow F^TF=I\Rightarrow F\in O(3), i.e. FF is an orthogonal matrix.
For common 3D transformations, rotation matrices and reflection matrices are both orthogonal. However, reflection will have the object's topology passing through itself, which is not rigid.

Therefore, we limit FSO(3)F\in SO(3), i.e. rigid bodies can only rotate, and SO(3)SO(3) is called the special orthogonal Group,

RSO(3).RT=R1det(R)>0\forall R\in SO(3). R^T = R^{-1}\land \det(R) > 0

In addition, rigid bodies can translate. Therefore, we have the rigid transformation from reference to world space as

x(X,t)=R(t)X+p(t)\mathbf x(\mathbf X, t) = R(t)\mathbf X + \mathbf p(t)

x,XR3\mathbf x, \mathbf X \in \mathbb R^3, RSO(3)R\in SO(3) is the rotation and pR3\mathbf p\in \mathbb R^3 is the translation (often written in t\mathbf t but we don't want this to be confused with time tt).

Generalized Coordinates

Consider the mapping

x(X,t)=R(t)X+p(t)\mathbf x(\mathbf X, t) = R(t)\mathbf X + \mathbf p(t)

clearly, only R,pR, \mathbf p is about time, so that our generalized coordinates is

q={R,p}\mathbf q = \{R, \mathbf p\}

Note that q\mathbf q is a set instead of a unified stacked vector.

Generalized Velocity

Consider the time derivative of the transformation on a single point

v(X,t)=x˙(X,t)=ddtR(t)X+ddtp(t)=R(t)˙X+p(t)˙\begin{align*} \mathbf v(\mathbf X,t ) &= \dot{\mathbf x}(\mathbf X, t) \\ &= \frac{d}{dt}R(t)\mathbf X + \frac{d}{dt}\mathbf p(t)\\ &= \dot{R(t)}\mathbf X + \dot{\mathbf p(t)} \end{align*}

p(t)˙\dot{\mathbf p(t)} is called the linear velocity since it's just the velocity of the origin of the object moves,
R(t)˙\dot{R(t)} is the time derivative of rotation matrix, and it's more complex, we can represent it with

R(t)˙=R[X]×TRTω\dot{R(t)} = R{[X]_{\times}}^TR^T\omega

where []×[\cdot]_\times is the cross product matrix,

XR3.[X]×=[0XzXyXz0XzXyXx0],a×b=[a]×bX\in\mathbb R^3. [X]_{\times} = \begin{bmatrix}0&-X_z& X_y\\X_z&0&-X_z\\-X_y&X_x&0\end{bmatrix}, a\times b = [a]_\times b

ωR3\omega\in\mathbb R^3 is the angular velocity, ω\|\omega\| (magnitude) encodes the speed or rotation and angle encodes the axis of rotation.

v(X,t)=R(t)˙X+p(t)˙=R[X]×TRTω+p(t)˙=R[[X]×TI][RT00RT][ωp˙]=R[[X]×TI]Aq˙\begin{align*} \mathbf v(\mathbf X,t ) &= \dot{R(t)}\mathbf X + \dot{\mathbf p(t)}\\ &= R{[X]_{\times}}^TR^T\omega + \dot{\mathbf p(t)}\\ &= R \begin{bmatrix}{[X]_{\times}}^T&I\end{bmatrix}\begin{bmatrix}R^T&0\\0&R^T\end{bmatrix}\begin{bmatrix}\omega\\\dot{\mathbf p}\end{bmatrix}\\ &= R \begin{bmatrix}{[X]_{\times}}^T&I\end{bmatrix}A\dot{\mathbf q} \end{align*}

Since R,XR, \mathbf X is constant about tt, we can define our generalized velocity q˙R6=[ω,p˙]T\dot{\mathbf q}\in\mathbb R^6 = [\omega, \dot{\mathbf p}]^T. and we define A=[RT00RT]A = \begin{bmatrix}R^T&0\\0&R^T\end{bmatrix} to simplify the expression.

Kinetic Energy

As usual, we can find kinetic energy as

12Ωρv(X)22dΩ=12ΩρvTvdΩ=12q˙TAT(ρΩ[[X]×[X]×T[X]×[X]×I]dΩ)Aq˙=12q˙TMq˙\begin{align*} \frac12\int_\Omega\rho\|\mathbf v(\mathbf X )\|_2^2 d\Omega &= \frac12\int_\Omega \rho \mathbf v^T\mathbf v d\Omega\\ &= \frac12\dot{\mathbf q}^TA^T \big(\rho\int_\Omega \begin{bmatrix}[\mathbf X]_\times {[\mathbf X]_\times}^T &[\mathbf X]_\times \\ [\mathbf X]_\times&I\end{bmatrix} d\Omega\big)A\dot{\mathbf q}\\ &= \frac12\dot{\mathbf q}^TM\dot{\mathbf q} \end{align*}

where MM is the mass matrix, M=AT(M0)AM = A^T \big(M_0\big)A and M0R6×6M_0 \in \mathbb R^{6\times 6} is a constant matrix integrating over the whole volume.

Center of Mass

The center of mass is defined as the weighted average of points of the volume,

C=1mΩρXdΩ,m=ΩρdΩ\mathbf C = \frac1m \int_\Omega\rho\mathbf X d\Omega, m = \int_\Omega \rho d\Omega

With the center of Mass, we can transfer our reference space to the center-of-mass coordinates system

Xˉ=XC\bar{\mathbf X} = \mathbf X - \mathbf C

This is convenient for the computation of M0M_0,

M0=ρΩ[[X]×[X]×T[X]×[X]×I]dΩ=ρΩ[[Xˉ]×[Xˉ]×T[Xˉ]×[Xˉ]×I]dΩM_0 =\rho\int_\Omega \begin{bmatrix}[\mathbf X]_\times {[\mathbf X]_\times}^T &[\mathbf X]_\times \\ [\mathbf X]_\times&I\end{bmatrix} d\Omega = \rho\int_\Omega \begin{bmatrix}[\bar{\mathbf X}]_\times {[\bar{\mathbf X}]_\times}^T &[\bar{\mathbf X}]_\times \\ [\bar{\mathbf X}]_\times&I\end{bmatrix} d\Omega

Consider the off-diagonal entries,

Ω[Xˉ]×dΩ=Ω[XC]×dΩ=Ω[X]×dΩΩ[C]×dΩ=Ω[X]×dΩΩ1dΩ[C]×=0\begin{align*} \int_\Omega [\bar{\mathbf X}]_\times d\Omega &= \int_\Omega [\mathbf X -\mathbf C]_\times d\Omega\\ &= \int_\Omega [\mathbf X]_\times d\Omega - \int_\Omega [\mathbf C]_\times d\Omega\\ &= \int_\Omega [\mathbf X]_\times d\Omega - \int_\Omega 1d\Omega \cdot [\mathbf C]_\times\\ &= 0 \end{align*}

So that we only need to compute the diagonal

Surface Integral via Divergence Theorem

Note that M0M_0 is an integral over the volume, naturally, we can use tetrahedral mesh so that we can do this integral. However, a more efficient method will use only the surface mesh.

Divergence Theorem states that

Vf(X)dV=Af(X)ndA\int_V \nabla \mathbf f (\mathbf X) dV = \int_A \mathbf f(\mathbf X) \cdot \mathbf ndA

so that we can turn a volume integral to a surface integral Therefore, consider the entries of M0M_0

M0=ρΩ[[Xˉ]×[Xˉ]×T00I]dΩ[Xˉ]×[Xˉ]×T=[Xyˉ2+Xzˉ2XxˉXyˉXxˉXzˉXxˉXyˉXxˉ2+Xzˉ2XyˉXzˉXxˉXzˉXyˉXzˉXxˉ2+Xyˉ2]\begin{align*} M_0 &= \rho \int_\Omega \begin{bmatrix}[\bar{\mathbf X}]_\times {[\bar{\mathbf X}]_\times}^T &0\\0&I\end{bmatrix}d\Omega \\ [\bar{\mathbf X}]_\times {[\bar{\mathbf X}]_\times}^T &= \begin{bmatrix} \bar{X_y}^2 + \bar{X_z}^2 & -\bar{X_x}\bar{X_y}&-\bar{X_x}\bar{X_z} \\ -\bar{X_x}\bar{X_y} &\bar{X_x}^2 + \bar{X_z}^2&-\bar{X_y}\bar{X_z}\\ -\bar{X_x}\bar{X_z}&-\bar{X_y}\bar{X_z}&\bar{X_x}^2 + \bar{X_y}^2\end{bmatrix} \end{align*}

I\mathcal I is called the inertia matrix and we can integrate M0M_0 over each entry.

Since each entry is simple enough, we can pick some integral of it and then apply divergence theorem. For example

ΩXyˉ2+Xzˉ2dΩ=T13(0Xyˉ3Xzˉ3)NdTΩXyˉXzˉdΩ=T(XxˉXyˉXzˉ00)NdTΩ1dΩ=T13(XxˉXyˉXzˉ)NdT\begin{align*} \int_\Omega \bar{X_y}^2 + \bar{X_z}^2d\Omega &= \int_T\frac13\begin{pmatrix}0\\\bar{X_y}^3\\\bar{X_z}^3\end{pmatrix}\cdot \mathbf N dT \\ \int_\Omega \bar{X_y}\bar{X_z}d\Omega &= \int_T\begin{pmatrix}\bar{X_x}\bar{X_y}\bar{X_z}\\0\\0\end{pmatrix}\cdot \mathbf N dT \\ \int_\Omega 1d\Omega &= \int_T\frac13\begin{pmatrix}\bar{X_x}\\\bar{X_y}\\\bar{X_z}\end{pmatrix}\cdot \mathbf N dT \end{align*}

Then, because we are using triangle meshes, we can integrate over TT through the barycentric coordinates via

Xˉ=Xˉ0ϕ0+Xˉ1ϕ1+Xˉ2ϕ2\bar{\mathbf X} = \bar{\mathbf X}_0 \phi_0 + \bar{\mathbf X}_1 \phi_1+\bar{\mathbf X}_2 \phi_2

Then, define h:(ϕ0,ϕ1,ϕ2)Rh: (\phi_0, \phi_1, \phi_2)\rightarrow \mathbb R and integrate as

2Area0101ϕh(1ϕ1ϕ2,ϕ1,ϕ2)dϕ2dϕ12\text{Area}\int_0^1\int_0^{1-\phi}h(1-\phi_1-\phi_2, \phi_1, \phi_2)d\phi_2d\phi_1

Therefore, we can obtain

M0=[ρΩ[Xˉ]×[Xˉ]×TdΩ00mI]=[I00mI]M_0 = \begin{bmatrix}\rho \int_\Omega[\bar{\mathbf X}]_\times {[\bar{\mathbf X}]_\times}^Td\Omega &0\\0&mI\end{bmatrix} = \begin{bmatrix}\mathcal I&0\\0&mI\end{bmatrix}

where I\mathcal I is the inertia matrix

Final Form of Kinetic Energy

Put every thing together

T=12q˙TATM0Aq˙=12q˙T[RT00RT]T[I00mI][RT00RT]q˙=12q˙T[RIRT00mI]q˙\begin{align*} T &= \frac12 \dot{\mathbf q}^TA^T M_0 A\dot{\mathbf q}\\ &= \frac 12\dot{\mathbf q}^T \begin{bmatrix}R^T&0\\0&R^T\end{bmatrix}^T \begin{bmatrix}\mathcal I&0\\0&mI\end{bmatrix} \begin{bmatrix}R^T&0\\0&R^T\end{bmatrix}\dot{\mathbf q} \\ &= \frac 12\dot{\mathbf q}^T \begin{bmatrix}R\mathcal I R^T&0\\0&mI\end{bmatrix}\dot{\mathbf q} \end{align*}

Potential Energy

Since the object is rigid and cannot deform,

V=0V = 0

Newton-Euler Equations (The Lagrangian)

Since V=0V = 0, we have L=TL = T so that the Euler-Lagrange Equation gives Newton-Euler Equations

  1. Conservation of Angular Momentum

    (RIRT)ω˙=ω×((RIRT)ω)+τext(R\mathcal I R^T)\dot\omega = \omega \times ((R\mathcal IR^T)\omega) + \tau_{ext}

    Where ω˙\dot\omega is the angular acceleration, τext\tau_{ext} is the external torque and ω×((RIRT)ω)\omega \times ((R\mathcal IR^T)\omega) is the quadratic velocity vector 2. Conservation of Linear Momentum

    mIp¨=fextmI\ddot{\mathbf p} = \mathbf f_{ext}

    where p¨\ddot{\mathbf p} is the acceleration and fext\mathbf f_{ext} is the external force.

External Torques and Forces

For a force applied to some position x\mathbf x as f(x)\mathbf f(\mathbf x), this will result in a rotation as well as a translation. Note that

v(Xˉ,t)=R[[X]×TI][RT00RT][ωp˙]\mathbf v(\bar{\mathbf X}, t) = R \begin{bmatrix}{[X]_{\times}}^T&I\end{bmatrix} \begin{bmatrix}R^T&0\\0&R^T\end{bmatrix} \begin{bmatrix}\omega\\\dot{\mathbf p}\end{bmatrix}

Let JR3×6=R[[X]×TI][RT00RT]J \in\mathbb R^{3\times 6} = R \begin{bmatrix}{[X]_{\times}}^T&I\end{bmatrix}\begin{bmatrix}R^T&0\\0&R^T\end{bmatrix} be the rigid body Jacobian matrix, which converts the generalized velocity to the velocity on some point

v(Xˉ,t)=Jq˙\mathbf v(\bar{\mathbf X}, t) = J \dot{\mathbf q}

Therefore, given x\mathbf x in the world space, we first inverse it back to reference space w.r.t. the center-of-mass

Xˉ=RT(xp)C\bar{\mathbf X} = R^T(\mathbf x - p)-\mathbf C

and then we can compute f(Xˉ)f(\bar{\mathbf X}) and J(Xˉ)J(\bar{\mathbf X}) so that

(τextfext)=J(Xˉ)Tf(Xˉ)\begin{pmatrix}\tau_{ext}\\\mathbf f_{ext}\end{pmatrix} = J(\bar{\mathbf X})^T\mathbf f(\bar{\mathbf X})

Time Integration

The NE Equations are 2 second order ODEs, and they are not inter-related. Therefore, we can derive update equations independently.

Update on Translation

The update for linear momentum is simple. Note that this is just the translation of the center of mass, we can think this as the update of the point mass, the rule is simply

mp˙t+1=mp˙t+Δtfextpt+1=pt+Δtp˙t\begin{align*} m\dot{\mathbf p}^{t+1} &= m\dot{\mathbf p}^t + \Delta t \mathbf f_{ext}\\ \mathbf p^{t+1} &= \mathbf p^t + \Delta t\dot{\mathbf p}^t \end{align*}

Update on Rotation

For the angular velocity update, we can simply using the Forward Euler, i.e. approximate

ω˙=Δt1(ωt+1ωt)\dot\omega = \Delta t^{-1}(\omega^{t+1}-\omega^t)

so that the update is

(RIRT)ωt+1=(RIRT)ωt+Δtωt×((RIRT)ωt)+Δtτextt(R\mathcal I R^T)\omega^{t+1} = (R\mathcal I R^T)\omega^t + \Delta t \omega^t\times ((R\mathcal I R^T)\omega^t)+\Delta t\tau^t_{ext}

However, when we update the rotation, we cannot simply use the xt+1=xt+ΔtvTx^{t+1} = x^t + \Delta t v^T, since this will destroy the orthogonality.
Alternatively, we can view this as solving an initial value problem.

dtx=vt,x(t0)=xtd_t\mathbf x = \mathbf v^t, \mathbf x(t_0) = \mathbf x^t

integrating yields that

dtx=vtt0t1dx=t0t1vtdtx(t)=vtΔt+xt\begin{align*} d_t\mathbf x &= \mathbf v^t\\ \int_{t0}^{t1}d\mathbf x &= \int_{t0}^{t1}\mathbf v^t dt\\ \mathbf x(t) &= \mathbf v^t \Delta t + \mathbf x^t \end{align*}

With this, we can view vt\mathbf v^t as the velocity of rotation around p\mathbf p so that

vt=ωt×(RXˉpp)y:=RXˉvt=dty=ωt×yy(t0)=yty(t1)=exp([ω]×tΔt)y(t0)Rt+1=exp([ω]×tΔt)Rt\begin{align*} \mathbf v^t &= \omega^t \times (R\bar{\mathbf X} - \mathbf p -\mathbf p)\\ \mathbf y&:= R\bar{\mathbf X}\\ \mathbf v^t &= d_t\mathbf y = \omega^t\times \mathbf y\\ \mathbf y(t_0) &= \mathbf y^t \mathbf y(t_1) = \exp([\omega]_\times^t\Delta t)\mathbf y(t_0)\\ R^{t+1} &= \exp([\omega]_\times^t\Delta t)R^t \end{align*}