#include "Sundance.h"


/*

	Heat equation with Robin boundary conditions.
 
 */

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

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

      int nx = 2;
      int ny = 2;
      double L = 1.0;
      Mesh mesh = quadRectMesh(0.0, L, nx, 0.0, L, ny);

      CellSet left(1, "left");
      CellSet right(1, "right");
      CellSet top(1, "top");
      CellSet bottom(1, "bottom");
			
      // 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);

      Expr grad = List(dx, dy);
			
      Expr x = new CoordExpr(0, "x");
      Expr y = new CoordExpr(1, "y");

      Expr e = -(grad*delU)*(grad*U);
      Expr u0 = (x + y - 1.0/3.0);
      Expr u1 = (x + y + 1.0);

      Expr eqn = Integral(e) 
        //				+ Integral(bottom, delU*(U-u1))
        //+ Integral(left, delU*(U-u1))
        + Integral(top, 3.0*delU*(U-u0))
        + Integral(right, 3.0*delU*(U-u0));
			
			
      EssentialBC bc = EssentialBC(left, delU*(U-y)) 
        && EssentialBC(bottom, delU*(U-x));
      //				&& EssentialBC(right, delU*(U-x-y)) 
      //	&& EssentialBC(top, delU*(U-x-y));
			


      DistributedMatrixBuilder::verbose();

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

      prob.buildMatrixAndVector();

      DistributedMatrixBuilder::silent();

      MatrixBase* A = prob.getMatrix();
      ofstream blah("matrix.dat");
      A->print(blah);
      Expr soln;
      prob.solve(soln);

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

      Expr exactSoln = x + y;
      Expr errorExpr = new DiscreteFunction(mesh, exactSoln-soln, 
                                            Lagrange(1), "errorExpr");
      ofstream ef("e.dat");
      errorExpr.matlabDump(ef);

      // compute the norm of the error
      double errorNorm = definiteIntegral(mesh, errorExpr*errorExpr);
      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__);
    }

}