+ +

+

Internal solver function/algorithm documentation¶

+

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.

+

+ +

+test togglecontent

+

Overview¶

+

Overview of the general algorithm goes here.

+

Integrate particles (integrateParticles())
Simulate wind/drag/lift (applyWind())
Distance constraints (constrainMotion())
Tether constraints (constrainTether())
Solve edge constraints (solveFabric())
Solve separation constraints (constrainSeparation())
Collision (collideParticles())
Self collision (selfCollideParticles())
Update sleep state (updateSleepState())

+

Particle invMass w component¶

+

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
+

+

+

Slack¶

+

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|
+

+

+

Integration¶

+

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?

+

Wind simulation¶

+

Drag and lift simulation is done in the applyWind() function. +The following algorithm is used:

+

for each triangle
+       \(p_j\) is the previous and \(c_j\) is the current position and \(m^{-1}_j\) is the inverse mass of the \(j\)th particle in of 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)
+       \(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
+       \(a = \left|n\right| \)
+
+       \(\Lrg{ \cos\left(\theta\right) = \frac{n \cdot d}{\left|n\right| \cdot \left|d \right|} }\)
+       \(\Lrg{ \sin\left(\theta\right) =  \sqrt{\max(0, 1 - \cos\left(\theta\right)^2)}}\)
+
+       //Lift direction is orthogonal to delta, in the delta-normal plane. Its length is \(\left|n\right|\cdot\left|d\right|\)
+       \(\Lrg{ l_{dir} = \frac{\left( d \times 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} \)
+       \(i_{drag} = c_{drag} \cdot \left|\cos\left(\theta\right)\right| \cdot a \cdot d \)
+       //\(c_{lift} \) and \(c_{drag}\) are the lift and drag coefficients
+       //\(c_{drag}\) should be set to compensate for the double area in \(a\)
+
+       //combined impulse
+       \(\Lrg{ i =
+               \begin{cases}
+               \epsilon < \left|d\right|^2,& 0\\
+               \text{else},& i_{lift} + i_{drag}
+               \end{cases}
+       }\)
+
+       //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.

+

Distance constraints¶

+

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¶

+

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
+

+

+

Edge constraints¶

+

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¶

+

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
+

+

+

Fabric data structure¶

+

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())

+

+ + +

Navigation