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.

a

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.comsfepy-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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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<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<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<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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

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

...

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:

...

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<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<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<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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

-- Andre

-- Andre

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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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.

On Tue, 25 Jan 2011, Andre Smit wrote:

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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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

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.comsfepy-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.comsfepy-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.comsfepy-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.comsfepy-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.comsfepy-devel%...@googlegroups.com . For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

-- Andre

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:

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

r.

On Wed, 26 Jan 2011, Andre Smit wrote:

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.

-- 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 - my assumption regarding Hooke's law was perhaps invalid.

On Mon, Feb 7, 2011 at 3:37 PM, Andre Smit freev...@gmail.com wrote: 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.cz wrote:

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

r.

On Wed, 26 Jan 2011, Andre Smit wrote:

  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.

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

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

On Tue, 8 Feb 2011, Andre Smit wrote:

Your reasoning regarding the time step size makes sense. Yep, I changed to a full-scale model instead of the 1/8th just to make sense of the strange behavior. I'll change back.

Or fix the two nodes as I said, that should give the same results.

The main change in the script is the calculation of stress, which I now do in the time loop:

line 387: stress = pb.evaluate('de_cauchy_stress.2.Omega(quasi.D, u)', copy_materials=False)

This is then passed to update_stiffness and repeated to get the stresses at the qp's:

ￜￜￜￜￜￜￜ stress=nm.repeat(stress,4) ￜￜￜￜￜￜￜ stress.shape=strain.shape ￜￜￜￜￜￜￜ szz1 = avgqp(stress[..., 2:3, :], problem) ￜ I'm doing this instead of calculating stress via: szz = ezz * smix

If you need the Cauchy stress at the qp's, use dq_cauchy_stress instead of de_cauchy_stress :)

Not sure of the relevance of "copy_materials=False" in the stress function.

You may try to remove this argument to be sure it does not cause some side-effects. The purpose of setting it to False is to speed things up in case one needs to evaluate a term with a material parameter, that was already used during solution. The material parameters then can be reused, instead of being evaluated again.

I also opted to only calculate the D matrix instead of using Lame's parameters, which seems to speed things up.

Yes, if you need that anyway, it's better.

r.

On Tue, Feb 8, 2011 at 2:41 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote: 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 -
  my assumption regarding Hooke's law was perhaps
  invalid.



  On Mon, Feb 7, 2011 at 3:37 PM, Andre Smit
  <freev...@gmail.com> wrote:
  ￜ ￜ ￜ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.cz> wrote:

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


  r.

  On Wed, 26 Jan 2011, Andre Smit wrote:

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



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




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