Thanks Robert!

This is brilliant - and a lot to digest.  

On Fri, Dec 17, 2010 at 12:24 PM, Robert Cimrman <> wrote:
I think I have something that might help, see the repository [1].

There is a new example examples/linear_elasticity/ which illustrates some sort of dependence of 'mu' on the strain rate magnitude. If it is similar to what you need, it would be cool if you modified it so that it uses a real constitutive relation just like you described below.

I have not updated the documentation yet to reflect the change of arguments passed to the material function (see commit 'rewrite updating of material parameters')...

Note that there is a hack to compute the strain rate, as the variables with history need to be updated w.r.t the new State class - the current implementation is not good.

To test it, just do:

./ examples/linear_elasticity/
./ cylinder.h5 -b --vector-mode=warp_norm -s 1 --layout=colrow


On 12/17/10 14:38, Andre Smit wrote:
Yes, that's it. I'm thinking to use an anisotropic model perhaps with time
dependency and then adjust the material properties based on the strain rates
in each of the elements per time step at which point the normal and shear
strengths of the elements (also strain rate dependent) are calculated and if
the stresses in the elements exceed the strengths then the element fails and
is given a very low modulus to indicate failure. If the element prevails,
the modulus of the element is adjusted based on the strain rate in the
element and the process moves to the next time step.

On Fri, Dec 17, 2010 at 3:01 AM, Robert Cimrman<>  wrote:

OK, I can see now what are you trying to do. You would like to implement a
material nonlinearity (material parameters depending on state) and linearize
it by taking state from the previous time step in each time step, right?

I guess it can be done the way you try, i.e. each element has its own
region/material/term, but that is going to be utterly slow, as it would
completely negate the fact that the terms are evaluate in C for a number of
elements simultaneously.

So again, it's a good catch, as it is a common need and the code should
support that in a clean and easy-to-use manner. The solution is to pass the
state to the material update function, that a user could write. I am working
on this now.

Then you would evaluate all just over Omega (strain, stress, ...).


On 12/16/10 18:18, Andre Smit wrote:


Yes, the materials are defined as:

# Dynamically assign materials
for i in range(NELS):
    exec "material_%s = { \
              'name' : 'mat_%s', \
            'values' : { \
                'lam' : youngpoisson_to_lame(young, poisson)[0], \
                 'mu' : youngpoisson_to_lame(young, poisson)[1], \
                  'D' : stiffness_tensor_youngpoisson(2, young,

So I have NELS number of materials in the model.
young will change depending on the strain rate in the individual elements
(I'm assuming the material is strain rate dependent) - so in effect the
material parameters change every time step when they are updated based on
young=f(strain rate). Each element is assigned a region and the equation

# Fixed regions
regions = {
     'Omega' : ('all', {}),
      'Left' : ('nodes in (x<   0.001)', {}),
    'Bottom' : ('nodes in (y<   0.001)', {}),
       'Top' : ('node 2', {}),

# Dynamically assigned regions
for i in range(NELS):
    exec "region_%s =  {'name' : 'ele_%s', 'select' : 'element

# Dynamically assigned equation:
for i in range(NELS):

equations = {
   'balance_of_forces' : Eq

This brings up the issue of having different materials defined in the
- I can calculate the strain over Omega as before:
strain=pb.evaluate('de_cauchy_strain.2.Omega(u)') but stress needs to be
done for the individual elements (and materials) separately, right?.

s='de_cauchy_stress.2.ele_%s(mat_%s.D, u)' %i
stress = pb.evaluate(s)

So I don't think 'de_volume_average_mat.2.Omega( material, parameter )'
work as it needs a material parameter - but in this case the material is
defined over Omega. Am I missing something?


On Thu, Dec 16, 2010 at 3:49 AM, Robert Cimrman<>

 Hi Andre,

I have some questions: how do you plan to use this D matrix? Do you
your material by a function? Note that to use a material parameter, it
to be evaluated in quadrature points, i.e. the coordinates passed into a
material function.

This code


assumes, that your state contains a single variable, which has the shape
like Dmatrix.ravel(). I doubt this is the case. SfePy unfortunately did
check this - I have just made a fix, so get the fresh git version.

Otherwise I would just use the term de_volume_average_mat, something like
the following code I use for the permeability tensor:

perm = pb.evaluate('de_volume_average_mat.2.Omega( m.K, pc )',
                  copy_materials=False, verbose=False)
out['pperm'] = Struct(name='K', mode='cell', data=perm, dofs=None)

Otherwise, it's great you play/work with the code, it really makes it
better, if it is actively used by other people than me. :)


On 12/16/10 06:27, Andre Smit wrote:

 Hi guys

I've generated an array comprising material D matrix or stiffness tensor
values. It looks like this - where each block represents the material
for a cell or element (yes, they're all the same but will change in
subsequent time steps).

[[[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]

 [[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]

 [[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]

 [[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]

 [[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]

 [[ 1604.9382716    864.19753086     0.        ]
  [  864.19753086  1604.9382716      0.        ]
  [    0.             0.           370.37037037]]]

I'd like to append this array to the vtk file generated by SfePy which
contains displacements, stress and strains, etc.
I was using the following to generate the array:

    for i in range(NELS):
    print "~~~~~~~~~~~~~~~~~~~~~~~~~~~"
    print Dmatrix


but get the following *numpy* error:

Traceback (most recent call last):
  File "/home/grassy/sfepy/", line 120, in<module>
  File "/home/grassy/sfepy/", line 117, in main
  File "/home/grassy/sfepy/sfepy/applications/", line 29,
    return **kwargs )
  File "/home/grassy/sfepy/sfepy/applications/", line 112,
  File "/home/grassy/sfepy/sfepy/solvers/", line 216, in
  File "/home/grassy/sfepy/sfepy/solvers/", line 148, in
  File "/home/grassy/sfepy/sfepy/fem/", line 535, in
    out = post_process_hook( out, self, state, extend = extend )
  File "", line 102, in strain_rate
  File "/home/grassy/sfepy/sfepy/fem/", line 188, in
    var_info, extend)
  File "/home/grassy/sfepy/sfepy/fem/", line 587, in
  File "/home/grassy/sfepy/sfepy/fem/", line 1468, in
    (self.n_dof / self.n_components, self.n_components))
  File "/usr/lib/python2.7/site-packages/numpy/core/",
170, in reshape
    return reshape(newshape, order=order)
ValueError: total size of new array must be unchanged

Can I use existing SfePy methods to write this array to the vtk or
rather write my own?

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