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