testFE-P3.edp

load "Element_P3" 
macro DD(f,hx,hy) ( (f(x1+hx,y1+hy)-f(x1-hx,y1-hy))/(2*(hx+hy))) // 
macro DD2(f,hx,hy) ( (-2*f(x1,y1)+f(x1+hx,y1+hy)+f(x1-hx,y1-hy))/(square(hx+hy))) // 
macro dn(f) ( N.x*dx(f)+N.y*dy(f)) // 
macro dnn(f) ( N.x*f#2+N.y*f#3) // 
mesh Th=square(1,1,[10*(x+y/3),10*(y-x/3)]); 

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

fespace Vh(Th,P3); 
fespace Eh(Th,P0edge); 

Eh  edges; 

Vh a1,b1,c1; 

varf vFlux([a],[e]) = intalledges(Th)( dn(a1)*e*(jump(real(nuTriangle))<= 0)); 
varf vMean([a],[e]) = intalledges(Th)( (a1)*e*(jump(real(nuTriangle))<= 0)/lenEdge); 


for (int i=0;i<Vh.ndofK;++i) 
	cout << i << " -> " << Vh(0,i) << endl; 
for (int i=0;i<Vh.ndofK;++i) 
{ 
  cout << " ***  node " << i << " of Traingle " << it1 << endl; 
  a1[]=0; 
  int j=Vh(it1,i); 
  a1[][j]=1; 
  edges[]=vFlux(0,Eh); 
  cout << "Flux  edges = " << edges[] << endl; 
  edges[]=vMean(0,Eh); 
  cout << " Mean   edges = " << edges[] << endl; 

  plot(a1, wait=1,cmm="w_"+i); 
  b1=a1; 

  plot(a1,b1,cmm="w"+i, wait=1); 

  c1[] = a1[] - b1[]; 

  cout << " ---------" << i << " " << c1[].max << " " << c1[].min << endl; 
  cout << " a = " << a1[] <<endl; 
  cout << " b = " << b1[] <<endl; 

  assert(c1[].max < 1e-5 && c1[].min > -1e-9); 

  cout << " dx(a1)(x1,y1) = " << dx(a1)(x1,y1) << " == " << DD(a1,h,0) << endl; 
  cout << " dy(a1)(x1,y1) = " << dy(a1)(x1,y1) << " == " << DD(a1,0,h)  << endl; 
  cout << " dxx(a1)(x1,y1) = " << dxx(a1)(x1,y1) << " == " << DD2(a1,h,0) << endl; 
  cout << " dyy(a1)(x1,y1) = " << dyy(a1)(x1,y1) << " == " << DD2(a1,0,h)  << endl; 

  assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-4); 
  assert( abs(dxx(a1)(x1,y1)-DD2(a1,h,0) ) < 1e-4); 
  assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-4); 
  assert( abs(dyy(a1)(x1,y1)-DD2(a1,0,h) ) < 1e-4); 



}