#include "Sundance.h"

#include "TSFTimer.h"

#include "ILUKPreconditionerFactory.h"
#include "TSFPreconditionerFactory.h"

using namespace TSF;
using namespace Sundance;


int main(int argc, void** argv)
{
  try
    {
      Timer::start ();
      MPIComm::init(&argc, &argv);

      TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(2);
      TSFLinearSolver solver = new BICGSTABSolver(precond, 1.e-14, 500);

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

      mesh.printCells(cerr);

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

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

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

			
      Expr x = new CoordExpr(0, "x");
      Expr exactSoln = x ;// -0.5*(x+2.0)*(x-2.0);

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

      EssentialBC bc = EssentialBC(left, delU*(U-exactSoln)) 
        && EssentialBC(right, delU*(U-exactSoln))
        && EssentialBC(top, delU*(U-exactSoln))
        && EssentialBC(bottom, delU*(U-exactSoln));





      StaticLinearProblem prob(mesh, eqn, bc, delU, U, new PetraMatrix());

      cerr << "operator: " << prob.getOperator() << endl;

      cerr << "rhs: " << prob.getRHS() << endl;

      Expr soln = prob.solve(solver);



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

      Expr errorExpr = new DiscreteFunction(discreteSpace, soln-exactSoln);
			
      // 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__);

      Timer::report();
    }

  catch(exception& e)
    {
      TSFOut::println(e.what());
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  MPIComm::finalize();
}