Hi there
I try to solve a very simple problem with getfem: a bended bar.
In fact I consider a bar fixed on one side and exerted to a force (in
negative z direction) on the other one. Clearly on the fixed end I have
Dirichlet boundary conditios with \vec{u}=0 and some force on the other
side. Therefore, I expect that the maximum value in uz would be zero
(due to the Dirichlet boundary condition).
However, in the solution I get the following:
...
solve done!
get displacement field
ux min: -1.090228e-08, max: 1.090042e-08
uy min: -4.643458e-07, max: 4.642960e-07
uz min: -5.047865e-06, max: 1.396889e-09
...
Obviously, the maximum in uz direction is not zero. Attached to this
mail you will find the small python program performing the calculation.
Is there any obvious mistake I did in the program?
best regards
Eugen Wintersberger
PS: I'm quite new to getfem and to FEM in general ;) so I hope my
question is not too stupid ;))
------------------------------------------------------------------------
#perform a simple calculation - interactive namespace from
#the ipython shell
import salome_mesh as smesh
import getfem
import numarray
[mesh,groups] = smesh.unv2getfem("simplebar.unv",3);
print groups
#some general nodes
degree = 1.0;
NEUMANN_BOUNDARY = 1;
DIRICHLET_BOUNDARY = 2;
#have to assemble an elasticity problem now
mfu = getfem.MeshFem(mesh,3);
mfe = getfem.MeshFem(mesh,1);
mfdata = getfem.MeshFem(mesh,1);
mim = getfem.MeshIm(mesh,getfem.Integ("IM_TETRAHEDRON(5)"));
mfu.set_fem(getfem.Fem("FEM_PK(3,%i)" %(degree)));
mfe.set_fem(getfem.Fem("FEM_PK(3,%i)" %(degree)));
#set the regions for the boundary conditions
fixed_faces = mesh.faces_from_pid(groups["FixedBoundary"]);
force_faces = mesh.faces_from_pid(groups["ForceBoundary"]);
print fixed_faces
print fixed_faces.shape
mesh.set_region(DIRICHLET_BOUNDARY,fixed_faces);
mesh.set_region(NEUMANN_BOUNDARY,force_faces);
#material parameters
c44 = 79.62e+9; c12 = 63.93e+9;
lam = c12;
mu = c44;
#calculate the source term
source_term = [0.0,0.0,-10.0];
b0 = getfem.MdBrick("isotropic_linearized_elasticity",mim,mfu);
b0.set_param("lambda",lam);
b0.set_param("mu",mu);
#set the source term
b1 = getfem.MdBrick("source term",b0,1);
b1.set_param("source_term",source_term);
b2 = getfem.MdBrick("dirichlet",b1,2,mfu,"penalized");
b2.set_param("eps",1.e-12)
mds = getfem.MdState(b2);
print "running solver ..."
b2.solve(mds,"noisy","max_res",1.e-10,"lsolver","cg/ildlt");
print "solve done!"
print "get displacement field"
U = mds.state()[0:mfu.nbdof()];
ux = U[0::3];
uy = U[1::3];
uz = U[2::3];
print "ux min: %e, max: %e" %(min(ux),max(ux))
print "uy min: %e, max: %e" %(min(uy),max(uy))
print "uz min: %e, max: %e" %(min(uz),max(uz))
print mfu.nbdof()
#vm = b0.von_mises(mds,mfe);
#save some data
print "export mesh to DX file"
mesh.export_to_dx("mesh.dx");
#slice the data and export it
print "export slice to DX file"
sl = getfem.Slice(("boundary",),mfu,degree);
slui = getfem.Slice(("boundary",),mfe,degree);
slui.export_to_dx("displacement.dx","ascii",mfe,uz,"Displacement Field uz");
slui.export_to_dx("displacement.dx","ascii","append",mfe,ux,"Displacement Field
ux");
slui.export_to_dx("displacement.dx","ascii","append",mfe,uy,"Displacement Field
uy");
sl.export_to_dx("displacement.dx","ascii","append",mfu,U,"Displacement Field
u");
------------------------------------------------------------------------
_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users