// macro the get the current mesh size 
// parameter 
//  in:  Th the mesh
//       Vh  P1 fespace on Th
//  out : 
// h: the Vh finite element finite set to the current mesh size 
macro  MeshSizecomputation(Th,Vh,h){
   /* Th mesh
	 Vh  P1 finite element space 
	 h   the P1 mesh size value */
	real[int]  count(Th.nv);
	/* mesh size  (lenEdge = $\int_e 1$)  */
	varf vmeshsizen(u,v)=intalledges(Th,qfnbpE=1)(v); 
	/* number of edge / par vertex */ 
	varf vedgecount(u,v)=intalledges(Th,qfnbpE=1)(v/lenEdge); 
   /*
	  computation of the mesh size
	  ----------------------------- */ 
	count=vedgecount(0,Vh);
	h[]=0.;
	h[]=vmeshsizen(0,Vh);
	cout << " count min = "<< count.min << " " << count.max << endl;
	h[]=h[]./count;
    cout << " -- bound meshsize = " <<h[].min << " " << h[].max << endl;
} // end of macro MeshSizecomputation

// macro to remesh according the residual indicator 
// in: 
//     Th the mesh
//     Ph  P0 fespace on Th
//     Vh  P1 fespace on Th
//     vindicator the varf of to evaluate the indicator to ^2
//     coef on etameam ..
// ------
macro  ReMeshIndicator(Th,Ph,Vh,vindicator,coef){
	Vh h=0;
	/*evalutate the mesh size  */
	MeshSizecomputation(Th,Vh,h); 
	Ph etak;  
	etak[]=vindicator(0,Ph);
	cout << " global  Eta  : " << sqrt(etak[].sum) << "  ......... " <<  Th.nv<< endl;
	etak[]=sqrt(etak[]); 
        plot(etak,cmm="arei-etak",fill=1,value=1);
	real etastar= coef*(etak[].sum/etak[].n);
	cout << " etastar = " << etastar << " sum=" << etak[].sum << " " << endl;

	/* here etaK is discontinous 
	   we use the P1 L2 projection with mass lumping . */
	
	Vh fn,sigma;
	varf veta(unused,v)=int2d(Th)(etak*v);
	varf vun(unused,v)=int2d(Th)(1*v);
	fn[]  = veta(0,Vh);
	sigma[]= vun(0,Vh);
	fn[]= fn[]./ sigma[];
	fn =  max(min(fn/etastar,3.),0.3333) ; 
	
	/* new mesh size */ 
	h = h / fn ; 
	/* plot(h,wait=1); */
	Th=adaptmesh(Th,IsMetric=1,h,splitpbedge=1,nbvx=10000);
} // EOM

// our L space problem. 

// mesh definition 
border ba(t=0,1.0){x=t;   y=0;  label=1;}
border bb(t=0,1.0){x=0;   y=-1+t;  label=2;}
border bc(t=0,3*pi/2){x=cos(t); y=sin(t); label=3;}
int n=2;
mesh Th = buildmesh (ba(4*n) + bb(4*n) + bc(30*n));

// FE space definition ---
fespace Vh(Th,P1); // for the mesh size 
fespace Ph(Th,P0); // for the indicator   

real hinit=0.2; // 
Vh   h=hinit; // the FE fonction  for the mesh size 
// to build a mesh with a given mesh size  : meshsize
Th=adaptmesh(Th,h,IsMetric=1,splitpbedge=1,nbvx=10000);
plot(Th,wait=1);

Vh u,v; 

func f=0;
func exactf=(y<-1.e-10? (x^2+y^2)^(1./3.)*sin( 2.*(2*pi+atan2(y,x))/3. ):
                        (x^2+y^2)^(1./3.)*sin( 2.*(     atan2(y,x))/3. ));

Vh truerror, exactu=exactf;

problem Poisson(u,v,solver=UMFPACK) =
    int2d(Th,qforder=5)( dx(u)*dx(v) + dy(u)*dy(v)) 
  + on(1,2,u=exactf) + on(3,u=exactf);

 varf indicator2(unused,chiK) = 
     intalledges(Th)(chiK*lenEdge*square(jump(N.x*dx(u)+N.y*dy(u))))
    +int2d(Th)(chiK*square(hTriangle*(f+dxx(u)+dyy(u))) );




for (int i=0;i< 10;i++)
{
	u=u;
	Poisson;
	plot(Th,u,wait=1);
        cout << u[].min << " " << u[].max << endl; 
        plot(u,cmm="solution",wait=1);
        truerror = u - exactu;
        real normerror = int2d(Th)( truerror^2 );
        real normsoln = int2d(Th)( u^2 );
        plot(truerror,cmm="true error",wait=true,value=true,fill=true);
        cout << " true rel error=" << normerror/normsoln << endl;

	real cc=0.7;
	if(i>5) cc=1;
        if(i<9)
        ReMeshIndicator(Th,Ph,Vh,indicator2,cc);
	plot(u,Th,cmm="remeshed solution",wait=1,value=1);
}