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


On Tue, 25 Jan 2011, Andre Smit wrote:

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

[grassy@myhost CFRAC]$ ~/sfepy/
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,, 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/", line 113, in <module>
  File "/home/grassy/sfepy/", line 110, in main
  File "/home/grassy/sfepy/sfepy/applications/", line 29, in call_basic
    return **kwargs )
  File "/home/grassy/sfepy/sfepy/applications/", line 121, in call
  File "/home/grassy/sfepy/sfepy/solvers/", line 217, in solve_direct
  File "/home/grassy/sfepy/sfepy/solvers/", line 137, in solve_evolutionary_op
    for ts, state in time_solver( state0 ):
  File "/home/grassy/sfepy/sfepy/solvers/", line 128, in __call__
    state = step_fun( self.ts, state0, *step_args )
  File "/home/grassy/sfepy/sfepy/solvers/", line 71, in time_step_function
    problem.time_update( ts )
  File "/home/grassy/sfepy/sfepy/fem/", line 460, in time_update
    functions, create_matrix)
  File "/home/grassy/sfepy/sfepy/fem/", line 422, in update_equations
    self.equations.time_update(self.ts, ebcs, epbcs, lcbcs, functions)
  File "/home/grassy/sfepy/sfepy/fem/", line 215, in time_update
    self.variables.equation_mapping(ebcs, epbcs, ts, functions)
  File "/home/grassy/sfepy/sfepy/fem/", line 239, in equation_mapping
    var.equation_mapping(bcs, var_di, ts, functions)
  File "/home/grassy/sfepy/sfepy/fem/", line 1372, in equation_mapping
    self.eq_map.map_equations(bcs, self.field, ts, functions, warn=warn)
  File "/home/grassy/sfepy/sfepy/fem/", 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?


On Tue, Jan 25, 2011 at 2:52 AM, Robert Cimrman <> 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

     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

     I think that makes sense.

     So good luck with the actual data fitting :)


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.


     On Mon, Jan 24, 2011 at 12:38 PM, Andre Smit <>
          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.


     On Mon, Jan 24, 2011 at 11:26 AM, Andre Smit
     <> 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!


     On Mon, Jan 24, 2011 at 9:08 AM, Robert Cimrman
     <> wrote:
          On Mon, 24 Jan 2011, Andre Smit wrote:


                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
                e *= 2 # With 1/8 model the applied
                deformation rate is half of actual

     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

     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?


          On Mon, Jan 24, 2011 at 2:05 AM, Robert
          Cimrman <>
               ok, so the patch is in... the script
          should work


          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

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


               ----- Reply message -----
               From: "Andre Smit"
               To: <>
               Subject: Conference
               Date: Fri, Jan 21, 2011 02:47


               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

               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
               c) I iterate the problem until the
          force converges and
               then increment the
               d) The element stiffness matrix (Ele)
          is saved to file
               and loaded at the
               start of the subsequent displacement

               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



               On Wed, Jan 19, 2011 at 10:22 AM,
          Robert Cimrman
                    On Wed, 19 Jan 2011, Andre Smit

                    Robert - I'm trying to change the
                    in the traction
                    function when the force
          stabilizes but for some
                    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()
                        f =
                        f.shape =
          (pb.domain.mesh.n_nod, 3)
                        p = []
                        for n in
                        p = nm.array(p)
                        if abs(old_force-p.sum()) <
                            disp -= ddisp

                    def traction(ts, coors, bc=None):
                        global disp
                        #disp = -(ts.dt*(ts.step+1))
                        print "traction disp: ", disp
                        val =
                        return val

               I would use brute force:

               globals()['disp'] = ...

               You cannot really change value of a
          global that is
               A common idiom is to use list instead,

               disp = [0]

               def traction(ts, coors, bc=None):
                  disp[0] = ...

                         What is the meaning of
          'smix' stiffness, for
                    strain rate

                    I'm not sure.

               That might be the problem causing the


You received this message because you are subscribed to the Google 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



You received this message because you are subscribed to the Google 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

You received this message because you are subscribed to the Google 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