1load "ff-Ipopt"
2int bfgs=0,constraint =1; // aglo ..
3real x0=1.5,y0=0.5;
4mesh Th=square(10,10,[x+1,y]);
5fespace Vh(Th,P2);
6func g = cos(pi*x)*cos(pi*y);; //acosh(sqrt(x*x+y*y));//cos(pi*x)*cos(pi*y);
7Vh ue= g;
8macro grad(u) [dx(u),dy(u)]//EOM
9func real J(real[int]& X)
10{
11    Vh u; u[]=X;
12    return int2d(Th)(sqrt(1+grad(u)'*grad(u))) ;
13}
14
15func real[int]  DJ(real[int]& X)
16{
17    Vh u; u[]=X;
18
19    varf vg(uu,v) = int2d(Th)((grad(u)'*grad(v)) / sqrt(1+grad(u)'*grad(u))) ;
20   real[int] G= vg(0,Vh);
21    return G;}
22matrix H;//global vairable for Hessien matrix Overwise  => seg fault in Ipopt
23func matrix  HJ(real[int]& X)
24{
25    Vh u; u[]=X;
26    varf vH(v,w) = int2d(Th)( (grad(w)'*grad(v)) / sqrt(1+grad(u)'*grad(u))
27     - (grad(w)'*grad(u))*(grad(v)'*grad(u)) *(1+grad(u)'*grad(u))^-1.5 ) ;
28    H = vH(Vh,Vh);
29    return H;}
30
31
32varf OnGamma(u,v) = on(1,2,3,4,u=1);
33Vh OnG;
34OnG[]=OnGamma(0,Vh,tgv=1); // 1 on Gamma
35Vh lb = OnG!=0 ? g : -1e19 ; //
36Vh ub = OnG!=0 ? g :  1e19 ; //
37Vh u = OnG!=0  ? g : 0 ; //  initial guest ..
38Vh clb = 3-square(10*(square(x-x0)+square(y-y0))); //  constraint ..
39if(constraint) lb = max(lb,clb);
40int ret;
41if(bfgs)
42  ret = IPOPT(J,DJ,u[],lb=lb[],ub=ub[],bfgs=1);
43else
44  ret = IPOPT(J,DJ,HJ,u[],lb=lb[],ub=ub[]);
45
46cout << " min = "  << J(u[]) << "  ~ " << J(ue[]) << " ret = " << ret << endl;
47plot(u, wait=1);
48