load "ff-Ipopt";
int nn=10;
mesh Th=square(nn,nn);
fespace Vh(Th,[P1,P1] );
fespace Wh(Th,[P1] );
int iter=0;

func f1 =  10;//right hand sides
func f2 = -15;
func g1 = -0.1;//Boundary conditions functions
func g2 =  0.1;

Vh [uz,uz2]=[1,1],[lz,lz2]=[1,1];
Wh lm=1.;
Vh [u1,u2]=[0,0];//starting point


while(++iter)
{
macro Grad(u) [dx(u),dy(u)]//gradient macro
varf vP([u1,u2],[v1,v2]) = int2d(Th)(Grad(u1)'*Grad(v1)+ Grad(u2)'*Grad(v2)) 
- int2d(Th)(f1*v1+f2*v2);

matrix A = vP(Vh,Vh);//Fitness function matrix...
real[int] b = vP(0,Vh);//and linear form

int[int] II1=[0],II2=[1];//Constraints matrix
matrix C1 =  interpolate (Wh,Vh, U2Vc=II1);
matrix C2 =  interpolate (Wh,Vh, U2Vc=II2);
matrix CC = -1*C1 + C2; // u2 - u1 >0
Wh cl=0;//constraints lower bounds (no upper bounds)

//Boundary conditions
varf vGamma([u1,u2],[v1,v2]) = on(1,2,3,4,u1=1,u2=1);
real[int] onGamma=vGamma(0,Vh);
Vh [ub1,ub2]=[g1,g2];
Vh [lb1,lb2]=[g1,g2];
ub1[] = onGamma ? ub1[] : 1e19  ; //Unbounded in interior
lb1[] = onGamma ? lb1[] : -1e19  ;


Vh [uzi,uzi2]=[uz,uz2],[lzi,lzi2]=[lz,lz2];
Wh lmi=lm;
Vh [ui1,ui2]=[u1,u2];

IPOPT([b,A],CC,ui1[],lb=lb1[],clb=cl[],ub=ub1[],warmstart=iter>1,uz=uzi[],lz=lzi[],lm=lmi[]);

//cout << "ADAPTMESH ITERATION " << iter << endl << "UZ = " << uzi[] << endl << "LZ = " << lzi[] << endl << "LM = " << lmi[] << endl;

plot(ui1,ui2,wait=1,nbiso=60,dim=3);
if(iter > 1) break;
Th= adaptmesh(Th,[ui1,ui2],err=0.004,nbvx=100000);
[uz,uz2]=[uzi,uzi2];
[lz,lz2]=[lzi,lzi2];
[u1,u2]=[ui1,ui2];
lm=lmi;
}
