// file: laplace.edp 
border C(t=0,2*pi){x=cos(t); y=sin(t); label=1;}; // defines the boundary 
mesh Th = buildmesh (C(100)); 
   //    the triangulated domain Th is on the left side of its boundary 
plot(Th, wait=1, ps = "laplacemeshplot.ps");
fespace Vh(Th,P1); // defines FE space as C^0 piecewise P1 functions 
Vh u,v;  //  the finite element space defined over Th is called Vh
real L2error, graderror;
func f= 8*(1-2*x^2-2*y^2);   //    definition of a called f function 
func uexact = (1- x^2-y^2)^2;
func dxuexact= -4*x*(1- x^2-y^2);
func dyuexact= -4*y*(1- x^2-y^2);
solve Poisson(u,v,solver=LU) =               //    defines the PDE 
     int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))      //    bilinear form 
    - int2d(Th)( f*v)                         //    right hand side 
    + on(1,u=0) ;                  //  Dirichlet boundary condition 
L2error= sqrt(int2d(Th)((u-uexact)^2));
cout << " L2error  = "<<  L2error <<endl;
graderror = sqrt(int2d(Th)((dx(u)-dxuexact)^2 + (dy(u)-dyuexact)^2)); 
cout << " graderror  = "<<  graderror <<endl;