LaplaceP3.edp

load "Element_P3" 

 mesh Th=square(10,10); 
 fespace Vh(Th,P3);     // P1 FE space 
 Vh uh,vh;              // unkown and test function. 
 func f=1;                 //  right hand side function 
 func g=0;                 //  boundary condition function 

 problem laplace(uh,vh) =                    //  Definion of  the problem 
    int2d(Th,qft=qf7pT)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) //  bilinear form 
  - int2d(Th,qft=qf7pT)( f*vh )                          //  linear form 
  + on(1,2,3,4,uh=g) ;                      //  boundary condition form 

  laplace; // solve the problem plot(uh); // to see the result 
  plot(uh,ps="LaplaceP3.eps",value=true);