[SOLVED] How to create bilateral constraint

[Nguyen]

Hi @cforest,

Could you please show me how to create bilateral constraint between the haptic tool and the object ? I have tried but it crashes.

Thanks,
Thien.

[Guillaume]

Hi Nguyen,

I just removed your reply from the topic grasping tool and created a topic for you.
Please create your own topic to ask a new question.

Regards,
Guillaume.

[Hugo]

Hi @tuanthienbk,

What you’re asking is pretty advanced and therefore complicated to support remotely. The best (but obvious) advice would be to get inspired from the original BilateralInteractionConstraint component. You can send us the error you have or give us your code to review, but as stressed before, this will be very complicated to support.

Best regards,

Hugo

[Nguyen]

Thank you all,

I was succeeded in creating a spring between the haptic tool and the collided point on the object. I used ContactMapper to create MechanicalState, then create a spring, then add the spring to the corresponding MechanicalState.

Now, instead of creating the spring, I create the BilateralInteractionConstraint, but it doesnt work.

[Hugo]

Dear Thien,

I am no expert in contraint resolution. Here is what I can explain. In a BilateralInteractionConstraint-like component, you need to implement the functions described below. Assuming your unconstrained linear system is A.x=b (where the solution is noted y) and your constraint problem is written trans(H).A.H λ = trans(H).y – d. Then:

• buildConstraint: this build the H matrix,
• getConstraintViolation: adds the -d to the right hand side,
• getConstraintResolution: implements a the diagonal resolution of the Gauss Seidel,

Hope this helps.

Hugo

[Nguyen]

Hi Hugo,

Thank you for your kind support.

Now, I can create the BilateralInteractionConstraint. I follow the code in the class StickContactConstraint. However, I sill have a problem that the haptic tool tip is not sticked to the object even though there is a interaction. What I dont understand is the function createContact in the class BilateralInteractionConstraint with many parameters and I dont understand the meaning of some parameters.

I am wondering if the BilateralInteractionConstraint only works with RigidMapping and IdentityMapping and does not work with BarycentricMapping ?

Any help is so much appreciated !

Thanks,
Thien.

[Hugo]

Hi Thien,

First, I can’t find the function createContact you are talking about in BilateralInteractionConstraint. So I cannot really help.

Second, the constraints are not related to the mappings. Therefore, I don’t see why this constraint would only work with some mappings and not others.

Hugo

[Nguyen]

I made a mistake. It’s addContact, not createContact.

Now, I solved the problem. The solution is to use GenericSolver (or similar name) instead of LCPSolver. The interesting thing is when I create a scene with BilateralInteractionConstraint statically, it works with LCPSolver but does not work dynamically i.e. added and removed by C++ code. I guess, LCPSolver needs to be refreshed whenever constraints are changed.

By the way, could you explain for me the differences between Generic and LCP ?

Thanks a lot !

[Hugo]

Dear ,

I am no expert in constraint solvers in SOFA. Nevertheless, I asked for some explanations, here are my understanding:

• The GenericConstraintSolver is a generic way of solving the constraints. It calls the constraint resolution and therefore allows to use any law (e.g. any contact law) that you want.
• LCPSolver is not really used anymore, but it corresponds to a pure Gauss Seidel resolution. This would be equivalent to a GenericConstraintSolver + BilateralConstraint

Hope this helps,
Cheers,

Hugo

[Wong]

Hello Hugo,

I want to know more about BilateralInteractionConstraint and I want to comprehend the concrete meaning of the equation trans(H).A.H λ = trans(H).y – d.

What are the concrete forms of H, A, λ and d?

What kind of reference can I find some content about constraint or constraint solver in?

Who is expert at constraint solvers?

Thanks,

Wong

[max]

Hi Wong,

If you want a nice intro to constrained dynamics, I suggest you start with David Baraff’s SIGGRAPH courses: https://www.cs.cmu.edu/~baraff/sigcourse/. A bit old but still relevant.

Generally, constrained dynamics involve solving a system of equations for constraint forces (usually noted lambda) so that the dynamics satisfy some constraints. In most cases, after proper time discretization the problem to solve is a quadratic program (QP), for which different (read: a whole lot of) solvers exist based on the properties of the QP.

In the case of Implicit Euler time stepping, the matrix A will be the inverse of the integration matrix A = inv(M – h^2 K) and matrix H will be the Jacobian matrix for the constraint mapping.

