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

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

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.

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

 r.

 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<sfepy-devel%...@googlegroups.com>

 --
 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<sfepy-devel%...@googlegroups.com>

 --
 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<sfepy-devel%...@googlegroups.com>

-- Andre

Think I've found the error - we need to increase fvar to 100 say and increase the number of n_steps = 100.

I believe the problem is now solved. Next step is to compare the output to actual data and code the tensile failure checks.

a

On Mon, Jan 24, 2011 at 12:38 PM, Andre Smit freev...@gmail.com wrote: Robert - sorry for all these posts but I'm making corrections as I go along.

  Please disregard the previous code link, I've corrected as
  shown at [1]. The force-displacement curve is now non-linear
  as one would expect. Still bothering me is the change in
  moduli for the same rate. If the displacement interval (dt)
  is -0.1 and the rate is abs(dt), one should get the same
  strain rate if dt is -0.01, but it increases. Still trying to
  track down this error.

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

On Mon, Jan 24, 2011 at 11:26 AM, Andre Smit freev...@gmail.com wrote: 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.

For comparison I was trying to run simple.py on the attached. It fails with the following error:

[grassy@myhost CFRAC]$ ~/sfepy/simple.py one.py sfepy: left over: ['lam', 'Dijkl', '__builtins__', '_filename', '__file__', 'stiffness_tensor_lame', '__name__', 'verbose', 'nm', 'mesh_file', 'youngpoisson_to_lame', 'young', 'mu', '__package__', 'traction', 'stress_strain', 'output_dir', 'poisson', '__doc__'] sfepy: guessing abaqus sfepy: reading mesh (/home/grassy/Dropbox/sfepy/msh/cyl_02.inp)... sfepy: ...done in 0.01 s sfepy: setting up domain edges... sfepy: ...done in 0.01 s sfepy: setting up domain faces... sfepy: ...done in 0.01 s sfepy: creating regions... sfepy:���� Omega sfepy:���� YZ-Plane sfepy:���� Top sfepy:���� XZ-Plane sfepy:���� XY-Plane sfepy: ...done in 0.01 s sfepy: equation "equilibrium": sfepy: dw_lin_elastic_iso.2.Omega(Asphalt.lam, Asphalt.mu, v, u) = 0 sfepy: setting up dof connectivities... sfepy: ...done in 0.00 s sfepy: using solvers: ��������������� ts: ts �������������� nls: newton ��������������� ls: ls sfepy: ====== time 0.000000e+00 (step�� 1 of 101) ===== sfepy: updating variables... sfepy: ...done Traceback (most recent call last): � File "/home/grassy/sfepy/simple.py", line 113, in <module> ��� main() � File "/home/grassy/sfepy/simple.py", line 110, in main ��� app() � File "/home/grassy/sfepy/sfepy/applications/application.py", line 29, in call_basic ��� return self.call( *kwargs ) � File "/home/grassy/sfepy/sfepy/applications/simple_app.py", line 121, in call ��� pre_process_hook=self.pre_process_hook) � File "/home/grassy/sfepy/sfepy/solvers/generic.py", line 217, in solve_direct ��� nls_status=nls_status) � File "/home/grassy/sfepy/sfepy/solvers/generic.py", line 137, in solve_evolutionary_op ��� for ts, state in time_solver( state0 ): � File "/home/grassy/sfepy/sfepy/solvers/ts.py", line 128, in __call__ ��� state = step_fun( self.ts, state0, step_args ) � File "/home/grassy/sfepy/sfepy/solvers/generic.py", line 71, in time_step_function ��� problem.time_update( ts ) � File "/home/grassy/sfepy/sfepy/fem/problemDef.py", line 460, in time_update ��� functions, create_matrix) � File "/home/grassy/sfepy/sfepy/fem/problemDef.py", line 422, in update_equations ��� self.equations.time_update(self.ts, ebcs, epbcs, lcbcs, functions) � File "/home/grassy/sfepy/sfepy/fem/equations.py", line 215, in time_update ��� self.variables.equation_mapping(ebcs, epbcs, ts, functions) � File "/home/grassy/sfepy/sfepy/fem/variables.py", line 239, in equation_mapping ��� var.equation_mapping(bcs, var_di, ts, functions) � File "/home/grassy/sfepy/sfepy/fem/variables.py", line 1372, in equation_mapping ��� self.eq_map.map_equations(bcs, self.field, ts, functions, warn=warn) � File "/home/grassy/sfepy/sfepy/fem/dof_info.py", line 316, in map_equations ��� if vv is not None: val_ebc[eq] = vv TypeError: array cannot be safely cast to required type

