#include "Sundance.h"

#include "TSFTimer.h"

#include "FieldWriter.h"
#include "VTKWriter.h"
#include "ILUKPreconditionerFactory.h"
#include "TSFPreconditionerFactory.h"

using namespace TSF;
using namespace Sundance;


bool onOuter(const Cell& cell)
{
  // x == -10
  const Point& a = cell.facet(0,0).point(0);
  const Point& b = cell.facet(0,1).point(0);
  double ra = sqrt(a*a);
  double rb = sqrt(b*b);
		
  if (ra > 199.0 && rb > 199.0)
    {
      return true;
    }
  return false;
}

int main(int argc, void** argv)
{
  try
    {
      Timer::start ();
      MPIComm::init(&argc, &argv);
			
      MeshReader reader = new ShewchukMeshReader("naca0012.1");
      Mesh mesh = reader.getMesh();
      mesh.labelCells(1, "outer", onOuter);
      CellSet outer(1, "outer");

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

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

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

      Expr x = new CoordExpr(0, "x");
			
      Expr eqn = Integral(e, new GaussLegendre(4));

      EssentialBC bc = EssentialBC(outer, delU*(U-x/200.0)) ;


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

      Expr soln = prob.solve(solver);



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

      FieldWriter writer = new VTKWriter("airfoil.vtk");
      writer.writeField(mesh, "velocityPotential", soln);

      Timer::report();
    }

  catch(exception& e)
    {
      TSFOut::println(e.what());
    }
  MPIComm::finalize();
}