// file CC-CG-OS.edp
//
// OPTIMAL CONTROL OF A CAPACITOR -- CONJUGATE GRADIENT WITH OPTIMAL STEP (AND PROJECTION)
//

verbosity=0;

//
// NUMERICAL VALUES TO PROVIDE
//
int C0=100, C1=99, C2=98; // the labels
int itermax=500;

real umin=0., umax=2.;
real acoef=.99, one=1., yd0=1.; // good values of acoef: .2, .3, .4
real gamma, rhofix=1.;
real[int] Uerror(1000);
real errGO=1.e-5, errmesh=1.e-1;

//
// AUXILIAR VARIABLES
//
int kiter=0, jiter=0, nmesh=1, itermm1=itermax-1;

real bcoef=1.-acoef, abc=rhofix*one, rho=rhofix, error, aux0, aux1, aux2,
    xnor2, xnorn2;
// real abcoef=acoef/bcoef;

//
// DESCRIPTION OF THE BOUNDARIES
//
// outer boundary ($\Gamma_0$):
border C00(t=0,2*pi){x=5*cos(t); y=5*sin(t); label=C0;}

// inner boundary ($\Gamma_1$):
border C11(t=0.,1.){ x=2.-t;  y=3.;      label=C1;}
border C12(t=0.,1.){ x=1.;    y=3.-6.*t;  label=C1;}
border C13(t=0.,1.){ x=1.+t;  y=-3.;     label=C1;}
border C14(t=0.,1.){ x=2.;    y=-3.+6.*t; label=C1;}

//
// FIGURES (I)
//
// boundaries:
plot (                C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20),
                    wait=true);

//
// MESH
//
mesh Th=buildmesh(    C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20));

//
// FIGURES (II)
//
// initial mesh:
plot(Th, cmm=" FIRST MESH", wait=true);

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//
//fespace Vh(Th,P1);
fespace Vh(Th,P2);
// Vh u=0.,v,fff,f0=f00,u0=u00,uold,udif;
Vh ystate, pstate, zeta, sol, phi, yd;
Vh un=0., u, des, desn, gr, grn, udif, fff;

int iref0=Th(2.1,3.1).region;
int iref1=Th(1.9,2.9).region;

//
// SOME FUNCTIONS
//
func unom=one*(region==iref1);
func abom=abc*(region==iref1);
func ydfun=yd0*(x*x-y*y);

yd=ydfun;
plot(yd, cmm="TARGET STATE", value=true,fill=true, wait=true);

// -------------------------------------
// THE SYSTEMS
// -------------------------------------
problem states(sol,phi,init=jiter)=int2d(Th)(dx(sol)*dx(phi)+dy(sol)*dy(phi))
        - int2d(Th)(fff*phi)
        + on(C0,sol=0.);
// -------------------------------------


// -------------------------------------
// THE FIRST PROJECTOR
// -------------------------------------
// gvec=max(fvec,umin)
//
func real[int] PUad(real[int] & fvec)
{
	real[int] gvec(fvec.n);
	for (int ifv=0;ifv<fvec.n;ifv++)
	{
		gvec[ifv]=max(fvec[ifv],umin);
	}
	
	return gvec;
}
// -------------------------------------

/*
// -------------------------------------
// THE SECOND PROJECTOR
// -------------------------------------
// gvec=max(min(fvec,umax),umin)
//
func real[int] PUad(real[int] & fvec)
{
	real[int] gvec(fvec.n);
	for (int ifv=0;ifv<fvec.n;ifv++)
	{
		gvec[ifv]=min(max(fvec[ifv],umin),umax);
	}
	
	return gvec;
}
// -------------------------------------
*/

//
// SOLVING THE PROBLEM
//

/*
fff=yd;
states;
yd=sol;
plot(yd, cmm="TARGET STATE", value=true,fill=true, wait=true);
*/

plot(un, value=true, fill=true, cmm="STARTING CONTROL", wait=1);

