// Stokes Quasi static avec tension superficiel
if (HaveUMFPACK) //  need UMPACK  computation 
{
verbosity=0; 
// voir https://fr.wikipedia.org/wiki/Pression_de_Laplace
load "isoline"
load "Curvature"

real X0 = 0.1, Y0 =1*X0,Y1=0, CC0=Y0; 
real meshsize=X0/10;
real CFL = 0.1; 
real dt = 0.005;
int niter = 100;
real Kc =10; 
real rhog = 0; 
real nu=1;
real Vin =1; 
int nlissage=2;
real X02=X0*X0; 
func u2in=(-1+x*x/X02)*Vin;
real fluxin;
border ctop(t=X0,0){x=t; y=Y0;label=1;};
border cleft(t=Y0,Y1){x=0; y=t;label=3;};
border cright(t=0,Y0){x=X0; y=t;label=0;}
int nsl=100; 
real[int,int] SL(3,nsl) ;
for(int i=0; i< nsl; ++i)
  { 
    real ll = i/(nsl-1.); 
    SL(0,i) = 0*(1-ll)+X0*ll;
    real l0=ll;
    real RR=1.2*1.2, YY=sqrt(RR-1);
    SL(1,i) = -CC0*(sqrt(RR-l0*l0)-YY);
    SL(2,i) = ll; 
  }
real lgc = setcurveabcisse(SL); 
 {  nsl = lgc/meshsize; 
   real[int,int] SL1=equiparameter(SL,nsl); 
   SL.resize(SL1.n,SL1.m);
   SL = SL1;  
 }
border bSL(t=0,1) 
{
    P= Curve(SL,t); 
    label=2; 
} 
Y1= SL(1,0);  
int nn = 20;
func cfix = ctop(X0/meshsize)+cleft(Y0/meshsize)+cright(Y0/meshsize); 
plot( cfix+ bSL(lgc/meshsize) , wait=1); 
func Bmesh = buildmesh( cfix+ bSL(lgc/meshsize) );
mesh Th= Bmesh;

plot(Th, wait=1);
// Formule (2.9) of  http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch02.d/AFEM.Ch02.pdf 
//  macro Eps(u,v) [ dx(u),dy(v), (dx(v)+dy(u))*0.5, (dx(v)+dy(u))*0.5 ] // 
// here r==x, z==y
macro Eps(u1,u2)  [dx(u1),dy(u2),u1/x, (dx(u2)+dy(u1))*0.5,(dx(u1)+dy(u2))*.5 ] // EOM

macro Grad(u,v) [ dx(u),dy(v), dx(v), dy(u) ] // 
macro div(u1,u2) (dx(u1)+dy(u2) ) // 
macro rdiv(ur,uy) (x*dx(ur) + ur + x*dy(uy)) //

func mN = [ N.x,N.y]*[ N.x,N.y]'; 
func mT = [ N.y,-N.x]*[ N.y,-N.x]'; 
real bN = 1e6, bT = 1e6, bT1 =1e6; // Coef de 
real bT0 = 1; // frottement a la paroi

fespace Mh(Th,P1);
fespace Wh(Th,P2);
Mh rfcourbure;

rfcourbure[]=0;
int[int] llc=[2];
//plot(frcourbure,wait=1,fill=1); 
//real fcourbure =1;
fespace Vh(Th,[P2,P2,P1]); 
Vh [u1,u2,p],[v1,v2,q];
Wh U;
real temps =0; 
for(int iter=0; iter<niter;++iter)
{
    rfcourbure=0;
    rfcourbure[]=raxicurvature(Th,2);
    u1[]=0;
solve PMeca2( [u1,u2,p],[v1,v2,q] ) =
int2d(Th) ( x*2*nu*Eps(u1,u2)'*Eps(v1,v2) - rdiv(v1,v2)*p -rdiv(u1,u2)*q)
 + on(1,u1=0,u2=u2in)
 + on(0,u1=0,u2=0)
 + on(3,u1=0)
 + int1d(Th,2) ( Kc*rfcourbure*[N.x,N.y]'*[v1,v2]) 
 - int2d(Th)( -x*rhog*v2)
;
//plot(rfcourbure,WindowIndex=1);
fluxin = int1d(Th,1)(-u2);

real[int,int] UU(2,nsl);
real Umax=0;
for(int i=0; i<nsl-1;++i)
  {
    x= SL(0,i);
    y= SL(1,i); 
    real Tx=1,Ty=0;
    if(i>0)
    {
    Tx=SL(0,i+1)-SL(0,i-1);
    Ty=SL(1,i+1)-SL(1,i-1);
    }
    
    real ll = sqrt(Tx*Tx+Ty*Ty); 
    real Nx = -Ty/ll;
    real Ny =  Tx/ll;
    real cN= Nx*u1+Ny*u2;
    UU(0,i)=cN*Nx;
    UU(1,i)=cN*Ny;
    Umax=max(Umax,abs(cN)); 
  }
  //  lissage de UU
  real bb= 0.9, aa = (1-bb)/2.; 
  for(int j=0; j<nlissage;++j)
  for(int i=1; i<nsl-2;++i)
    { 
        UU(0,i) = aa*UU(0,i-1)+bb*UU(0,i)+aa*UU(0,i+1);
        UU(1,i) = aa*UU(1,i-1)+bb*UU(1,i)+aa*UU(1,i+1);
    }

  dt = meshsize*CFL/Umax;
  real Vadd = fluxin*dt;
  real volUU=0;
  
  cout << "nsl =" << nsl << " "<<UU.n << " " << UU.m << " " << SL.m << endl; 
  UU(0,nsl-1)=0;
  UU(1,nsl-1)=0;
  for(int i=0; i<nsl-1;++i)
   { 
     SL(0,i) += dt*UU(0,i);
     SL(1,i) += dt*UU(1,i);
   }
temps += dt;

 cout << iter << temps << " " << " dt = " << dt << " " << Umax << endl; 
lgc = setcurveabcisse(SL); 
if( iter% 10 ==9) 
{
   nsl = lgc/meshsize; 
   cout << "nsl =" << nsl << endl; 
   real[int,int] SL1=equiparameter(SL,nsl); 
   SL.resize(SL1.n,SL1.m);
   SL = SL1;  
   if(verbosity>2) 
   for(int i=1; i<nsl-1;++i)
     cout << i << " " << SL(0,i) << " " << SL(1,i) << " " <<  SL(2,i)  <<endl; 
}
plot([u1,u2],p,bSL(lgc/meshsize),coef=dt,cmm=" temps " + temps+ " / " + dt + " / " +Umax );
 Y1= SL(1,0);
Th= Bmesh;
}
cout << " End loop \n"; 
}
