/* *****************************************************
							      *
								*   *
							  *       *
					        *       *
                          *    2  *
            			*       * 
            			  *   *	
* * * * * * * * *  *        *
*		1		   *
*                  *
* * * * * * * * *  *      *
					    *   *
					  *		  *
					    *	3	*
					      *		  *
					        *	*
					          *		 * * * * * *
					    			 *	  4	   *
					    			 * * * * * *
					    	* * *
					    	*	*
					    	* 5 *
					    	*	*
					    	*	*
					    	* * *
***************************************************** */ 
border down1(t=0,1){x=t;y=-0.5;label=1;};
border lowerRight1(t=-0.5,0){x=1;y=t;label=2;};
border upperRight1(t=0,0.5){x=1;y=t;label=3;};
border up1(t=1,0){x=t;y=0.5;label=4;};
border left1(t=0.5,-0.5){x=0.0;y=t;label=5;};
int nm=32;

mesh Th1=buildmesh(down1(nm)+lowerRight1(nm)+upperRight1(nm)+up1(nm)+left1(nm)); // first domain
plot(Th1,wait=1);

ofstream merror("materrors.m");

merror << "error=[ ";

fespace X1h(Th1,P2); 
X1h u11,u12,u11old,u12old,v11,v12,k11,k12;
fespace Q1h(Th1,P1); //2
Q1h p1,q1;

//------------------------------------------------------------------
border down2(t=1.5,2.5){x=t;y=t-1.5;label=1;};
border right2(t=2.5,2){x=t;y=-t+3.5;label=2;};
border up2(t=2,1){x=t;y=t-0.5;label=3;};
border left2(t=1,1.5){x=t;y=-t+1.5;label=4;};

mesh Th2=buildmesh(down2(nm)+right2(nm)+up2(nm)+left2(nm)); // second upper-right domain
plot(Th2,wait=1);

fespace X2h(Th2,P2); 
X2h u21,u22,u21old,u22old,v21,v22,k21,k22;
fespace Q2h(Th2,P1); //2
Q2h p2,q2,p2old;

//------------------------------------------------------------------

border down3(t=1,2){x=t;y=-t+0.5;label=1;};
border downright3(t=2,2.25){x=t;y=t-3.5;label=2;};
border upright3(t=2.25,2.5){x=t;y=t-3.5;label=3;};
border up3(t=2.5,1.5){x=t;y=-t+1.5;label=4;};
border left3(t=1.5,1){x=t;y=t-1.5;label=5;};

mesh Th3=buildmesh(up3(nm)+downright3(nm)+upright3(nm)+down3(nm)+left3(nm)); //third lower-right domain
plot(Th3,wait=1);


fespace X3h(Th3,P2); 
X3h u31,u32,u31old,u32old,v31,v32,k31,k32;
fespace Q3h(Th3,P1); //2
Q3h p3,q3,p3old;

//------------------------------------------------------------------

border down4(t=2.5,3.5){x=t;y=-1.5;label=1;};
border right4(t=-1.5,-1){x=3.5;y=t;label=2;};
border up4(t=3.5,2.5){x=t;y=-1;label=3;};
border left4(t=-1,-1.5){x=2.5;y=t;label=4;};

mesh Th4=buildmesh(up4(nm)+right4(nm)+down4(nm)+left4(nm)); //right child domain
plot(Th4,wait=1);

fespace X4h(Th4,P2); 
X4h u41,u42,u41old,u42old,v41,v42,k41,k42;
fespace Q4h(Th4,P1); //2
Q4h p4,q4,p4old;

//------------------------------------------------------------------
border down5(t=2,2.5){x=t;y=-2.5;label=1;};
border right5(t=-2.5,-1.5){x=2.5;y=t;label=2;};
border up5(t=2.5,2){x=t;y=-1.5;label=3;};
border left5(t=-1.5,-2.5){x=2;y=t;label=4;};

mesh Th5=buildmesh(up5(nm)+right5(nm)+down5(nm)+left5(nm)); //down child domain
plot(Th5,wait=1);

fespace X5h(Th5,P2); 
X5h u51,u52,u51old,u52old,v51,v52,k51,k52;
fespace Q5h(Th5,P1); //2
Q5h p5,q5,p5old;

