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') 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?
./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.