7 November 2020 at 02:31 #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:09 #17629
Have a look at the
StaticSolverwhich does a series of Newton-Raphson 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.
Jean-Nicolas11 November 2020 at 15:10 #17631
Great, 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 using
5) Recompute forces, now in the new state
Somewhere between 4 and 5, the mappings have to be updated. Who is in charge of this?
I will hack something fast later and push it to github so you can have a look.11 November 2020 at 20:15 #17636
You understood correctly.
mop.computeForces()is only responsible to propagate forces from the mapped mechanical object to their parent.
MechanicalPropagateOnlyPositionAndVelocityVisitoris 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:19 #17641
I think I still have one doubt.
// 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 08:57 #17643
Yeah this part isn’t very straightforward to understand.
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 the
sofa::core::MultiVecCoordId xResultwhich is defaulted to
VecCoordId::position(), hence the MechanicalObject.position datafield. Now, the
x.peq(dx);operation starts a visitor that stops at every mapped mechanical objects, i.e. any mechanical object where a mechanical mapping is found.
MechanicalPropagateOnlyPositionAndVelocityVisitorwill go through all mapped mechanical objects and compute their new position (MechanicalObject.position datafield) using
mapped_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:
root.object1.mo.position += root.object1.mo.dx
root.object2.mo.position += root.object2.mo.dx
root.object1.surface.position = root.object1.surface.mapping.J * root.object1.mo.position
Hope that cleared things a little bit.
J-N12 November 2020 at 13:26 #17646
I 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:40 #17647
Yes sorry, I’ve corrected it.
- You must be logged in to reply to this topic.