I believe the input is correct. Any ideas?

a

On Tue, Jan 25, 2011 at 2:52 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: It works for me *with a tweak, see below) even for fvar = 0.1, the other settings left as they were, for both compression and extension.

  One thing I noted is, that you compare floating point numbers using the '=='
operator - this is not safe due to floating point
  approximation. Use

  is_min(stiffness0, smin)

  instead of

  stiffness0 == smin

  otherwise you may get wrong indices...

  So the current stiffness update logic is:
  1. mark elements where stress > strengh using

  stiffness = nm.where(abs(szz) > abs(fc), smin, smix)

  2. using

  stiffness = nm.where(is_min(stiffness0, smin), smin, stiffness)

  ensures, that already failed elements remain failed. The "already failed" are taken from
the previous load step, not the
  sub-iterations.

  I think that makes sense.

  So good luck with the actual data fitting :)

  r.

On Mon, 24 Jan 2011, Andre Smit wrote:

  Think I've found the error - we need to increase fvar to 100 say and
  increase the number of n_steps = 100.

  I believe the problem is now solved. Next step is to compare the output
  to actual data and code the tensile failure checks.


  a

  On Mon, Jan 24, 2011 at 12:38 PM, Andre Smit <freev...@gmail.com>
  wrote:
  � � �Robert - sorry for all these posts but I'm making corrections
  � � �as I go along.

  � � �Please disregard the previous code link, I've corrected as
  � � �shown at [1]. The force-displacement curve is now non-linear
  � � �as one would expect. Still bothering me is the change in
  � � �moduli for the same rate. If the displacement interval (dt)
  � � �is -0.1 and the rate is abs(dt), one should get the same
  � � �strain rate if dt is -0.01, but it increases. Still trying to
  � � �track down this error.

  � � �[1] http://paste.pocoo.org/show/326129/


  On Mon, Jan 24, 2011 at 11:26 AM, Andre Smit
  <freev...@gmail.com> wrote:
  � � �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.

Thanks Robert!

I've attached a modification that shows the strain/stress output for the elastic case of the cylinder under compression. The slope of the line is Young's modulus as you'd expect. As with our non-linear analyses, the forces calculated after the first time step appear to be too high - not sure what the reason for this is.�

On Wed, Jan 26, 2011 at 2:25 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: There was some old code in the stress_strain() function. The exception was caused by putting directly the traction function into the EBC definition, instead of its name. The attached file should work.

  r.

