//=========================================================================================
// PDE Coupling Newtonian and Viscoelastic fluid using Cahn-Hilliard equation
//=========================================================================================
 // this code is modifyed to involve the effect of red blood cell 

// phi3 is the volume fraction of clot
//1-phi3 is the volume fraction of plasma
// phi2*phi3 is the volume fraction of fibrin and platelet 
//(1-phi2)*phi3 is the volume fraction of red blood cell
// phi1*phi2*phi3 is the volume fraction of fibrin 
//(1-phi1)*phi2*phi3 is the volume fraction of platelet 

  




int     meshsize=128;
int     meshcircle = 300;
real    dt=0.001;//1./real(meshsize)/10.;

real    Re=2.5e-2;
real    etaf=1e-1,etap=5e1,etan=1e1,etar=2.5e0; // viscosity 
real    beta1=0.2,beta2=0.2,beta3=0.1;// cohesion energy 
real 	betae=0.1;// eleasticity 
real    s1,s2,s3; // convex splitting 
 

real    tau1=1e-4,tau2=1e-4,tau3=1e-4;// mobility 
real    alphanp=1e1,alphanr=1.0; // elasticity diff
real 	alpha=1e0,beta=0.5/3.0,gamma=0.0; // G_1 coffecient 
real 	alpha22=1e1,beta22=0.2/3.0;
real    eps1=0.01,eps2=0.01,eps3=0.01;// surface width
real    thetas=pi/2.0;
//real 	a=abs(3*alpha-6*alpha*beta)/2.0;
		s1=	abs(3*alpha-6*alpha*beta);
		s2=5.0;
		s3=5.0;
//real 	a2=abs(3*alpha22-6*alpha22*beta22)/2.0;

//real    b=7.0;
real 	eps; // used in SUPG should be update in each time setp
real 	umax=40.0;
real 	dh=0.5,dl=2.0;
func 	uin=umax*(dh-y)*y;
func 	psi1in=-y;
real    t;
func 	psi2in=x-umax*y*(dh-y)*t;
real 	ul2;

cout<<"dt is "<< dt<<endl;
cout<<"Re is "<< Re<<endl;
cout<<"beta1 is "<< beta1<<" beta2 is "<<beta2<<" beta3 is "<<beta3<<endl;
cout<<"eps1 is "<< eps1<<" eps2 is "<<eps2<<endl;
cout<<"tau1 is "<< tau1<<" tau2 is "<<tau2<<endl;
cout<<"alpha is "<<alpha<<" beta is "<<beta<<" s1 is "<< s1 <<endl;
cout<<" cos(contact angle) is " << cos(thetas)<<endl;
real elastic;
real kinetic;
real mix;
real ewall;
real etot;
real    phi3max=0.3;
real    keff=5.0e-3;
real    kperS=1e4;
real 	kperC=1e6;
real    knt=1e-1;
 

 

// *************** Mesh**********************************
 
border 	bb1(t=0,1.0){x=t; y=0; label=1;};
border 	bb2(t=1.0,3.0){x=t; y=0; label=1;};
border 	br(t=0,dh){x=3.0; y=t; label =2;};
border 	bt1(t=1.0,0){x=t; y=dh; label =3;};
border	bt2(t=3.0,1.0){x=t;y=dh; label =3;};
border 	bl(t=dh,0){x=0; y=t; label =4;};
 

mesh Th= buildmesh (bt1(meshsize)+bt2(meshsize*2)+bb1(meshsize)+bb2(meshsize*2)
		+ bl(meshsize*dh) + br(meshsize*dh));
 




//plot(Th,ps="meshe.eps");
savemesh(Th,"umax50.msh");
savemesh(Th,"meshe.txt");

 
// ********* 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.
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  phi1,mu1,phi01,phi1h,mu1h;
Vh  phi2,mu2,phi02,phi2h,mu2h;
Vh  phi3,mu3,phi03,phi3h,mu3h;
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;
Vh  growth,delta,kper;
Wh  stress,stressn,stresstnorm,dxu1,dyu1,dxu2,dyu2,nx,ny,dxphi3,dyphi3,length, stresst;
//***********Define the Problem ********************

