Home › Forum › SOFA › Programming with SOFA › Non linear implicit solver
Tagged: 64_bits, Linux_other, newton, nonlinear, SOFA_2006, Solver
 This topic has 7 replies, 2 voices, and was last updated 1 year, 2 months ago by jnbrunet.

AuthorPosts

7 November 2020 at 2 h 31 min #17605
I have looked at the EulerImplicitSolver and it only does a Newton step integration. I am trying to implement a version of it doing several Newton steps for integrating more complex energies and systems but I am stuck on the design of it.
Normally, when I do this kind of solvers outside of SOFA, I implement a solve very similar to what EulerImplicitSolver has, do a lineSearch using the direction provided by the solve and take the step when the error is smaller than the error I had.
Then I repeat this linearizing in the new state until the system error is smaller than a cerain value.
In SOFA is not possible to implement this directly as I need to update the state of the system, which is not done directly from the solver from what I have seen.
Could someone guide me on how to implement a non linear newton solver?
11 November 2020 at 14 h 09 min #17629Hey @jjcasmar,
Have a look at the
StaticSolver
which does a series of NewtonRaphson iterations. It should give you hints on how to achieve such state updates within the same time step.Let us know if you hit specific difficulties and don’t hesitate to show your work to the community, this is a component that could be very useful to others.
JeanNicolas
11 November 2020 at 15 h 10 min #17631Great, that is what I was looking for.
So if I understand correctly,
1) Compute forces and current error to bootstrap the solver
2) Compute a dx
3) Compute x as x0 + dx
4) Propagate this changes through the graph usingMechanicalPropagateOnlyPositionAndVelocityVisitor
5) Recompute forces, now in the new state
6) RepeaSomewhere between 4 and 5, the mappings have to be updated. Who is in charge of this?
MechanicalPropagateOnlyPositionAndVelocityVisitor
ormop.computeForces()
I will hack something fast later and push it to github so you can have a look.
11 November 2020 at 20 h 15 min #17636Hey @jjcasmar,
You understood correctly.
mop.computeForces()
is only responsible to propagate forces from the mapped mechanical object to their parent.MechanicalPropagateOnlyPositionAndVelocityVisitor
is responsible to propagate the position and velocity vectors from parent (top level, or system) vectors to mapped vectors.You can have a look at this StaticSolver from the Caribou plugin. It is a little bit more advanced compared to the SOFA one and its documentation is also richer.
11 November 2020 at 23 h 19 min #17641I think I still have one doubt.
From Caribou
// Updating the geometry x.peq(dx); // x := x + dx // Solving constraints // Calls "solveConstraint" method of every ConstraintSolver objects found in the current context tree. // todo(jnbrunet): Shouldn't this be done AFTER the position propagation of the mapped nodes? mop.solveConstraint(x, sofa::core::ConstraintParams::POS); // Propagate positions to mapped mechanical objects, for example, identity mappings, barycentric mappings, ... // This will call the methods apply and applyJ on every mechanical mappings. sofa::simulation::MechanicalPropagateOnlyPositionAndVelocityVisitor(&mparams).execute(this>getContext());
You set x as x+dx. Is that operation already triggering to store the new values in the mechanical objects or when is that done? I understand that the mappings are update using the visitor, but the DoF values should be set somewhere, no?
12 November 2020 at 8 h 57 min #17643Yeah this part isn’t very straightforward to understand.
Here,
x.peq(dx);
will store the new positions inside the mechanical object’s vector pointed by x. Now, which one is x? If you follow its definition, it is set to thesolve()
method’s parametersofa::core::MultiVecCoordId xResult
which is defaulted toVecCoordId::position()
, hence the MechanicalObject.position datafield. Now, thex.peq(dx);
operation starts a visitor that stops at every mapped mechanical objects, i.e. any mechanical object where a mechanical mapping is found.Next,
MechanicalPropagateOnlyPositionAndVelocityVisitor
will go through all mapped mechanical objects and compute their new position (MechanicalObject.position datafield) usingmapped_x = mapping.J*parent.x
.Let’s say you have this scene:
<Node name="root"> <StaticSolver /> <LinearSolver /> <Node name="object1"> <MechanicalObject name="mo" /> <Forcefield name="volume_ff" /> <Node name="surface"> <MechanicalObject name="mo" /> <Forcefield name="surface_ff" /> <Mapping name='mapping'/> </Node> </Node> <Node name="object2"> <MechanicalObject name="mo" /> <Forcefield name="volume_ff" /> </Node> </Node>
It goes like this:
1.
x.peq(dx);
root.object1.mo.position += root.object1.mo.dx
root.object2.mo.position += root.object2.mo.dx2.
MechanicalPropagateOnlyPositionAndVelocityVisitor
root.object1.surface.position = root.object1.surface.mapping.J * root.object1.mo.positionHope that cleared things a little bit.
JN12 November 2020 at 13 h 26 min #17646I think this solves all the issue.
I think you did a small typo in
root.object2.mo.position += root.object1.mo.dx
right?12 November 2020 at 13 h 40 min #17647Yes sorry, I’ve corrected it.

AuthorPosts
 You must be logged in to reply to this topic.