Because there are some many possible combination of (constrained) time-stepping schemes and constraint solvers, there are many ways of formulating constrained dynamics in SOFA. Unfortunately, I do not know what is the preferred/official way of dealing with this.

However, I can point you to the Compliant plugin that has a simple approach to it: the Compliant plugin defines Compliance components that act like the inverse of Stiffness components (i.e. ForceFields). In some sense, kinematic constraints are like infinitely stiff materials, which correspond to zero compliance.

For instance, you can consider a distance constraint between two points as an infinitely stiff spring between the points. So the Compliant plugin approach is this: you just put Compliance components at the end of mappings, the output of which should be constrained to zero, then using the adequate solvers you get constrained dynamics. Using non-zero compliance gives you soft constraints. For more information, you may check the examples.

Hope this helps,

[Hugo]

Dear Wong,

y is the solution (x_sol) of the linear system Ax=b, without constraints
H is the (holonomic) constraint matrix
λ are your Lagrange multipliers
and d is the constraint violation.
Correct me if I am wrong Max.

You can find information in the paper about SOFA, from page 21.

There is only few expert in constraint solvers, we do have some of them in Strasbourg.
The Consortium may provide you training if needed, if you become a Consortium member.

Best regards,

Hugo

[Wong]

Hello @hugo,

I have read the page 21 of the SOFA PAPER literally but it does not mention the violation and the equation trans(H).A.H λ = trans(H).y – d.

The bilateral constraint law is Ф(x1,x2)=0 but what is the concrete expression of Ф ?

Could you give me more information about this?

Thanks,
Wong

[Wong]

Hello max,

Thank you. The courses are fantastic and it is actually a nice introduce of dynamics simulation.

And I have also read about the theory of the compliant approach. It converts a linear system of equations into a quadratic program (QP) using KKT conditions. So what the solver does is to solve the QP not the linear system ?

And from the sofa documantion, when a collision occurs, the contact force must be computed. So how to extract the value of the contact force when the simulation is running ?

Wong

[max]

Hello Wong,

You are correct, when bilateral/unilateral constraints are defined on the system, the compliant plugin builds the KKT system corresponding to the QP to be solved.

This is the job of the OdeSolver/integrator component: it takes an integration scheme and the system state, and builds a QP for the numerical solver to solve. The standard one in Compliant is CompliantImplicitSolver, which is Implicit Euler time integration.

Now, there are various choices for solving the QP: you could work from the large indefinite KKT system directly, or you could compute the Schur complement then solve a smaller LCP instead. The standard for Compliant would be SequentialSolver, which is a Projected Gauss-Seidel-like LCP solver. It can also handle frictional constraints, but there may be anisotropy in the result so use it with caution.

If you only have bilateral constraints though, you might prefer iterative solvers like MinresSolver, or direct solvers like LDLTSolver depending on your requirements for speed/precision.

So the numerical solver solves the QP and computes Lagrange multipliers corresponding to contraint forces, then feed them back to the ode solver, which steps the system forward in time.

IIRC, most (all?) ode solvers in Compliant have a “propagate_lambda” flag that will cause the constraint forces to be written at the end of the time step in the “force” state vector of the degrees of freedom where the corresponding compliance component is defined, and these forces will be propagated all the way up to the independent dofs through the mapping Jacobian transposes.

This means that if you enable the flag on a scene with contacts, at the end of any time step the force vector for independent dofs will containt the **net** constraint forces (the sum of all the constraints acting on this object), which you can further process using e.g. a python script controller.

Let me know if you need help with this, I will try to add a simple example to demonstrate it.

Best,

[Wong]

Hello max,

After searching from the Sofa source code, a “propagate_lambda” flag is indeed in Compliant. However, I cannot find this flag in other components but I can find a “lambda” flag in some Correction and ForceField components.

As far as I know, a contact constraint is an unilateral constraint. Sofa should handle the collision by constructing such a constraint. And I do not know how and in which component Sofa constructing and solving this constraint.I have tried to read some relevant source code of Sofa, but the internal mechanism of sofa is complex so I cannot understand it.

Sofa can still handle the collision even if not using the Compliant. So could you provide me some documents and some simple examples about how Sofa constructing and solving the contact constraint(involving specific mesh vertexes and pertinent formulas)?

In the component LCPForcefeedback, the method computeForce can compute the contact force. Do you know how it computing the contact force and could you provide me a simple example to demonstrate it?

[Wong]

Hi @mtournie

Max,could you give me a reply?

[Author]

Posts