
load "Element_P3dc"
load "Element_P2pnc"
func real cc(real aa) {real a=aa; if(abs(a)<1e-10) a=0; return a;} 
int[int] ne1=[1,2,0];
int[int] ne2=[2,0,1];

real[int] xx=[0,1,0,0],yy=[0,0,1,0],zz=[0,0,0,1];
 
int[int] ll=[2,0,0,1];
mesh Th=square(1,1,flags=2,label=ll);

Th = trunc(Th,x<0.5,label=0);


//Th = movemesh(Th,[x+0.1*y,y-0.2*x]);
//Th = movemesh(Th,[x*2,y*2]);

mesh Thg = trunc(Th,1,split=2,label=-1);
plot(Th, wait=1);

int it0=0; 

fespace Lh(Th,P1);
Lh[int] l(3);
l[0][][0]=1;
l[1][][1]=1;
l[2][][2]=1;


fespace Wh(Th,P3dc);
fespace Vh(Th,P2pnc);
Vh u,v;
for (int j=0;j<Vh.ndofK;++j)
{
	int i = Vh(0,j);
	cout << j << " -> " <<  i << " : " ;
    u=0;
	u[][i]=1;
	for(int k=0; k<3;k++)
	{
		int i1= (k+1)%3, i2=(k+2)%3;
		cout << cc(int1d(Th,k)(u*l[i1]/lenEdge)) << " " << cc(int1d(Th,k)(u*l[i2]/lenEdge)) << " "; 
	}
	cout << cc(int2d(Th)( u/area)) << endl; 
	//plot(u,wait=1);
	v=u;
	for (int j=0;j<Vh.ndofK;++j)
	cout << cc(v[][Vh(0,j)]) << " " ;
	cout << endl;

}
 Th=square(10,10);
 u=x*x*y*y;
 
 solve PL2(u,v)=int2d(Th,qforder=7)(u*v)-int2d(Th,qforder=7)(x*x*y*y*v);

 plot(u,wait=1);
 real err=0,err1=0,err2=0; 
for (int i=0;i<Vh.ndof;++i)
{
    u=0;

	u[][i]=1;
	Wh u3=u;
	v=u;
	v[]-=u[];
	err += v[].linfty;
	err1 += sqrt(int2d(Th)(square(dx(u)-dx(u3))));
	err2 += sqrt(int2d(Th)(square(dy(u)-dy(u3))));
}
cout << " err=" << err << " " << err1 << " " << err2 << endl; 
assert( err+err1+err2 < 1e-10);
