Hi Robert,

Thanks for your reply. I tried your suggestions with the assistance of http://sfepy.org/doc-devel/users_guide.html#defining-material-parameters but seems not working. There are some error massages which I cannot fix. Could you please provide more help. I additionally have some explanations and questions regarding my wield codes.

1. eta is dynamic scalar field, I solve the PDE with respect to eta by fipy.because I am familiar with fipy. But the elastic equation is stationary, and eta is as a known input for elastic equation.

2. I generated my own mesh(as attached) because, a) I am not able to use gmsh although learned hard, b) I want the mesh numbering or ordering agrees with fipy mesh. If the generated mesh is applicable for sfepy?

3. why when printing coors, values are [[0.21132487 0.78867513]........not real coordinates of mesh vertices?

4. I used integral = Integral('i', order=3), but what is the most appropriate integral for my problem ?

5. for stationary problem, if pb.time_update() and pb.update_materials() are necessary before pb.solve()?

My complete codes are shown below. also attached. Thank you.

Regards

Ronghai """

A 2D-stationary-inclusion-elastic problem.

equation is written in Latex

\begin{equation}

\begin{aligned}

\sigma_{ij,j} = 0\\

\text{in which}\\

\sigma_{ij} &= C_{ijkl} \epsilon_{kl}^{elastic} \\

&= C_{ijkl} (\epsilon_{kl}^{total} - \epsilon_{kl}^{pre}) \\

&= C_{ijkl} \big[\frac{1}{2}(\frac{\partial u_{k}}{\partial x_{l}} + \frac{\partial u_{l}}{\partial x_{k}})-\eta \epsilon_{kl}^{00} \big]

\end{aligned}

\end{equation}

"""

import fipy as fp

import numpy as np

from sfepy.discrete.fem import (Mesh, FEDomain, Field)

from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem)

from sfepy.solvers.ls import ScipyDirect

from sfepy.base.base import IndexedStruct

from sfepy.solvers.nls import Newton

from sfepy.terms import Term

def fem_mesh_generation(nx, ny, dx, dy, path='/home/ronghaiwu/permanent', name='/mesh.mesh'):

outfile = open(path + name,'w')

outfile.write('MeshVersionFormatted 1\t\t\t\t\t\t\nDimension\t2\t\t\t\t\t\t\nVertices\t\t\t\t\t\t\n')

outfile.write('%4d'%((nx+1)*(ny+1)))

outfile.write('\n')

for j in range(ny+1):

for i in range(nx+1):

outfile.write('\t\t')

outfile.write('%8f'%(i*dx))

outfile.write('\t\t')

outfile.write('%8f'%(j*dy))

outfile.write('\t\t')

outfile.write('%8i'%(0))

outfile.write('\n')

outfile.write('\t\t\t\t\t\t\n\t\t\t\t\t\t\nQuadrilaterals\t\t\t\t\t\t\n')

outfile.write('%4d'%(nx*ny))

outfile.write('\n')

for q in range(1,nx*ny + ny+1):

if (q % (nx+1)) != 0:

outfile.write('%4d'%(q))

outfile.write('\t')

outfile.write('%4d'%(q+1))

outfile.write('\t')

outfile.write('%4d'%(q+nx+2))

outfile.write('\t')

outfile.write('%4d'%(q+nx+1))

outfile.write('\t')

outfile.write('%4d'%(2))

outfile.write('\n')

outfile.close()

return path + name

dx = 1.

dy = dx

nx = 3 #number of elements in x direction

ny = nx #number of elements in y direction

E = 198.3 # young's modulus

nu = 0.33 #posisson's ratio

mu = 0.5*E/(1. + nu) # Lame second parameter

lambda_ = nu*E/((1. + nu)*(1.-2.*nu)) # Lame first parameter

stress00 = np.array([30.03881026, 30.03881026, 43.6856391])

################################################################################################

#fipy part, just for controling scalar field eta, therefore defining the prestress tensor field

#eta.value.shape is (nx*ny,) , eta.value is array([0., 0., 0., 0., 1., 0., 0., 0., 0.])

################################################################################################

fvm_mesh = fp.Grid2D(dx=dx,dy=dy,nx=nx,ny=ny)

