gistfile1.txt

real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1;
real Outside = -1;
real Pe= 62;
real Ch=0.01;

func real initialCondition(real x, real y)
{
    if (sqrt(pow(x - centerX, 2) + pow(y - centerY, 2)) <= radius || 
       sqrt(pow(x - centerX2, 2) + pow(y - centerY, 2)) <= radius || 
       sqrt(pow(x - centerX, 2) + pow(y - 0.8, 2)) <= radius || 
       sqrt(pow(x - centerX2, 2) + pow(y - 0.8, 2)) <= radius)
    return Inside;
  else
    return Outside;
}

mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,P1);

real dt=0.01;
int i,k;  

Vh2 u,w,phi,oldv,oldz,psi;

u=initialCondition(x,y);
plot (u,wait=1,cmm="Cahn-Hilliard",value=true);


for (int i=0;i<30;i++)
      {oldv[]=u[];

    solve CahnHill1(w,phi) = 
     int2d(Th)(Ch*Ch*(dx(oldv)*dx(phi)+dy(oldv)*dy(phi)) )+ int2d(Th)((oldv^3-oldv)*phi) -int2d(Th)(w*phi);

    solve CahnHill2(u,phi)= int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);

      plot (u,wait=1,fill=true, cmm="Cahn-Hilliard_"+i,value=true);

      }