I'm trying to implement a frictional constraint using position based dynamics, but it is only half working. I am using the formulation in the paper "Unified Particle Physics for Real-Time Applications".

Here is my implementation:

	Particle *np = particle->nbs[j].neighbour; 

    vec3  r = particle->x - np->x;
    float r_length = glm::length(r);

    float distDiff = r_length - restDistance;

    if(distDiff >= 0 ) continue;

    //Frictional Constraint
    vec3 n   = r/r_length;
    vec3 xxi = particle->x - xi[particle->getIndex()];
    vec3 xxj = np->x - xi[np->getIndex()];

    vec3 tangentialDisplacement = (xxi - xxj) - glm::dot(xxi - xxj, n) * n;

    float td_length = glm::length(tangentialDisplacement);

    float genMass = ( imass / (imass + imass) );


    if(td_length < (staticFriciton * distDiff)){
        particle->x += genMass * tangentialDisplacement;

        np->x += -genMass * tangentialDisplacement; 
    }else{
        float upper = kineticFriction * distDiff;
        particle->x += genMass * tangentialDisplacement * std::min(upper/td_length, 1.f); 

        np->x += -genMass * tangentialDisplacement * std::min(upper/td_length, 1.f); 
    }

Half working how? What problems are you having, what is it doing that is different from what you think the code should be doing?

Here is a youtube link to the working penetration constraint:

And a link tot the failed friction constraint (Friction constraint applied after penetration constraint):

And I just realized that only the else statement gets executed in this part of the code:

  if(td_length < (staticFriciton * distDiff)){
        particle->x += genMass * tangentialDisplacement;

        np->x += -genMass * tangentialDisplacement; 
    }else{
        float upper = kineticFriction * distDiff;
        particle->x += genMass * tangentialDisplacement * std::min(upper/td_length, 1.f); 

        np->x += -genMass * tangentialDisplacement * std::min(upper/td_length, 1.f); 
    }

Can (upper/td_length) ever be negative?

distDiff will always be negative, otherwise there's no interpenetration. So the answer is yes. 

I'm thinking the tangentialDisplacement is being calculated wrongly because everything largely depends on it. It is perpendicular to the contact normal, a dot product gives zero.

The other thing I am not sure about is the which interpenetration distance to use, the one from the penetration constraint or the one from the friction. The paper does not make this clear.

Hoping someone has implemented one before, I can't find a reference anywhere.

If (upper/td_length) is always negative then those std::min calls are redundant because they would only take effect if the value is > 1.  Mainly what I was concerned about was:  Often when min or max calls are done, someone is thinking of magnitudes (i.e.  "make sure the magnitude is between 0 and 1") but accidentally misses the fact that the value is not a magnitude, but can be negative (in which case you would either want to get the magnitude using the absolute value or clamp between -1 and +1).

I'm not familiar with the system you're trying to implement; I'm just looking at "suspicious" code based on my own experience with bugs

Thanks for the help, I understand what you're saying. I will test the values thoroughly and see if they are producing negative values in the wrong places.

Finally got this to work.