//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "mpi-cmaes"


int NN = 6;
mesh Th = square(NN,NN);
func f1=1.;
func f2=-1.;
func g1=0.;
func g2=0.1;
int iter=0;
int nadapt=1;
real starttol=1e-10,bctol=6.e-12;
fespace Vh(Th,P1);
Vh ou1,ou2;
real pena=1000;


for(int al=0;al<nadapt;++al)
{
	
varf BVF(v,w) = int2d(Th)(0.5*dx(v)*dx(w) + 0.5*dy(v)*dy(w));
varf LVF1(v,w) = int2d(Th)(f1*w);
varf LVF2(v,w) = int2d(Th)(f2*w);
matrix A =  BVF(Vh,Vh);
real[int] b1 = LVF1(0,Vh) , b2 = LVF2(0,Vh); 

varf Vbord(v,w) = on(1,2,3,4,v=1);
//real[int] bord = Vbord(0,Vh);
//real[int] in = bord ? 0 : 1;
Vh In,Bord;
Bord[] = Vbord(0,Vh,tgv=1);
In[] = Bord[] ? 0:1;
Vh gh1=Bord*g1,gh2=Bord*g2;



//Function which create a vector of the search space type from
//two finite element functions
func int FEFToSSP(real[int] &fef1,real[int] &fef2,real[int] &ssp)
{
	int kX=0;
	for(int i=0;i<Vh.ndof;++i)
	{
		if(In[][i])
		{
			ssp[kX] = fef1[i];
			ssp[kX+In[].sum] = fef2[i];
			++kX;
		}
	}
	return 1;
}
//Function spliting a vector from the search space and fills
//two finite element functions with it
func int SSPToFEF(real[int] &fef1,real[int] &fef2,real[int] &ssp)
{
	int kX=0;
	for(int i=0;i<Vh.ndof;++i)
	{
		if(In[][i])
		{
			fef1[i] = ssp[kX];
			fef2[i] = ssp[kX+In[].sum];
			++kX;
		}
		else
		{
			fef1[i] = gh1[][i];
			fef2[i] = gh2[][i];
		}
	}
	return 1;
}


func real IneqC(real[int] &X)
{
	real[int] constraints(In[].sum);
	for(int i=0;i<In[].sum;++i)
	{
		constraints[i] = X[i] - X[i+In[].sum];
		constraints[i] = constraints[i] <= 0 ? 0. : constraints[i];
	}
	return constraints.l2;
}


func real J(real[int] &X)
{
	Vh u1,u2;
	SSPToFEF(u1[],u2[],X);
	iter++;
	real[int] Au1 = A*u1[], Au2 = A*u2[];
	Au1 -= b1;
	Au2 -= b2;
	real val = u1[]'*Au1 + u2[]'*Au2;
	val +=  pena * IneqC(X);
	if(iter%100==0) plot(u1,u2,nbiso=30,fill=1,dim=3,cmm="adapt level " + al + " - iteration " + iter + " - J = " + val,value=1);
	return val ;
}



real[int] start(2*In[].sum);

if(al==0)
{
	start(0:In[].sum-1) = 0.;
	start(In[].sum:2*In[].sum-1) = 0.1;
}
else FEFToSSP(ou1[],ou2[],start);


real mini = cmaesMPI(J,start,stopMaxFunEval=10000*(al+1),stopTolX=1.e-4/(10*(al+1)),initialStdDev=(0.025/(pow(100.,al))));
Vh best1,best2;
SSPToFEF(best1[],best2[],start);

Th = adaptmesh(Th,best1,best2);
ou1 = best1;
ou2 = best2;

}
