testFE.edp

load "BernadiRaugel" 
// 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,[10*(x+y/3),10*(y-x/3)]); 

real x1=0.7,y1=0.9, h=1e-7; 
int it1=Th(x1,y1).nuTriangle; 

fespace Vh(Th,P2BR); 

Vh [a1,a2],[b1,b2],[c1,c2]; 


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,a2], wait=1); 
  [b1,b2]=[a1,a2];// 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 


  cout << " dx(a1)(x1,y1) = " << dx(a1)(x1,y1) << " == " << DD(a1,h,0) << endl; 
  assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-5); 
  assert( abs(dx(a2)(x1,y1)-DD(a2,h,0) ) < 1e-5); 
  assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-5); 
  assert( abs(dy(a2)(x1,y1)-DD(a2,0,h) ) < 1e-5); 


}