//------------------------------------------------------------------
u11=0.; u12=0.; u11old = 0.; u12old=0.; 
u21=0.; u22=0.; u21old = 0.; u22old=0.; 
u31=0.; u32=0.; u31old = 0.; u32old=0.; 
u41=0.; u42=0.; u41old = 0.; u42old=0.;
u51=0.; u52=0.; u51old = 0.; u52old=0.;
p1=0.; p2=0.;p2old=0.; p3=0.;p3old=0.;p4=0.;p4old=0.; p5=0.;p5old=0.;
q1=0.; q2=0.;q3=0.;q4=0.;q5=0.;
k11=u11;k12=u12;
k21=u21;k22=u22;
k31=u31;k32=u32;
k41=u41;k42=u42;
k51=u51;k52=u52;

func g=(0.4-y)*(y+0.4)*(y<0.4)*(y>-0.4); 
real  nu=1.e-1;
real toll=1.e-4;
real gamma=.3;
real P11= 0.35;
real errL2=0.,errL2max=0.;

problem Stokes1 ([u11,u12,p1],[v11,v12,q1],solver=UMFPACK) =
    int2d(Th1)(//dti*u1*v1 + dti*u2*v2 +
            nu * ( dx(u11)*dx(v11) + dy(u11)*dy(v11)
            +     dx(u12)*dx(v12) + dy(u12)*dy(v12) )) 

  - int2d(Th1)(p1*dx(v11) + p1*dy(v12)) 
  - int2d(Th1)(dx(u11)*q1 + dy(u12)*q1)

      + int2d(Th1)( (k11*dx(u11)+k12*dy(u11))*v11 
            +   (k11*dx(u12)+k12*dy(u12))*v12 )

  + int1d(Th1,3)(((p2-nu*dx(u21))*N.x-nu*dy(u21)*N.y)*v11+
             ((p2-nu*dy(u22))*N.y-nu*dx(u22)*N.x)*v12)
  
  + int1d(Th1,2)(((p3-nu*dx(u31))*N.x-nu*dy(u31)*N.y)*v11+
             ((p3-nu*dy(u32))*N.y-nu*dx(u32)*N.x)*v12)
  
  + on(5,u11=g,u12=0.)
  + on(1,4,u11=0.,u12=0.); 


problem Stokes2 ([u21,u22,p2],[v21,v22,q2],solver=UMFPACK) =
    int2d(Th2)(//dti*u1*v1 + dti*u2*v2 +
            nu * ( dx(u21)*dx(v21) + dy(u21)*dy(v21)
            +        dx(u22)*dx(v22) + dy(u22)*dy(v22) )) 
  - int2d(Th2)(p2*dx(v21) + p2*dy(v22)) 
  - int2d(Th2)(dx(u21)*q2 + dy(u22)*q2)  
    + int2d(Th2)( (k21*dx(u21)+k22*dy(u21))*v21 
            +   (k21*dx(u22)+k22*dy(u22))*v22 )
  + on(4,u21=-u11*N.x*3.0/4.0,u22=-u11*N.y*3.0/4.0)
  + on(1,3,u21=0.,u22=0.); 

problem Stokes3 ([u31,u32,p3],[v31,v32,q3],solver=UMFPACK) =
    int2d(Th3)(//dti*u1*v1 + dti*u2*v2 +
            nu * ( dx(u31)*dx(v31) + dy(u31)*dy(v31)
            +        dx(u32)*dx(v32) + dy(u32)*dy(v32) )) 
  - int2d(Th3)(p3*dx(v31) + p3*dy(v32)) 
  - int2d(Th3)(dx(u31)*q3 + dy(u32)*q3)
    + int2d(Th3)( (k31*dx(u31)+k32*dy(u31))*v31 
            +   (k31*dx(u32)+k32*dy(u32))*v32 )
  + int1d(Th3,3)(((p4-nu*dx(u41))*N.x-nu*dy(u41)*N.y)*v31+
             ((p4-nu*dy(u42))*N.y-nu*dx(u42)*N.x)*v32) 
  + int1d(Th3,2)(((p5-nu*dx(u51))*N.x-nu*dy(u51)*N.y)*v31+
             ((p5-nu*dy(u52))*N.y-nu*dx(u52)*N.x)*v32)
  + on(5,u31=-u11*N.x*1.0/4.0,u32=-u11*N.y*1.0/4.0)
  + on(1,4,u31=0.,u32=0.); 

