[SOLVED] Floor as a projective constraint

[jjcasmar]

Hello!

I would like to implement a floor as a projective constraint. If I understand correctly, a projective constraint would project the newton step of the newton-rapson method to the feasible space (this means, it would generate a velocity which is tangent to the constraint) and then updates the position, is it like this?

Could you please give some more info (mathematical background) about what the project methods should do?

virtual void projectResponse(const MechanicalParams* mparams, MultiVecDerivId dxId) = 0;
virtual void projectJacobianMatrix(const MechanicalParams* mparams, MultiMatrixDerivId cId) = 0;
virtual void projectVelocity(const MechanicalParams* mparams, MultiVecDerivId vId) = 0;
virtual void projectPosition(const MechanicalParams* mparams, MultiVecCoordId xId) = 0;
virtual void projectMatrix( sofa::defaulttype::BaseMatrix* /*M*/, unsigned /*offset*/ );

It would also be interesting if you give a clue of when is better to use a projective constraint and a LM constraint.

Thanks!!

[Matthieu Nesme]

Hi,

projective constraints can only be applied on independent dofs while Lagrange Multipliers based constraints can be applied to any dofs (even mapped ones).

Projective constraints are way cheaper (almost free) compared to LM that increases the system size. They also are generally numerically more stable but can only represent a limited range of constraints (representable by a projection matrix, constant for a time step). In my knowledge, they were first introduced in http://dl.acm.org/citation.cfm?id=280821 sec 5.3.

projectPosition and projectVelocity are called at the end of the time step to project the current state in the feasible space.
projectResponse implements P*dx (with P the projection matrix and dx is relative to a position variation = position difference, velocity, acceleration…) Used for iterative, non-assembled solvers.
projectJacobianMatrix implements P*J where J is a matrix. I guess it is used in the core LM constraint solver, but has to be checked!
projectMatrix is, for now, used as a way to get the projection matrix from several projective constraints (by doing the first call with the identity matrix). It is used in the Compliant plugin (another LM solver implementation).

Implementing a floor as a projective constraint is not a good idea as it would only be a post-treatment. Nothing would prevent a particle to penetrate the floor during the solve.
The idea would be to move the violated particles in projectPosition, and to adjust their velocity in projectVelocity (inelastic or elastic contacts).
The particles should stay free to move and the projection matrix should be the identity matrix (ie projectResponse does not do anything).

While clean unilateral constraints would prevent a particle to penetrate the floor during the solve but it would be way more costly.

A third way is to use penalities (see PlaneForceField). Note that penalities could be mixed with the projective constraint (used here as a kind of post-stabilization).

I hope this helps.

Matthieu

[jjcasmar]

Hi Matthieu,

thanks for your answer. I have been some time doing some other stuffs but now I want to come back to constraints.

Thanks for the BW98, it’s an interesting paper, I will have a closer look to projective constraints defined there.

Could you please explain how to use LM constraints in SOFA? I know that I have to add a LM Solver, but, as I’m coding my own ODESolver, I don’t know if I need to add something in the ODESolver code or something. How the pipe for LM constraints works in SOFA? If there is any doc about LM constraints in SOFA, please, link it, as I havent found it in the doc section :S

Thanks!

[Eulalie]

Hi,

To use LM constraints you need a ConstraintSolver in your scene. I’m used to the GenericConstraintSolver which is based on a Gauss Seidel algorithm. Also add a FreeMotionAnimationLoop and a ConstraintCorrection.

You can have a look at these examples:
examples/Components/constraint/SlidingConstraint.scn
examples/Components/constraint/BilateralInteractionConstraint.scn

Here is a link https://hal.archives-ouvertes.fr/tel-00785118/document to help you understand how the constraint problem is set and solved (see section 2.2.3 Constraint-based response).

