// ------------------------------------------------------------
// main.cc
//
// ------------------------------------------------------------


#include "source/ugblock.h"

using namespace ::_COLSAMM_;



          /* Poisson on a cylinder */


int main()
{
  //  cout.precision(10);
  //  cout.setf(ios::fixed,ios::floatfield);
  int iteration, i, n;
  double delta, delta_prime, beta, tau, eps;              

  n = 10;                                      // number of grid points
  iteration = 100;
  eps = 0.0000000001;

  cout << "\n Solving Poisson's equation \n";
  cout << " with Dirichlet boundary conditions on a cylinder by cg: ";
  cout << "\n Solving Poisson's equation on a cylinder by cg: ";
  cout << "\n Number of grid points in each direction: " << n;
  cout << "\n Maximal number of iterations: " << iteration;
  cout << "\n Stopping parameter for cg: " << eps << endl;

  // Part 1: Construction of a domain and markers:
  //-----------------------------------------------
  Cylinder cylinder(2.0,1.0,4.0);

  Boundary_Marker boundary(&cylinder);
  Unstructured_grid_Marker inner_points(&cylinder);
  inner_points.complement(&boundary);


  // Part 2: Construction of a grid:
  //--------------------------------
  Blockgrid grid(&cylinder,n);                    // Construction of grid

  // Part 3: Definition Variables and other object which need storage:
  //------------------------------------------------------------------
  Variable<double> f(grid);
  Variable<double> g(grid);
  Variable<double> d(grid);               

  Variable<double> r(grid);
  Variable<double> e(grid);
  Variable<double> u(grid);
  Variable<double> u_exakt(grid);


  // Part 4: The mathematical code:
  //-------------------------------
  Function1<double> Sin(sin);                   // definition of functions
  Function1<double> Sinh(sinh);

  X_coordinate X(grid);
  Y_coordinate Y(grid);
  Z_coordinate Z(grid);

  Local_stiffness_matrix<double> Laplace(grid);
  Local_stiffness_matrix<double> Helm(grid);
		  
  Laplace.Calculate(grad(v_())*grad(w_()));
  Helm.Calculate(v_()*w_());

  u_exakt = Sin(X)*Sinh(Y)*Z;     // exact solution
  
  u = 0.0;
  u = u_exakt | boundary;         // setting Dirichlet boundary values
  
  r = 0.0;                        // right hand side is zero
  f = Helm(f)                          | inner_points;
  
  cout << "\n cg: " << endl;
  cout << "-------------------- " << endl;
  
  r = Laplace(u) - f                   | inner_points;
  d = -1.0 * r                         | inner_points;
  
  delta = product(r,r,inner_points);
  
  for(i=1;i<=iteration && delta > eps;++i) {
    cout << i << ". ";
    
    g = Laplace(d)                     | inner_points;
    
    tau = delta / product(d,g,inner_points);
    
    r = r + tau * g                    | inner_points;
    u = u + tau * d                    | inner_points;

    delta_prime = product(r,r,inner_points);
    beta = delta_prime / delta;   
    delta = delta_prime;
    d = beta*d - r                     | inner_points;
    e = u - u_exakt;

    cout << "Iteration: Error: Max:" << L_infty(e) << endl;
  }

  cout << "\n end. " << endl;
}                                                 






    /*


#define U u[level] 
#define F f[level] 
#define R r[level] 
#define Laplace laplace[level]

double func (double x){
  return ABS(x-0.5);
}

double sinus (double x){
  return sin(x);
}

double cosinus (double x){
  return cos(x);
}

double sinush (double x){
  return sinh(x);
}



template <typename DTyp>
//void Gauss_Seidel(Variable<DTyp>& u, Variable<DTyp>& f, int iter,
//		  Local_stiffness_matrix<DTyp>&  lsm,
void Gauss_Seidel(Variable<DTyp>& u, Variable<DTyp>& f, int iter,
		  Local_stiffness_matrix<DTyp>&  lsm,
		  Marker& marker) {

  for(int i=0;i<iter;++i) {
    u = u - (lsm(u) -f) / lsm.diag() | marker;
  }
}

template <typename DTyp>
void MG( MG_array<Variable<DTyp> > u, 
	 MG_array<Variable<DTyp> > f,
	 MG_array<Variable<DTyp> > r,
	 int level, int max_level, int iter,
	 MG_array<Local_stiffness_matrix<DTyp> >  lsm,
	 MG_array<Multigrid_operator > MGop,
	 Marker& marker) {
  if(level > 0) {
  //  if(false) {
    // if(level > max_level-1) {
    // pre-smoothing
    Gauss_Seidel(U,F,iter,lsm[level],marker);

    // restriction
    R = F - lsm[level](U) | marker;


    f[level-1] = MGop[level-1].Restrict(r[level]);

    
    level = level-1;

    U = 0.0;

    // coarse-grid correction
    MG(u,f,r,level,max_level,iter,lsm,MGop,marker);

    // prolongation
    level = level+1;

    r[level] = MGop[level-1].Prolongate(u[level-1]);

    U = U + R;

    // post-smoothing
    Gauss_Seidel(U,F,iter,lsm[level],marker);
  }
  else {
    // smoothing on coarsest grid
    Gauss_Seidel(U,F,iter,lsm[level],marker);
  }
}



using namespace ::_COLSAMM_;

int main(int argc, char** argv)
{
  cout.precision(10);
  cout.setf(std::ios::fixed,std::ios::floatfield);

  std::ifstream PARAMETER;
  int iteration, i, n, m, m_max, iter_relax;
  int level;

  double normi;  

  ofstream DATEI;

  PARAMETER.open("para.dat",std::ios::in);
  PARAMETER >> n >> iteration >> iter_relax;
  PARAMETER.close();

  cout << "\n Solving Poisson's equation on a ball \n";
  cout << " with Dirichlet boundary conditions by multigrid: ";
  cout << "\n Depth of grid: " << n;
  cout << "\n Maximal number of iterations: " << iteration;
  cout << "\n Number of smoothing steps: " 
       << iter_relax << endl;     
  std::complex<double> I(0.0,1.0);

  //Hexahedron2_test_mod hexa(1.0, 1.0, 1.0, 1.0);
  //Hexahedron2_test hexa(2.0, 1.0, 2.0);
  //Hexahedron8 hexa(1.0,1.0,1.0);
  //Hexahedron4_periodic hexa(1.0,1.0,1.0);


  Hexahedron hexa(1.0, 1.0, 1.0);

  Boundary_Marker boundary(&hexa);
  Unstructured_grid_Marker inner_points(&hexa);
  inner_points.complement(&boundary);

  DATEI.open("hex.dx",std::ios :: out);
  hexa.Print_Dx(&DATEI);
  DATEI.close();

  DATEI.open("hex.dat",std::ios :: out);
  hexa.Print_AVS(&DATEI);
  DATEI.close();

  cout << " Die Anzahl der Freiheitsgrade sind "
       << hexa.degree_of_freedom() << endl;


  //  block_grids
  ////////////////

  MG_array<Blockgrid> blockgrids(n);

  m_max = 2;
  for(i=1;i<n;++i) m_max = m_max * 2;

  m = 2; 
  for(i=0;i<n;++i) {
    blockgrids(i) = new Blockgrid(&hexa,m,m,m_max);
    //blockgrids(i) = new Blockgrid(&hexa,m_max,m_max,m);
    //blockgrids(i) = new Blockgrid(&hexa,m);
    m = m * 2;
  }
  m = m/2;

  cout << " Blockgrids_hexa erstellt! " << endl;


  
  //  variables
  ////////////////
  MG_array<Variable<double> > u(n);
  MG_array<Variable<double> > r(n);
  MG_array<Variable<double> > f(n);


  for(i=0;i<n;++i) {
    u(i) = new Variable<double>(blockgrids[i]);
    r(i) = new Variable<double>(blockgrids[i]);
    f(i) = new Variable<double>(blockgrids[i]);

    u[i] = 0.0;
    f[i] = 0.0;
    r[i] = 0.0;
  }

  
  Variable<double> U_exakt(blockgrids[n-1]);


  X_coordinate X(blockgrids[n-1]);
  Y_coordinate Y(blockgrids[n-1]);
  Z_coordinate Z(blockgrids[n-1]);

  Function1 Func(func);
  Function1 Sin(sinus);
  Function1 Cos(cosinus);
  Function1 Sinh(sinush);


  //  multigrid
  ////////////////
  MG_array<Multigrid_operator> MGop(n-1);

  for(i=0;i<n-1;++i) {
    MGop(i) = new Multigrid_operator(blockgrids[i],
				     blockgrids[i+1]);
  }


  //  local stiffness matrices
  ////////////////
  MG_array<Local_stiffness_matrix<double> > laplace(n);
  MG_array<Local_stiffness_matrix<double> > helm(n);

  for(i=0;i<n;++i) {
    laplace(i) = new Local_stiffness_matrix<double>(blockgrids[i]);
    helm(i)    = new Local_stiffness_matrix<double>(blockgrids[i]);

    //laplace[i].Calculate(grad(v_())*grad(w_()));
    //helm[i].Calculate(v_()*w_());
  }


  laplace[n-1].Calculate(grad(v_())*grad(w_()));
  helm[n-1].Calculate(v_()*w_());
  for(i=n-2;i>=0;--i) {
    laplace[i]=MGop[i].Restrict(laplace[i+1]);
  }


  Local_stiffness_matrix<double> system(blockgrids[n-1]);
  system.Calculate(grad(v_())*grad(w_()));


  //  start
  ////////////////
  level = n-1;

  u[level] = 1.0 | inner_points;
  //normi = product(U,U,inner_points);
  normi = product(u[level],u[level]);


  cout << "\n Solving equation! \n";

  U_exakt = Sin( X * M_PI ) * Sinh( Y * M_PI);   // exact solution
  //U_exakt = X;   // exact solution

  u[level] = U_exakt | boundary;
  u[level] = 0.0     | inner_points;

  f[level] = 0.0;
  r[level] = 0.0;


  
  cout << "\n MG: " << endl;
  cout << "-------------------- " << endl;

  for(i=1;i<=iteration;++i) {
    MG(u,f,r,level,level,iter_relax,laplace,MGop,inner_points);
    cout << "Iteration: " << i << "  Error: Max:" 
	 << L_infty(u[level]-U_exakt,inner_points) << endl;
    //    cout << ",  L2: " << sqrt(product(r,r)/normi) << endl;
  }

  DATEI.open("u_exakt.dx",std::ios :: out);
  U_exakt.Print_Dx(&DATEI);
  DATEI.close();

  DATEI.open("u.dx",std::ios :: out);
  u[level].Print_Dx(&DATEI);
  DATEI.close();


  cout << " Variable ausgegeben! " << endl;

  //Fernziel:
  //Blockgrid blockgrid(&cylinder,10);
  //Blockgrid blockgrid_c(blockgrid.coarse_x());
  //Blockgrid blockgrid_c(&cylinder,5);

}


    */