#ifndef CLASS_PDE
#define CLASS_PDE

#include <string>
#include <iostream>
#include <iomanip>

#include "Field.h"
#include "TriMatrix.h"

using namespace std;

class PDE {
 public:
  PDE(int N, double dt, double a, double b, double theta);
  ~PDE();
  
  void        timestep();
  void        construct();
  double      xv(int i) const;
  
  const Field getU() const;
  void        print();
  
  friend ostream &operator<<(ostream &stream, PDE &pde);
  
  virtual void setOperator() =0;
  virtual void setBC() =0;
  virtual void setIC() =0;
  
 protected:
  int        N;              // Number of spatial grid points.
  Field      u;              // Field to store the function u at each grid point.

  TriMatrix  lhs, rhs;       // Left-hand side and right-hand side operators which 
                             // are created by the construct() function.

  TriMatrix  le, li;         // Explicit and implicit differential operators for this 
                             // PDE, which will be set by the subclass implementing
                             // setOperator().

  double     theta;          // Value of theta to use for the theta method.

  double     dx, dt;         // Grid spacing in space and time. 
  double     a, b;           // Left-hand and right-hand endpoints of the domain.
  double     *x;             // Vector of length N storing the spatial grid points; i.e.
                             // x[i] = a + i*dx.

  int        ts;             // Timestep number.
  double     t;              // Current time (i.e. dt * ts)

  bool       rebuild;        // If rebuild is set, then the construct() function should 
                             // be called each time we take a timestep.
  bool       built;          // Flag indicating whether lhs and rhs have ever been 
                             // constructed.
};

#endif