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?