#include "Sundance.h"


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

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

      int nx = 4;
      int ny = 4;
      Mesh mesh = rectMesh(0.0, 1.0, nx, 0.0, 1.0, ny, false);

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

      Expr grad = List(dx, dy);
			
      Expr delQ = List(delQx, delQy);
      Expr Q = List(Qx, Qy);
			
      Expr e = delQ*Q + U*(grad*delQ) + delU*(grad*Q - 1.0);

			
      Expr eqn = Integral(e, GaussLegendre(4));

      EssentialBC bc = EssentialBC(top, delQy*Qy)
        && EssentialBC(bottom, delQy*Qy);

      StaticLinearProblem prob(mesh, eqn, bc, List(delU, delQx, delQy),
                               List(U, Qx, Qy), solver);

      Expr soln;
      prob.solve(soln);

      ofstream of("qx.dat");
      soln[1].matlabDump(of);

      Expr x = new CoordExpr(0, "x");
      Expr exactSoln = x - 0.5;
      Expr errorExpr = new DiscreteFunction(mesh, exactSoln-soln[1], 
                                            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__);
    }

}