fvm_x = fvm_mesh.cellCenters[0]

fvm_y = fvm_mesh.cellCenters[1]

eta = fp.CellVariable(mesh=fvm_mesh)

eta.setValue(1., where= (fvm_x > 0.5*nx*dx-0.6*dx) & (fvm_x < 0.5*nx*dx+0.6* dx) & (fvm_y > 0.5*ny*dy-0.6*dy) & (fvm_y < 0.5*ny*dy+0.6*dy))

eta_value = np.zeros((nx*ny*4))

for i in range(nx*ny):

eta_value[i*4] = eta.value[i]

eta_value[i*4+1] = eta.value[i]

eta_value[i*4+2] = eta.value[i]

eta_value[i*4+3] = eta.value[i]

eta_value.shape = (nx*ny*4,1,)

prestress_value = eta_value*stress00

###################################################################################################

fem_mesh = Mesh.from_file(fem_mesh_generation(nx, ny, dx, dy))

domain = FEDomain('domain', fem_mesh)

omega = domain.create_region('Omega', 'all')

field = Field.from_args('displacement', np.float64, 'vector', omega, 1)

u = FieldVariable('u', 'unknown', field)

v = FieldVariable('v', 'test', field, primary_var_name='u')

integral = Integral('i', order=3)

m = Material('m', lam=lambda_, mu=mu)

t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u )', integral, omega, m=m , v=v, u=u)

######################################################################################################

#two ways to define prestress according to your suggestion

######################################################################################################

#1. by a function

def prestress_function(ts, coors, mode, **kwargs):

val = prestress_value

val.shape = (nx*ny*4,1,3)

print coors

return {'prestress_value':val}

prestress = Material('prestress', function=prestress_function)

t2 = Term.new('dw_lin_prestress(prestress.function, v )',integral, omega, prestress=prestress, v=v)

#2. by a place holder constant

#prestress = Material('prestress', val=0.)

#prestress.set_all_data(prestress_value)

#t2 = Term.new('dw_lin_prestress(prestress.val, v )',integral, omega, prestress=prestress, v=v)

######################################################################################################

eq = Equation('balance', t1-t2 )

eqs = Equations([eq])

ls = ScipyDirect({})

nls_status = IndexedStruct()

nls = Newton({'i_max' : 1,'eps_a' : 1e-6,'problem' : 'nonlinear'}, lin_solver=ls, status=nls_status)

pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)

pb.time_update()

pb.update_materials()

displacement = pb.solve()

pb.save_state('/home/ronghaiwu/python_package/sfepy/linear_elasticity.vtk', displacement)

strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg')

在 2014年11月6日星期四UTC+1下午3时17分50秒，Robert Cimrman写道：

Hi Ronghai,

...

Dear sfepy users and developers,

I am using sfepy for solving an 2D inclusion elastic problem. So, I implemented based on "linear_elasticity/prestress_fibres" example in interactive way. However, the prestress in my problem is not constant, so I manually assign values. I have two questions:

1. why when the integral order is 3, the array shape of

On 11/06/2014 01:16 PM, Ronghai Wu wrote: prestress.datas[(

...

'Omega', 'i')][0]['val'] is (nx*ny,4,3,1), so when I assign by prestress_value, I have to make it the same shape. If change integral order to 1, the shape changes to (nx*ny,8,3,1), why?

This seems strange. You have quadrilateral elements, right?

./script/plot_quadratures.py -g 2_4 -n 3 sfepy: geometry: 2_4 order: 3 num. points: 4 true_order: 3

./script/plot_quadratures.py -g 2_4 -n 1 sfepy: geometry: 2_4 order: 1 num. points: 3 true_order: 2

So for the order one (or two) there should be 3 points, not 8.

...

2. When I use mesh smaller than 70*70, no problem. When 90*90 very slow. When 120*120 "MemoryError" comes out. Thanks a lot. My codes and errors are as follow:

field = Field.from_args('displacement', np.float64, 'vector', omega, 1) u = FieldVariable('u', 'unknown', field)

v = FieldVariable('v', 'test', field, primary_var_name='u') prestress_value = np.zeros((nx*ny,3,1)) m = Material('m', lam=lambda_, mu=mu) prestress = Material('prestress', val=prestress_value)

