// this file is used to test the two-phase case 



// parameter //// 
int     meshsize=128;
real    dt=0.001;//1./real(meshsize)/10.;
real 	dh=0.5,dl=2.0;


real    etaf=1e-1,etap=1e0;
real    beta2=10.0,betae=1.0;
real    tau2=1e-4;
real    eps2=0.01;

real    Re=1e-1;
real    thetas=pi/4.0;
real 	eps; // used in SUPG should be update in each time setp
real 	umax=0.0;
func 	uin=umax*(dh-y)*y;
func 	psi1in=-y;
real    t;
func 	psi2in=x-umax*y*(dh-y)*t;
real 	ul2;
real    ls=0.002;
real    gamma=0.01;

real    rhobar;
real    rho12=0.001;
real    Bo=20;
real    theta=0.0;


real elastic;
real kinetic;
real mix;
real ewall;
real etot;

// Mesh 
border 	bb(t=0,dl){x=t; y=0; label=1;};
border 	br(t=0,dh){x=dl; y=t; label =2;};
border	bt(t=dl,0){x=t; y=dh; label =3;};
border 	bl(t=dh,0){x=0; y=t; label =4;};
 

mesh Th= buildmesh (bt(meshsize*dl)+bb(meshsize*dl)+ bl(meshsize*dh) + br(meshsize*dh));
savemesh(Th,"mesh.msh");
savemesh(Th,"meshe.txt");

plot(Th,ps="meshe.eps");

// ********* Define the Finite Element Space
fespace     Vh(Th,P2);
fespace     Wh(Th,P1);

Vh  u1,u2,u01,u02,ut01,ut02,us1,us2,v1,v2;
// femspace of the velocity, step n+1: u1,u2; step n: u01,u01; step half ut01,ut02; ustart: us1,us2,test:v1,v2.
// femspace of the velocity, step n+1: u1,u2; step n: u01,u01; step half ut01,ut02; ustart: us1,us2,test:v1,v2.

Wh  w11,w12,w21,w22,p,p01,p02,q,dux; 
// fespace of the pressure  step n+1, p; step n: p01; test:q.
// w11=dx(psi1),w12=dy(psi1),w21=dx(psi2),w22=dy(psi2);

Vh  psi1,psi2,psi01,psi02,psih1,psih2;

Vh  phi2,mu2,phi02,phi2h,mu2h;
Vh 	lambdae; //  which is the effective elactis in thrombus and equal to (phi_1+(1-phi)alaph_np)phi_2. it need to be update in each time step.
Vh 	eta,rho;

//***********Define the Problem ********************

// ++++++++++step 1:The problem of phi2, mu2+++++++++
 
//     energy split method 
problem step1([phi2,mu2],[phi2h,mu2h])=
int2d(Th)(1./dt*phi2*mu2h+tau2*(dx(mu2)*dx(mu2h)+dy(mu2)*dy(mu2h))
		-mu2*phi2h
		-phi2*(u01*dx(mu2h)+u02*dy(mu2h))
		+beta2/eps2*phi2*phi2h+eps2*beta2*(dx(phi2)*dx(phi2h)+dy(phi2)*dy(phi2h))
		+betae*phi2*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))*phi2h
)
+int1d(Th,1)(beta2*gamma/dt*phi2*phi2h+beta2*gamma*u01*dx(phi2)*phi2h)

+int2d(Th)(-1./dt*phi02*mu2h
//		-phi02*(u01*dx(mu2h)+u02*dy(mu2h))
		-beta2/eps2*phi02*phi2h+beta2/eps2*(phi02^3-1.5*(phi02)^2+0.5*phi02)*phi2h
		+Re*0.5*(1-rho12)*(u01^2+u02^2)*phi2h
//+beta3*(phi01+(1.0-phi01)*alphanp)*0.5
//*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))*phi2h
)   
+int1d(Th,2)(u01*phi02*mu2h)
+int1d(Th,4)(-u01*phi02*mu2h)
+int1d(Th,1)(beta2*sqrt(2.0)/12.0*cos(thetas)/2.0*(3.*(2.*phi02-1.0)^2-3.0)*phi2h
			-beta2*gamma/dt*phi02*phi2h);

// ++++++++++++step3: update the Psi++++++++++++++++


problem step3([psi1,psi2],[psih1,psih2])=
int2d(Th)(psi1*(psih1+hTriangle*eps*(u01*dx(psih1)+u02*dy(psih1)))
		+dt*(u01*dx(psi1)+u02*dy(psi1))*(psih1+hTriangle*eps*(u01*dx(psih1)+u02*dy(psih1)))  
		+psi2*(psih2+hTriangle*eps*(u01*dx(psih2)+u02*dy(psih2)))
		+dt*(u01*dx(psi2)+u02*dy(psi2))*(psih2+hTriangle*eps*(u01*dx(psih2)+u02*dy(psih2)))  
)
+int1d(Th,4)(-(1.0+eps)*(psi1*psih1)-(1.0+eps)*(psi2*psih2))
+int2d(Th)(
-psi01*(psih1+hTriangle*eps*(u01*dx(psih1)+u02*dy(psih1)))
-psi02*(psih2+hTriangle*eps*(u01*dx(psih2)+u02*dy(psih2)))
//+dt*(u01*dx(psi01)+u02*dy(psi01))*(psih1+hTriangle*eps*(u01*dx(psih1)+u02*dy(psih1)))
//+dt*(u01*dx(psi02)+u02*dy(psi02))*(psih2+hTriangle*eps*(u01*dx(psih2)+u02*dy(psih2)))  
)
+int1d(Th,4)((1.0+eps)*(psi1in*psih1)+(1.0+eps)*(psi2in*psih2));

