getfem-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] problem with dirichlet boundary conditions


From: julien pommier
Subject: Re: [Getfem-users] problem with dirichlet boundary conditions
Date: Tue, 23 Jan 2007 14:38:44 +0100
User-agent: Thunderbird 1.5.0.8 (X11/20061117)

Eugen Wintersberger wrote:
Thanks for your fast response
 The results are getting better using "augmented" or "eliminated" for
the Dirichlet brick. However, the convergence of the solver (cg - in
this case) is extremely bad. I used superlu instead for this simple
example. But for larger problems (the real ones) a direct solver would
not be appropriate. Is it somehow possible to change the penalization
weight from the python interface?

Well, no it was not possible. Now it is (I just commited the change). You can now use
brick.penalization_epsilon(1e-12)

or whatever value of epsilon (default is 1e-9)

--
Julien
Eugen

On Tue, 2007-01-23 at 10:38 +0100, julien pommier wrote:
Hi Eugen,
Your Dirichlet condition is imposed by penalization, which may be the source of your problem. Try with 'eliminated' or 'augmented' , I think the error that you observe should improve

(especially since we use a 1e9 weight for the penalization and you have very large coefficients which are the same order of magnitude)

--
Julien

Eugen Wintersberger wrote:
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



--
Julien Pommier, bureau 111
GMM, INSA Toulouse, tél:05 61 55 93 42




reply via email to

[Prev in Thread] Current Thread [Next in Thread]