/* 
   
   Copyright (C) Bernard Piette 
                 University of Durham (UK) 
                 Department of Mathematical Sciences
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   You can modify or use any part of the source code but only under the  
   condition that this Copyright notice is kept in the resulting source file.
*/
#include <pde.h>
#include <iostream>

using namespace std;

// pde constructor
// xmin, xmax : lattice range
// nx : number of lattice points
//*******************************
pde::pde(double xmin,double xmax,int nx): 
  t_(0.0),Xmin_(xmin),Xmax_(xmax),nx_(nx)
{ dx_ = (Xmax_-Xmin_)/(nx_-1);
  dt_ = dx_*dx_*0.25;
  grf_index_ = 0;

  try 
  { Fcts_  = new double[nx_];
  } 
  catch(...) { throw((string)"No more memory!"); }
}

// Destructor
//*************************
pde::~pde()
{ delete [] Fcts_;
}

// Set f_, fl_, fll_, fr_, frr_, x_ and i_. 
// f : fct array.                        
// i : index position of the grid          
//*****************************************
void pde::set_neighbour_fcts(double *f,int i)
{ 
  x_ = Xmin_+i*dx_;
  i_ = i;

  f_   = &f[i];
  fl_  = &f[i-1];
  fll_ = &f[i-2];
  fr_  = &f[i+1];
  frr_ = &f[i+2];

  if(i==0)           { fl_ =  0; fll_ = 0; }
  else if(i==1)      { fll_ = 0; }
  else if(i==(nx_-2)){ frr_ = 0; }
  else if(i==(nx_-1)){ fr_ =  0; frr_ = 0; }

}

// Evaluate the derivatives of the function.
// f_, fl_, fr_, point respectively to current lattice point, 
//               the point on left and the one on the right.
// The result is put into fx_ and fxx_.
// The method is called by the integration routine before calling
//   the equation method.
//***********************************************************
void pde::eval_derivatives()
{ if(!fl_) // left edge of the grid
  { fx_ = (fr_[0]-f_[0])/dx_;
    fxx_ = (f_[0]+fr_[1]-2*fr_[0])/(dx_*dx_); 
    return; 
  }
  if(!fr_) // right edge of the grid 
  { fx_ = (f_[0]-fl_[0])/dx_;
    fxx_ = (fl_[-1]+f_[0]-2*fl_[0])/(dx_*dx_); 
    return; 
  }

  fx_ = (fr_[0]-fl_[0])/(2*dx_);
  fxx_ = (fl_[0]+fr_[0]-2*f_[0])/(dx_*dx_);
}

// Initialise the Fcts_
// Calls Init() for every value x_i
//*********************************
void pde::initial_value()
{ int i;
  double x;

  x = Xmin_;

  for(i=0; i < nx_; i++)
  { Init(x,Fcts_[i]);
     x += dx_;
  }
}

// Create a graphic file for the selected fcts and densities
// Throw a string on error                                  
//***********************************************************
void pde::write_grf()
{ grfp_->write_1_grf(grf_index_,Fcts_,t_);

  grf_index_++;
}

// Perform n integration steps on the grid 
// By calling  integrate_1step several times.  
// Return : the updated time
//********************************************
double pde::integrate_until(double t)
{ 
  while(integrate_1step() < t);

  return(t_);
}

// Perform n integration step on the grid with graphics         
// Output graphs every grf_dt                                  
// Update the time t_                                          
// tmax : integration bound.                                   
// grf_dt : time interval between graphics                     
// Return : the updated time                                   
//***************************************************************
double pde::integrate_grf_until(double tmax,double grf_dt)
{ int i;
  string grf_name;
  double grf_t;

  grf_t = t_ + grf_dt;

  while(t_ < tmax)
  { if(t_ > grf_t) { write_grf(); grf_t += grf_dt; }
    integrate_1step();
    i++;
  }

  if(t_ > grf_t) write_grf();
 
  return(t_);
}