Sorry - forgot to attach the figure - your point is taken re the nodal residual. Forces are calculated and summed on each of the nodes in the Top region. I checked this a while back and plotted in Paraview - the forces are zero at the other nodes within the model but equal and opposite in direction to the Top at the corresponding Bottom nodes of the model (as you'd expect).

a

On Wed, Jan 26, 2011 at 9:26 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: Where can I see the line? maybe you forgot to attach a figure?

  As for the forces, it might be better to compute them using a
  (inner) surface integral of stress instead of the nodal
  rezidual. Btw. you look at the nodal force in the first node
  of the Top region only, right? What are the values in the
  other nodes?

  r.

  On Wed, 26 Jan 2011, Andre Smit wrote:

        Thanks Robert!

        I've attached a modification that shows the
        strain/stress output for the
        elastic case of the cylinder under compression.
        The slope of the line is
        Young's modulus as you'd expect. As with our
        non-linear analyses, the
        forces calculated after the first time step
        appear to be too high - not
        sure what the reason for this is.ￜ

        On Wed, Jan 26, 2011 at 2:25 AM, Robert Cimrman
        <cimr...@ntc.zcu.cz>
        wrote:
        ￜ ￜ ￜThere was some old code in the
        stress_strain() function. The
        ￜ ￜ ￜexception was caused by putting directly the
        traction
        ￜ ￜ ￜfunction into the EBC definition, instead of
        its name. The
        ￜ ￜ ￜattached file should work.

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

Alright, that fixes it! Now I understand the reason for quasistatic! In the non-linear analyses, the time stepper is set programatically - is is possible to force solving for the first time step here as well as in:

��� for ii, disp in ds: ������� #my_output(dformat % (ii + 1, ds.n_step, disp))

������� force0 = 0.0

������� pb.ebcs['Load'].dofs['u.2'] = disp

������� pb.ts.is_quasistatic = True� # <<<=========== Force quasistatic ������� pb.time_update(ts)

On Wed, Jan 26, 2011 at 9:57 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote:

  So the problem is the single dot outside the linear curve,
  right?

  What if you add

  � � � �'quasistatic' : True,

  to the time-stepper options? As it is now, the first time
  step is not solved but taken from the initial conditions
  (unspecified = zero) and boundary conditions (nonzero!), so
  it is not in equlibrium!

r.

On Wed, 26 Jan 2011, Andre Smit wrote:

Sorry - forgot to attach the figure - your point is taken re the nodal residual. Forces are calculated and summed on each of the nodes in the Top region. I checked this a while back and plotted in Paraview - the forces are zero at the other nodes within the model but equal and opposite in direction to the Top at the corresponding Bottom nodes of the model (as you'd expect).

a

On Wed, Jan 26, 2011 at 9:26 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: � � �Where can I see the line? maybe you forgot to attach a � � �figure?

� � �As for the forces, it might be better to compute them using a � � �(inner) surface integral of stress instead of the nodal � � �rezidual. Btw. you look at the nodal force in the first node � � �of the Top region only, right? What are the values in the � � �other nodes?

� � �r.

� � �On Wed, 26 Jan 2011, Andre Smit wrote:

� � � � � �Thanks Robert!

� � � � � �I've attached a modification that shows the � � � � � �strain/stress output for the � � � � � �elastic case of the cylinder under compression. � � � � � �The slope of the line is � � � � � �Young's modulus as you'd expect. As with our � � � � � �non-linear analyses, the � � � � � �forces calculated after the first time step � � � � � �appear to be too high - not � � � � � �sure what the reason for this is.�

� � � � � �On Wed, Jan 26, 2011 at 2:25 AM, Robert Cimrman � � � � � �cimr...@ntc.zcu.cz � � � � � �wrote: � � � � � �� � �There was some old code in the � � � � � �stress_strain() function. The � � � � � �� � �exception was caused by putting directly the � � � � � �traction � � � � � �� � �function into the EBC definition, instead of � � � � � �its name. The � � � � � �� � �attached file should work.

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

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

I think the iterations are indirectly sorting this out. In fact I found that the initial stiffness(0) doesn't influence the results at all. I've attached results of compression (comp) and tension (tens)� tests at rate of 1mm/sec as well as the script in its current state. Based on lab tests, comparative compression failure at this rate occurs at 6.6 N/mm2 whereas tension failure occurs at 3.0 N/mm2. So the script is accurately predicting compression failure but under-estimating tension failure. Still exploring ...

a

On Wed, Jan 26, 2011 at 10:19 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: Well, this flag is used essentially only in the generic solver used by simple.py. In the strain rate solver script we solve every time step regardless this flag.

  The first time step is different, however, as strain equals
  strain0 (no previous strain defined), and so initial dstrain
  is zero. Maybe this causes some problems?

r.

On Wed, 26 Jan 2011, Andre Smit wrote:

  Alright, that fixes it! Now I understand the reason for
  quasistatic! In
  the non-linear analyses, the time stepper is set
  programatically - is is
  possible to force solving for the first time step here
  as well as in:


  ��� for ii, disp in ds:
  ������� #my_output(dformat % (ii + 1, ds.n_step,
disp))

  ������� force0 = 0.0

  ������� pb.ebcs['Load'].dofs['u.2'] = disp

  ������� pb.ts.is_quasistatic = True� #
