load "Element_P2pnc"
load "Element_P3dc"

 mesh Th=square(10,10);
 fespace Vh(Th,P2pnc);     // P2+ nc  in P FE space
 cout << " Vh.ndof = "<< Vh.ndof << endl;
 Vh uh,vh;              // unkown and test function.
 func f=1;                 //  right hand side function
 func g=0;                 //  boundary condition function

 varf vLap(uh,vh) =                    //  Definion of  the problem
    int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) //  bilinear form
   - int2d(Th)( f*vh )                          //  linear form
   + on(1,2,3,4,uh=g) ;                      //  boundary condition form
  ;
  matrix A = vLap(Vh,Vh,solver="SPARSESOLVER");
  real[int] b= vLap(0,Vh);
  //cout << A << endl;
  uh[]=A^-1*b ;// solve the problem plot(uh); // to see the result
  plot(uh,ps="o/LaplaceP2pnc.eps",value=true);
