About the convergence problems:
initialguess will only matter if you prescribe unknowns (it basically smears out the change in Dirichlet b.c.s over the entire domain). You claim that initialguess makes a difference when you change the force magnitude, and that I cannot explain, nor can I repeat the simulations.
Are you sure you are not mixing up the results from simulations where you make use of "boundarycondition 5" ? Here, initialguess makes a difference. I cannot get unloading to work at all (initialguess=0 and initialguess=1 does not matter).
Two possibilities:
1. Another bug somewhere (if there is, I would guess in the layered cross-section code).
2. It is simply a very nonlinear problem and the standard Newton-Raphson solver is no good (which is common).
I tried using a modified NR by forcing the elastic tangent, and using linesearch, but neither helped.
From the output
NRSolver: 98 D_u: *1.297e-11 D_w: 4.277e-02 R_v: 3.154e-02
NRSolver: 99 D_u: *5.302e-11 D_w: 6.547e-02 R_v: 3.080e-01
NRSolver: 100 D_u: *7.860e-12 D_w: 4.277e-02 R_v: 3.154e-02
it seems clear that the solution is flip-flop'ing around the equilibrium. Adding "constrainednrminiter 5" seems like a good choice.
With this, everything converge.
Note: Since the constrained parameter hasn't been used for a long time, there was a bug making it so that it was never used, which needs to be fixed first
Around line 180 in NRSolver, it should set the "constrainedNRFlag":
this->constrainedNRminiter = 0;
IR_GIVE_OPTIONAL_FIELD(ir, this->constrainedNRminiter, _IFT_NRSolver_constrainedNRminiter);
this->constrainedNRFlag = this->constrainedNRminiter != 0;
(I'm considering trying to modify NRSolver so that it detects problems like this by default, but this is a much bigger problem to solve, general nonlinear solvers are actually quite difficult)
--------------------------------------
Now, a lengthy reply (sorry) on that second question:
Now, if I understand the problem correctly:
For metal forming in general, there would be no way to solve this problem more efficiently than to try different the loads, and re-compute the problem (different load gives a completely different result). Since you have a simple load shape (a single nodal load) there is only a single variable (the load amplitude) that you care about, and the optimization problem becomes much simpler.
Unfortunately, OOFEM doesn't have any built in support for any type of inverse/optimization problem.
Spring-back is not necessarily a completely linear problem; and here you will branch out and create, e.g. 20, different (possibly nonlinear, plastic) solutions. All the solvers in OOFEM only store a single (current) internal state (stress, strain, plastic strain), and it can't be easily changed to allow this type of branching (or going back to try again).
There is really nothing built into OOFEM that will help you solve this type of problem, except for possibly the save & restore functionality and modify the analysis to perform the multiple different spring-back simulations. This is a bit messy and unsafe, since the save and & restore functionality wasn't intended for this purpose (a high risk for segmentation faults). A long and thorny road ahead.
The only benefit for all of this work would be *performance* (and you must determine if this mess is worth the effort).
Or, as you suggest, have an external script that solves this optimization problem (by whatever means it can). Doing 20 solutions with different load levels and check the closest, or perhaps do some response-surface-optimization or bisection-search (if you really only need to adjust the load amplitude).
The reason why we have the Python interface is that specialized, conditional problem setups like these are hard to support with a simple input file.
You can of course code your solver in C++ and add it as a special EngineeringModel as well (though, this is not trivial to set up, and requires quite detailed knowledge about OOFEM).
