#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include "dumka3.cpp"
#include "euler.cpp"
#include <vector>
#include "Point2D.h"
//CC Brusselator.cpp Point2D.cpp -lm -LANG:std -64 -apo
using namespace std;
ofstream  solution( "solDumkaCPP8" ); // open output file for solution
class Stencil {
  public:
    int ID;
    int r;
    int l;
    int t;
    int b;
    Stencil() {
    }
    Stencil(int pointId, int right, int left,int top,int bottom) {
      this->ID = pointId;
      this->r=right;
      this->l=left;
      this->t=top;
      this->b=bottom;
    }
};

template <class Vector, class Element>
class VectorFunction: public VectorFunctionBase<Vector,Element> {
protected:
  double alf;
  int nx;
  int ny;
  Stencil* stencils;
  Point2D* points;
public:
  VectorFunction() {
  }
  VectorFunction ( double alf, int nx, int ny, Stencil* stencils, Point2D* points ) {
    this->alf = alf;
    this->nx = nx;
    this->ny = ny;
    this->stencils = stencils;
    this->points = points;
  }
  bool isStepAccepted;
  void compute(Vector& y, double x, Vector& result) {
    long size = y.size();
    if ( isStepAccepted ) 
    {
      //isStepAccepted = false;
      //output(solution, ns, y, x);
    }
    double radsq=0.01;
    double bet;
    if(x >= 1.1) {
      bet=5.00e0;
    } else {
      bet=0.00e0;
    };
// ----- big loop -----
    for ( int i = 0; i < size; i++ ) {
        double uleft,vleft,uright,vright,ulow,vlow,uup,vup,uij,vij;
// ----- the derivative -----
        uij       = y[i].x;
        vij       = y[i].y;
        uleft     = y[stencils[i].l].x;
        vleft     = y[stencils[i].l].y;
        uright    = y[stencils[i].r].x;
        vright    = y[stencils[i].r].y;
        ulow      = y[stencils[i].b].x;
        vlow      = y[stencils[i].b].y;
        uup       = y[stencils[i].t].x;
        vup       = y[stencils[i].t].y;
        result[i].x=1.e0+uij*uij*vij-4.4e0*uij +alf*nx*ny*(uleft+uright+ulow+uup-4.e0*uij);
        result[i].y=3.4e0*uij - uij*uij*vij +alf*nx*ny*(vleft+vright+vlow+vup-4.e0*vij);
// ----- inhomogenity -----    
        double yy=points[i].x;
        double xx=points[i].y;
        if(((xx-0.3e0)*(xx-0.3e0)+(yy-0.6e0)*(yy-0.6e0)) <= radsq) {
          result[i].x=result[i].x + bet;
          result[i].y=result[i].y + bet;
        }
      }
  };
  Element GetNonZeroElement() {     
     Element p(0.01*(rand()%100),0.01*(rand()%100));
     return p;
     //return 0.01*(rand()%100);
   }
  double cour(Vector& y, double t) {
    return 0.0340264;//3.53695e-05;//0.019;//0.038948
  };
};
template<class Vector> 
void output(ofstream& solution, int nx, int ny, Vector& y, double t) { 
   solution.setf(ios::fixed,ios::floatfield);
   solution.precision(15);
   for ( int j = 0; j < ny; j++ ) {  
//     double yy = ((double)j)/(ny); 
     for ( int i = 0; i < nx; i++ ) {
//       double xx = ((double)i)/(nx);
       int ID = j*ny+i;
//       solution << xx << " " << yy << "  " << y[ID].x << " " << y[ID].y << endl; 
       solution << y[ID].x << endl;
    }
  }
   for ( int j = 0; j < ny; j++ ) {  
     for ( int i = 0; i < nx; i++ ) {
       int ID = j*ny+i;
       solution << y[ID].y << endl;
     }
   }
}
int main ( int argc, char** argv ) {
  // y - solution of ODE 
  int nx=128; 
  int ny=128;
  
  double alf=1.0e-1;
  vector<Point2D>  y(nx*ny);
// ----- Set up configuration -----
  Stencil* stencils = new Stencil[nx*ny];
  Point2D* points = new Point2D[nx*ny];
  for ( int j = 0; j < ny; j++ ) 
  {   
    double yy=(double)j/(double)ny; 
    for ( int i=0; i < nx; i++ ) 
    {
      int ID = j*ny+i;
      stencils[ID].ID = ID;
      stencils[ID].r = (i==nx-1)?(j*ny):(ID+1);
      stencils[ID].l = (i==0)?(j*ny+(nx-1)):(ID-1);
      stencils[ID].t = (j==ny-1)?(i):(ID+ny);
      stencils[ID].b = (j==0)?(ny*(ny-1)+i):(ID-ny);

      double xx=(double)i/(double)nx; 
      points[ID].x = xx;
      points[ID].y = yy;   
    }
  }

// ----- initial values -----;
  for ( int j = 0; j < ny; j++ ) 
  {
    for ( int i=0; i < nx; i++ ) 
    {
      int ID = j*ny+i;
      y[ID].x=22.e0*points[ID].y*pow(1.e0-points[ID].y,1.5e0);
      y[ID].y=27.e0*points[ID].x*pow(1.e0-points[ID].x,1.5e0);
    }        
  }


  VectorFunction<vector<Point2D>,Point2D> f( alf, nx, ny, stencils, points );


  double tinitial=0.;
  double t=tinitial;    // t - initial time
  double tend=1.5;       // tend - interval of integration
  double rtol = 0.00000001; // rtol - relative error (should be strictly positive) 
  double atol=rtol;      // tolerance of computation(should be strictly positive)
  double h0 = 0.000001;  // h0 - initial step size
                      // f - class VectorFunction which has 
  //output(solution, nx, ny, y, t);
  dumka3(t, tend, h0, atol, rtol, f, y, true );
  //euler( t, tend, f, y );                    //     method Vector  
  output(solution, nx, ny, y, t);
  solution.close();

  return 1;
}