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