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:
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((nxny,3,1)) m = Material('m', lam=lambda_, mu=mu) prestress = Material('prestress', val=prestress_value) 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((nxny,4,3,1))
###### #
this is how I assign values to prestress, eta is a known scalar field.
###### #
for i in range(nxny): 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()
###### #
following is error: sfepy: updating materials... sfepy: m sfepy: prestress Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/local/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py" , line 580, in runfile execfile(filename, namespace) File "/home/ronghaiwu/temporary/test.py", line 304, in <module> main() File "/home/ronghaiwu/temporary/test.py", line 179, in main pb.update_materials() File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/problem.py" , line 522, in update_materials problem=self, verbose=verbose) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/equations.py" , line 313, in time_update_materials verbose=verbose) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/materials.py" , line 70, in time_update mat.time_update(ts, equations, mode=mode, problem=problem) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/materials.py" , line 339, in time_update self.update_data(key, ts, equations, term, problem=problem) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/materials.py" , line 251, in update_data self.extra_args) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/functions.py" , line 34, in __call__ return self.function(*args, _kwargs) File "/home/ronghaiwu/mypylibs/lib/python2.7/site-packages/sfepy/discrete/functions.py" , line 66, in get_constants out[key] = nm.tile(val, (coors.shape[0], 1, 1)) File "/usr/local/lib/python2.7/dist-packages/numpy/lib/shape_base.py", line 860, in tile c = c.reshape(-1, n).repeat(nrep, 0) MemoryError
Regards Ronghai
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:
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.
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(nxny): 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,
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.
I used integral = Integral('i', order=3), but what is the most appropriate integral for my problem ?
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_ = nuE/((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.5nxdx-0.6dx) & (fvm_x < 0.5nxdx+0.6 dx) & (fvm_y > 0.5nydy-0.6dy) & (fvm_y < 0.5nydy+0.6dy))
eta_value = np.zeros((nxny4))
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 = (nxny4,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 = (nxny4,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,
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:
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.
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(nxny): 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 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.
So eta is then defined in the mesh vertices? Or elsewhere?
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".
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.
Your field is bi-linear, so 3 should be ok.
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 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.
So eta is then defined in the mesh vertices? Or elsewhere?
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".
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.
Your field is bi-linear, so 3 should be ok.
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 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.
On 11/14/2014 01:03 PM, Ronghai Wu wrote:
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.
So eta is then defined in the mesh vertices? Or elsewhere?
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".
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.
Your field is bi-linear, so 3 should be ok.
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.
On 11/14/2014 01:03 PM, Ronghai Wu wrote:
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.
So eta is then defined in the mesh vertices? Or elsewhere?
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".
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.
Your field is bi-linear, so 3 should be ok.
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.
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,
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.
On 11/17/2014 11:04 AM, Ronghai Wu wrote:
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
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,
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:
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写道： >
On 11/17/2014 11:04 AM, Ronghai Wu wrote:
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
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...
On 11/21/2014 11:29 AM, Ronghai Wu wrote:
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:
Maybe you have not the umfpack solver installed. In that case superlu (which is budled in scipy) is used.
With no boundary conditions, the solutions can contain _arbitrary_ translations and rotations of the whole domain as a rigid body. So the agreement of comsol and sfepy is a coincidence, maybe the linear solver uses the same strategy when dealing with singular matrices.
Yes.
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?
It is an average in the elements, so yes, it can be said that it is a strain at cell centers.
r.
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.
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.
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.
On 11/26/2014 07:10 PM, Ronghai Wu wrote:
Hi Robert,
Should I prescribe Dijkl to each vertex like the prestress, or only need to prescribe to each element?
Ideally, it would be to evaluate/prescribe Dijkl in the coordinates that are passed to the material function - those coordinates are the quadrature points, where the Dijkl is needed. In the last version of your script, the prestress is not prescibed in vertices. Instead, it is prescribed to each element - it is just repeated to all the points that fall to each element. So if having Dijkl piece-wise constant (over elements) is ok for you, you can define in exactly the same way as the prestress - define it in elements (centres) and repeat n_qp times using np.repeat().
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?
See [1].
[1] http://sfepy.org/doc-devel/src/sfepy/terms/terms_elastic.html#sfepy.terms.te...
Cheers, r.
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