//=========================================================================================
// 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.01;//1./real(meshsize)/10.;

real    Re=2.5e-2;
real    etaf=1e-1,etap=1e0,etan=1e0,etar=1.e0; // viscosity 
real    beta1=0.2,beta2=0.2,beta3=0.1;// cohesion energy 
real 	betae=0.0;// 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=8.88889;
real 	dh=0.3,dl=2.0;
func 	uin=umax*(dh-y)*y;
func 	psi1in=-y;
real    t;
func 	psi2in=x;
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=0.0;//5.0e-3;
real    kperS=2e4;//1e4;
real 	kperC=5e4;//1e6;
real    knt=1e-1;
 
real    theta=1.0,DCW=6e-5,DCS=1e-5,DCC=4e-6;
 

// *************** 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.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,"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)))  
)
+on(4,psi1=psi1in,psi2=psi2in);

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);



Wh C,C0,Ch,Ct;

Wh DC;
 
problem diff(C,Ch)=
int2d(Th)(C*Ch
+theta*(u01*dx(C)+u02*dy(C))*Ch*dt
+theta*DC*(dx(C)*dx(Ch)+dy(C)*dy(Ch))*dt
+theta*(C+u01*dx(C)*dt+u02*dy(C)*dt)*hTriangle*(u01*dx(Ch)+u02*dy(Ch))
)
+int2d(Th)(-C0*(Ch+ (hTriangle)*(u01*dx(Ch)+u02*dy(Ch)))
+(1-theta)*(u01*dx(C0)+u02*dy(C0))*Ch*dt
+(1.0-theta)*DC*(dx(C0)*dx(Ch)+dy(C0)*dy(Ch))*dt
+(1.0-theta)*(C0+u01*dx(C0)*dt+u02*dy(C0)*dt)*hTriangle*(u01*dx(Ch)+u02*dy(Ch))
) 
//+int1d(Th,4)((1.0+hTriangle)*(0.0*Ch));
+on(4,C=0);



// ****************initial value *****************
 
 
Vh phia1,phia2,phia3;
 

phia1=1./2.+tanh((-sqrt(((x-0.6)/0.4)^2+((y/0.15)^2))+1.0)/0.05/sqrt(8.0))/2.0;
phia2=1./2.+ (tanh((-sqrt(((x-0.4)/0.2)^2+(y/0.1)^2)+1.0)/0.05/sqrt(8.0)))/2.0;


// phi3 initial

phi3=0.6*phia1+0.2*phia2;
phi3=max(min(phi3,1.0),0.0);
 


// phi2 initial 

phi2=phia2;//-0.3*phia4;
phi2=max(min(phi2,1.0),0.0);
phi1=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;

C=tanh(x/0.01);
C=phia1;
//C=1.0;
real VolC=int2d(Th)(C);
 
cout<<"Vol of C " <<VolC<<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(C);
//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=700.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;
C0=C;
//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);
//};*/
 
growth = 0.0;
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;
}

eps=1.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;


diff;
C=max(C,0.0);
C=min(C,1);
 VolC=int2d(Th)(C);
//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%100==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;}
}


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

}

}
}
};

 
 
 
//**********************The End of Process************************