If you want to implement your own constraint, you can have a look at the SlidingConstraint and BilateralInteractionConstraint implementation.
The three main functions are buildConstraintMatrix, getConstraintViolation and getConstraintResolution.
These functions are called by the GenericConstraintSolver to respectively, build the constraint matrix H (in the link I gave you), get the constraint violation delta_free which is computed from the free motion x_free and finaly you will need to provide a class for the resolution of your constraint. This is specific to the Gauss Seidel resolution which iterates and converges by solving each constraint while freezing the others.

No changes in your ODESolver is required. The correction is computed and applied by the ConstraintCorrection component.

Hope this helps.

Eulalie

[jjcasmar]

I’m a bit confused about your post Eulalie.

I have been playing around with the “examples/Components/constraint/DistanceLMConstraint.scn” and it doesn’t have the FreeMotionAnimationLoop nor the ConstraintCorrection components.

I have tried to simulate that scene using my Optimizator ODESolver (there is a post in this forum about it) and the constraints didn’t work.

I have also tried to add to that scene the FreeAnimationLoop and the ConstraintCorrection component but it’s still not working. The ConstraintCorrection components says that it can’t find a linear solver (althought a CG component is in the scene)

[Eulalie]

I see the confusion. I’m not used to the LMConstraintSolver component and LMConstraint class. As I said, I’m used to the GenericConstraintSolver which is different in the way of solving the Lagrange Multipliers (and thus in the way of building the scene).

I think that LMConstraintSolver actually increases the size of the system (as Matthieu mentioned). Unlike the GenericConstraintSolver which projects the mechanics into the constraint space and solves the reduced problem.

Maybe I should let someone else help you with LMConstraintSolver.
But if you want to try the approach described in https://hal.archives-ouvertes.fr/tel-00785118/document you should have a look at those examples:
examples/Components/constraint/SlidingConstraint.scn
examples/Components/constraint/BilateralInteractionConstraint.scn

Eulalie

[jjcasmar]

I have also tried the GenericConstraintSolver with the FreeAnimationLoop and the UncoupledConstraintCorrection for the scene I mentioned in my last post, but the constraint are not being solved. The console gives me the following output:

Bad size of constraintsResolutions in GenericConstraintProblem

The visual output is this

May the error be because the DistanceLMConstraint is not compatible with the GenericConstraintSolver?

[Eulalie]

Yes it is.

Only the constraints that implement buildConstraintMatrix(), getConstraintViolation() and getConstraintResolution() can be solved by the GenericConstraintSolver.

[jjcasmar]

Okay, I understand how the GenericConstraintSolver works, but I’m still confused about the ConstraintCorrection components.

Also, I would like to know how the LMConstraintSolver, LMConstraintDirectSolver and the LCPConstraintSolver work.

Is a bit strange that you have code several ways of solving constraints. How does each approach work and which are the benefits and drawbacks for each of them?

[Eulalie]

The ConstraintCorrection component computes and applies the correction:

dx_correction = A^{-1}*Ht*lambda
dx = dx_free + dx_correction

With lambda given by the Gauss Seidel algorithm implemented in the GenericConstraintSolver. Note that we do not compute the inverse of A. In practice we use a LDLt factorization of the matrix A.

I think that the main benefit of this approach is the computation time.
About the three other solvers, I don’t know. Maybe someone else can help you.

[Hugo]

Hi Juan Jo,

Does the explanations of Eulalie help you in solving your issue?
Cheers,

Hugo

[Hugo]

Hi @jjcasmar,

Did the previous answers help you in solving your issue ?
If yes, I would then close the topic.

Best wishes,

Hugo

[jjcasmar]

Sorry Hugo. Right now I’m doing some development in other techniques so I haven’t really try the solutions that appear in this topic yet.

[François Faure]

Hi,
if you need a bilateral constraint to maintain particles on a plane, try ProjectToPlaneConstraint. There is an example scn file.
best,

[Hugo]

Hi Juan Jo,

Did Francois’ hint help you?
Cheers,

Hugo

[Author]

Posts