1load "ff-Ipopt";
2int nn=10;
3mesh Th=square(nn,nn);
4fespace Vh(Th,[P1,P1] );
5fespace Wh(Th,[P1] );
6int iter=0;
7
8func f1 =  10;//right hand sides
9func f2 = -15;
10func g1 = -0.1;//Boundary conditions functions
11func g2 =  0.1;
12
13Vh [uz,uz2]=[1,1],[lz,lz2]=[1,1];
14Wh lm=1.;
15Vh [u1,u2]=[0,0];//starting point
16
17
18while(++iter)
19{
20macro Grad(u) [dx(u),dy(u)]//gradient macro
21varf vP([u1,u2],[v1,v2]) = int2d(Th)(Grad(u1)'*Grad(v1)+ Grad(u2)'*Grad(v2))
22- int2d(Th)(f1*v1+f2*v2);
23
24matrix A = vP(Vh,Vh);//Fitness function matrix...
25real[int] b = vP(0,Vh);//and linear form
26
27int[int] II1=[0],II2=[1];//Constraints matrix
28matrix C1 =  interpolate (Wh,Vh, U2Vc=II1);
29matrix C2 =  interpolate (Wh,Vh, U2Vc=II2);
30matrix CC = -1*C1 + C2; // u2 - u1 >0
31Wh cl=0;//constraints lower bounds (no upper bounds)
32
33//Boundary conditions
34varf vGamma([u1,u2],[v1,v2]) = on(1,2,3,4,u1=1,u2=1);
35real[int] onGamma=vGamma(0,Vh);
36Vh [ub1,ub2]=[g1,g2];
37Vh [lb1,lb2]=[g1,g2];
38ub1[] = onGamma ? ub1[] : 1e19  ; //Unbounded in interior
39lb1[] = onGamma ? lb1[] : -1e19  ;
40
41
42Vh [uzi,uzi2]=[uz,uz2],[lzi,lzi2]=[lz,lz2];
43Wh lmi=lm;
44Vh [ui1,ui2]=[u1,u2];
45
46IPOPT([b,A],CC,ui1[],lb=lb1[],clb=cl[],ub=ub1[],warmstart=iter>1,uz=uzi[],lz=lzi[],lm=lmi[]);
47
48//cout << "ADAPTMESH ITERATION " << iter << endl << "UZ = " << uzi[] << endl << "LZ = " << lzi[] << endl << "LM = " << lmi[] << endl;
49
50plot(ui1,ui2,wait=1,nbiso=60,dim=3);
51if(iter > 1) break;
52Th= adaptmesh(Th,[ui1,ui2],err=0.004,nbvx=100000);
53[uz,uz2]=[uzi,uzi2];
54[lz,lz2]=[lzi,lzi2];
55[u1,u2]=[ui1,ui2];
56lm=lmi;
57}
58