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.

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!

You should have one lambda for each constraint, not one lambda for each particle. So, get rid of the lambdasA and lambdasB, and instead have a single lamdbas array of size cnstrsCount.

Then, fix up your lambda calculation, using equation 26 from the paper. So, something like:

float num = -C - alphaHat * lambdas[i] - gamma * (math.dot(gradA, velA) + math.dot(gradB, velB));
float den = (1 + gamma) * (math.dot(gradA * wA, gradA) + math.dot(gradB * wB, gradB)) + alphaHat;
float deltaLambda = num / den;
lambdas[i] += deltaLambda;

predictedPositions[cnstr.idA] +=  deltaLambda * gradA * wA;
predictedPositions[cnstr.idB] +=  deltaLambda * gradB * wB;

Hope that helps!

Hey Chris,

thank you very much for taking the time to write me back!

Of course, there should be a single lambda per-constraint, stupid me☹ I have been pulling my hair out on this for a while new. These dual lambdas was a leftover of some crazy brute-force techniques…

I have found another problem related with wrong graph-colouring of lambdas ids. It is fixed now and things kind of work now (attached). Below corrected code:

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
	float3 posA = positions[cnstr.idA];
	float3 posB = positions[cnstr.idB];

	//currently projected postions
	float3 predPosA = predictedPositions[cnstr.idA];
	float3 predPosB = predictedPositions[cnstr.idB];

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

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

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

	float C = len - cnstr.restLength;

	//normalize
	N = N / len;

	float3 gradA = -N;
	float3 gradB =  N;

	float wSum = dot(gradA * wA, gradA) + dot(gradB * wB, gradB) 
	
	float vA = dot(gradA, predPosA - posA);
	float vB = dot(gradB, predPosB - posB);

	float num = C + alphaHat * springLagrangeBuffer[id.x] + gamma * (vA + vB);
	float den = (1.0f + gamma) * wSum + alphaHat;

	float lambda = -num / den;

	springLagrangeBuffer[id.x] += lambda;

	predPosBuffer[idA] += float4(lambda * gradA * wA, 0);
	predPosBuffer[idB] += float4(lambda * gradB * wB, 0);
}

However, there is a weird stretch going on when I tune in the damping in constraints projection. I have even recorded a quick 3min YT video demonstrating the problem https://youtu.be/tvD7JNxFWMM if you fancy having a look.

However, my main worry is that, as long as such damping is clearly noticeable on simple chains, it seems that it does not have any noticeable impact on more complex objects such as cloths or soft-bodies.

I am really not sure whether I am still doing something wrong or it is just the way it is. If the latter, maybe you could please recommend some other damping method worth trying out to make cloths more realistic?

Cheers

Kay

Glad you got it working - sometimes implementing these papers can be challenging.

Adding position based damping will increase the stretching of the constraints, so what you are seeing is not unreasonable. You can somewhat compensate for the increased stretching by lowering compliance, or increasing iterations or substeps.

Some recent research suggests that increasing substeps is better than increasing iterations: http://mmacklin.com/smallsteps.pdf

In my own implementation, I also use a variant of hierarchical PBD which drastically reduces the number of iterations / substeps needed: https://matthias-research.github.io/pages/publications/hpbd.pdf​

However, this paper was written in the context of normal PBD, not XPBD. So, it does not consider compliance or damping values. I modified the algorithm by using the XPBD solver instead of PBD for each level. This modification does not converge to the exact correct result - it tends to be stiffer than expected, and stiffness is increased as more hierarchy levels are added. However, in my application I just manually tune compliance and damping to compensate for the extra stiffness.

You might also want to look more into your velocity damping - there might be an issue there. It should work for the cloth… Are you using the method from the original PBD paper?

Thanks Chris,

I will have a look at HPBD paper.

I am using velocity level damping taken from Matthis Mueller git (the author behind PBD papers). However, it requires additional iteration over all constraints compared to embedding it directly into positions projections.

https://github.com/matthias-research/pages/blob/master/challenges/pendulum.html

Below is my port to HLSL. I am using Gauss-Seidel appraoch on the GPU with graph-coloring. Hence, no atomics.

for (int i = 0; i < cnstrsCount; i++)
{
	DistanceConstraint cnstr = springsBuffer[id.x];

	//particles state after positions and velocities update
	float3 posA = posBuffer[cnstr.idA];
	float3 posB = posBuffer[cnstr.idB];

	float3 velA = velBuffer[cnstr.idA];
	float3 velB = velBuffer[cnstr.idB];

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

	float3 N = normalize(posB - posA);

	float v0 = dot(N, velA);
	float v1 = dot(N, velB);

	float dv0 = (v1 - v0) * min(0.5, springsVelocityDamping * dt * wA);
	float dv1 = (v0 - v1) * min(0.5, springsVelocityDamping * dt * wB);

	if (abs(dv0) > float.Epsilon &amp;&amp; abs(dv1) > float.Epsilon)
	{
		velBuffer[idA] += N * dv0;
		velBuffer[idB] += N * dv1;
	}
}

That looks reasonable for a distance constraint. Does your cloth use bending constraints? I think for cloth, dampening the bending motion is important for a realistic appearance.

At the end of the day, I suspect that the velocity based damping may actually perform better than position based - once you factor in the additional iterations / substeps required to make the position based damping stiff enough. You could profile it both ways to check.

Or if you go with hierarchical PBD, it should be stiff enough even with just a couple substeps, so you could use the position based damping there. HPBD is more complicated and managing the constraint hierarchy can be a little tricky, so you'll have to decide if its worth the effort or not for your application.