problem Stokes4 ([u41,u42,p4],[v41,v42,q4],solver=UMFPACK) =
    int2d(Th4)(//dti*u1*v1 + dti*u2*v2 +
            nu * ( dx(u41)*dx(v41) + dy(u41)*dy(v41)
            +        dx(u42)*dx(v42) + dy(u42)*dy(v42) )) 
  - int2d(Th4)(p4*dx(v41) + p4*dy(v42)) 
  - int2d(Th4)(dx(u41)*q4 + dy(u42)*q4) 
    + int2d(Th4)( (k41*dx(u41)+k42*dy(u41))*v41 
            +   (k41*dx(u42)+k42*dy(u42))*v42 )
  + on(4,u41=-u31*N.x*0.5,u42=-u31*N.y*0.5)
  + on(1,3,u41=0.,u42=0.); 

problem Stokes5 ([u51,u52,p5],[v51,v52,q5],solver=UMFPACK) =
    int2d(Th5)(//dti*u1*v1 + dti*u2*v2 +
            nu * ( dx(u51)*dx(v51) + dy(u51)*dy(v51)
            +        dx(u52)*dx(v52) + dy(u52)*dy(v52) )) 
  - int2d(Th5)(p5*dx(v51) + p5*dy(v52)) 
  - int2d(Th5)(dx(u51)*q5 + dy(u52)*q5)
   + int2d(Th5)( (k51*dx(u51)+k52*dy(u51))*v51 
   +   (k51*dx(u52)+k52*dy(u52))*v52 )
  + on(3,u51=-u31*N.x*0.5,u52=-u31*N.y*0.5)
  + on(4,2,u51=0.,u52=0.); 

real[int]xx(50); real[int]yy(50);
real thr=0.001;
real Error=10^5;
int nmax = 200;
 int j=0;
while(Error>thr && j<nmax )
{
 cout << "Iteration " << j << endl;
 u11old=u11; u12old=u12;
 u21old=u21; u22old=u22; p2old=p2;
 u31old=u31; u32old=u32; p3old=p3;
 u41old=u41; u42old=u42; p4old=p4;
 u51old=u51; u52old=u52; p5old=p5;

 Stokes3;
 plot([u31,u32],value=true,nbiso=10,wait=false);//,ps="es10_vel_a.eps");
 plot(p3,fill=1,value=true,wait=false);//,ps="es10_pre_a.eps");

 Stokes4;
plot([u41,u42],value=true,nbiso=10,wait=false);//,ps="es10_vel_a.eps");
plot(p4,fill=1,value=true,wait=false);//,ps="es10_pre_a.eps");

 Stokes5;
plot([u51,u52],value=true,nbiso=10,wait=false);//,ps="es10_vel_a.eps");
plot(p5,fill=1,value=true,wait=false);//,ps="es10_pre_a.eps");

 Stokes1;
plot([u11,u12],value=true,nbiso=10,wait=false);//,ps="es10_vel_a.eps");
plot(p1,fill=1,value=true,wait=false);//,ps="es10_pre_a.eps");

 Stokes2;
plot([u21,u22],value=true,nbiso=10,wait=false);//,ps="es10_vel_a.eps");
plot(p2,fill=1,value=true,wait=false);//,ps="es10_pre_a.eps");

real Error1 = int2d(Th1)((u11-u11old)^2+(u12-u12old)^2);
 Error1 = sqrt(Error1);
cout << "Error1 " << Error1 << endl; 
 real Error2 = int2d(Th2)((u21-u21old)^2+(u22-u22old)^2);
 Error2 = sqrt(Error2);
cout << "Error2 " << Error2 << endl; 
 real Error3 = int2d(Th3)((u31-u31old)^2+(u32-u32old)^2);
 Error3 = sqrt(Error3);
cout << "Error3 " << Error3 << endl; 
real Error4 = int2d(Th4)((u41-u41old)^2+(u42-u42old)^2);
 Error3 = sqrt(Error4);
cout << "Error4 " << Error4 << endl; 
real Error5 = int2d(Th5)((u51-u51old)^2+(u52-u52old)^2);
 Error5 = sqrt(Error5);
cout << "Error5 " << Error5 << endl; 
 Error = int2d(Th1)((u11-u11old)^2+(u12-u12old)^2)+int2d(Th2)((u21-u21old)^2+(u22-u22old)^2) + int2d(Th3)((u31-u31old)^2+(u32-u32old)^2)+ 
              int2d(Th4)((u41-u41old)^2+(u42-u42old)^2)+ int2d(Th5)((u51-u51old)^2+(u52-u52old)^2);
 Error = sqrt(Error);
 cout << "Convergence " << Error << endl; 

u11 = gamma*u11 + (1-gamma)*u11old;
u12 = gamma*u12 + (1-gamma)*u12old;

u21 = gamma*u21 + (1-gamma)*u21old;
u22 = gamma*u22 + (1-gamma)*u22old;
p2 = gamma*p2 + (1-gamma)*p2old;

u31 = gamma*u31 + (1-gamma)*u31old;
u32 = gamma*u32 + (1-gamma)*u32old;
p3 = gamma*p3 + (1-gamma)*p3old;

u41 = gamma*u41 + (1-gamma)*u41old;
u42 = gamma*u42 + (1-gamma)*u42old;
p4 = gamma*p4 + (1-gamma)*p4old;

u51 = gamma*u51 + (1-gamma)*u51old;
u52 = gamma*u52 + (1-gamma)*u52old;
p5 = gamma*p5 + (1-gamma)*p5old;

k11=u11;k12=u12;
k21=u21;k22=u22;
k31=u31;k32=u32;
k41=u41;k42=u42;
k51=u51;k52=u52;
 xx[j]=Error;
 yy[j]=j;
merror << xx[j] << "; \n";


j++;
}

