Hi Robert,

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

"""

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

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

Hi Ronghai,

On 11/06/2014 01:16 PM, Ronghai Wu wrote:
> 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 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?

sfepy: geometry: 2_4 order: 3 num. points: 4 true_order: 3

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.