Viscous friction ?

[aleblanc30]

Hi SOFA community !

I'm looking to add a viscous friction force, i.e. a friction force proportionnal to the slipping velocity.

Is there a component that can help me do this ?

Otherwise, how can I get started to contribute one ?

[aleblanc30]

I read some of the code for FrictionContact. It seems to me that when a collision is detected, a lagrangian constraint is created that ensures that there is no interpenetration or slipping. Then the constraint is resolved in UnilateralConstraintResolutionWithFriction which represents a Gauss Seidel algorithm. The resolutionmethod is one iteration of the Gauss Seidel algorithm. The computed force is then stored to be used later ?

[aleblanc30]

Still in UnilateralConstraintResolutionWithFriction:resolution, the following code :

    if(normFt > fN)
    {
        const SReal factor = fN / normFt;
        force[line+1] *= factor;
        force[line+2] *= factor;
    }

returns a tangential force that either blocks tangential movement if normFt < fN, or is clamped to fN and represents the dry friction force that opposes tangential movement.

So I need to find a way to bring velocity information to this function and replace force[line+1] *= factor by something like force[line+1] *= factor*(1+alpha*slipping_velocity).

I guess this means I may have to write new classes similar to UnilateralLagrangianConstraint, UnilateralConstraintResolutionWithFriction, and FrictionContact.

[bakpaul]

Hello !

Yes to define new constraint behavior you need to write a class that inherits from ConstraintResolution. Those objects are the one defining how the constraint violation is linked to the constraint multiplier (force). Then you'll need to create an new UnilateralLagrangianConstraintWithViscousFriction that inherits from UnilateralLagrangianConstraint and just instead of using the UnilateralConstraintResolutionWithFriction uses your new constraint resolution.

The only difficulty I foresee is that delta is currently homogeneous to a position. You'll need to be able to access to the velocity from the input of the solvemethod. Dividing by dt might work but I would need more time to really understand what this means and how to relate it the the actual slipping velocity.

Maybe @ChristianDuriez might help you on this.

[aleblanc30]

Hi @bakpaul, thank you for your answer ! Can you be more specific about the solve method ? Which class does it belong to ? I'll keep working on this and report my progress here.

[hugtalbot]

I am not sure what you mean by the solve function @bakpaul do you mean to access the velocity data within the resolution function of the UnilateralConstraintResolutionWithViscousFriction ?

[aleblanc30]

I think this is what @bakpaul meant. I don't think dividing delta (by which I assume he means d, the third argument of the method) by dt would work. From what I understand dt would be the constraint violation, so here the algebraic distance that must remain positive ? Dividing it by a time would be homogeneous to a velocity but would definitely not be the slipping velocity. @hugtalbot @bakpaul any further input would be very welcome.

[aleblanc30]

I worked on figuring out where the constraints are instantitated. CollisionResponse implements ContactManager and thus determines the contact response, in the method createNewContacts.

In CollisionResponse.cpp, l. 168 :
auto contact = core::collision::Contact::Create(responseUsed, model1, model2, intersectionMethod,notMuted());

responseUsed is a string that is passed to a factory in Contact::Create. Contact.cpp, l.45:
return Factory::CreateObject(type,std::make_pair(std::make_pair(model1,model2),intersectionMethod));
where type is responseUsed.

This means you have to register the collision response model with the contact manager. This sends us back to the FrictionContact class, where the registration must happen in l.39 of FrictionContact.h :
SOFA_CLASS(SOFA_TEMPLATE3(FrictionContact, TCollisionModel1, TCollisionModel2, ResponseDataTypes), core::collision::Contact);

From this, a Contact object must receive all its information about the scene from two CollisionModels and an IntersectionMethod.

In order to implement the force I'm looking for directly with a contact, I must be able to get the velocity information from the collision models and the intersection method. If the data I need is accessible from those classes, I may not be able to implement viscous friction from the contact response. I may need to define a force field that takes into account collision information, instead of a collision response taking into account velocity.

edit : goodish news, TetrahedronModel has methods v1, v2, v3, v4 that returns velocities of the vertices, and TriangleModel has a method velocity that returns the average of the velocities of the vertices. So, at least in the particular case of triangular mesh collision models, which is the one I'm interested in, this should work.