#include "Sundance.h"


int main()
{
  try
    {
      Profiler::start("main");
      int n = 10;
      double h = 1.0/((double) n);
      Mesh mesh = rectMesh(0.0, 1.0, n, 0.0, 1.0, n);
      CellSet left(1, "left");
      CellSet right(1, "right");

			
      // create an expr
      Expr varU = TestFunction(Lagrange(1), "varU");
      Expr U(Lagrange(1), "U");
      Expr lambda(Lagrange(1), "lambda");
      Expr varLambda = TestFunction(Lagrange(1), "varLambda");
      Derivative dx(0,1);
      Derivative dy(1,1);

      Expr grad = List(dx, dy);
			
      Expr e = -1.0*(grad*varU)*(grad*U);

			
      Expr eqn = Integral(e, GaussLegendre(4))
        + Integral(left, -1.0*varLambda*U + varU*lambda, GaussLegendre(4))
        + Integral(right, -1.0*varLambda*(U-1.0) + varU*lambda, GaussLegendre(4));


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

      StaticLinearProblem prob(mesh, eqn, List(varU, varLambda),
                               List(U, lambda), solver);

      Expr solnU;

			
      prob.solve(solnU);


      Expr x = CoordExpr(0, "x");
      Expr exactSoln = x;

      Expr residFunc(Lagrange(1), "residFunc");
      Expr varResidFunc = TestFunction(Lagrange(1), "varResidFunc");
			
      Expr errorEqn = Integral(varResidFunc*(residFunc-(solnU[0]-exactSoln)));
      LinearSolver s2 = new BICGSTABSolver(1.e-12, 5000);

      StaticLinearProblem errorProb(mesh, errorEqn, 
                                    varResidFunc, residFunc, s2);

      Expr errorExpr;
			
      errorProb.solve(errorExpr);

      ofstream of("tmp.dat");
      errorExpr.matlabDump(of);
			
      mesh.diagnostics();

      Profiler::stop("main");
      Profiler::report(cerr);
    }
  catch(exception& e)
    {
      e.print();
      exit(1);
    }
				

}