one problem in ''biot.py"

[Kid Guo]

Hi,

I have a problem based on  "../sfepy/examples/multi_physics/biot.py".

In the attached scripy 'block2.py', I have calculated the stress and strain according to the displacement, at the 111 and 112 lines you can see it.

Now I want to output the K of every cell, I gave a value of k at the 64 line in this script, it is a 3*3 array.

I don't know which function or equation to use, so I was stucked at the 113 line. 

if I change the K value,  will where be some difference results?

Thank you advance

K.

[Kid Guo]

在 2016年12月26日星期一 UTC+8下午8:52:44，Kid Guo写道：

Hi,

I have a problem based on  "../sfepy/examples/multi_physics/biot.py".

In the attached scripy 'block_hc.py', I have calculated the stress and strain according to the displacement, at the 111 and 112 lines you can see it.

[Kid Guo]

Where Ki is the permeability of the i-th compartment and Kij is the general diffusion permeability.

how to ouput the Kij of each cell.

在 2016年12月26日星期一 UTC+8下午8:52:44，Kid Guo写道：

Hi,

[Robert Cimrman]

Hi Kid,

you can use

hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m.k, p)', m=m,
mode='el_avg')

...

out['hydraulic_conductivity'] = Struct(name='output_data', mode='cell',
data=hydraulic_conductivity)


Essentially, ev_integrate_mat can be used to integrate any material parameter,
check [1] to see what the mode argument means.

r.

[1]
http://sfepy.org/doc-devel/src/sfepy/terms/terms_basic.html#sfepy.terms.terms_basic.IntegrateMatTerm

[Kid Guo]

Hi Robert,

Thank you very much, I have output the value of Kij in my script.

In the script, I have give a initial value of K, of course , I assumed it is a homogeneous material, so it has the same value of every cell.

Now I have new thinking, if it is a heterogeneous material, that Kij value are different of every cell.

I want to give different K value in every cell, for example,

in the cell 1, the Kij value is K1 = nm.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], dtype = nm.float64)

in the cell 2, the Kij value is K2 = nm.array([[1.1, 0.0, 0.0], [0.0, 1.2, 0.0], [0.0, 0.0, 10.1]], dtype = nm.float64)

etc.

The new Kij like this:

K2 = nm.array([[[[ 1.,  0.,  0.],

                 [ 0,  1.,  0.],

                 [ 0.,  0.,  1.]]],

   

               [[[ 1.1,  0.,  0.],

                 [ 0,  1.2,  0.],

                 [ 0.,  0.,  10.1]]]], dtype = nm.float64)

when I give the new value to the problem,it has a error.

"ValueError: material parameter array must have three dimensions! ('k' has 4)"

or the K is

K3 = nm.array([[[1.0, 0.0, 0.0], 

                [1.0, 0.0, 0.0], 

                [0.0, 0.0, 1.0]],

               [[1.0, 0.0, 0.0], 

                [0.0, 2.0, 0.0], 

                [0.0, 0.0, 3.0]]], dtype = nm.float64)

it said:

"ValueError: incompatible shapes! ((27, 8, 3, 8), out: (27, 1, 3, 3), arr: (54, 8, 3, 3))"

what should I do?

Any idea?

Thank you advance.

K.

在 2016年12月28日星期三 UTC+8下午9:14:36，Robert Cimrman写道：

[Robert Cimrman]

Hi Kid,

If you want to have the permeability varying with position, you have to define
the material parameter using a function [1]- it seems that you assigned your
new K2 directly to the material constructor, right?

