#include "Sundance.h"

/** \example exp2D.cpp
 * Solve Poisson's equation with an exponential source term on the 
 * rectangle [0,1] x [0,2] with the following boundary conditions:
 *
 * Left:   Neumann, du/dx = e^y
 * Bottom: Dirichlet, u= 0.5 x^2 + e^x
 * Right:  Neumann, du/dx = 1 + e^(y+1)
 * Top:    Dirichlet, u=0.5 x^2 + 2/3 + e^(x+2)
 *
 * Let the source term be f(x,y) = 1 + 2 e^(x+y). 
 *
 * The solution is u(x,y) = 0.5*x^2 + y/3 + e^(x+y).
 *
 * This problem cannot be solved exactly in the space of second order 
 * polynomials. Provided that a sufficiently accurate quadrature 
 * rule is used, the error norm should go as h^2. 
 */

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

      /* create a simple mesh on the rectangle */
      int nx = 64;
      int ny = 64;
      MeshGenerator mesher = new RectangleMesher(0.0, 1.0, nx, 0.0, 2.0, 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 == 0.0 );
      CellSet right = boundary.subset( x == 1.0 );
      CellSet bottom = boundary.subset( y == 0.0 );
      CellSet top = boundary.subset( y == 2.0 );
			
      /* Create a vector space factory, used to 
       * specify the low-level linear algebra representation */
      TSFVectorType petra = new PetraVectorType();
			
      /* create a discrete space on the mesh */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(1), petra);


      /* Define source expression */
      Expr f = 1.0 + 2.0*exp(x+y);

      /* 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);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);

      /* create a high-order quadrature rule for accuracy on the 
       * exponential terms */
      QuadratureFamily quad = new GaussianQuadrature(6);

      /* Write symbolic weak equation and natural BCs */
      /* Flux on left = -exp(y) */
      /* Flux on right = 1.0 + exp(1+y) */
      Expr poisson = Integral(-(grad*v)*(grad*u)) 
        + Integral(-f*v)
        + Integral(left, -v*exp(y), quad) 
        + Integral(right, v*(1.0 + exp(1.0+y)), quad);


      /* Write essential BCs: 
       * Bottom: u= 0.5 x^2 + e^x
       * Top:  u= 0.5 x^2 + 2/3 + e^(x+2)
       */
      Expr uBottom = 0.5*x*x + exp(x);
      Expr uTop = 0.5*x*x + 2.0/3.0 + exp(x+2.0);
      EssentialBC bc = EssentialBC(bottom, v*(u - uBottom), quad)
        && EssentialBC(top, v*(u - uTop), quad);

      /* Assemble everything into a problem object, with a specification that
       * Petra be used as the low-level linear algebra representation */
      StaticLinearProblem prob(mesh, poisson, bc, v, u, petra);

      /* create a preconditioner and solver */
      TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(2);
      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("exp2D.dat");
      writer.writeField(soln);

      /* compare to known solution */
      Expr exactSoln = 0.5*x*x + y/3.0 + exp(x+y);
      Expr errorExpr = soln - exactSoln;
			
      double errorNorm = errorExpr.norm(2, quad);
      double tolerance = 1.0e-5;
				
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }

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