kiter=0;
fff=un;
states;
ystate=sol;
plot(ystate, cmm="STATE - ITERATE "+kiter, value=true,
    fill=true, wait=true);

jiter=1;	
fff=ystate-yd;
states;
pstate=sol;
plot(pstate, cmm="ADJOINT STATE - ITERATE "+kiter, value=true,
    fill=true, wait=true);
	
des=unom*(acoef*pstate+bcoef*un);
fff=des;
states;
zeta=sol;
aux1=int2d(Th)(des*des);
aux2=int2d(Th)(zeta*zeta);
rho=aux1/(acoef*aux2+bcoef*aux1);
		
u=un-rho*des;
//	u[]=PUad(u[]);
plot(u, cmm="CONTROL - ITERATE "+kiter, value=true,
    fill=true, wait=true);

un=u;
desn=des;
grn=des;
xnorn2=int2d(Th)(grn*grn);

// ITERATES

for (int k=1;k<itermax;k++)
{
	kiter=k;
	
	fff=un;
	states;
	ystate=sol;
	plot(ystate, cmm="STATE - ITERATE "+kiter, value=true,
	    fill=true, wait=true);
	
	jiter=1;
	fff=ystate-yd;
	states;
	pstate=sol;
	plot(pstate, cmm="ADJOINT STATE - ITERATE "+kiter, value=true,
	    fill=true, wait=true);
	
	gr=unom*(acoef*pstate+bcoef*un);
	xnor2=int2d(Th)(gr*gr);
	gamma=xnor2/xnorn2;
	
	des=gr+gamma*desn;
	
// ????
	fff=des;
	states;
	zeta=sol;
	aux0=int2d(Th)(gr*des);
	aux1=int2d(Th)(des*des);
	aux2=int2d(Th)(zeta*zeta);
	rho=aux0/(acoef*aux2+bcoef*aux1);
// ????

	u=un-rho*des;
//	u[]=PUad(u[]);
	plot(u, cmm="CONTROL - ITERATE "+kiter, value=true,
	    fill=true, wait=true);

	udif=u-un;
	error=udif[].l2;
	Uerror[kiter]=error;
	cout << " ITER=" << kiter << " ERROR=" << error << endl;
	if (error <= errGO) break;
/*	
	if ((kiter+1) % 100 == 0)
	{
		nmesh=nmesh+1;
		jiter=0;
		Th=adaptmesh(Th,u,err=errmesh);
		plot(Th, cmm=" MESH NO. "+nmesh, wait=true);
		u=u;
	};
*/
	if (kiter == itermm1) cout << " ATTENTION: TOO MANY ITERATES" << endl;
	un=u;
	desn=des;
	xnorn2=xnor2;
}

//
// DISPLAY, RECORD AND PRINT
//

plot(ystate, value=true, fill=true, cmm="STATE - FINAL ITERATE",
    wait=1);
plot(pstate, value=true, fill=true, cmm="ADJOINT STATE - FINAL ITERATE",
    wait=1);
plot(u, value=true, fill=true, cmm="CONTROL - FINAL ITERATE",
    wait=1);

// WRITING

/*
{ // save solution
	ofstream file("sol.txt");
	file << u[] << endl;
} // close the file (end block)
// { // read
// ifstream file("sol.txt");
// file >> u[] ;
// } // close reading file (end block)

kiter1=kiter+1;
for (int k=0;k<kiter1;k++)
{
	cout << " FIXED-POINT ITERATE AND ERROR: " << k << " - "<< Uerror[k] << endl;
}

if (kiter > 1)
{
	aux=log(Uerror[kiter-1]/Uerror[kiter])/log(2.);
    cout << " FIXED-POINT RATE=            " << aux << endl;
}

//
// FIGURES (III)
//
// the solution:
plot(u, wait=true, value=true, fill=true);
// plot(u, wait=true, value=true, fill=true, ps="nonsymsol.eps");
*/