I've revised the code as shown in [1] to illustrate what I mean by assigning smin at convergence. It is interesting to see that the problem still converges when the elements start to fail!



[1] http://paste.pocoo.org/show/326090/

On Mon, Jan 24, 2011 at 9:08 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
On Mon, 24 Jan 2011, Andre Smit wrote:

Thanks

works like a charm! Two comments:

1) Since the model is 1/8 of actual the deformation rate is half of
actual test. So the strain rate in the Fc, Ft and Ec functions should be
doubled:
e *= 2 # With 1/8 model the applied deformation rate is half of actual
test

Yes, modify it as needed :)


2) Don't you think failed elements (smin) should only be assigned at the
end of an iteration run and not during sub-iterations.

I do not know - if you do not change the stiffness in sub-iterations, then nothing changes after the first one, and there is no need for them, right?

r.


On Mon, Jan 24, 2011 at 2:05 AM, Robert Cimrman <cimr...@ntc.zcu.cz>
wrote:
     ok, so the patch is in... the script should work
     out-of-the-box.

     r.


On Sat, 22 Jan 2011, Robert Cimrman wrote:

     I have overhauled the script, so that there are no
     globals, subiterations
     can be ended as needed etc. The traction() function is
     removed, the
     loading displacement is given simply by value. I hope
     it does what you
     intended.

     It needs a patch to TimeStepper, that I will be able to
     upload first on
     Monday, so either just look at it, or fix it
     yourself... :)

     r.

     ----- Reply message -----
     From: "Andre Smit" <freev...@gmail.com>
     To: <sfepy...@googlegroups.com>
     Subject: Conference
     Date: Fri, Jan 21, 2011 02:47


     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, Jan 19, 2011 at 10:22 AM, Robert Cimrman
     <cimr...@ntc.zcu.cz>
     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()
              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


     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
               zero?


          I'm not sure.


     That might be the problem causing the oscillations...

     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.




     --
     Andre

     --
     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.

     --
     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.



--
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.




--
Andre

--
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.



--
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.




--
Andre