load "Element_PkEdge"
// a macro the compute numerical derivative
macro DD(f,hx,hy) ( (f(x1+hx,y1+hy)-f(x1-hx,y1-hy))/(2*(hx+hy))) //
mesh Th=square(1,1,flags=0,[10*(x+y/3),10*(y-x/3)]);
Th=square(1,1,flags=0   );
plot(Th,wait=1);
//savemesh(Th,"Th.msh");
real x1=0.4,y1=0.4;//10*(x+y/3),y1=10*(y-x/3), 
real h=1e-5 ;
int it1=Th(x1,y1).nuTriangle; 

// FFCS: store reference value
real regtest;

 macro Check(PkEdge)
{ 
  cout << " Test Finite Element " <<endl << endl;
  fespace Vh(Th,PkEdge);

  Vh a1,b1,c1;
  c1=x+y;
  real[int] viso(100);
  real v0= -100;
  real dv = 200./(viso.n+1);
  for(int i=0;i<viso.n;++i)
     viso[i]=v0+i*dv;
  plot(c1,fill=1,wait=1,viso=viso);
  for (int i=0;i<Vh.ndofK;++i)
    cout << i << " " << Vh(0,i) << endl;
  for (int i=0;i<Vh.ndofK;++i)
    {
      a1[]=0;	
      int j=Vh(it1,i);
      a1[][j]=1;
      plot(a1, wait=1); 
      b1=a1;/* do the interpolation */
      
      
      c1[] = a1[] - b1[];
      cout << " ---------" << i << " " << c1[].max << " " << c1[].min << endl;	
      cout << " a = " << a1[] <<endl;
      cout << " b = " << b1[] <<endl;
      assert(c1[].max < 1e-9 && c1[].min > -1e-9); /* check if the interpolation is correct */
      
      /*  check the derivative and numerical derivative */
      /* no derivative in this element.  */
      
      
	cout << " dx(a1)(x1,y1) =  " << dx(a1)(x1,y1) << "  == " 
	   << DD(a1,h,0) << " DD " << " /  " 
	   << (a1(0.001,0.001) -  a1(0.999,0.999))/2 <<endl; 
	assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-4);
	assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-4);
	/*sassert( abs(dy(a2)(x1,y1)-DD(a2,0,h) ) < 1e-5);*/
     
      
    } 
  regtest=a1'*a1;
}  //EOM

Check(P1edge)
Check(P2edge)
Check(P3edge)
Check(P4edge)
Check(P5edge)