This causes the memory error - you set the (nx*ny,3,1) as a constant value parameter to the material, so it is then repeated in pb.update_materials() below into all quadrature points in the whole mesh.

I do not recommend to work with the datas attribute directly - instead, use a function:

def eval_prestress(ts, coors, mode, **kwargs): """Evaluate/return the prestress here"""

m = Material('m', function=eval_prestress)

or, alternatively, use a placeholder constant above (for example prestress = Material('prestress', val=0.0)) and then modify the datas by using Material.set_all_data(). The former way (using a user-defined function) is cleaner, though.

...

integral = Integral('i', order=3) t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u )', integral, omega, m=m , v=v, u=u) t2 = Term.new('dw_lin_prestress(prestress.val, v )',integral, omega, prestress=prestress, v=v) eq = Equation('balance', t1-t2 ) eqs = Equations([eq]) ls = ScipyDirect({}) nls_status = IndexedStruct() nls = Newton({'i_max' : 1,'eps_a' : 1e-6,'problem' : 'nonlinear'}, lin_solver=ls, status=nls_status) pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls) pb.update_materials() prestress_value = np.zeros((nx*ny,4,3,1)) ####################################################################### this is how I assign values to prestress, eta is a known scalar field. ####################################################################### for i in range(nx*ny): for j in range(4): prestress_value[i][j] = np.array([[stress00_11], [stress00_22], [ stress00_12]]) * eta.value[i] prestress.datas[('Omega', 'i')][0]['val'] = prestress_value pb.update_materials() ########################################################################

Let us know if the above helps. If not, send here the complete script, so that we can try it.

Cheers, r.

Hi Robert,

I still cannot fix the bug in my codes nor find other better ways. I am wondering have you tried my codes found the problem?

Regards Ronghai

在 2014年11月11日星期二UTC+1下午4时45分43秒，Robert Cimrman写道：

Hi Ronghai,

...

Hi Robert,

Thanks for your reply. I tried your suggestions with the assistance of http://sfepy.org/doc-devel/users_guide.html#defining-material-parameters but seems not working. There are some error massages which I cannot fix. Could you please

On 11/11/2014 03:03 PM, Ronghai Wu wrote: provide

...

more help. I additionally have some explanations and questions regarding my wield codes.

1. eta is dynamic scalar field, I solve the PDE with respect to eta by fipy.because I am familiar with fipy. But the elastic equation is stationary, and eta is as a known input for elastic equation.

So eta is then defined in the mesh vertices? Or elsewhere?

...

2. I generated my own mesh(as attached) because, a) I am not able to use gmsh although learned hard, b) I want the mesh numbering or ordering agrees with fipy mesh. If the generated mesh is applicable for sfepy?

Yes, it seems fine. The same mesh (only "transposed" - the numbering of nodes goes y-axis first) can be obtained by:

./script/blockgen.py --2d -s [4,4,4] -o out.mesh

that uses "from sfepy.mesh.mesh_generators import gen_block_mesh".

...

