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 <cimr...@ntc.zcu.cz> 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, ...).

r.



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

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,
poisson),},}"%(i+1,i)

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

# 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
%s'}"%(i+5,i,i)

# Dynamically assigned equation:
Eq=""
for i in range(NELS):
    Eq+="dw_lin_elastic_iso.2.ele_%s(mat_%s.lam,mat_0.mu,v,u)+"%(i,i)
Eq=Eq.rstrip("+")
Eq+="=0"

equations = {
   'balance_of_forces' : Eq
}


This brings up the issue of having different materials defined in the model
- 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 )' will
work as it needs a material parameter - but in this case the material is not
defined over Omega. Am I missing something?


Andre




On Thu, Dec 16, 2010 at 3:49 AM, Robert Cimrman<cimr...@ntc.zcu.cz>  wrote:

Hi Andre,

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

This code


     Dstate=pb.create_state()
     Dstate.set_full(Dmatrix)

assumes, that your state contains a single variable, which has the shape
like Dmatrix.ravel(). I doubt this is the case. SfePy unfortunately did not
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 )',
mode='el_avg',
                  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. :)

r.


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

    Dlist=[]
    for i in range(NELS):
        Dlist.append(pb.conf.materials.values()[i].values['D'])
    Dmatrix=array(Dlist)
    print "~~~~~~~~~~~~~~~~~~~~~~~~~~~"
    print Dmatrix
    Dstate=pb.create_state()
    Dstate.set_full(Dmatrix)
    outD=Dstate.create_output_dict()

out['D_matrix']=Struct(name='output_data',mode='cell',data=outD,dofs=None)


but get the following *numpy* error:

Traceback (most recent call last):
  File "/home/grassy/sfepy/simple.py", line 120, in<module>
    main()
  File "/home/grassy/sfepy/simple.py", line 117, in main
    app()
  File "/home/grassy/sfepy/sfepy/applications/application.py", line 29, in
call_basic
    return self.call( **kwargs )
  File "/home/grassy/sfepy/sfepy/applications/simple_app.py", line 112, in
call
    post_process_hook_final=self.post_process_hook_final)
  File "/home/grassy/sfepy/sfepy/solvers/generic.py", line 216, in
solve_direct
    nls_status=nls_status)
  File "/home/grassy/sfepy/sfepy/solvers/generic.py", line 148, in
solve_evolutionary_op
    ts=ts)
  File "/home/grassy/sfepy/sfepy/fem/problemDef.py", line 535, in
save_state
    out = post_process_hook( out, self, state, extend = extend )
  File "its.py", line 102, in strain_rate
    outD=Dstate.create_output_dict()
  File "/home/grassy/sfepy/sfepy/fem/state.py", line 188, in
create_output_dict
    var_info, extend)
  File "/home/grassy/sfepy/sfepy/fem/variables.py", line 587, in
state_to_output
    fill_value=fill_value))
  File "/home/grassy/sfepy/sfepy/fem/variables.py", line 1468, in
create_output
    (self.n_dof / self.n_components, self.n_components))
  File "/usr/lib/python2.7/site-packages/numpy/core/fromnumeric.py", line
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 should
I
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 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