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


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