3. why when printing coors, values are [[0.21132487 0.78867513]........not real coordinates of mesh vertices?

The materials are evaluated in quadrature points for numerical integration and not in the mesh vertices. If your values are in the mesh vertices, you would need to interpolate them to the quadrature points, for example using a field.

...

4. I used integral = Integral('i', order=3), but what is the most appropriate integral for my problem ?

Your field is bi-linear, so 3 should be ok.

...

5. for stationary problem, if pb.time_update() and pb.update_materials() are necessary before pb.solve()?

Yes, those set up the boundary conditions and materials. A stationary problem is solved as a single step of a quasi-static solver.

...

My complete codes are shown below. also attached. Thank you.

I will look at it.

r.

Hi Robert,

I have some explanations of fipy and my idea. Please see attached PDF. If you need more information, please just ask.

Regards Ronghai

在 2014年11月14日星期五UTC+1下午2时54分04秒，Robert Cimrman写道：

Hi Ronghai,

I have just installed fipy and tried it. Could you explain more the shape of the prestress array? It is (nx * ny * 4, 3), where nx * ny is the number of cells (elements of the mesh). Why the factor 4 there? The values seem to be constant in each cell.

r.

...

Hi Robert,

I still cannot fix the bug in my codes nor find other better ways. I am wondering have you tried my codes found the problem?

Regards Ronghai

在 2014年11月11日星期二UTC+1下午4时45分43秒，Robert Cimrman写道：

...

Hi Ronghai,

On 11/11/2014 03:03 PM, Ronghai Wu wrote:

...

Hi Robert,

Thanks for your reply. I tried your suggestions with the assistance of

http://sfepy.org/doc-devel/users_guide.html#defining-material-parameters

...

...

but seems not working. There are some error massages which I cannot fix. Could you please provide more help. I additionally have some explanations and questions regarding my wield codes.

1. eta is dynamic scalar field, I solve the PDE with respect to eta by fipy.because I am familiar with fipy. But the elastic equation is stationary, and eta is as a known input for elastic equation.

So eta is then defined in the mesh vertices? Or elsewhere?

...

2. I generated my own mesh(as attached) because, a) I am not able to use gmsh although learned hard, b) I want the mesh numbering or ordering agrees with fipy mesh. If the generated mesh is applicable for sfepy?

Yes, it seems fine. The same mesh (only "transposed" - the numbering of nodes goes y-axis first) can be obtained by:

./script/blockgen.py --2d -s [4,4,4] -o out.mesh

that uses "from sfepy.mesh.mesh_generators import gen_block_mesh".

...

