NOTE: This documentation is work in progress and may contain errors and unfinished sections.
This document describes the internal workings of the solver. We describe how the different parts work independent from platform specific features. Although we might occasionally comment on particular implementation details.
Todo: give this solver a name, so that we can differentiate when a new/different solver is added.
Overview of the general algorithm goes here.
Todo: check and rewrite:
// note on invMass (stored in current/previous positions.w):
// integrateParticles()
// - if(current.w == 0) current.w = previous.w
// constraintMotion()
// - if(constraint.radius <= 0) current.w = 0
// computeBounds()
// - if(current.w > 0) current.w = previous.w
// collideParticles()
// - if(collides) current.w *= 1/massScale
// after simulate()
// - previous.w: original invMass as set by user
// - current.w: zeroed by motion constraints and mass-scaled by collision
The position constraint for keeping distance between points a and b is usually written as:
error = |b-a| - restLength
offset = error * (b-a)/|b-a|
= (b-a) - restLength*(b-a)/|b-a|
The equation calculating slack pops up often in the solver code. This does exactly the same as the above:
slack = 1 - restLength / |b-a|
offset = (b-a) * slack
= (b-a) - restLength*(b-a)/|b-a|
Many constraints have a stiffness parameter that can be used to influence the rate at which the constraint recovers from errors. Stiffness can be applied, in the most basic forms, as follows: $$ p_1 = p_0 + \Delta p k $$ where \(p_0\) and \(p_1\) are the current and next particle positions, \(\Delta p\) is the constraint correction offset (as seen in the Slack section), and \(k\) is the stiffness constant.
The constraint error will be reduced by a factor of \(k\) per iteration, making it solver frequency dependent. The remaining constraint error is \(\Delta p(1-k)^n\) after \(n\) iterations.
We adjust the stiffness constant based on the delta time to get framerate independence. We want to pick \(k_{\Delta t} \) so that the error after one second at reference frequency \(f\), \(\Delta p(1-k)^f\), and the error after one second with time steps of \(\Delta t\), \(\Delta p(1-k_{\Delta t})^{\frac{1}{\Delta t}}\), are equal: \begin{align} (1-k)^{f} &= (1-k_{\Delta t})^{\frac{1}{\Delta t}}\\ (1-k)^{f\Delta t} &= 1-k_{\Delta t}\\ k_{\Delta t} &= 1 - (1-k)^{f\Delta t}\\ \end{align}
This solution is most likely based on page 8 of this paper.
The user specifies the stiffness \(k\) (using the constraint phase configs and functions like Cloth::setTetherConstraintStiffness()) for the reference frequency \(f\) (set using Cloth::setStiffnessFrequency()). Instead of storing \(k\) as is we store it in logarithmic space: $$ k_{\log} = \begin{cases} \log_2(1-k),& 1-k>0\\ \text{-FLT_MAX_EXP},&\text{otherwise} \end{cases} $$
This way \(k_{\Delta t}\) can be calculated without using the powf() function: $$ k_{\Delta t} = 1 - \mathrm{e}^{\left(f \Delta t \log(2) k_{\log}\right)} $$
Note: this still uses the expf() function. We also need the extra \(\log(2)\) constant as for some reason \(k_{\log}\) is in the base 2 logarithm. This is probably to make the condition work with FLT_MAX_EXP.
Also note that even though \(k_{\Delta t}\) should be timestep independent the time step still affects the stiffness of the simulation, both because of error propagation and because the integrator is not timestep independent.
The first step of the cloth simulation is the integration. Explicit Euler integration is used to calculate the new position of the particles. The velocity of the particles is not stored, instead we use the position from the previous frame to calculate the velocity. We need to compensate for differences in delta time. Large fluctuations can cause problems, so the delta time is damped.
The following pseudo code describes how this is implemented:
//calculate damping constants from user setting 'damping'
logDamping = log_2(1-damping) //from ClothImpl<T>::setDamping
dampExponent = stiffnessFrequency * dt1 //from IterationState::create
//calculate delta time ajustment
scale = e^{logDamping * dampExponent} * {dt1/dt0} //from IterationState::create
//Do the integration
delta = (particle_position1-particle_position0) * scale + acceleration
particle_position1 = particle_position1 + delta
The integration code also replaces the current inverse mass with the previous inverse mass if the current is zero.
Todo: rotating reference frame.
Todo: how does damping work?
Drag and lift simulation is done in the applyWind() function. We use the following equations to model air drag and air lift:
$$ F_d = \frac{1}{2} \rho v^2 c_{drag} A $$
$$ F_l = c_{lift}\frac{1}{2} \rho v^2 A $$
where \(F_d\) and \(F_l\) are the drag and lift forces, \(\rho\) is the fluid density, \(v\) is the flow speed, \(A\) is the surface area and \(c_{drag}\) and \(c_{lift}\) are the drag and lift coefficients. We use different symbols and notation in the explanation below to match closer to the source implementation. Note that the equations above work with velocity in \(\mathrm{m.s^{-1}}\), while we work mostly with the offset per frame in \(\mathrm{m}\) (meanning velocity multiplied by the iteration delta time).
The following algorithm is used:
for each triangle
\(p_j\), \(c_j\) and \(m^{-1}_j\) are the previous position, current position and is the inverse mass of the \(j\)th particle in the triangle.
//Calculate current and previous center of the triangle
\(c_t = \frac{1}{3} \cdot \left( c_0 + c_1 + c_2 \right)\)
\(p_t = \frac{1}{3} \cdot \left( p_0 + p_1 + p_2 \right)\)
//Position difference including wind (same as flow speed times dt; in \(\mathrm{m}\))
\(d = c_t - p_t + wind\)
if rotating reference frame
\(d = c_t + R\left(d-c_t\right)\) //where \(R\) is the inverse local space rotation for one solver iteration
//Todo check/explain this
//Unnormalized normal of the triangle
\(n = \left(c_2 - c_0\right) \times \left(c_1 - c_0\right) \)
//Double the area of the triangle in \(\mathrm{m^2}\)
\(a = \left|n\right| \)
//Normalized triangle normal
\(\hat{n} = \frac{n}{a}\)
//Calculate the cos and sin of the angle \(\theta\) between the \(d\) and \(n\)
\(\Lrg{ \cos\left(\theta\right) = \frac{\hat{n} \cdot d}{\left|d \right|} }\)
\(\Lrg{ \sin\left(\theta\right) = \sqrt{\max(0, 1 - \cos\left(\theta\right)^2)}}\)
//Lift direction is orthogonal to \(d\) and in the \(d\) \(n\) plane. Its length is \(\left|d\right|\)
\(\Lrg{ l_{dir} = \frac{\left( d \times \hat{n}\right) \times d}{\left|d \right|} }\)
//Calculate the lift and drag impulse
//The lift is at its maximum when \(d\) is at an \(45^\circ\) angle to the triangle
//We use \(\sin\left(\theta\right)\cdot\cos\left(\theta\right) = \frac{1}{2}\sin\left(2\cdot \theta\right)\) to calculate this
\(i_{lift} = c_{lift} \cdot \cos\left(\theta\right) \cdot \sin\left(\theta\right) \cdot l_{dir} \cdot \left|d\right| \cdot \Delta t^{-1}\)
\(i_{drag} = c_{drag} \cdot \left|\cos\left(\theta\right)\right| \cdot d \cdot \left|d\right| \cdot \Delta t^{-1} \)
//\(c_{lift} \) and \(c_{drag}\) are the lift and drag coefficients
//combined impulse
\(\Lrg{ i =
\begin{cases}
\epsilon < \left|d\right|^2,& 0\\
\text{else},& \left(i_{lift} + i_{drag}\right) \cdot \rho \cdot a
\end{cases}
}\)
//where \(\rho\) is the fluid density
//\(\rho\) should be set to compensate for the double area in \(a\)
//apply impulses to the particle positions
for each particle \(j = \left\{ 0, 1, 2 \right\} \)
\(c_j = c_j - i \cdot m^{-1}_j \)
Note that when we combine the impulses we check if \(\left|d\right|\) is too small as this causes instabilities due to float division inaccuracies. This can cause different drag and lift behavior depending on the time step size (or solver frequency). Smaller time steps result in smaller position deltas which reach the \(\epsilon\) threshold sooner. This results in less drag/lift at small time steps (high solver frequencies) for slow moving particles.
A distance constraint can limit the movement of a particle to the volume of a sphere. The constraints are specified with an array of 4 dimensional vectors (physx::PxVec4) where the x, y, z elements define the center and w the radius of the sphere.
The constraints are interpolated between the start and target spheres if both are given.
The constrainMotion() function in the solver is responsible for enforcing these constraints. The following pseudo code describes how this is implemented:
delta = sphere_center - particle_position
sqrLength = epsilon + |delta|^2
radius = max(0, sphere_radius * scale + bias)
slack = 1 - radius / sqrt{sqrLength}
if(slack>0)
{
if(radius <= 0)
particle_invMass.w = 0 //Set the inverse mass to 0 if we are constrained to a zero radius sphere
slack = slack * stiffness
particle.xyz += delta * slack
}
Tether constraints help with reducing the stretchiness of the cloth without increasing the solver iteration count. This is done by adding additional max distance constraints connecting simulated particles with their nearest fixed particle(s). The tether constraints are stored as an index & length pair. The array of constraints has a multiple of particle count elements. The constraints in the array are stored in order so that the first particle of the constraint is element%numParticles. The index stored in the constraint defines the other (anchor) particle of the constraint.
The constraint for one particle is solved as follows:
offset = 0
for each tether connecting pb
//pa is the anchor particle
delta = pa - pb
radius = tetherLength * scale
slack = 1 - radius / (|delta| + epsilon)
offset += delta * max(slack, 0)
pb += offset * stiffness
The stiffness is calculated as follows:
stiffness = tetherConstraintStiffness * numParticles / numTethers
Algorithm:
r = restlength
pi = particle i
piw = inv weight of particle i
h = pb-pa
e2 = epsilon + |h|^2
er = r>0 ? (1 - r / sqrt{e2}) : 0
if(useMultiplier)
{
//from PhaseConfig.h cloth::transform
multiplierC = log2(stiffnessMultiplier)
compressionLimitC = 1 - 1 / compressionLimit
strechLimitC = 1 - 1 / stretchLimit
er -= multiplierC * max(compressionLimitC, min(er, stretchLimitC))
}
stiffnessExponent = stiffnessFrequency * dt/iterations
stiffness = log2(1-stiffness) //check when this happens //from PhaseConfig.h cloth::transform
stiffnessC = 1-2^{stiffness * stiffnessExponent}
ex = er * stiffnessC / (pbw+paw)
pa = pa + h*ex * paw
pb = pb - h*ex * pbw
Separation constraints do exactly the opposite of distance constraints. These constraints ensure that the particles remain outside the defined spheres.
The constraints are interpolated between the start and target spheres if both are given.
The constrainSeparation() function in the solver is responsible for enforcing these constraints. Solving a single constraint is done as follows:
//p is the particle position
//c is the sphere center
//r is the sphere radius
d = c-p
l = |d|
slack = 1 - r/l
if(slack < 0)
p += d * slack
The fabric data structure contains the constraint information that can be reused by multiple cloth instances. The constraints are divided in different phases. Each phase usually contains a different type of constraints such as (horizontal/vertical) stretching, searing, bending constraints.
mPhases contains indices indicating which set from mSets belongs to a phase.
mSets contains a list of start (and end) indices for mRestvalues and mIndices (if multiplied by 2). The first rest value of set s is mRestvalues[mSets[s]] and the last is mRestvalues[mSets[s+1]-1].
mRestvalues contains the rest lengths / edge lengths of the constraints.
mIndices contains pairs of indices connected by a constraint. (mIndices.size()*2 == mRestvalues.size())