#include "Sundance.h"


static Mesh mesh ;

static CellSet left(1, "left");
static CellSet right(1, "right");

static Expr delU = new TestFunction(Lagrange(2), "delU");

static Expr U = new UnknownFunction(Lagrange(2), "U");
static Expr dx = new Derivative(0,1);
static Expr dy = new Derivative(1,1);

static Integral eqn;
static EssentialBC bc;

int main(int argc, void** argv)
{
  try
    {
      PMachine::init(&argc, &argv);

      LinearSolver solver = new BICGSTABSolver(1.e-14, 5000);

      int nx = 10;
      int ny = 10;
      mesh = rectMesh(-2.0, 2.0, nx, -1.0, 1.0, ny, false);

      //CellSet left(1, "left");
      //CellSet right(1, "right");
			
      // create an expr
      //Expr delU = new TestFunction(Lagrange(2), "delU");
      //Expr U = new UnknownFunction(Lagrange(2), "U");
      //Expr dx = new Derivative(0,1);
      //Expr dy = new Derivative(1,1);

      double fValue = 1.0;			
      Expr f = new DiscreteFunction(mesh, fValue, Lagrange(1), "f");
      double bcValue = 0.0;			
      Expr bcFunc = new DiscreteFunction(mesh, bcValue, Lagrange(1), "bc");

      Expr grad = List(dx, dy);
			
      Expr e = -(grad*delU)*(grad*U)  + f*delU;

			
      /* Integral */ eqn = Integral(e, GaussLegendre(4));

      /* EssentialBC */ bc = EssentialBC(left, delU*(U-bcFunc)) 
                          && EssentialBC(right, delU*(U-bcFunc));





      StaticLinearProblem prob(mesh, eqn, bc, delU, U, solver);

      Expr soln;
      prob.solve(soln);

      ofstream of("u.dat");
      soln.matlabDump(of);

      Expr x = new CoordExpr(0, "x");
      Expr exactSoln = -0.5*(x+2.0)*(x-2.0);
      Expr errorExpr = new DiscreteFunction(mesh, exactSoln-soln, 
                                            Lagrange(1), "errorExpr");
			
      // compute the norm of the error
      double errorNorm = errorExpr.norm();
      double tolerance = 1.0e-9;
				
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);

    }
  catch(exception& e)
    {
      e.print();
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }

}