#include "Sundance.h"
#include "NewtonSolver.h"

/** \example inlinePoissonBoltzmann1D.cpp
 * Solve the Poisson-Boltzmann equation \f$\nabla^2 u = e^-u$ on the unit
 * rectangle with boundary conditions:
 * Left: Natural, du/dx=0
 * Right: Dirichlet u = 2 log(cosh(1/sqrt(2)))
 * Top and Bottom: Natural
 *
 * The solution is 2 log(cosh(x/sqrt(2))).
 *
 * The problem is nonlinear, so we use Newton's method to iterate 
 * towards a solution. 
 *
 */

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

      /* create a simple mesh on the unit line */
      double L=1.0;
      int nx = 16;
      int ny = 16;
      int npx = 1;
      int npy = 1;
      MeshGenerator mesher 
        = new PartitionedRectangleMesher(0.0, L, nx, npx, 0.0, L, ny, npy);
      Mesh mesh = mesher.getMesh();

      /* define an expression representing the x-coordinate function */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      /* create a cell set representing the right boundary */
      CellSet boundary = new BoundaryCellSet();
      CellSet top = boundary.subset( y == L );

      /* create a discrete space on the mesh */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(2));

      /* create an expression for the initial guess. This will be reused as the
       * starting point for each newton step. Assume u(x)=x as an initial 
       * guess, and discretize it.
       */
      Expr exactSoln = 2.0*log(cosh(y/sqrt(2.0)));
      Expr uBC = 2.0*log(cosh(L/sqrt(2.0)));
      Expr u0 = new DiscreteFunction(discreteSpace, x+y);

      /* create symbolic objects for test and unknown functions. */
      Expr u = new UnknownFunction(new Lagrange(2), "u");
      Expr v = u.variation();
      Expr du = u.differential();

      /* create a differential operator representing the x-derivative. */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);
			
      /* write the nonlinear equation and BC */
      Expr eqn = Integral((grad*u)*(grad*v) + exp(-u)*v);
      EssentialBC bc = EssentialBC(top, (u - exactSoln)*v);

      Expr linearizedEqn = eqn.linearization(u, u0);
      EssentialBC linearizedBC = bc.linearization(u, u0);

      NewtonLinearization::verbose() = true;

      /* create linear problem */
      StaticLinearProblem prob(mesh, linearizedEqn, linearizedBC, v, du); 

      /* create linear solver */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(prec, 1.0e-14, 1000);
			
      NewtonLinearization newton(prob, u0, solver);
      Expr soln = newton.solve(NewtonSolver(solver, 15, 1.0e-6, 1.0e-6));

      double errorNorm = (exactSoln - soln).norm(2);			
      /* write to matlab */
      string filename = "pb2D." + TSF::toString(MPIComm::world().getRank())
        + ".dat";
      FieldWriter writer = new MatlabWriter(filename);
      writer.writeField(soln);	

      /* compute the norm of the error */

      double tolerance = 1.0e-4;
      TSFOut::printf("error = %g\n", errorNorm);

      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();

}