yuki-koyama / elasty

A research-oriented elastic body simulator
MIT License
363 stars 29 forks source link

Damping in XPBD #34

Open korzen opened 4 years ago

korzen commented 4 years ago

Hi Yuki,

great repo, thanks!

I was wondering whether you could add a proper damping support to Distance Constraints as described in the original XPBD paper (Eq. 26)

Thanks!

yuki-koyama commented 4 years ago

Thanks!

Yes, I recognize supporting the constraint damping in XPBD would be an important future TODO.

However, I'm not sure I can work on it shortly because many other higher-priority features are missing such as collisions.

korzen commented 4 years ago

I have come up with the following code but something is still not right. The motion is dampened also when cloth rotates, similarly to viscous dampening:

`DistanceConstraint spring = springsBuffer[id.x];

    int idA = spring.idA;
    int idB = spring.idB;

//pos at start of the sub-step float4 posA = posBuffer[idA]; float4 posB = posBuffer[idB];

    float4 particelA = predPosBuffer[idA];
    float4 particelB = predPosBuffer[idB];

//predicted position during constraints solve float3 predPosA = particelA.xyz; float3 predPosB = particelB.xyz;

//inv mass float wA = particelA.w; float wB = particelB.w;

    float restLength = spring.restLength;
    float stiff = spring.stiffness;

    float3 dir = predPosB - predPosA;

    float dist = length(dir);

    if (dist > EPS)
    {
        float wSum = wA + wB;
        if (wSum == 0)
            return;

        float C = dist - restLength;
        float3 N = dir / dist;

        float3 gradA = -N;
        float3 gradB =  N;

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

        float alphaHat = _compliance/ (dt * dt);
        float betaHat = _damping * (dt * dt);
        float gamma = (alphaHat * betaHat) / dt;

        float vA = dot(gradA, velA);
        float vB = dot(gradB, velB);

        float numA = -C - gamma * vA;
        float numB = -C - gamma * vB;

        float den = (1.0 + gamma) * wSum + alphaHat;

        float lambdaA = (numA / den);
        float lambdaB = (numB / den);

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

    }`

Have a try if you can.

Cheers