load "MetricPk"

int kDeg=4;
int nIter=5;
mesh Th=square(20,20,[2*x-1,2*y-1]);
real c=5;
func f = x^2*y + y^3 + tanh(c*(-2*x + sin(5*y)));
real ratio = 1.3;
//real massRef=40*Th.nt; //40 : constante pour passer de la masse au nb de triangles (?).
real massRef=Th.nt/2; //The number of triangles is approx twice the mass

fespace Metric(Th,[P1,P1,P1]);
Metric [m11,m12,m22];
for(int i=0; i<nIter; i++){
	plot(Th,wait=true);
	[m11,m12,m22]=[0,0,0];
	m11[]=MetricPk(Th,f,kDeg=kDeg,rDeg=1,mass=massRef,MetricType=1,TriangulationType=1);
	Th = adaptmesh(Th,m11,m12,m22,IsMetric=true);
	massRef=massRef*ratio;
}
plot(Th,wait=true,ps="o/lastMesh.ps");
