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

[Eugen Wintersberger]

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?

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