load "ff-Ipopt"
int bfgs=0,constraint =1; // aglo .. 
real x0=1.5,y0=0.5; 
mesh Th=square(10,10,[x+1,y]);
fespace Vh(Th,P2);
func g = cos(pi*x)*cos(pi*y);; //acosh(sqrt(x*x+y*y));//cos(pi*x)*cos(pi*y);
Vh ue= g;
macro grad(u) [dx(u),dy(u)]//EOM
func real J(real[int]& X)
{
    Vh u; u[]=X;
    return int2d(Th)(sqrt(1+grad(u)'*grad(u))) ;
}

func real[int]  DJ(real[int]& X)
{
    Vh u; u[]=X;
   
    varf vg(uu,v) = int2d(Th)((grad(u)'*grad(v)) / sqrt(1+grad(u)'*grad(u))) ;
   real[int] G= vg(0,Vh); 
    return G;}
matrix H;//global vairable for Hessien matrix Overwise  => seg fault in Ipopt     
func matrix  HJ(real[int]& X)
{
    Vh u; u[]=X;   
    varf vH(v,w) = int2d(Th)( (grad(w)'*grad(v)) / sqrt(1+grad(u)'*grad(u)) 
     - (grad(w)'*grad(u))*(grad(v)'*grad(u)) *(1+grad(u)'*grad(u))^-1.5 ) ;
    H = vH(Vh,Vh); 
    return H;}  
       
    
varf OnGamma(u,v) = on(1,2,3,4,u=1); 
Vh OnG;
OnG[]=OnGamma(0,Vh,tgv=1); // 1 on Gamma 
Vh lb = OnG!=0 ? g : -1e19 ; // 
Vh ub = OnG!=0 ? g :  1e19 ; //  
Vh u = OnG!=0  ? g : 0 ; //  initial guest ..
Vh clb = 3-square(10*(square(x-x0)+square(y-y0))); //  constraint .. 
if(constraint) lb = max(lb,clb); 
int ret;
if(bfgs) 
  ret = IPOPT(J,DJ,u[],lb=lb[],ub=ub[],bfgs=1); 
else 
  ret = IPOPT(J,DJ,HJ,u[],lb=lb[],ub=ub[]); 

cout << " min = "  << J(u[]) << "  ~ " << J(ue[]) << " ret = " << ret << endl; 
plot(u, wait=1);
