#include "Sundance.h"

/** \example heat2D.cpp

 */

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

      /* create a simple mesh on the rectangle */
      int nx = 41;
      int ny = 41;
      double L = 4.0;
      MeshGenerator mesher = new RectangleMesher(-L, L, nx, -L, L, ny);
      Mesh mesh = mesher.getMesh();

      /* 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 == -L );
      CellSet right = boundary.subset( x == L );
      CellSet bottom = boundary.subset( y == -L );
      CellSet top = boundary.subset( y == L );

      /* create symbolic objects for test and unknown functions */
      Expr v = new TestFunction(new Lagrange(2));
      Expr u = new UnknownFunction(new Lagrange(2));

      /* create symbolic differential operators */
      Expr dx = new Derivative(0,1);
      Expr dy = new Derivative(1,1);
      Expr grad = List(dx, dy);

      double peclet = 1.0;

      /* */
      Expr vx = (L+y)*(L-y);
      Expr vy = 0.0;
      Expr velocity = List(vx, vy);

      /* Define a Gaussian source centered at (x0,y0) */
      double x0 = -2.0;
      double y0 = 1.0;
      Expr deltaX = hold(x-x0);
      Expr deltaY = hold(y-y0);
      Expr d2 = deltaX*deltaX + deltaY*deltaY;
      Expr source = exp(-d2);

      /* Write symbolic weak equation and Neumann and Robin BCs */
      Expr eqn = Integral(-(grad*v)*(grad*u)/peclet 
                          - v*(velocity*grad)*u
                          + v*source);


      /* 
       * Boundary conditions:
       * u=0 on left, Neumann elsewhere
       */
      EssentialBC bc = EssentialBC(left, u*v);

      /* 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, v, u);
			
      /* 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("cd2D.dat");
      writer.writeField(soln);

      /* collect sensor data and write to file */
      int nsx = 3;
      int nsy = 3;
      ofstream of("sensors.dat");
      of << (nsx+1)*(nsy+1) << endl;
      for (int i=0; i<=nsx; i++)
        {
          for (int j=0; j<=nsy; j++)
            {
              double sx = -3.0 + 6.0*((double) i)/((double) nsx);
              double sy = -3.0 + 6.0*((double) j)/((double) nsy);
              Cell sensor = mesh.cellNearestToPoint(Point(sx, sy));
              of << sx << " " << sy << " " << soln.average(sensor).value()
                 << endl;
            }
        }


    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}