I believe I've sorted out the oscillation problem, which was a bug in the code. I've pasted the latest version at [1] that now converges tightly after 3 or 4 iterations. I have a few comments.

As you will see:

a) I'm runnng the solver repeatedly with increasing displacements. This was mainly because I'm unable to change the displacement variable (disp) in traction() outside of the function. Globals and lists didn't work for me here.
b) To overcome this obstacle, the displacement for each solver/displacement cycle is saved to file and read from disk in traction().
c) I iterate the problem until the force converges and then increment the displacement.
d) The element stiffness matrix (Ele) is saved to file and loaded at the start of the subsequent displacement interval.

1) I considered adding a variable to ts, say ts.disp and adding the displacement to that, making it available in traction. I'm sure there is a better way though.  
2) I also cannot terminate the ts loop and am forced to run through all ts.n_steps. This is a problem since the forces converge at different number of steps at which point I would like the loop to end.



On Wed, Jan 19, 2011 at 10:22 AM, Robert Cimrman <> wrote:
On Wed, 19 Jan 2011, Andre Smit wrote:

Robert - I'm trying to change the displacement (disp) in the traction
function when the force stabilizes but for some reason the scope of disp
within traction() appears different to that within calc_force(), even
though disp is defined as global. Any ideas?

old_force = 0
def calc_force(pb, ts, state):
    global disp,old_force
    d = state()
    f = pb.evaluator.eval_residual(d)
    f.shape = (pb.domain.mesh.n_nod, 3)
    p = []
    for n in pb.domain.regions['Top'].vertices[0]:
    p = nm.array(p)
    if abs(old_force-p.sum()) < fvar:
        disp -= ddisp

def traction(ts, coors, bc=None):
    global disp
    #disp = -(ts.dt*(ts.step+1))
    print "traction disp: ", disp
    val = nm.empty_like(coors[:,0])
    return val

I would use brute force:

globals()['disp'] = ...

You cannot really change value of a global that is immutable (IMHO)...
A common idiom is to use list instead, like:

disp = [0]

def traction(ts, coors, bc=None):
   disp[0] = ...

     What is the meaning of 'smix' stiffness, for strain rate

I'm not sure.

That might be the problem causing the oscillations...


You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email to
For more options, visit this group at