//(FreeFem++)
//The following code allows to see some R-D equations on the Road-Field space R^{2}_{+} for the reaction f(u)=u^(1+p)(1-u)

verbosity=0;	//To not show the text in the terminal...

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

int n = 30;	//Number of vertical subdivisions****

real H = 40.;		//Height****
real W = 200.;	//Width****

//PARAMETERS OF THE ROAD-FIELD SPACE
real d = 1.;  //Diffusion intensity in the Field****
real D = 10.;  //Diffusion intensity on the Road****
real mu = 1.;  //Migratory equilibrium between Road and Field****
real nu = 1.;  //Migratory equilibrium between Road and Field****

//REACTION FUNCTION
// 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 = 2;	//Allee Effect intensity****

macro reactionF(xx,vv)vv*(1-vv)//OEM reaction function (logistic (KPP))
//macro reactionF(xx,vv)vv*(sin(xx)-vv)//OEM reaction function (logistic (KPP))
//macro reactionF(xx,vv)vv^(1+p)*(1-vv)//OEM reaction function (monostable degenerate (Allee Effect))

//COMPACTLY SUPPORTED INITIAL DATUM ON FIELD (Indicator functions)

//Indicator functions of discs
macro indicatriceDF(x0,y0,radius)0.5*(1+sign(radius - sqrt((x-x0)^2+(y-y0)^2)))//OEM indicator function of the ball B(x0,y0,radius)
macro indicatriceDR(x0,radius)0.5*(1+sign(radius - abs(x-x0)))//OEM indicator function of the ball B(x0,radius)
real x0 = 0.;	//x-coordinate of the disc center****
real y0 = 0;	//y-coordinate of the disc center****
real radius = 10;	//Radius of the disc****
real v0max = 0.9;	//Maximum of the initial datum v0****
real u0max = 0.9;	//Maximum of the initial datum u0****
func InitialDatumField = v0max*indicatriceDF(x0,y0,radius);  //Building of the initial datum v0
func InitialDatumRoad = u0max*indicatriceDR(x0,radius);  	//Building of the initial datum u0

/*
//Indicator functions of rectangles
macro indicatriceX(a,b)0.5*(1+sign(0.5*abs(a-b) - abs(x-0.5*(a+b))))//OEM indicator Field-function of [a,b]²
macro indicatriceY(a,b)0.5*(1+sign(0.5*abs(a-b) - abs(y-0.5*(a+b))))//OEM indicator Road-function of [a,b]
real v0max = 0.8;	//Maximum of the initial datum u0****
real u0max = 0.8;	//Maximum of the initial datum v0****
func InitialDatumField = v0max*indicatriceX(-4,4)*indicatriceY(-4,4);	//Building of the initial datum v0
func InitialDatumRoad = u0max*indicatriceX(-4,4);				//Building of the initial datum u0
*/

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

mesh ThF=square(m, n, [-W/2+W*x, H*y]);		//To build a rectangular mesh for the Field
mesh ThR=square(m, 1, [-W/2+W*x, H*y/3]);		//To build a rectangular mesh for the Road

mesh ThRbis=square(m, 1, [-W/2+W*x, H*y/3-H/2]);	//To build the same mesh as ThR shifted downward (for ploting!)

plot(ThRbis,ThF, wait=true, cmm="(Meshes)		Press ENTER");	//Display the meshes of the Field and the Road

//DEFINITION OF FUNCTIONAL SPACES AND FUNCTIONS
fespace VhF(ThF,P1);		//Field functions 
VhF v0=InitialDatumField;	//Initial datum
VhF v=v0, vold=v0;		//Unknowns
VhF wF;				//Test function

fespace VhR(ThR,P1,periodic=[[1,x],[3,x]]);		//Road functions 
VhR u0=InitialDatumRoad;					//Initial datum
VhR u=u0, uold=u0;						//Unknowns
VhR wR;								//Test function

fespace VhRbis(ThRbis,P1,periodic=[[1,x],[3,x]]);	//Road solution shifted downward  (for ploting!)
VhRbis ubis=u(x,y+H/2);

//PLOT THE INITIAL DATUM
real[int] colorhsv=[	// To plot in black(0) & white(1)
	0, 0 , 0,			// min is in black (note min is here zero)
	0, 0 , v0max		// u0max is in grey (1 is in white)
	];
plot(ubis,v, nbiso=255, value = 0, dim = 2, fill = 1, wait=1, hsv=colorhsv, cmm="(Initial data)		Press ENTER");

//DEFINE THE WEAK PROBLEMS (FIELD & ROAD)
problem field(v, wF)	//wF is the test function
	= int2d(ThF)(v*wF/dt)-int2d(ThF)(vold*wF/dt)	//Estimate the time derivative

	+ int2d(ThF)(d* (dx(v)*dx(wF)+dy(v)*dy(wF)))	//Diffusion

	+ int1d(ThF, 1)(nu*v*wF)	//Migration Field>>Road
	- int1d(ThF, 1)(mu*u*wF)	//Migration Road>>Field

	- int2d(ThF)(reactionF(x,vold)*wF)	//Reaction

	//+ on(2,3,4, v=0)		//Dirichlet on the frontier boundaries (if commented : Neumann)
;

//Note the u is considered as a function whom the domain is the same as v but it does not matter because we only consider u at the bottom frontier...
problem road(u, wR)		//wR is the test function
	= int2d(ThR)(u*wR/dt)-int2d(ThR)(uold*wR/dt)	//Estimate the time derivative

	+ int2d(ThR)(D* (dx(u)*dx(wR)))	//Diffusion

	+ int2d(ThR)(mu*u*wR)	//Migration Road>>Field
	- int2d(ThR)(nu*v(x,0)*wR)		//Migration Field>>Road

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

{
ofstream usavemin("umin.dat");
ofstream usavemax("umax.dat");
ofstream vsavemin("vmin.dat");
ofstream vsavemax("vmax.dat");

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

	vold = v;	//ie. u^{n-1} = u^n
	field;	//Solve field problem

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

	real vmax=vold[].max;  //Get the min and the max of v
	real vmin=vold[].min;
	real vminRounded = round(abs(vmin)*100)/100;	//The rounded min and max of v
	real vmaxRounded = round(abs(vmax)*100)/100;

	real[int] colorhsv=[	// To plot in black(0) & white(1) (assert that (u0,v0)<1...)
		0, 0 , vmin,
		0, 0 , vmax
	];

	real umax=uold[].max;  //Get the min and the max of u
	real umin=uold[].min;
	real uminRounded = round(abs(umin)*100)/100;	//The rounded min and max of u
	real umaxRounded = round(abs(umax)*100)/100;

	ubis=u(x,y+H/2);
	plot(ubis,v, nbiso=255, cmm="(Time = "+t+" on "+T+")       (Min = "+uminRounded+" | "+vminRounded+")     (Max = "+umaxRounded+"  |  "+vmaxRounded+")       (Diffusion : (d,D) = ("+d+","+D+"))       (Migration equilibrium : (mu,nu) = ("+mu+","+nu+"))", value = 0, dim = 2, fill = 1, wait=0, hsv=colorhsv);

	usavemin << t << "	" << umin << endl;	//save the solution's min
	usavemax << t << "	" << umax << endl;	//save the solution's max
	vsavemin << t << "	" << vmin << endl;	//save the solution's min
	vsavemax << t << "	" << vmax << endl;	//save the solution's max
}//END SOLVING PROBLEMS

}//END OFSTREAM