testfe.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

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


}