3. why when printing coors, values are [[0.21132487 0.78867513]........not real coordinates of mesh vertices?

The materials are evaluated in quadrature points for numerical integration and not in the mesh vertices. If your values are in the mesh vertices, you would need to interpolate them to the quadrature points, for example using a field.

...

4. I used integral = Integral('i', order=3), but what is the most appropriate integral for my problem ?

Your field is bi-linear, so 3 should be ok.

...

5. for stationary problem, if pb.time_update() and

On 11/14/2014 01:03 PM, Ronghai Wu wrote: pb.update_materials()

...

...

...

are necessary before pb.solve()?

Yes, those set up the boundary conditions and materials. A stationary problem is solved as a single step of a quasi-static solver.

...

My complete codes are shown below. also attached. Thank you.

I will look at it.

r.

Hi Robert,

Thank you very much. But I do not understand what you mean by "Note that you still need to specify some boundary conditions so that the matrix is not singular". Because it seems no boundary condition is constrained in the modified codes, but it still runs, and the result looks correct. Would you please explain more about your saying?

Regards Ronghai

在 2014年11月15日星期六UTC+1下午7时02分31秒，Robert Cimrman写道：

On 11/14/2014 06:02 PM, Ronghai Wu wrote:

...

Hi Robert,

I have some explanations of fipy and my idea. Please see attached PDF. If you need more information, please just ask.

That's fine, thanks. I have modified your file and it runs. Note that you still need to specify some boundary conditions so that the matrix is not singular.

r.

Hi Robert,

I compared the displacement calculated by sfepy and by Comsol (fine mesh, both for free boundary condition and fixed bottom boundary condition). No matter what boundary condition it is, sfepy result and Comsol result agree with each other perfectly. I mean sfepy free boundary result == Comsol free boundary result, sfepy fixed bottom result == Comsol fixed bottom result. But free boundary result =! fixed bottom result. So, I still have three questions, please forgive my ignorance on FEM: then free boundary condition can be solved. So, the question is, if it is

1. you mentioned, with free boundary condition, there will be UmfpackWarning: (almost) singular matrix!. However, I do not have, why ?
2. In Comsol, the general solver's relative tolerance is 0.0001, in this case, fixed bottom boundary condition can be solved successfully, but free boundary condition cannot. Only if I set relative tolerance to be over 1.5,

the case that, the free boundary condition actually cannot be solved theoretically, but if we do some tricks to the solver, it might be solved. So, can I say the free boundary condition result is kind of 'artificial', therefore not very physically meaningful? 3. Assuming I have a 3*3 quad mesh with following vertex coordinates (0,3) (1,3) (2,3) (3,3) (0,2) (1,2) (2,2) (3,2) (0,0) (1,1) (2,1) (3,1) (0,0) (1,0) (2,0) (3,0) if the displacement result is the displacement at these vertices?

Assuming I use strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg'), if the strain result is the strain at cell centers?

Regards Ronghai

在 2014年11月17日星期一UTC+1下午3时47分20秒，Robert Cimrman写道：

...

Hi Robert,

Thank you very much. But I do not understand what you mean by "Note that you still need to specify some boundary conditions so that the matrix is not singular". Because it seems no boundary condition is constrained in

On 11/17/2014 11:04 AM, Ronghai Wu wrote: the

...

modified codes, but it still runs, and the result looks correct. Would you please explain more about your saying?

Regards Ronghai

sfepy: nls: iter: 0, residual: 7.497687e+01 (rel: 1.000000e+00) /usr/local/lib/python/dist-packages/scikits/umfpack/umfpack.py:720: UmfpackWarning: (almost) singular matrix! (estimated cond. number: 5.91e+15) warnings.warn(msg, UmfpackWarning) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: nls: iter: 1, residual: 4.099472e-14 (rel: 5.467649e-16)

Yes, it is solved, because the direct solver somehow deals with the singular matrix (see the warning above). But without any boundary conditions the matrix coming from the linear elasticity is singular - in 2D, it has three zero eigenvalues for the two rigid body displacements and one rotation. These need to be fixed to have a unique solution. So, for example, fix the bottom left corner vertex completely ({'u.all' : 0.0}), and fix the 'y' displacement component of the bottom right vertex ({'u.1' : 0.0}). See [1], where two ways (a constant and by a function) of specifying the Dirichlet boundary conditions (EssentialBC) are shown for the interactive use. Let us know if you need more help.

r. [1]

http://sfepy.org/doc-devel/examples/standalone/interactive/linear_elasticity...

Hi Robert,

Now I encounter a more complicated case. What if the stiffness tensor is not constant, but a function of eta, as shown in PDF "non-constant elastic tensor.pdf". Is it possible to implement it based on the modified code "Ronghai.py"?

Regards Ronghai

在 2014年11月15日星期六UTC+1下午7时02分31秒，Robert Cimrman写道：

On 11/14/2014 06:02 PM, Ronghai Wu wrote:

...

Hi Robert,

I have some explanations of fipy and my idea. Please see attached PDF. If you need more information, please just ask.

That's fine, thanks. I have modified your file and it runs. Note that you still need to specify some boundary conditions so that the matrix is not singular.

r.

Hi Robert,

Should I prescribe Dijkl to each vertex like the prestress, or only need to prescribe to each element? Besides, if using strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg'), the shear value in strain is sigma_xy or 2*sigma_xy?

Regards Ronghai

在 2014年11月26日星期三UTC+1下午3时14分44秒，Robert Cimrman写道：

On 11/26/2014 02:39 PM, Ronghai Wu wrote:

...

Hi Robert,

Now I encounter a more complicated case. What if the stiffness tensor is not constant, but a function of eta, as shown in PDF "non-constant elastic tensor.pdf". Is it possible to implement it based on the modified code "Ronghai.py"?

So eta is a known function of position? This can be done by a material parameter given by a user-defined function again [1], just like you do the prestress.

The function would be something like:

def get_d(ts, coors, mode=None, **kwargs): if mode == 'qp': x = coors[:, 0] y = coors[:, 1]

     value = ... 

     return {'D' : value} 

and in the code, use something like:

d_fun = Function('d_fun', get_d) m = Material('m', function=d_fun)

Also, if you need to work with the elastic tensors, use dw_lin_elastic term instead of dw_lin_elastic_iso.

r. [1] http://sfepy.org/doc-devel/users_guide.html#defining-material-parameters

...

Regards Ronghai

在 2014年11月15日星期六UTC+1下午7时02分31秒，Robert Cimrman写道：

...

On 11/14/2014 06:02 PM, Ronghai Wu wrote:

...

Hi Robert,

I have some explanations of fipy and my idea. Please see attached PDF. If you need more information, please just ask.

That's fine, thanks. I have modified your file and it runs. Note that

you

...

still need to specify some boundary conditions so that the matrix is not singular.

r.