Wh w1,w2,w1h,w2h;

// +++++++++++step4: Fluid ++++++++++++

problem NS([u1,u2],[v1,v2])=
int2d(Th)(Re/dt*rho*(u1*v1+u2*v2)+Re*rho*(u01*dx(u1)+u02*dy(u1))*v1+Re*rho*(u01*dx(u2)+u02*dy(u2))*v2
		+0.5*Re*rho*(dx(u01)+dy(u02))*(u1*v1)+0.5*Re*rho*(dx(u01)+dy(u02))*u2*v2
		 +eta*(2.0*dx(u1)*dx(v1)+(dy(u1)+dx(u2))*dy(v1)+(dx(u2)+dy(u1))*dx(v2)+2.0*dy(u2)*dy(v2))
	 	 )

+int1d(Th,1)(u1*v1/ls) 
+int2d(Th)(-Re/dt*rho*(u01*v1+u02*v2)-(2.0*p01-p02)*(dx(v1)+dy(v2))
		 	+(phi2*dx(mu2))*v1
			+(phi2*dy(mu2))*v2
			-betae*(dx(psi1)*w1*v1+dx(psi2)*w2*v1)
			-betae*(dy(psi1)*w1*v2+dy(psi2)*w2*v2)
			+Re*0.5*(1.0-rho12)*(dx(phi2)*v1+dy(phi2)*v2)*(u01^2+u02^2)
//			-Bo*rho*sin(theta)*v1
			+Bo*rho*v2	
		 )
+int1d(Th,1)(beta2*eps2*dy(phi2)*dx(phi2)*v1+betae*lambdae*(w11*w12+w22*w21)*v1
			-beta2*sqrt(2.0)/12.0*cos(thetas)/2.0*(3.*(2.*phi2-1.0)^2-3.0)*dx(phi2)*v1)
+on(3,u1=0.0,u2=0.0)
+on(4,u1=uin,u2=0.0)
+on(1,2,u2=0.0);

problem press(p,q)=
    int2d(Th)(dx(p)*dx(q)+dy(p)*dy(q))
    +int2d(Th)(Re/dt*rhobar*(dx(u1)+dy(u2))*q-dx(p01)*dx(q)-dy(p01)*dy(q))
   +on(2,p=0.0);

// ****************initial value *****************
//phi2=1./2.+tanh((-sqrt(((x-2.0))^2+(y)^2)+0.5)/2/sqrt(2)/eps2)/2.0;
Vh phi21,phi22,phi23;
phi23=1./2.+(tanh((-sqrt((x-1.0)^2+(y)^2)+0.25)/eps2/sqrt(8)))/2;
 
phi2=phi23;
phi2=max(min(phi2,1.0),0.0);
real phimax=phi2[].max;
real phimin=phi2[].min;
real  vophi2 ,vophi20, verr2,vol;
 
vophi2=int2d(Th)(phi2);
vol=int2d(Th)(1.0);
cout<<"  the volume of phi2 is "<<vophi2<<endl;

 
ofstream fphi2("phi2.txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{fphi2<<phi2[][Vh(i,j)]<<endl;}
}

//plot(phi2,cmm= "  min=" + phi2[].min + ", max=" + phi2[].max,fill=true,wait=true,value=true);
//plot(phi1,cmm= "  min=" + phi1[].min + ", max=" + phi1[].max,fill=true,wait=true,value=true);
cout<<"the min value of phi2 is "<<  phimin<<endl;
cout<<"the max value of phi2 is "<< phimax<<endl;

psi1=-y;
psi2=x;

u1=0.0;
u2=0.0;
p=0.0;
p01=0.0;


u01=0.0;
u02=0.0;



//int Tstep=0;
real T=5.0;
real u1max;
real u2max;

// ******************Total Time iteration **********************

 
vophi20=vophi2;
int tstep=0;
Vh tmp;
real voltmp;






