#include <iostream>
#include <fstream>
#include <set>
#include <boost/shared_ptr.hpp>
#include "include/linalg.hpp"
#include "include/rbf.hpp"
#include "include/bvp.hpp"
#include "include/error.hpp"
#include "include/utils.hpp"
#include "include/diff_op.hpp"
#include "include/ddm.hpp"
#include "include/func.hpp"

using namespace kwantix;

class conv_diff : public linear_diff_op2{
public:
  conv_diff(double u_in) : u(u_in) {} ;
  double at(const realfunc& f, const point &p) const;
private:
  Laplacian L;
  del1 D;
  double u;
};

double conv_diff::at(const realfunc& f, const point &p) const{
  return L.at(f,p) - u*D.at(f,p);
}

class bdry_conv_diff : public bdry_diff_op{
public:
  double at(const realfunc& f, const point &p) const{
    f(p);
    throw "this function shouldn't be called";
  };
  double at(const realfunc& f, const point &p, 
	    const vector &n) const{
    if(p(1) == 0 or p(1) == 1)
      return D(f,p,n);
    else
      return N(f,p,n);
  };
private:
  neumann_op N;
  dirichlet_op D;
};

double f(const kwantix::point& p){
  p.size();
  return 0;
}

double g(const kwantix::point& p){
  if(p(1) == 0)
    return 1;
  if(p(1) == 1)
    return 2;
  return 0;  
}

using boost::shared_ptr;

int main(){
  
  gsl_set_error_handler(&kwantix::errorHandler);

  try{
    using namespace std;
    using kwantix::vector;
    shared_ptr<domain> 
      Omega(new domain("data/intr.v", "data/bdry.v", "data/nrml.v"));

    shared_ptr<overlapping_domain>
      O1(new overlapping_domain("data/intr1.v", 
				"data/bdry1.v", "data/nrml1.v")),
      O2(new overlapping_domain("data/intr2.v", 
				"data/bdry2.v", "data/nrml2.v")),
      O3(new overlapping_domain( "data/intr3.v", 
				 "data/bdry3.v","data/nrml3.v")),
      O4(new overlapping_domain("data/intr4.v", 
				"data/bdry4.v", "data/nrml4.v"));
    set<shared_ptr<overlapping_domain> > dd;
    dd.insert(O1);
    dd.insert(O2);
    dd.insert(O3);
    dd.insert(O4);
    set_overlapper_info(dd);
    
    shared_ptr<conv_diff> L(new conv_diff(10));
    shared_ptr<bdry_conv_diff> B(new bdry_conv_diff);
    shared_ptr<linear_BVP2> thebvp (new linear_BVP2(Omega,L,B,f,g));
    
    set<shared_ptr<const domain> > dd_in;
    dd_in.insert(O1);
    dd_in.insert(O2);
    dd_in.insert(O3);
    dd_in.insert(O4);

    piecewise_polynomial::set_n(7);
    thin_plate_spline::set_n(8);
    c_infty_kwantix::set_epsilon(7.234);

    additive_schwarz_ddm<thin_plate_spline> solver(dd_in, thebvp);
    
    matrix out(Omega->get_interior().size() + Omega ->get_boundary().size(),
	       3);
    slice s(2,3);
    
    size_t i = 1;
    set<point>::iterator I;
    for(I = Omega -> get_interior().begin(); I != Omega->get_interior().end();
	I++){
      out(i,1) = solver(*I);
      out(i,s) = *I;
      i++;
    }
    for(I = Omega -> get_boundary().begin(); I != Omega->get_boundary().end();
	I++){
      out(i,1) = solver(*I);
      out(i,s) = *I;
      i++;
    }

    matrix::set_precision(30);
    ofstream ofs("results/result.matrix");
    ofs << out;
    ofs.close();

    return 0;
    
  }
  catch(error exc){
    kwantix::show_exception(exc);
  }
  return 0;
}