Hth!
r.
On 12/17/10 20:12, Andre Smit wrote:
Thanks Robert!
This is brilliant - and a lot to digest.
On Fri, Dec 17, 2010 at 12:24 PM, Robert Cimrman<cimr...@ntc.zcu.cz>wrote:
I think I have something that might help, see the repository [1].
There is a new example examples/linear_elasticity/material_nonlinearity.py 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:
./simple.py examples/linear_elasticity/material_nonlinearity.py ./postproc.py cylinder.h5 -b --vector-mode=warp_norm -s 1 --layout=colrow
r. [1] https://github.com/rc/sfepy
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<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?