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

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.

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
rate, then the results change significantly, both the maximum
unexpected. I can't think what the problem can be - any ideas

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

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

'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
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,
its name. The
attached file should work.

r.

--
You received this message because you are
subscribed to
the
Groups "sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email
to
For more options, visit this group at

--
Andre

--
You received this message because you are
subscribed to
the
"sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email
to
For more options, visit this group at

--
You received this message because you are
subscribed to
Groups "sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email
to
For more options, visit this group at

--
Andre

--
You received this message because you are
subscribed to
"sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email
to
For more options, visit this group at

--
You received this message because you are
Groups "sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email to
For more options, visit this group at

--
Andre

--
You received this message because you are
"sfepy-devel" group.
To post to this group, send email to
To unsubscribe from this group, send email to
For more options, visit this group at

--
You received this message because you are subscribed to the
To post to this group, send email to
To unsubscribe from this group, send email to
For more options, visit this group at

--
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
For more options, visit this group at

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