load "MetricPk"
load "Element_P3";
load "Element_P4";

int nIter=4;
mesh Th=square(20,20,[x,y]);
//fespace Vh(Th,P1); int kDeg=1;  //degr� des FE utilis�s
//fespace Vh(Th,P2); int kDeg=2;  //degr� des FE utilis�s
fespace Vh(Th,P3); int kDeg=3;  //degr� des FE utilis�s
int rdeg=1;	 //on approche le gradient ou des d�riv�es d'ordre sup�rieur

real ratio = 2.;//number of triangles is multiplied by ratio at each step
real rmax = 1; 	//max triangle size desired.

Vh u,v;
fespace Metric(Th,[P1,P1,P1]);
Metric [m11,m12,m22];
////pour de l'adaptation isotrope
//problem EpsDelta(u,v) =
//		int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))+
//		on(1,u=1)+on(2,3,4, u=0);
//avec une forte anisotropie, et des singularit�s aux coins.
real eps=0.01;
problem EpsDelta(u,v) =
		int2d(Th)(eps*dx(u)*dx(v)+eps*dy(u)*dy(v)+(dx(u)+dy(u))*v)+
		on(1,u=atan(10*x*(1-x)))+on(2,3,4, u=0);

//real massRef = 1000*Th.nt;
real massRef = Th.nt/2;
EpsDelta;
plot(u, wait=1);
for(int i=0; i<nIter; i++){
	plot(Th, wait=1); //on affiche avant de tout changer
	plot(u, wait=1);
	[m11,m12,m22]=[0,0,0];	//absolument crucial...
	m11[]=MetricPk(Th,u,kDeg=kDeg,rDeg=rdeg, mass=massRef);
	Th = adaptmesh(Th,m11,m12,m22,IsMetric=true);
	cout << "expected number of triangles : " << massRef/20 << "; obtained " << Th.nt << "\n";
	EpsDelta;
	massRef=ratio*massRef;
}

plot(Th, wait=1, ps="o/triangulationEpsDelta.ps");
