Damping in XPBD

Please don't post Bullet support questions here, use the above forums instead.
Post Reply
korzen303
Posts: 29
Joined: Thu Dec 19, 2013 12:13 pm

Damping in XPBD

Post by korzen303 »

Hello, I am working on a surgical simulator using XPBD. I am trying to implement a better damping as currently with simple viscous damping the virtual anatomy is either unresponsive or oscillates too much. I have adapted the solveDistVel() method from this repo https://github.com/matthias-research/pa ... dulum.html which corrects the velocities after the position projections but it does not make any noticeable difference for volumetric meshes or cloths (it works fine however for things like pendulums or hanging chains ).

So I came up with the following implementation embedded directly in distance constraints projection (https://matthias-research.github.io/pag ... s/XPBD.pdf, equation 26) but there is something not quite right in its behavior. Would you be so kind to have a quick look. Maybe you would spot an error in code.

Code: Select all

float dtSq = dt * dt;
float alphaHat = compliance / dtSq;
float betaHat = damping * dtSq;
float gamma = (alphaHat * betaHat) / dt;

for (int i = 0; i < cnstrsCount; i++)
{
DistanceConstraint cnstr = cnstrs[i];

//positions at start of the substep
float4 posA = positions[cnstr.idA];
float4 posB = positions[cnstr.idB];

float4 predPosA = predictedPositions[cnstr.idA];
float4 predPosB = predictedPositions[cnstr.idB];

float wA = massesInv[cnstr.idA];
float wB = massesInv[cnstr.idB];

float wSum = wA + wB;

if (cnstr.restLength == 0 || wSum == 0)
    continue;

float4 N= predPosB - predPosA;
float len = math.length(N);
if (len <= float.Epsilon)
    continue;

float C = len - cnstr.restLength;

N = N/ len;

float4 gradA = -N;
float4 gradB = N;

float4 velA = predPosA - posA;
float4 velB = predPosB - posB;

float numA = C + (alphaHat * lambdasA[i]) + (gamma * (math.dot(gradA, velA)));
float denA = ((1.0f + gamma) * wSum) + alphaHat;
float lambdaA = -numA / denA;

float numB = C + (alphaHat * lambdasB[i]) + (gamma * (math.dot(gradB, velB)));
float denB = ((1.0f + gamma) * wSum) + alphaHat;
float lambdaB = -numB / denB;

//these lambda arrays are zeroed at start of each substep
lambdasA[i] += lambdaA;
lambdasB[i] += lambdaB;

predictedPositions[cnstr.idA] +=  lambdaA * gradA * wA;
predictedPositions[cnstr.idB] +=  lambdaB * gradB * wB;
}
Thanks for your help!
Post Reply