// bs2 option
int s=10; // y-scale
int m=50;
int S1max=200; // x: variable S1 
int S2max=S1max; // y: variable s2;
border aa(t=0,S1max){x=t;y=0;};
border bb(t=0,S2max){x=S1max;y=t;}; // bord droit
border cc(t=S1max,0){x=t;y=S2max;}; // bord haut
border dd(t=S2max,0){x = 0; y = t;};

mesh th = buildmesh(aa(m)+bb(m)+cc(m)+dd(m));
fespace Vh(th,P1);

real sigmax=0.2;
real sigmay=0.2;
real rho=0.3; 
real r=0.05;
real K=100;
real T=1;
real N=10; real nstop=1;
real dt=T/N;

real eps=0.3;

func f = max(K-max(x,y),0.);
//func f = max(K-x,0.);

Vh u=f,v,w;
Vh g=0; //- pour la condition limite du/dn=g ou u = g . 

plot(u,th,wait=1);

th = adaptmesh(th,u,abserror=1,nbjacoby=2,
   err=0.004, nbvx=5000, omega=1.8,ratio=1.8, nbsmooth=3,
   splitpbedge=1, maxsubdiv=5,rescaling=1 ); 
plot(th,wait=1);
u=u;

Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
Vh yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;


int j=0;
int n;
int jstop=N/nstop;
cout << "N=" << N << endl;
cout << "jstop=" << jstop << endl;

//-----------------
//- Utiliser l'option qft=qf1pTlump pour le "Mass Lumping"
//-----------------
problem eq1(u,v,init=j,solver=CG) = int2d(th)(
//problem eq1(u,v,init=j,solver=CG) = int2d(th,qft=qf1pTlump)( 
    u*v*(r+1/dt) 
    + dx(u)*dx(v)*(x*sigmax)^2/2. 
    + dy(u)*dy(v)*(y*sigmay)^2/2.
    + dy(u)*dx(v)*rho*sigmax*sigmay*x*y/2.
    + dx(u)*dy(v)*rho*sigmax*sigmay*x*y/2.   )
  + int2d(th)( -convect([xveloc,yveloc],-dt,w)/dt * v ) 
//  + int2d(th,qft=qf1pTlump)( -convect([xveloc,yveloc],-dt,w)/dt * v ) 
  + on(bb,cc,u=0)
;


int ww=1;
for ( n=0; n*dt < T; n++) 
{
  cout <<" iteration " << n   <<  " j=" << j << endl; 
  w=u; 
  eq1; 
  // u=max(u,f); //- Pour l'aspect Americain (splitting); 
  plot(u,wait=ww,cmm="Isovaleurs de u");
  ww=0;
  j=j+1; 
  if(j==jstop)  { cout << " adaptmesh " << endl;
    th = adaptmesh(th,u,verbosity=1,abserror=1,nbjacoby=2,
    err=0.001, nbvx=5000, omega=1.8, ratio=1.8, nbsmooth=3,
    splitpbedge=1, maxsubdiv=5,rescaling=1) ; 
    plot(th,wait=1);
    j=0;
    xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
    yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;

    u=u;
    ww=1;
    };
  //j=j+1; 
  cout << " j = " <<  j << endl;
};
//v = max(u-f,0.);
//plot(v,wait=1,value=1);
plot(u,wait=0,value=1);
//real normvL2=sqrt(int2d(th)(v*v));
//cout << " ||v||_L2 = " << normvL2 << endl;
real s1=K,s2=K;
cout << " P(s1,s2)= " << u(s1,s2) << endl;
cout << " programme end normaly " << endl;