/*! \file diff_op.hpp 
 *
 *  \brief Differential operators for boundary-value problems (BVP)
 *  are declared here.
 *  \file diff_op.cpp
 *  \brief Implementations and instantiations of the functions and
 *  classes declared in diff_op.hpp
 */

#pragma once

#include "linalg.hpp"
#include "func.hpp"

namespace kwantix{
/*
 *Class graph: (other leaves may also be present)
 *
 *              diff_op----------                 Dirichlet
 *             /       \         \               /
 *linear_diff_op       diff_op2   bdry_diff_op---
 *            \        /                         \
 *           linear_diff_op2                      Neumann
 *                 |
 *                 |        
 *                 |        
 *                 |        
 *                 +-- del1   
 *                 |
 *                 +-- del2
 *                 |
 *                 +-- Laplacian
 *                 |
 *                 ... etc
 *
 */

/*! \brief A general differential operator. 
 *
 * All this class is evaluate the operator applied at a specified
 * point for a real-valued function that take a vector
 * argument. This is a pure abstract class.
 */
class diff_op{ 
public:
  ///Evaluate the function \f$f\f$ at the point \f$p\f$.
  /** All this does is call the pure virtual at(realfunc, point)
   * function below.
   */
  double operator()(const realfunc &f, const point &p) const;
  /** \brief Pure virtual function. 
   *
   * Actual implementation of the differential operator goes here.
   */
  virtual double at(const realfunc &f, const point &p) const = 0;
  ///Nothing to destroy.
  virtual ~diff_op(){};
};

/// A linear diff_op. Also pure abstract class for now.  
class linear_diff_op : virtual public diff_op{  };

/// An at most 2nd-order differential operator. 
class diff_op2 : virtual public diff_op{ };

/*! \brief The heat, wave, and Laplace's equation use this kind of
 * operators: they're both linear and at most 2nd order.
 */
class linear_diff_op2 : public linear_diff_op, public diff_op2{
public:
  double operator()(const realfunc &f, const point &p) const;
};
  

/// An operator for specifying boundary conditions
class bdry_diff_op : public diff_op{
public:
    
  double operator()(const realfunc &f, const point &p) const;
  /// Only calls at(realfunc, point, vector) below.
  double operator()(const realfunc &f, const point &p, 
                    const vector &n) const;
  virtual double at(const realfunc &f, 
                    const point &p) const =0;
  /*! \brief A boundary differential operator may need to know about
   *  the normal vector at the boundary.
   *
   *  A boundary differential operator must provide evaluation at
   *  the boundary, or explicitly ignore it. Robin and Neumann
   *  boundary operators both make use of this normal vector.
   *  \param f - The function to which the operator is applied.
   *  \param p - The point at which the operator is evaluated
   *  \param n - The normal vector for this point on the boundary.
   */
  virtual double at(const realfunc &f, const point &p, 
                    const vector &n) const =0;
};

/// Dirichlet boundary conditions
class dirichlet_op : public bdry_diff_op{
public:
  double at(const realfunc &f, const point &p) const;
  double at(const realfunc &f, const point &p, const vector &n) const;
};


/// Neumann boundary conditions
class neumann_op : public bdry_diff_op{
public:
  double at(const realfunc &f, const point &p, const vector &n) const;    
private:
  double at(const realfunc &f, const point &p) const {return f(p);};
};

/// Identity operator
class Id_op : public linear_diff_op2{
  double at(const realfunc &f, const point &p) const;
}; 

/*! \brief Partial wrt to some direction. 
 *
 * By default (i.e with the overloaded at() function), wrt first
 * coordinate. Else, last argument gives the direction to evaluate
 * the partial.
 */
class del1 : public linear_diff_op2{
public:
  double at(const realfunc &f, const point &p) const;
  double at(const realfunc &f, const point &p, size_t i) const;
};

/*! \brief Second partials wrt some direction. 
 *
 * By default (i.e with the overloaded at() function), wrt first
 * coordinate, twice. Else, last two arguments give the directions
 * to evaluate the partials.
 */
class del2 : public linear_diff_op2{
public:    
  double at(const realfunc &f, const point &p) const;
  double at(const realfunc &f, const point &p, size_t i, size_t j) const;
};

///The Laplacian
class Laplacian : public linear_diff_op2{ 
public:
  double at(const realfunc &f, const point &p) const;
};

}//namespace kwantix