Robert
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.
Comments:
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.
[1] http://paste.pocoo.org/show/324169/
Andre
On Wed, 19 Jan 2011, Andre Smit wrote:
I would use brute force: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()
pb.remove_bcs()
f = pb.evaluator.eval_residual(d)
pb.time_update()
f.shape = (pb.domain.mesh.n_nod, 3)
p = []
for n in pb.domain.regions['Top'].vertices[0]:
p.append(f[n][2])
p = nm.array(p)
if abs(old_force-p.sum()) < fvar:
my_output(disp,",",p.sum())
disp -= ddisp
old_force=p.sum()
def traction(ts, coors, bc=None):
global disp
#disp = -(ts.dt*(ts.step+1))
print "traction disp: ", disp
val = nm.empty_like(coors[:,0])
val.fill(disp)
return val
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]disp[0] = ...
def traction(ts, coors, bc=None):That might be the problem causing the oscillations...
What is the meaning of 'smix' stiffness, for strain rate
zero?
I'm not sure.
r.
--You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to sfepy...@googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.