#include "Sundance.h"

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

      MeshReader reader = new ShewchukMeshReader("Tests/step.1");
      Mesh mesh = reader.getMesh();

      double H = 2.0;
      double xLeft = -10.0;
      double xRight = 20.0;

      /* define coordinate functions for x and y coordinates */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      /* define cells sets for each of the four sides of the rectangle */
      CellSet boundary = new BoundaryCellSet();
      CellSet left = boundary.subset( x == xLeft );
      CellSet right = boundary.subset( x == xRight );
      CellSet bottomRight = boundary.subset( y == -H );
      CellSet bottomLeft = boundary.subset( y == 0.0 );
      CellSet notch = boundary.subset( x == 0.0 );
      CellSet top = boundary.subset( y == H );

      Expr uxLeft = 0.5*y*(H-y);

      /* define variations and unknowns for U, V, and P
       * using the same basis functions for velocity and pressure */
      Expr delU = new TestFunction(new Lagrange(1), "delU");
      Expr U = new UnknownFunction(new Lagrange(1), "U");
      Expr delV = new TestFunction(new Lagrange(1), "delV");
      Expr V = new UnknownFunction(new Lagrange(1), "V");
      Expr delP = new TestFunction(new Lagrange(1), "delP");
      Expr P = new UnknownFunction(new Lagrange(1), "P");

      Expr delVelocity = List(delU, delV);
      Expr velocity = List(U, V);

      /* create differential operators for x and y directions, and
       * then form gradient operator. */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);
			
      /* momentum continuity equation */
			
      Expr momentumEqn  = -(grad*U)*(grad*delU) - (grad*V)*(grad*delV)
        + P*(dx*delU + dy*delV);

      /* incompressibility constraint equation with stabilization */

      Expr beta = new ParameterExpr(0.02);
      Expr h = new CellDiameterExpr();

      Expr continuityEqn = -delP*(dx*U + dy*V) 
        //				- beta*h*h*(grad*P)*(grad*delP);
        - beta*h*h*(grad*P)*(grad*delP);

      Expr eqn = Integral(momentumEqn) + Integral(continuityEqn);

      /* boundary conditions:
       * v=0 everywhere
       * u=1 on top, u=0 elsewhere
       * p floats everywhere
       **/
      EssentialBC bc = EssentialBC(left, delU*(U-uxLeft))
        && EssentialBC(bottomLeft, delU*U + delV*V)
        && EssentialBC(bottomRight, delU*U + delV*V)
        && EssentialBC(notch, delU*U + delV*V)
        && EssentialBC(top, delU*U + delV*V);


      /*
        Create a solver object: stablized biconjugate gradient solver
      */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(2);
      TSFLinearSolver solver = new BICGSTABSolver(prec, 1.0e-12, 4000);

      StaticLinearProblem prob(mesh, eqn, bc, List(delU, delV, delP),
                               List(U, V, P), new PetraVectorType());

      Expr soln = prob.solve(solver);
      Expr u0 = soln[0];
      Expr v0 = soln[1];

      /* We've now solved the problem for the primitive variables.
       * For visualization we next compute the streamfunction */
      //			Expr delPsi = new TestFunction(new Lagrange(1));
      //Expr psi = new UnknownFunction(new Lagrange(1));

      //Expr vorticity = dx*v0 - dy*u0;
      //Integral streamfunctionEqn(-(grad*delPsi)*(grad*psi) - delPsi*vorticity);

      /* streamfunction is zero along entire boundary */
      //EssentialBC streamfunctionBC(boundary, delPsi*psi);
			
      //StaticLinearProblem streamfunctionProb(mesh, streamfunctionEqn, 
      //																		 streamfunctionBC,
      //																		 delPsi, psi, 
      //																		 new PetraVectorType());

      //Expr psi0 = streamfunctionProb.solve(solver);

      /*
        write the streamfunction in a form readable by matlab
      */
      //FieldWriter psiWriter = new MatlabWriter("psi.dat");
      //psiWriter.writeField(psi0);
      FieldWriter vWriter = new MatlabVectorWriter("v.dat");
      vWriter.writeField("velocity", List(u0, v0));
			
			

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