//(FreeFem++)
//The following code allows to see some R-D equations in a rectangular domain of R2 for the reaction f(u)=u^(1+p)(1-u)

real T = 25;	//Final time****
real dt = 0.1;	//Time step****

int J=2;	//Number of adjustments of the mesh (>0)

int n = 100;	//Number of subdivisions****

real H = 3;	//Height****
real W = 3;	//Width****

real d = 0.01;  //Diffusion intensity****

//FUJITA'S EXPONENT
// p Fujita = pF = 2/N = 2/2 = 1
// p<pF => systematic HTE
// p>pF => non-systematic HTE
// p=0  => logistic
real p = 0;	//Allee Effect intensity****

//COMPACTLY SUPPORTED INITIAL DATUM (Indicator function of a disc)
real x0 = 0;	//x-coordinate of the disc****
real y0 = H/2;	//y-coordinate of the disc****
real radius = 0.2;	//Radius of the disc****
real u0max = 0.7;	//Maximum of the initial datum (<1)
func InitialDatum = u0max*0.5*(1+sign(radius - sqrt((x-x0)^2+(y-y0)^2)));	//Building of the initial datum

//COMPUTE TO GET A REGULAR MESH
real ratio = H/W;
int m = n*W/H;

//BUILDING OF THE RECTANGULAR MESH
border bottom(t = 0,1) {label = 1 ; x = -W/2 + t*W ; y = 0 ;};
border right(t = 0,1) {label = 1 ; x = W/2 ; y = 0 + t*H ;};
border top(t = 0,1) {label = 1 ; x = W/2 - t*W ; y = H ;};
border left(t = 0,1) {label = 1 ; x = -W/2 ; y = H - t*H ;};
mesh Th=buildmesh(bottom(m)+right(n)+top(m)+left(n));

//plot(Th, wait=false);	//Display the mesh

//DEFINITION OF FUNCTIONAL SPACES AND FUNCTIONS
fespace Vh(Th,P1);
Vh u0=InitialDatum;
Vh u=u0, v, uold;

//PLOT THE INITIAL DATUM
real[int] colorhsv=[	// To plot in black & white
	0, 0 , 0,			// min is in black (note min is here zero)
	0, 0 , u0max		// u0max is in grey (1 is in white)
	];
	
cout << endl<< endl <<"Norm L1 of u0 = " << int2d(Th)(u0) << endl<< endl<< endl;	//Display the L1-norm of u0

plot(u, nbiso=255, value = 0, dim = 2, fill = 1, wait=1, hsv=colorhsv, cmm="Press ENTER");

//DEFINE THE WEAK PROBLEM
problem RD(u, v, solver=UMFPACK)	//v is the test function
	= int2d(Th)(u*v/dt) - int2d(Th)(uold*v/dt)	//Estimate the time derivative

	+ int2d(Th)(d* (dx(u)*dx(v)+dy(u)*dy(v)))	//Diffusion

	- int2d(Th)(uold^(1+p)*(1-uold)*v)	//Reaction

    //+ on(1, u=0)		//Dirichlet on the frontier boundaries (if commented : Neumann)
;

{
ofstream savemin("min.dat");	//To save min and max of the solution for each time
ofstream savemax("max.dat");

//SOLVING THE PROBLEM AND SAVE MIN AND MAX
for(real t = 0; t <= T; t += dt){
	uold = u;	//ie. u^{n-1} = u^n

	for (int j=1;j<J+1;j++)
	{
	RD;		//Solve the problem RD
	Th=adaptmesh(Th,u);		//Adapt the mesh to the function u
	//plot(Th,wait=1);		//To see how adaptmesh works...
	}

    uold=abs(uold); //To get rid of potential negative values due to approximation
    u=abs(u);		//Idem

	real umax=uold[].max;  //Get the min and the max of the solution
	real umin=uold[].min;

	cout << "min = " << umin << endl;
	cout << "max = " << umax << endl;
    
    //real uminRounded = round(abs(umin)*100)/100;	//The rounded min and max
    //real umaxRounded = round(abs(umax)*100)/100;

	real[int] colorhsv=[	// To plot in black & white
		0, 0 , umin,	// min is in black
		0, 0 , umax		// max is in white
	];

    plot(u, nbiso=255, value = 0, cmm="Time="+t+"          min="+umin+"          max="+umax, dim = 2, fill = 1, wait=0, hsv=colorhsv);	//Plot the solution

    savemin << t << "	" << umin << endl;	//save the solution's min
    savemax << t << "	" << umax << endl;	//save the solution's max
}//END SOLVING PROBLEM

}//END OFSTREAM