load "BernardiRaugel"

// Parameters
real x1 = 0.7;
real y1 = 0.9;
real h = 1e-7;

// Mesh
mesh Th = square(1, 1, [10*(x+y/3), 10*(y-x/3)]);

// Macro
//a macro to compute numerical derivative
macro DD(f, hx, hy) ((f(x1+hx, y1+hy) - f(x1-hx, y1-hy))/(2*(hx+hy))) //

// Fespace
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;

// Test P2BR
int it1 = Th(x1, y1).nuTriangle;
for (int i = 0; i < Vh.ndofK; ++i){
	a1[] = 0;
	int j = Vh(it1, i);
	a1[][j] = 1;
	plot([a1, a2], wait=true);
	[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);
} 