<<<===========
  Force quasistatic
  ������� pb.time_update(ts)




  On Wed, Jan 26, 2011 at 9:57 AM, Robert Cimrman
  <cimr...@ntc.zcu.cz>
  wrote:

  � � �So the problem is the single dot outside the
  linear curve,
  � � �right?

  � � �What if you add

  � � �� � � �'quasistatic' : True,

  � � �to the time-stepper options? As it is now, the
  first time
  � � �step is not solved but taken from the initial
  conditions
  � � �(unspecified = zero) and boundary conditions
  (nonzero!), so
  � � �it is not in equlibrium!


  r.

  On Wed, 26 Jan 2011, Andre Smit wrote:

  Sorry - forgot to attach the figure - your point is
  taken re
  the nodal
  residual. Forces are calculated and summed on each of
  the
  nodes in the
  Top region. I checked this a while back and plotted in
  Paraview - the
  forces are zero at the other nodes within the model but
  equal
  and
  opposite in direction to the Top at the corresponding
  Bottom
  nodes of the
  model (as you'd expect).

  a

  On Wed, Jan 26, 2011 at 9:26 AM, Robert Cimrman
  <cimr...@ntc.zcu.cz>
  wrote:
  � � �Where can I see the line? maybe you forgot to
  attach a
  � � �figure?

  � � �As for the forces, it might be better to compute
  them
  using a
  � � �(inner) surface integral of stress instead of the
  nodal
  � � �rezidual. Btw. you look at the nodal force in the
  first
  node
  � � �of the Top region only, right? What are the values
  in
  the
  � � �other nodes?

  � � �r.

  � � �On Wed, 26 Jan 2011, Andre Smit wrote:

  � � � � � �Thanks Robert!

  � � � � � �I've attached a modification that shows the
  � � � � � �strain/stress output for the
  � � � � � �elastic case of the cylinder under
  compression.
  � � � � � �The slope of the line is
  � � � � � �Young's modulus as you'd expect. As with our
  � � � � � �non-linear analyses, the
  � � � � � �forces calculated after the first time step
  � � � � � �appear to be too high - not
  � � � � � �sure what the reason for this is.�

  � � � � � �On Wed, Jan 26, 2011 at 2:25 AM, Robert
  Cimrman
  � � � � � �<cimr...@ntc.zcu.cz>
  � � � � � �wrote:
  � � � � � �� � �There was some old code in
the
  � � � � � �stress_strain() function. The
  � � � � � �� � �exception was caused by
putting
  directly the
  � � � � � �traction
  � � � � � �� � �function into the EBC
definition,
  instead of
  � � � � � �its name. The
  � � � � � �� � �attached file should work.

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




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

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

Robert

I'm having problems relating the model's output to the actual test output. One issue that is confusing is that if I reduce the displacement interval from 0.001 to 0.0001, say and the time during this interval to ensure a consistent loading rate, then the results change significantly, both the maximum load at failure and the displacement at maximum load. This is unexpected. I can't think what the problem can be - any ideas on your side?

On Wed, Jan 26, 2011 at 10:50 AM, Robert Cimrman cimr...@ntc.zcu.czwrote:

>

 The first time step is different, however, as strain equals
 strain0 (no previous strain defined), and so initial dstrain
 is zero. Maybe this causes some problems?

 Alright, that fixes it! Now I understand the reason for
 quasistatic! In
 the non-linear analyses, the time stepper is set
 programatically - is is
 possible to force solving for the first time step here
 as well as in:


     for ii, disp in ds:
         #my_output(dformat % (ii + 1, ds.n_step, disp))

         force0 = 0.0

         pb.ebcs['Load'].dofs['u.2'] = disp

         pb.ts.is_quasistatic = True  # <<<===========
 Force quasistatic
         pb.time_update(ts)




 On Wed, Jan 26, 2011 at 9:57 AM, Robert Cimrman
 <cimr...@ntc.zcu.cz>
 wrote:

      So the problem is the single dot outside the
 linear curve,
      right?

      What if you add

             'quasistatic' : True,

      to the time-stepper options? As it is now, the
 first time
      step is not solved but taken from the initial
 conditions
      (unspecified = zero) and boundary conditions
 (nonzero!), so
      it is not in equlibrium!


 r.

 On Wed, 26 Jan 2011, Andre Smit wrote:

 Sorry - forgot to attach the figure - your point is
 taken re
 the nodal
 residual. Forces are calculated and summed on each of
 the
 nodes in the
 Top region. I checked this a while back and plotted in
 Paraview - the
 forces are zero at the other nodes within the model but
 equal
 and
 opposite in direction to the Top at the corresponding
 Bottom
 nodes of the
 model (as you'd expect).

 a

 On Wed, Jan 26, 2011 at 9:26 AM, Robert Cimrman
 <cimr...@ntc.zcu.cz>
 wrote:
      Where can I see the line? maybe you forgot to
 attach a
      figure?

      As for the forces, it might be better to compute
 them
 using a
      (inner) surface integral of stress instead of the
 nodal
      rezidual. Btw. you look at the nodal force in the
 first
 node
      of the Top region only, right? What are the values
 in
 the
      other nodes?

      r.

      On Wed, 26 Jan 2011, Andre Smit wrote:

            Thanks Robert!

            I've attached a modification that shows the
            strain/stress output for the
            elastic case of the cylinder under
 compression.
            The slope of the line is
            Young's modulus as you'd expect. As with our
            non-linear analyses, the
            forces calculated after the first time step
            appear to be too high - not
            sure what the reason for this is.

            On Wed, Jan 26, 2011 at 2:25 AM, Robert
 Cimrman
            <cimr...@ntc.zcu.cz>
            wrote:
                 There was some old code in the
            stress_strain() function. The
                 exception was caused by putting
 directly the
            traction
                 function into the EBC definition,
 instead of
            its name. The
                 attached file should work.

                 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.




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

-- Andre

I do not have a definite answer, only some remarks that might be related to the problem:

When I run the input srs_03.py, I get:

Loading rate (mm/s): -0.02117 warning: (almost) singular matrix! (estimated cond. number: 4.38e+14) Initing ... 6.66666666667 0.000905933418981

• the singular matrix means (provided the Lame's parameters are ok), that the EBCs are not sufficient to prevent rigid body motions. The current EBCs just fix the 'z' axis motion, but you need to fix at least one other node completely and another node to prevent rotations around the 'z' axis.

This alone might cause the problems, as who knows what the linear solver does when it solves with a singular matrix.

Otherwise to your problem - you keep the loading rate the same, but change the displacement increment size in a single time step, right?

There is always a trade-off when approximating (time) derivatives with differences - if the step is too large, the derivative is not accurate, while if the step is too small, one gets big rounding-off error due to floating point arithmetics. So there is an optimal step, but I cannot tell you what it is for your problem...

The scheme we have is sort of ad-hoc devised, maybe a mathematician should look at the equations in future...

As for the stress computation, how did you change the way of its calculation?

So you see, I have more questions than answers :)

r.

On Mon, 7 Feb 2011, Andre Smit wrote:

Changing the way the stress is calculated gives more realistic results -

 I'm having problems relating the model's output to the actual
 test output. One issue that is confusing is that if I reduce
 the displacement interval from 0.001 to 0.0001, say and the
 time during this interval to ensure a consistent loading
 rate, then the results change significantly, both the maximum
 load at failure and the displacement at maximum load. This is
 unexpected. I can't think what the problem can be - any ideas
 on your side?

 Ok, interesting. Well, I cannot help you much at this
 point :)

 You might want to try a finer mesh, just to compare
 with the current one...

 I think the iterations are indirectly sorting
 this out. In fact I found
 that the initial stiffness(0) doesn't influence
 the results at all. I've
 attached results of compression (comp) and
 tension (tens)  tests at rate
 of 1mm/sec as well as the script in its current
 state. Based on lab
 tests, comparative compression failure at this
 rate occurs at 6.6 N/mm2
 whereas tension failure occurs at 3.0 N/mm2. So
 the script is accurately
 predicting compression failure but
 under-estimating tension failure.
 Still exploring ...


 a

 On Wed, Jan 26, 2011 at 10:19 AM, Robert Cimrman
 <cimr...@ntc.zcu.cz>
 wrote:
      Well, this flag is used essentially only in
 the generic
      solver used by simple.py. In the strain rate
 solver script we
      solve every time step regardless this flag.

      The first time step is different, however,
 as strain equals
      strain0 (no previous strain defined), and so
 initial dstrain
      is zero. Maybe this causes some problems?


 r.

 On Wed, 26 Jan 2011, Andre Smit wrote:

      Alright, that fixes it! Now I understand the
 reason for
      quasistatic! In
      the non-linear analyses, the time stepper is
 set
      programatically - is is
      possible to force solving for the first time
 step here
      as well as in:


          for ii, disp in ds:
              #my_output(dformat % (ii + 1,
 ds.n_step, disp))

              force0 = 0.0

              pb.ebcs['Load'].dofs['u.2'] = disp

              pb.ts.is_quasistatic = True  #
 <<<===========
      Force quasistatic
              pb.time_update(ts)




      On Wed, Jan 26, 2011 at 9:57 AM, Robert
 Cimrman
      <cimr...@ntc.zcu.cz>
      wrote:

           So the problem is the single dot
 outside the
      linear curve,
           right?

           What if you add

                  'quasistatic' : True,

           to the time-stepper options? As it is
 now, the
      first time
           step is not solved but taken from the
 initial
      conditions
           (unspecified = zero) and boundary
 conditions
      (nonzero!), so
           it is not in equlibrium!


      r.

      On Wed, 26 Jan 2011, Andre Smit wrote:

      Sorry - forgot to attach the figure - your
 point is
      taken re
      the nodal
      residual. Forces are calculated and summed
 on each of
      the
      nodes in the
      Top region. I checked this a while back and
 plotted in
      Paraview - the
      forces are zero at the other nodes within
 the model but
      equal
      and
      opposite in direction to the Top at the
 corresponding
      Bottom
      nodes of the
      model (as you'd expect).

      a

      On Wed, Jan 26, 2011 at 9:26 AM, Robert
 Cimrman
      <cimr...@ntc.zcu.cz>
      wrote:
           Where can I see the line? maybe you
 forgot to
      attach a
           figure?

           As for the forces, it might be better
 to compute
      them
      using a
           (inner) surface integral of stress
 instead of the
      nodal
           rezidual. Btw. you look at the nodal
 force in the
      first
      node
           of the Top region only, right? What are
 the values
      in
      the
           other nodes?

           r.

           On Wed, 26 Jan 2011, Andre Smit wrote:

                 Thanks Robert!

                 I've attached a modification that
 shows the
                 strain/stress output for the
                 elastic case of the cylinder
 under
      compression.
                 The slope of the line is
                 Young's modulus as you'd expect.
 As with our
                 non-linear analyses, the
                 forces calculated after the first
 time step
                 appear to be too high - not
                 sure what the reason for this
 is.

                 On Wed, Jan 26, 2011 at 2:25 AM,
 Robert
      Cimrman
                 <cimr...@ntc.zcu.cz>
                 wrote:
                      There was some old code in
 the
                 stress_strain() function. The
                      exception was caused by
 putting
      directly the
                 traction
                      function into the EBC
 definition,
      instead of
                 its name. The
                      attached file should work.

                      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.




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

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