// ******************step  1 solve phi3 ********************

problem step1([phi3,mu3],[phi3h,mu3h])=
int2d(Th)(1./dt*phi3*mu3h+tau3*(dx(mu3)*dx(mu3h)+dy(mu3)*dy(mu3h))
			-mu3*phi3h
			-phi3*(u01*dx(mu3h)+u02*dy(mu3h))
			+beta3*s3/eps3*phi3*phi3h+beta3*eps3*(dx(phi3)*dx(phi3h)+dy(phi3)*dy(phi3h))
			+betae*((1.0-phi02)^2*alphanr+phi02^2*(phi01+(1-phi01)*alphanp))*phi3*phi3h
			*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))

)

+int2d(Th)(-1./dt*phi03*mu3h			
//			-phi03*(u01*dx(mu3h)+u02*dy(mu3h))
			-growth*mu3h
			-beta3*s3/eps3*phi03*phi3h+beta3/eps3*(phi03^3-1.5*(phi03)^2+0.5*phi03)*phi3h
			+beta2/eps2*(phi02^2*(1.0-phi02)^2/4.0)*phi3h
			+beta2*eps2*(dx(phi02)*dx(phi02)+dy(phi02)*dy(phi02))*phi3h*0.5
			+beta1/eps1*phi02*(alpha*phi01^4/4.0-alpha*beta*phi01^3)*phi3h
			+beta1*eps1*phi02*(dx(phi01)*dx(phi01)+dy(phi01)*dy(phi01))*phi3h*0.5
//			+betae*((1.0-phi02)*alphanr+phi02*(phi01+(1-phi01)*alphanp))*phi03*phi3h
//			*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))

)

+int1d(Th,2)(u01*phi03*mu3h)
+int1d(Th,4)(-u01*phi03*mu3h);




// ++++++++++step 2:The problem of phi2, mu2+++++++++
 
//     energy split method 
problem step2([phi2,mu2],[phi2h,mu2h])=
int2d(Th)(1./dt*phi2*mu2h+tau2*(dx(mu2)*dx(mu2h)+dy(mu2)*dy(mu2h))
// 		+dt/Re*phi02^2*(dx(mu2)*dx(mu2h)+dy(mu2)*dy(mu2h))
		-mu2*phi2h
		-phi2*(u01*dx(mu2h)+u02*dy(mu2h))
		+beta2/eps2*s2*phi2*phi2h
		+eps2*beta2*(dx(phi2)*dx(phi2h)+dy(phi2)*dy(phi2h))
//		+beta3*(phi01+(1.0-phi01)*alphanp)*phi2
//*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))*phi2h
		+betae*((phi01+(1.0-phi01)*alphanp)+alphanr)*phi3^2*phi2*phi2h
*(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))

)
//+int1d(Th,1)(//beta2*gamma/dt*phi2*phi2h-beta2*gamma*dt/Re*dx(mu2)*dx(phi02)*phi02*phi2h
//+beta2*sqrt(2.0)/12.0*phi2*phi2h
//)

