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

Not sure of the relevance of "copy_materials=False" in the stress function. I also opted to only calculate the D matrix instead of using Lame's parameters, which seems to speed things up.

thx


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