[SOLVED] How to create bilateral constraint

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

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

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

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

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,

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

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,

Hi @mtournie

Max,could you give me a reply?