+int2d(Th)(-1./dt*phi02*mu2h
//	-phi02*(u01*dx(mu2h)+u02*dy(mu2h))
-beta2/eps2*s2*phi02*phi2h+beta2/eps2*phi3*(phi02^3-1.5*(phi02)^2+0.5*phi02)*phi2h
-eps2*beta2*(dx(phi02)*dx(phi2h)+dy(phi02)*dy(phi2h))
+eps2*beta2*phi3*(dx(phi02)*dx(phi2h)+dy(phi02)*dy(phi2h))
+beta1/eps1*phi3*(alpha*phi01^4/4.0-alpha*beta*phi01^3)*phi2h
+beta1*eps1*phi3*(dx(phi01)*dx(phi01)+dy(phi01)*dy(phi01))*phi2h*0.5
+betae*(-alphanr)*phi3^2
*(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*(3.*(2.*phi02-1.0)^2-3.0)*phi2h
//-beta2*sqrt(2.0)/12.0*phi02*phi2h
//-beta2*gamma/dt*phi02*phi2h
//);


// ++++++++++step3: The problem of phi1, mu1+++++++++

problem step3([phi1,mu1],[phi1h,mu1h])=
int2d(Th)(phi1*mu1h/dt+tau1*(dx(mu1)*dx(mu1h)+dy(mu1)*dy(mu1h))
// 			+dt/Re*(phi01)^2*(dx(mu1)*dx(mu1h)+dy(mu1)*dy(mu1h))
 			-mu1*phi1h
			-phi1*(u01*dx(mu1h)+u02*dy(mu1h))
			+beta1/eps1*s1*phi1*phi1h
			+beta1*eps1*(dx(phi1)*dx(phi1h)+dy(phi1)*dy(phi1h))
)

+int2d(Th)(-phi01*mu1h/dt
//			-phi01*(u01*dx(mu1h)+u02*dy(mu1h))
			-beta1/eps1*s1*phi01*phi1h 
			-beta1*eps1*(dx(phi01)*dx(phi1h)+dy(phi01)*dy(phi1h))
			+beta1*eps1*phi2*phi3*(dx(phi01)*dx(phi1h)+dy(phi01)*dy(phi1h))
			+beta1/eps1*phi3*phi2*(alpha*(phi01)^3-3.0*alpha*beta*(phi01)^2)*phi1h
			+betae*(1.0-alphanp)*phi2^2*phi3^2*0.5*
			(dx(psi01)*dx(psi01)+dx(psi02)*dx(psi02)+dy(psi01)*dy(psi01)+dy(psi02)*dy(psi02))*phi1h
)
+int1d(Th,2)(u01*phi01*mu1h)
+int1d(Th,4)(-u01*phi01*mu1h);


// ++++++++++++step4: update the Psi++++++++++++++++


problem step4([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*(u1*v1+u2*v2)+Re*(u01*dx(u1)+u02*dy(u1))*v1+Re*(u01*dx(u2)+u02*dy(u2))*v2
		+0.5*Re*(dx(u01)+dy(u02))*(u1*v1)+0.5*Re*(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))
		+kper*phi3^2*eta*(u1*v1+u2*v2)

	//	-beta3*lambdae*(dx(psi1)*dx(w11*v1)+dy(psi1)*dy(w11*v1))
	//	-beta3*lambdae*(dx(psi2)*dx(w21*v1)+dy(psi2)*dy(w21*v1))
	//	-beta3*lambdae*(dx(psi1)*dx(w12*v2)+dy(psi1)*dy(w12*v2))
	//	-beta3*lambdae*(dx(psi2)*dx(w22*v2)+dy(psi2)*dy(w22*v2))
	 )
 
+int2d(Th)(-Re/dt*(u01*v1+u02*v2)-(2.0*p01-p02)*(dx(v1)+dy(v2))
		 	+(phi1*dx(mu1)+phi2*dx(mu2)+phi3*dx(mu3))*v1
			+(phi1*dy(mu1)+phi2*dy(mu2)+phi3*dy(mu3))*v2
		-betae*(dx(psi1)*w1*v1+dx(psi2)*w2*v1)
		-betae*(dy(psi1)*w1*v2+dy(psi2)*w2*v2)	
//		-beta3*lambdae*(dx(psi1)*(dx(w11)*v1+w11*dx(v1))+dy(psi1)*(dy(w11)*v1+w11*dy(v1)))
//		-beta3*lambdae*(dx(psi2)*(dx(w21)*v1+w21*dx(v1))+dy(psi2)*(dy(w21)*v1+w21*dy(v1)))
//	 	-beta3*lambdae*(dx(psi1)*(dx(w12)*v2+w12*dx(v2))+dy(psi1)*(dy(w12)*v2+w12*dy(v2)))
//		-beta3*lambdae*(dx(psi2)*(dx(w22)*v2+w22*dx(v2))+dy(psi2)*(dy(w22)*v2+w22*dy(v2)))
		 )
 
+on(1,3,u1=0.0,u2=0.0)
+on(4,u1=uin,u2=0.0)
+on(2,u2=0.0);

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






// ****************initial value *****************
 
 
Vh phia1,phia2,phia3;
//phia1=1./2.+tanh((-sqrt(((x-0.6))^2+((y+0.)^2))+0.35)/0.01/sqrt(8.0))/2.0;
//phia2=1./2.+ (tanh((-sqrt(((x-1.0)/0.5)^2+((1.4*y*x^(1./2.)-0.2)/0.1)^2)+1.0)/0.05/sqrt(8.0)))/2.0;
//phia3=1./2.+(tanh((-sqrt(((x-0.55)/0.2)^2+(y/0.1)^2)+1.0)/0.05/sqrt(8.0)))/2;
 
//phia1=1./2.+tanh((-sqrt(((x+0.3*y-0.6)/0.4)^2+((y/0.35)^2))+1.0)/0.02/sqrt(8.0))/2.0;
//phia2=1./2.+ (tanh((-sqrt(((x+y-1.0)/0.5)^2+((1.8*y*x^(1.5)-0.18)/0.12)^2)+1.0)/0.05/sqrt(8.0)))/2.0;
//phia3=1./2.+(tanh((-sqrt(((x-0.6)/0.2)^2+(y/0.1)^2)+1.0)/0.05/sqrt(8.0)))/2;


phia1=1./2.+tanh((-sqrt(((x+0.3*y-0.9)/0.4)^2+((y/0.35)^2))+1.0)/0.02/sqrt(8.0))/2.0;
phia2=1./2.+ (tanh((-sqrt(((x+y-1.3)/0.5)^2+((1.8*y*x^(1.2)-0.25)/0.15)^2)+1.0)/0.05/sqrt(8.0)))/2.0;
phia3=1./2.+(tanh((-sqrt(((x-0.9)/0.2)^2+(y/0.15)^2)+1.0)/0.05/sqrt(8.0)))/2;

// phi3 initial

phi3=phia1+(x>1.0)*phia2*0.8;
phi3=max(min(phi3,1.0),0.0);
phi3=0.3*phi3+0.4*phia3; 


// phi2 initial 

phi2=phia3;//-0.3*phia4;
phi2=max(min(phi2,1.0),0.0);
//phi1= (y-0.05)*(y>0.05&&phi2>0.1)*2+0.5*(phi2>=0.1);


// phi1 initial 


phi1=0.2*phia3;//+0.1*phi21;
phi1=max(min(phi1,1.0),0.0);
real phimax=phi2[].max;
real phimin=phi2[].min;

real vophi1,vophi2,vophi3, vophi10,vophi20,vophi30,verr1,verr2,verr3,vol;
vophi1=int2d(Th)(phi1);
vophi2=int2d(Th)(phi2);
vophi3=int2d(Th)(phi3);

vol=int2d(Th)(1.0);
cout<<"the volume of phi1 is "<< vophi1<<"  the volume of phi2 is "<<vophi2<<endl;

ofstream fphi1("phi1.txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{fphi1<<phi1[][Vh(i,j)]<<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;}
}

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

Vh core, shell, ple,fib;
core=phi2*phi3;
shell=(1.0-phi2)*phi3;
ple=(1.0-phi1)*phi2*phi3;
fib=phi1*phi2*phi3;
//plot(phi3,cmm= "  min=" + phi3[].min + ", max=" + phi3[].max,fill=true,wait=true,value=true);
//plot(core,fill=true,wait=true,value=true);
//plot(shell, fill=true,wait=true,value=true);
//plot(fib, fill=true,wait=true,value=true);
//plot(ple, 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=3.0;
real u1max;
real u2max;

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

vophi10=vophi1;
vophi20=vophi2;
vophi30=vophi3;
int tstep=0;

Vh tmp;
real voltmp,growthall;


growthall=0.0;
eta=(phi1*phi2*etan+(1.0-phi1)*phi2*etap+(1-phi2)*etar)*phi3+(1.0-phi3)*etaf;

{ofstream ff("energy.txt");
{ofstream fvol("vol.txt");
//{ofstream ft("time.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;

phi01=phi1;
phi02=phi2;
phi03=phi3;

//vophi10=vophi1;
//vophi20=vophi2;
//vophi30=vophi3;


//cout<<" the max eta is "<< eta[].max<<endl;

// solve the system phi3 


dxu1=dx(u01);
dxu2=dx(u02);
dxphi3=dx(phi03);
dyu1=dy(u01);
dyu2=dy(u02);
dyphi3=dy(phi03);
length=sqrt(dxphi3^2+dyphi3^2)+eps3;

nx =dxphi3 /length ;
ny =dyphi3 /length;
stressn=p-eta*(nx*(2.0*dxu1*nx+(dxu2+dyu1)*ny)+ny*((dxu2+dyu1)*nx+2.0*dyu2*ny));
stresst=eta*(ny*(2.0*dxu1*nx+(dxu2+dyu1)*ny)-nx*((dxu2+dyu1)*nx+2.0*dyu2*ny));
stresstnorm=abs(stresst);
cout<<"tangential stress"<<stresstnorm[].max<<"  "<<stresstnorm[].min<<endl;
cout<<"normal stress"<<stressn[].max<<"  "<<stressn[].min<<endl;
//stressnorm=sqrt((2.0*dxu1*nx+(dxu2+dyu1)*ny)^2+((dxu2+dyu1)*nx+2.0*dyu2*ny)^2);
//if(tstep>10)
//{//delta=6.0*sqrt(2.0)*eps3*(dx(phi03)^2+dy(phi03)^2);
delta=sqrt(dx(phi03)^2+dy(phi03)^2);
stress=knt*stressn-stresstnorm;
//growth=(phi03>=0.01&&phi03<phi3max&&stress>=0)*keff*(stress)*(1.0-phi03/phimax)*phi03;//*sqrt(dx(phi03)^2+dy(phi03)^2);
growth=(phi03>=0.1&&phi03<phi3max&& delta>=4.0)*keff*(stress)*delta*phi03;
cout<<growth[].max<<"  "<<growth[].min<<endl;
//plot(stress,fill=true,wait=true,value=true);
//plot(growth,fill=true,wait=true,value=true);
//};
 

growthall=growthall+int2d(Th)(growth); 


 

 
step1;
vophi30=int2d(Th)(phi3);
phi3=max(min(phi3,1.000001),-0.000001);
phi3=(phi3>0.0001)*phi3;
tmp=(phi3<0.999&phi3>0.01)*1.0;
voltmp=int2d(Th)(tmp);
//vophi3=int2d(Th)(tmp*phi3);
//vophi30=int2d(Th)(tmp*phi03);
//verr3=(vophi30-vophi3-(int1d(Th,2)(u01*phi03)
//+int1d(Th,4)(-u01*phi03)))/voltmp;
//cout<<"the modify error is : "<< verr3<<endl;
vophi3=int2d(Th)(phi3);
verr3=(vophi30-vophi3)/voltmp;
phi3=phi3+tmp*verr3;
//phi3=max(min(phi3,1.0),0.0);
vophi3=int2d(Th)(phi3);



cout<<"phi3 min is "<<  phi3[].min<<"phi3 max is "<<  phi3[].max<<endl;
//cout<<"the max value of phi2 is "<< phi2[].max<<endl;
// modefy the phi2
// solve the system phi2
step2;
//vophi20=int2d(Th)(phi2);
phi2=max(min(phi2,1.000001),-0.000001);
phi2=(phi2>0.0001)*phi2;
tmp=(phi2<0.999&phi2>0.01)*1.0;
voltmp=int2d(Th)(tmp);
vophi2=int2d(Th)(phi2);
//verr2=(vophi20-vophi2-(int1d(Th,2)(u01*phi02)
//+int1d(Th,4)(-u01*phi02)))/vol;
verr2=(vophi20-vophi2)/voltmp;
phi2=phi2+tmp*verr2;
phi2=max(min(phi2,1.0),0.0);
vophi2=int2d(Th)(phi2);

// solve the system phi118
step3;
//vophi10=int2d(Th)(phi1);
phi1=max(min(phi1,1.000001),-0.000001);
phi1=(phi1>0.0001)*phi1;
tmp=(phi1<0.999&phi1>0.01)*1.0;
vophi1=int2d(Th)(phi1);
//verr1=(vophi10-vophi1-(int1d(Th,2)(u01*phi01)
//+int1d(Th,4)(-u01*phi01)))/vol;
verr1=(vophi10-vophi1)/voltmp;
phi1=phi1+tmp*verr1;
phi1=max(min(phi1,1.0),0.0);
vophi1=int2d(Th)(phi1);


//cout<<"the volume of phi1 is "<< vophi1<<"  the volume of phi2 is "<<vophi2<<endl;
fvol<<vophi1<<"        "<<vophi2<<"      " << vophi3<<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);


//plot([u1,u2],p,wait=true);
//plot(psi1,wait=true);
//plot(psi2,wait=true);
//getus;
//us1=u01-dt/Re*(phi1*dx(mu1)+phi2*dx(mu2));
//us2=u02-dt/Re*(phi1*dy(mu1)+phi2*dy(mu2));


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

step4;

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

//getw;
lambdae=(phi1+(1.0-phi1)*alphanp)*phi2^2*phi3^2+(1-phi2)^2*phi3^2*alphanr;
eta=(phi1*phi2*etan+(1.0-phi1)*phi2*etap+(1-phi2)*etar)*phi3+(1.0-phi3)*etaf;
kper=kperS*(1.0-phi2)+kperC*phi2;
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=beta3/2.*int2d(Th)(lambdae*(w11^2+w12^2+w21^2+w22^2));
mix=int2d(Th)(beta1*phi2*phi3*((alpha*phi1^4/4.-alpha*beta*phi1^3-(alpha*(3.*beta)^4/4.-alpha*beta*(3.*beta)^3))/eps1
			+eps1/2.0*(dx(phi1)*dx(phi1)+dy(phi1)*dy(phi1)))
			+beta2*phi3*( (0.25*phi2^2*(1-phi2)^2)/eps2+eps2/2.0*(dx(phi2)*dx(phi2)+dy(phi2)*dy(phi2)))
			+beta3*( (0.25*phi3^2*(1-phi3)^2)/eps3+eps3/2.0*(dx(phi3)*dx(phi3)+dy(phi3)*dy(phi3)))
);
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 fphi1t ("phi1t.txt"); fphi1t << phi1 []; // to save the solution
ofstream fphi2t ("phi2t.txt"); fphi2t << phi2 [];
ofstream fphi3t ("phi3t.txt"); fphi3t << phi3 [];
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%10==0 ){

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

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

}

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

}

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

}

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

}

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

}

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

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

}

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

}

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 ff1("phi1"+tstep+".txt");
for (int i=0;i<Th.nt;i++)
{
for (int j=0;j<3;j++)
{ff1<<phi1[][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************************