// just a function which converts integers to strings
// and pads with zeros; for saving photos...
func string f(int i, int len){
   string ss = i;
   int nn = ss.length; 
   for(int i=0;i<len-nn;i++)
     ss = "0"+ss;
   return ss;
}

// if you set the following parameter to 1 then at each iteration 
// you'll save the image of the current shape
// BEWARE that this will fill your folder with images
int makemovie = 0; 

real r = 0.05;  // radii of the small disks
// coordinates of the initial disk centers
real[int] xx = [0.05, -0.06, 0.3, 0.6, 0.8, 0.9, -0.2, -0.5, -0.4, -0.8];
real[int] yy = [0.09, -0.1, -0.1, 0.5, -0.5, 0.7, -0.6, 0.2, -0.9, 0.3];
int iter = 1;   
real optv = 0;  // initialization maximal value
real epst = 0.01;   // step size
 
int nhole = 1;      // number of holes. If greater than 10
                    // you need to add more coordinates to xx and yy
int[int] one(nhole);
one = 1;
//verbosity = 0;

int[int] NC=-30*one;  // parameters for the mesh creation
                      // one number per small disk
// optimization loop
while (iter<200){

// border of the square
border C1(t=-1,1){x = -1; y = t; label = 1;};
border C2(t=-1,1){x = t; y = 1; label = 1;};
border C3(t=-1,1){x = 1; y = -t;label = 1;};
border C4(t=-1,1){x = -t; y = -1;label = 1;};

// multi-border for the small disks
border cc(t=0,2*pi;i){
    x=xx[i]+r*cos(t); 
    y=yy[i]+r*sin(t);
    label=i+2;
}

int numb = -20; // number of points on the sides of the square
                // sign chosen so that the mesh is on the good side
// build the mesh

mesh Th = buildmesh (C1(numb)+C2(numb)+C3(numb)+C4(numb)+
            cc(NC));

// If you want to adapt the mesh...
// Th= adaptmesh(Th,1./40,IsMetric=1,nbvx=15000);

// plot(Th);
// define the P1 finite element space defined over Th
fespace Vh(Th,P1);
// variables on the mesh
Vh u1,u2;

int[int] BC = 2:1:2+nhole; // label indexes for each disk
                           // used for the boundary condition

// Define the problem in weak form
varf a(u1,u2) = int2d(Th)
(dx(u1)*dx(u2) + dy(u1)*dy(u2))+on(1,u1=0)//on(C1,C2,C3,C4,u1=1)
                               +on(BC,u1=0);     
varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ; 
// define matrices for the eigenvalue problem
matrix A= a(Vh,Vh,solver=Crout,factorize=1); 
matrix B= b(Vh,Vh,solver=CG,eps=1e-20); 

// we are interested only in the first eigenvalue
int eigCount = 1;
real[int] ev(eigCount); // Holds eigenvalues
Vh[int] eV(eigCount);   // holds eigenfunctions
// Solve Ax=lBx
int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV);

// if the eigenfunction is negative
// then turn it positive
if(int2d(Th)(eV[0])<0){
  eV[0] = -eV[0];
}


// compute the direction of the gradient for each hole
real[int] ar1(nhole);
for(int i=0;i<nhole;i++) 
ar1[i]= -int1d(Th,i+2)((dx(eV[0])*N.x+dy(eV[0])*N.y)^2*N.x);
real[int] ar2(nhole);
for(int i=0;i<nhole;i++)
ar2[i] = -int1d(Th,i+2)((dx(eV[0])*N.x+dy(eV[0])*N.y)^2*N.y);

if(ev[0]>optv){
// if this configuration is better then advance
xx = xx+epst*ar1;
yy = yy+epst*ar2;
optv = ev[0];

if(makemovie==1){
// if makemovie then save each figure in png format
string nh    = nhole;
string photo = "holes"+nh+"_"+f(iter,4)+".png";

plot(eV[0],fill=true,wait=true,ShowAxes = 0,ps = photo);
}
else{
// else just plot the figure
 plot(eV[0],fill=true,wait=true,ShowAxes = 0,value =true );
}
cout << iter << " " << epst << " " << ev[0] << endl;
epst = 1.2*epst;
if(epst>0.01) epst = 0.01;
iter = iter+1;
}
else{
// else decrease the step size
 epst = epst*0.8;
 if(epst<1e-10)  break;
// perturb values a bit to avoid local minima
 xx = xx+epst/1000*ar1;
 yy = yy+epst/1000*ar2;
 optv = ev[0];
}

}