[SOLVED] Optimization integrator

[jjcasmar]

Hi!

I’m trying to develop an optimization integrator as defined in http://alexey.stomakhin.com/research/tvcg2015_integrator.pdf.

This integrator, instead of computing the velocity and then integrating it to compute the positions, tries to calculate a position state which minimizes the energy.

I have been playing with the integrators in SOFA but I’m finding some trouble with MultiVecCoord and MultiVecDeriv. I explain this problem:

First of all, the integrator calculates a position, and after that, I should compute the velocity that gives that position. So, given a newPos for the next frame (t+1), the velocity in the next frame (t+1) will be newVel = (newPos – pos)/dt where pos is the positions in t. The problem with this is that I can’t do a vector operation where a Coord type gives a Deriv type.

If I try to do

newVel.eq(newPos, pos, -1.0);
newVel.teq(1/dt);

it gives me a warning that the operation is not possible.

[WARNING] [MechanicalObject "mObject2"] Invalid vOp operation 3 (velocity(V_DERIV),position(V_COORD),position(V_COORD),-1)

I also found a similar problem when trying to multiply the M matrix with the (newPos – (pos+vel*dt)). As the second term should be a Coord type to be computed, newPos is also a Coord type. But I can only do mop.addMdx with Deriv types.

How could I avoid this?

[Hugo]

Hi Juan Jo,

As we discussed, in your linear system Ax=b, x corresponds to the Newton step Δx (the 4) in the section3.1 of your article). In SOFA, your unknown vector Δx is a VecDeriv. You can therefore compute:

newPos.eq(pos, DeltaX, alpha)
newVel.eq(alpha.DeltaX, 1/dt)

You can then propagate the constraints as follows:

mop.solveConstraint(newVel,core::ConstraintParams::VEL);
mop.solveConstraint(newPos,core::ConstraintParams::POS);

Moreover, you were wondering about how to make sure that the position pos at the end of the time step is well update, even through the mapping (to allow several iterations for a Newton-Raphson-like algorithm). Here is the solution: at the end of each iteration, you need to execute the UpdateMapping visitor as it is usually done at the end of one step.

sofa::simulation::UpdateMappingVisitor(&execparams).execute(this->getContext());

This way, you should be able to compute several iterations in the solve() function in your integration scheme and therefore compute your line search (Newton). No need to modify the AnimationLoop.

I never experienced it myself, let me know if this helps.
Cheers,

Hugo

[Matthieu Nesme]

Hi,

you should also have a look to the Compliant Plugin where I did implement a non-linear solver CompliantNLImplicitSolver.
Note that the implementation itself is crappy (playing with both flat vectors and graph-scattered multi-vec is not elegant and I did not clean anything).
But it can inspire you.

I would also be interested in such an integrator, please let me know if you do implement it!

Matthieu

[jjcasmar]

I have implement it but I still need to check that it works correctly in different scenarios.

I still want to check how it works with constraints and I also need to implement some forcefields which come from an energyfield. SOFA has a lot of forcefields which are great, but for this integrator I need the energyfield for the linesearch. The elastic forces defined in SOFA dont have the potential energy defined, so I’m trying to use the elastic models defined in VegaFEM for this purpose.

If you are interested in the component, write me an email and I can share it with you, but as said, a lot of forcefields are not supported and I havent check yet how to work with constraints.

[Hugo]

Hi Juan Jo,

Are you then implementing the potential energy for the elastic forces in SOFA (based on Vega)? That would be very nice to share!
I know this is already implemented in TetrahedronFEMForceField::getPotentialEnergy()

If you need help on this, I can put you in contact with researchers in the field.

Cheers,

Hugo

[jjcasmar]

No, I’m taking another approach. I dont want to touch the elastic forces in SOFA as I dont know if they are integrable.

Instead of that, I’m creating a component which wraps the elastic forces (in fact, only the force models that are compatible with isotropicHyperelasticityFEM) in Vega. Of course, I have to link against Vega to do this, so right now I’m doing this in my plugin.

[Hugo]

Alright,

If you change your mind do no hesitate to let us know so that we can help.
And, if you need any assistance in the linking process, keep us informed (in a new topic), we would be happy to help on this as well.

Cheers,

Hugo

[jjcasmar]

Hi,

long time has passed since the last message in this topic.

I’m going to reimplement the components talked in this topic (the integrator and the vega forcefields) but some guidelines would be helpful.

For the integrator, I’m planning to use NLOPT [1] for the optimization step. Should I add it in the extlibs folder?

Also, for the vega fem, its clear that I need VEGA [2]. Should I add it in the extlibs?

Should I implement this as a new plugin or as part of the core? I believe that as they are a simple integrator and some forcefields, there should be placed in the core.

I will still have to check how to enable constraints for this integrator. As it’s an optimization, constraints are well defined as part of the optimization and NLOPT will take care of them.

[Hugo]

Hi JuanJo,

Very good news! Nice to hear!
I think the easiest way to proceed would be to create your own plugin (in a separate repo github or gitlab) and start developing there your components.
You can find information about how to create a plugin. There, you will be able to add as an extlib of your plugin the different libraries.

Do you already know if you would like to release your code under opensource (e.g. LGPL) licence?

Let us know (creating a new topic) if you face any issue !
Cheers,

Hugo

[jjcasmar]

Okay, I will create a plugin under the 16.08 branch and when I have it ready I will do a merge request or something.

I feel very comfortable with the opensource manifest, so I will release my code under an opensource license (still want to check which one).

The only thing I asked is my name to appear in the code/plugin/component whatever so in the future I can say that I code that and be able to prove it 🙂

[Hugo]

hi JuanJo,

We are very happy to have active developers ans contributors like you in the community.
We will see how to integrate your work when you’ll release it. We will be very happy to promote it with you!

Cheers,

Hugo

[Author]

Posts