merror << "]; \n";
plot([xx,yy],value=true,nbiso=10,wait=true);
plot(p1,[u11,u12],value=true,wait=false);
plot(p2,[u21,u22],value=true,wait=false);
plot(p3,[u31,u32],value=true,wait=false);
plot(p4,[u41,u42],value=true,wait=false);
plot(p5,[u51,u52],value=true,wait=false);
plot([u11,u12],value=true,nbiso=10,wait=1);//,ps="es10_vel_a.eps");
plot(p1,fill=1,value=1,wait=1);//,ps="es10_pre_a.eps");
errL2 = int2d(Th1)((u11-g)*(u11-g)+u12*u12);
errL2 = sqrt(errL2);
cout << "Error vs Exact " << errL2 << endl;
plot([u21,u22],value=true,nbiso=10,wait=1);//,ps="es10_vel_a.eps");
plot(p2,fill=1,value=1,wait=1);//,ps="es10_pre_a.eps");
errL2 = int2d(Th2)((u21-g)*(u21-g)+u22*u22);
errL2 = sqrt(errL2);
cout << "Error vs Exact " << errL2 << endl;
plot([u31,u32],value=true,nbiso=10,wait=1);//,ps="es10_vel_a.eps");
plot(p3,fill=1,value=1,wait=1);//,ps="es10_pre_a.eps");
errL2 = int2d(Th3)((u31-g)*(u31-g)+u32*u32);
errL2 = sqrt(errL2);
cout << "Error vs Exact " << errL2 << endl;
plot([u41,u42],value=true,nbiso=10,wait=1);//,ps="es10_vel_a.eps");
plot(p4,fill=1,value=1,wait=1);//,ps="es10_pre_a.eps");
errL2 = int2d(Th4)((u41-g)*(u41-g)+u42*u42);
errL2 = sqrt(errL2);
cout << "Error vs Exact " << errL2 << endl;
plot([u51,u52],value=true,nbiso=10,wait=1);//,ps="es10_vel_a.eps");
plot(p5,fill=1,value=1,wait=1);//,ps="es10_pre_a.eps");
errL2 = int2d(Th5)((u51-g)*(u51-g)+u52*u52);
errL2 = sqrt(errL2);
cout << "Error vs Exact " << errL2 << endl;