//Solves 
//
//       nabla^2 u = f on  Omega
//               u = g on dOmega
//
//Input data is read from interior.matrix and boundary.matrix
//files. Each row in interior.matrix corresponds to an interior
//point of Omega, and similarly for boundary.matrix. Although it's not
//used, also need a normal vector at each point of the boundary as
//given by normals.matrix

#include <iostream>
#include <fstream>
#include <list>
#include <cmath>

#include <boost/shared_ptr.hpp>

#include "include/linalg.hpp"
#include "include/error.hpp"
#include "include/rbf.hpp"
#include "include/diff_op.hpp"
#include "include/bvp.hpp"
#include "include/interpolator.hpp"
#include "include/utils.hpp"
#include "include/func.hpp"

#define RBF_TYPE inverse_multiquadric //Replace by template later.

using namespace std;
using namespace kwantix;

//Some arbitrary interior function for Poisson  equation (change to
//"return 0*p(1);", for example, for Laplace equation).
class intr_func : public realfunc{
public:
  double at(const point& p) const{return 0*(p(1)+p(2));};
};

//Boundary function for Poisson equation
class bdry_func : public realfunc{
public:
  //sin(theta)^3
  double at(const point& p) const{return 
      std::pow(std::sin(std::atan2(p(1),p(2))),3);};
};

void output_result(const interpolator<RBF_TYPE>& u);

int main(){
  using kwantix::vector;

  try{
    gsl_set_error_handler(&kwantix::errorHandler);
    
    //Define the domain from given data
    shared_ptr<domain> Omega(new domain("data/interior.matrix",
                                        "data/boundary.matrix",
                                        "data/normals.matrix"));

    //Define the interior and boundary operators
    shared_ptr<Laplacian> L(new Laplacian);
    shared_ptr<dirichlet_op> B(new dirichlet_op);

    //rhs of BVP
    intr_func f;
    bdry_func g;
    
    //Define the Poisson boundary value problem.
    shared_ptr<linear_BVP2> Poisson_BVP (new linear_BVP2(Omega, L, B, f, g));

    //Define interpolator and solve the BVP.
    interpolator<RBF_TYPE> u(Poisson_BVP);

    //Write result to file
    output_result(u);

  }
  catch(error& exc){
    kwantix::show_exception(exc);
  }

  return 0;
}

void output_result(const interpolator<RBF_TYPE>& u)
{
  //Output result as a 41x41 meshgrid (to be plotted, say, with
  //Octave; show_sol script included to work with this data).
  matrix sol(41,41);    
    
  double r, th;
  const double pi = atan(1)*4;
  point p(2);
  th = 0;
  for(size_t i = 1; i <= 41; i++){
    r = 0;
    for(size_t j = 1; j <= 41; j++){
      p(1) = r*cos(th); 
      p(2) = r*sin(th);
      sol(i,j) = u(p); //Evaluate interpolator at point p
      r += 1.0/40;
    }

    th += 2.0*pi/40;
  }

    
  ofstream ofs("data/sol.matrix");
  ofs << sol;
}