// file nl-elliptic-fp.edp
//
// solution of $-a\Delta u + b\partial_1 u + cu = F(u) + f_0$
// together with $u = u_0$ on $\Gamma$
//
// fixed-point iterates
//
// $a$, $b$, $c$ are positive constants
// $F = F(s)$, $u_0 = u_0(x,y)$ are functions

//
// NUMERICAL VALUES AND FUNCTIONS TO PROVIDE
//
int C0=100, C1=99, C2=98; // the labels
real a=.1, b=5., bt2 = 2.*b, c=10.; // the values on the boundaries
real errmesh=1.e-2;

int itermax=100, itermm1=itermax-1, kiter, kiter1, nmesh=1, jiter;
real[int] Uerror(1000);
real errFP=1.e-15, error, aux;

func f00=0.;
func u00=x*x - y*y;                     // the function u_0

//
// DESCRIPTION OF THE BOUNDARY AND THE INTERFACE
//
// 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=1+t;  y=3;      label=C1;}
border C12(t=0,1){ x=2;    y=3-6*t;  label=C1;}
border C13(t=0,1){ x=2-t;  y=-3;     label=C1;}
border C14(t=0,1){ x=1;    y=-3+6*t; label=C1;}

//
// FIGURES (I)
//
// boundary and interface:
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)
//
// mesh:
plot(Th, cmm=" FIRST MESH", wait=true);

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//
fespace Vh(Th,P1); Vh u,v,fff,f0=f00,u0=u00,uold,udif;

// -------------------------------------
// THE NONLINEAR FUNCTION
// -------------------------------------
func real[int] FFF(real[int] & uvec)
{
	real[int] F(uvec.n);
	for(int i=0;i<uvec.n;i++)
	{
		real ui = uvec[i];
		F[i] = exp(1.-ui*ui)*sin(ui);
	}
	return F;
}
// -------------------------------------

// -------------------------------------
// THE LINEAR PROBLEMS
// -------------------------------------
problem nonsymk(u,v,init=jiter) = int2d(Th)(a*(dx(u)*dx(v)+dy(u)*dy(v))
        + b*dx(u)*v + c*u*v)
        - int2d(Th)(fff*v)
        + on(C0,C1, u=u0);
// -------------------------------------


//
// SOLVING THE PROBLEM
//

// STARTING

kiter=0;
fff = f0;
nonsymk;
uold  = u;
jiter = 1;

// ITERATES

for (int k=0;k<itermax;k++)
{
	kiter = k;
	fff[] = FFF(u[])+f0[];
	nonsymk;
	
	plot(u, cmm="ITERATE "+kiter, value=true, fill=true, wait=true);

	udif = u-uold; error = int2d(Th)(udif*udif);
	Uerror[kiter] = sqrt(error);
	cout << " ITER = " << kiter << " ERROR = " << Uerror[kiter] << endl;
	if (Uerror[kiter] <= errFP) break;
	
	if (kiter % 10 == 0)
	{
		nmesh = nmesh+1; jiter = -1;
		Th = adaptmesh(Th,u,err=errmesh);
		plot(Th, cmm=" MESH NO. "+nmesh, wait=true);
	};
	
	if (kiter == itermm1) cout << " ATTENTION: TOO MANY ITERATES" << endl;
	uold = u; jiter = jiter+1;
}

//
// DISPLAY, RECORD AND PRINT
//

plot(u, value=true, fill=true, cmm="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, ps="nonsymsol.eps");