r.
[1[
https://groups.google.com/forum/?_escaped_fragment_=topic/sfepy-devel/6QWAwCUNeWE#!topic/sfepy-devel/6QWAwCUNeWE

[Kid Guo]

Hi Robert,

Thank you for your answer.

The initial K value is 

K = nm.array([[2.0, 0.2, 0.0], [0.2, 1.0, 0.0], [0.0, 0.0, 0.5]], dtype = nm.float64).

Assuming my mesh file has 211  tetrahedra element,

Is it for the every element in the material, or the actual K total array like this

K = nm.array([[[2.0, 0.2, 0.0], 
              [0.2, 1.0, 0.0], 
              [0.0, 0.0, 0.5]], 

              [[2.0, 0.2, 0.0], 
              [0.2, 1.0, 0.0], 
              [0.0, 0.0, 0.5]],

              ......,
              [[2.0, 0.2, 0.0], 
              [0.2, 1.0, 0.0], 
              [0.0, 0.0, 0.5]]],
              dtype = nm.float64).
K.shape = (211,3,3)

Is it right?

If it is right, now I want to give a total K array. for example.

K = nm.array([[[2.0, 0.2, 0.0], 
              [0.2, 1.0, 0.0], 
              [0.0, 0.0, 0.5]], 

              [[3.0, 0.2, 0.0], 
              [0.2, 2.0, 0.0], 
              [0.0, 0.0, 1.0]],

              ......,
              [[10.0, 0.2, 0.0], 
              [0.2, 2.0, 0.0], 
              [0.0, 0.0, 0.4]]],
              dtype = nm.float64).
K.shape = (211,3,3)

I change the tensor k value in every tetrahedra element. 

then re-solve this problem. 

Is it possible?How can I do it?

Thank you advance.

Regards,

Kid

在 2016年12月30日星期五 UTC+8上午4:26:02，Robert Cimrman写道：

[Robert Cimrman]

Hi Kid,

first, if you do (e.g. before pb.solve() call in your block_hc.py script):

pb.update_materials()

you could see the following:

ipdb> pb.get_materials()['m'].datas[('Omega', 2)]['k'].shape
(211, 4, 3, 3)
ipdb> pb.get_materials()['m'].datas[('Omega', 1)]['k'].shape
(211, 1, 3, 3)

the first shape corresponds to the order 2 integral, the second shape to the
order 1 integral. -> You may want to use two materials instead, one for each
order of integration, so that the data are not uselessly duplicated.

This also means, that the material does not require a value per cell, but a
value per quadrature point. Put the following stub after your definition of 'm'

def get_pars(ts, coors, mode=None, **kwargs):
if mode == 'qp':
out = {}

print(coors.shape)
print(domain.shape.n_el)
n_qp = coors.shape[0] // domain.shape.n_el
print n_qp

k = nm.zeros((coors.shape[0], 3, 3))

# Fill your values here...

out['k'] = k
return out

m2 = Material('m2', function=get_pars)

edit it to use your values, and use m2.k instead of m.k everywhere (remove k
from m).

r.

[Kid Guo]

Hi Robert,

First thanks for your reply.

In my scripy, I have replaced the m.k with new m2.k, it worked.

At the 98line of the script, the new term is "t22 = Term.new('dw_diffusion(m2.k, q, p)', integral1, omega, m2=m2, q=q, p=p)"

Now I output the 'hydraulic_conductivity' and save it in the 'block_pointload_script_hc_strain.vtk'. I can view the result through the 142 and 143 line in the script.

I have a new question, if I output 'hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m.k, p)', m=m, mode='el_avg')' ,all k value is the same and the result can view it,the hydraulic_conductivity is 0.104, you can see it in the attached file.

and if I output 'hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m2.k, p)', m2=m2, mode='el_avg')' ,the output value should be different, but the hydraulic_conductivity is 0,you can see it in the attached file.

What is wrong?

Thank you advance. 
Regards
Kid 

在 2016年12月30日星期五 UTC+8下午6:14:38，Robert Cimrman写道：

[Kid Guo]

Hi, 

This is the supplement.

If I set 

't22 = Term.new('dw_diffusion(m.k, q, p)', integral1, omega, m=m, q=q, p=p) and

hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m.k, p)', m=m, mode='el_avg')',

the output result is the m_k.png file.

then 

if I set

't22 = Term.new('dw_diffusion(m2.k, q, p)', integral1, omega, m2=m2, q=q, p=p) and

hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m2.k, p)', m2=m2, mode='el_avg')

'

the output result is the m2_k.png file.

By comparing the two file, I can see the displacement,cauchy_strain and hydraulic_conductivity is different, that is to say, when I get the m2.k, the Kij is different of each element.it works.But the view result is 0.

Why?

Thank you advance. 
Regards, 

Kid 

在 2016年12月30日星期五 UTC+8下午6:14:38，Robert Cimrman写道：

Hi Kid,

[Robert Cimrman]

Hi Kid,

It is a mayavi tensor visualization issue. In paraview I can see the Kij
different in each cell (with the mySample.mesh you sent previously).

The problem is, that in mayavi, tensors as a whole can be displayed only with
the following three modes: component, determinant, effective_stress. SfePy uses
effective_stress, which is zero for a diagonal tensor.

You can workaround that, for example, by computing the matrix norm of Kij
yourself and store that in the results file, maybe alongside the full Kij:

out['hc_norm'] = Struct(name='output_data', mode='cell',
data=nm.linalg.norm(hydraulic_conductivity,
axis=(2,3), keepdims=True))

BTW. you can select data to display, as:

view(is_scalar_bar=True, is_wireframe=True, is_3d=True,
only_names=['hc_norm'])

r.

On 01/09/2017 07:12 AM, Kid Guo wrote:
> Hi,
>
> This is the supplement.
> If I set
> 't22 = Term.new('dw_diffusion(m.k, q, p)', integral1, omega, m=m, q=q, p=p)
> and
> hydraulic_conductivity = ev('ev_integrate_mat.3.Omega(m.k, p)', m=m, mode=
> 'el_avg')',
> the output result is the m_k.png file.
> then
> if I set
> 't22 = Term.new('dw_diffusion(m2.k, q, p)', integral1, omega, m2=m2, q=q, p=p)

> *and*

[Kid Guo]

在 2017年1月10日星期二 UTC+8下午4:52:28，Robert Cimrman写道：

Hi,

I have a new problem about the biot.py. it is a small model, but now if the model is very big, such as the z-coordinate is 100m or bigger. I have to consider the init stress in the model.
I try to load a stress(the weight of the material) in the top plane of the model, and get a result, but I don't know whether it is right?
And if I consider pore pressure and the the weight of the model,what should I do?

There are two files in the attachment, one is biot_add_gravity.py and the other is the mid_exca_z10.mesh. I set the max z coordinate is 10m.
I add a equation in the 105 line of the file and some additional terms.

Any suggestions?
Thank you advance.
Kid.

[Kid Guo]

在 2017年6月7日星期三 UTC+8下午12:28:31，Kid Guo写道：