testFE-PkEdge.edp

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,[10*(x+y/3),10*(y-x/3)]); 
// Th=square(1,1); 
//savemesh(Th,"Th.msh"); 
real x1=0.7,y1=0.9, h=1e-7; 
int it1=Th(x1,y1).nuTriangle; 


 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) << 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);       */ 

    } 

}  //EOM 

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