load "Element_P4"
load "Element_P3"
load "qf11to25"

mesh Th=square(10,10);
mesh3 Th3=cube(8,8,8);
verbosity=4;

for (int i=1;i<=26;i++)
{
  cout <<" int1d order= "<< i ;
  int1d(Th,qforder=i)(1.);
}

for (int i=1;i<=26;i++)
{
 cout <<" int2d order= "<< i ;
  int2d(Th,qforder=i)(1.);
}

for (int i=1;i<=14;i++)
{
 cout <<" int3d order= "<< i ;
  int3d(Th3,qforder=i)(1.);
}

verbosity=1;

 fespace Vh(Th,P4);     // 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=qf11pT)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) //  bilinear form
  - int2d(Th,qft=qf25pT)( 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="o/LaplaceP4.eps",value=true,wait=1);

  real[int,int] qq1=[ [0.5,0],
                      [0.5,1]];
  real[int,int] qq2=[[1./3,0,0],[1./3.,1,0],[1./3.,0,1]];
  real[int,int] qq3=[[1./4,0,0,0],[1./4.,1,0,0],[1./4.,0,1,0],[1./4.,0,0,1]];

  QF1 qf1(1,qq1);
  QF2 qf2(1,qq2);
  QF3 qf3(1,qq3);
  real I1 = int1d(Th,qfe=qf1)(x^2) ;
  real I1l = int1d(Th,qfe=qf1pElump)(x^2) ;

  real I2 = int2d(Th,qft=qf2)(x^2) ;
  real I2l = int2d(Th,qft=qf1pTlump)(x^2) ;

  cout << I1 << " == " << I1l << endl;
  cout << I2 << " == " << I2l << endl;
  assert( abs(I1-I1l) < 1e-10);
  assert( abs(I2-I2l) < 1e-10);

   fespace Vh3(Th3,P4);     // P1 FE space

Vh3 u3,v3;              // unkown and test function.

 solve laplace3(u3,v3) =                    //  Definion of  the problem
    int3d(Th3,qforder=9)( dx(u3)*dx(v3) + dy(u3)*dy(v3) +  dz(u3)*dz(v3) ) //  bilinear form
  - int3d(Th3,qforder=9)( v3 )                          //  linear form
  + on(1,2,3,4,5,6,u3=1) 
  ;                      //  boundary condition form

plot(u3,wait=1);
cout << " u3 max "<< u3[].max << " " << u3[].min  << endl; 