#include "Sundance.h"


int main(int argc, void** argv)
{
  try
    {
      Sundance::init(&argc, &argv);
			
      MeshReader reader = new ShewchukMeshReader("../Meshes/wing3");
      Mesh mesh = reader.getMesh();
      
      /* define coordinate functions for x and y coordinates */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);
      Expr z = new CoordExpr(2);


      /* define cells sets for each of the six sides of the box */
      CellSet boundary = new BoundaryCellSet();
      CellSet leading = boundary.subset( x == -2.0 ); 
      CellSet trailing = boundary.subset( x == 8.0 );
      CellSet bottom = boundary.subset( z == -2.0 );
      CellSet top = boundary.subset( z == 2.0 );
      CellSet inner = boundary.subset( y == 0.0 );
      CellSet outer = boundary.subset( y == 4.0 );
      CellSet walls = leading + trailing + top + bottom + inner + outer;
      CellSet wing = boundary - walls;
			
      /* Create a vector space factory, used to 
       * specify the low-level linear algebra representation */
      TSFVectorType petra = new PetraVectorType();
			

      Expr varU = new TestFunction(new Lagrange(1));
      Expr U = new UnknownFunction(new Lagrange(1));

      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr dz = new Derivative(2);

      Expr grad = List(dx, dy, dz);     

      Expr e = (grad*varU)*(grad*U);      

      Expr eqn = Integral(e, new GaussianQuadrature(2));
      
      EssentialBC bc = EssentialBC(leading, varU*U) && EssentialBC(trailing, varU*(U-1));
      

      /* Assemble everything into a problem object, with a specification that
       * Petra be used as the low-level linear algebra representation */
      StaticLinearProblem prob(mesh, eqn, bc, varU, U, petra);
			
      /* create a preconditioner and solver */
      TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(precond, 1.e-14, 500);

      /* solve the problem and return the solution as a symbolic object */
      Expr soln = prob.solve(solver);

      /* write to matlab */
      FieldWriter writer = new MatlabWriter("wing3.dat", wing+inner);
      writer.writeField(soln);

      /* compare to known solution */
      Expr exactSoln = z;
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}