{ofstream ff("energy.txt");
{ofstream fvol("vol.txt");
for (t=0.0;t<T;t+=dt){
//output the time
tstep+=1; 
cout<<" The time is "<<t<<endl;
// update data//
u01=u1;
u02=u2;
p02=p01;
p01=p;
psi01=psi1;
psi02=psi2;
 
phi02=phi2;
 
//vophi20=vophi2;



//cout<<" the max eta is "<< eta[].max<<endl;
// solve the system  
step1;
//vophi20=int2d(Th)(phi2);
// modefy the phi2

phi2=max(min(phi2,1.000001),-0.000001);
phi2=(phi2>0.001)*phi2;
vophi2=int2d(Th)(phi2);
tmp=(phi2<0.999&phi2>0.01)*1.0;
voltmp=int2d(Th)(tmp);
//verr2=(vophi20-vophi2-(int1d(Th,2)(u01*phi02)
//+int1d(Th,4)(-u01*phi02)))/voltmp;
verr2=(vophi20-vophi2)/voltmp;
phi2=phi2+tmp*verr2;
//phi2=max(min(phi2,1.0),0.0);
vophi2=int2d(Th)(phi2);
cout<<"ph2 max is "<<phi2[].max<<"phi2 min is "<<phi2[].min<<endl;
 
//cout<<"the volume of phi1 is "<< vophi1<<"  the volume of phi2 is "<<vophi2<<endl;
fvol<<vophi2<<endl;

 ul2=sqrt(int2d(Th)(u01^2+u02^2));
if (ul2==0.0) {eps=0.0;
}else{
eps=1.0/ul2/2.0;
}

step3;

w11=dx(psi1);
w12=dy(psi1);
w21=dx(psi2);
w22=dy(psi2);

//getw;
lambdae= phi2*phi2;
eta= phi2*etap+(1.0-phi2)*etaf;
rho= phi2+(1.0-phi2)*rho12;
rhobar= rho12;
w1=-dx(lambdae)*w11-dy(lambdae)*w12-lambdae*(dx(w11)+dy(w12));
w2=-dx(lambdae)*w21-dy(lambdae)*w22-lambdae*(dx(w21)+dy(w22));

NS;
//dux=2.0*dx(u1);
press;

//plot(u2,cmm= "  min=" + u2[].min + ", max=" + u2[].max,fill=true,wait=true,value=true);
//plot(u1,cmm= "  min=" + u1[].min + ", max=" + u1[].max,fill=true,wait=true,value=true);
// *************************Check the Energy*********************
kinetic=Re/2.*int2d(Th)(u1^2+u2^2);
elastic=betae/2.*int2d(Th)(lambdae*(w11^2+w12^2+w21^2+w22^2));
mix=int2d(Th)(beta2*( (0.25*phi2^2*(1-phi2)^2)/eps2+eps2/2.0*(dx(phi2)*dx(phi2)+dy(phi2)*dy(phi2)))
);
ewall=-int1d(Th,1,2,3)(beta2*eps2*sqrt(2.0)/12.0*cos(thetas)*(2.0*phi2-1.0)*(3.0-(2.0*phi2-1)^2))/4;

etot=kinetic+elastic+mix+ewall;

cout<<"the kinetic energy is "<<kinetic<<" the mix energy is "<<mix <<" the elastic energy is "<<elastic<<" the wall energy is "<<ewall<<endl;
cout<<"the total energy is "<< etot<<endl;
ff<<"the kinetic energy is "<<kinetic<<" the mix energy is "<<mix <<" the elastic energy is "<<elastic<<" the wall energy is "<<ewall<<endl;
ff<<"the total energy is "<< etot<<endl;
//ft<<t<<endl;
//u1max=max(u1[].max, abs(u1[].min));
//u2max=max(u2[].max,abs(u2[].min));
//if (max(u1max,u2max)<10){
//dt=min(1e-4, min(0.001/u1max,0.001/u2max));
//}
//else
//{dt=1e-5;
//}

cout<<"the time step is"<< dt<<endl;
// remeber  the temporary  data

ofstream fphi2t ("phi2t.txt"); fphi2t << phi2 [];
ofstream fu1t ("fu1t.txt"); fu1t << u1 [];
ofstream fu2t ("fu2t.txt"); fu2t << u2 [];
ofstream fpsi1t ("fpsi1t.txt"); fpsi1t << psi1 [];
ofstream fpsi2t ("fpsi2t.txt"); fpsi2t << psi2 [];
ofstream fpt ("fpt.txt"); fpt << p [];
ofstream fp01t ("fp01t.txt"); fp01t << p01 [];

ofstream ftimet ("ftime.txt"); ftimet << t;





// output the data 

if(tstep%100==0 ){

ofstream ff2("phi2"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ff2<<phi2[][Vh(i,j)]<<endl;}

}

ofstream ffu1("u1"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ffu1<<u1[][Vh(i,j)]<<endl;}
}

ofstream ffu2("u2"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ffu2<<u2[][Vh(i,j)]<<endl;}
}

ofstream ffp("press"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ffp<<p[][Wh(i,j)]<<endl;}
}

ofstream ffpsi1("psi1"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ffpsi1<<psi1[][Vh(i,j)]<<endl;}
}

ofstream ffpsi2("psi2"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ffpsi2<<psi2[][Vh(i,j)]<<endl;}
}




}

}
}
};
 
 
 
//**********************The End of Process************************