#include "Sundance.h"


/** \example inlinePoissonBoltzmann1D.cpp
 * Solve the Poisson-Boltzmann equation \f$\nabla^2 u = e^-u$ on the unit
 * line with boundary conditions:
 * Left: Natural, du/dx=0
 * Right: Dirichlet u = 2 log(cosh(1/sqrt(2)))
 *
 * The solution is 2 log(cosh(x/sqrt(2))).
 *
 * We solve using fixed-point iteration
 */

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

      /* create a simple mesh on the unit line */
      double L=1.0;
      int n = 80;
      MeshGenerator mesher = new LineMesher(0.0, L, n);
      Mesh mesh = mesher.getMesh();

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

      /* create a cell set representing the right boundary */
      CellSet boundary = new BoundaryCellSet();
      CellSet right = boundary.subset( x == 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 u0 = new DiscreteFunction(discreteSpace, x);

      /* create symbolic objects for test and unknown functions. At each newton
       * step we will solve a linearized equation for a step du, so our
       * unknown is du. */
      Expr v = new TestFunction(new Lagrange(2));
      Expr u = new UnknownFunction(new Lagrange(2));

      /* create a differential operator representing the x-derivative. */
      Expr dx = new Derivative(0);
			
      /* linearized weak equation  */
      Expr newton = (dx*u)*(dx*v) + v*exp(-u0) ;
      Expr eqn = Integral(newton);

      /* Dirichlet boundary condition */
      double uBC = 2.0*log(cosh(L/sqrt(2.0)));					
      EssentialBC bc = EssentialBC(right, (u-uBC)*v) ;
			
      /* linear problem for the improved solution u */
      StaticLinearProblem prob(mesh, eqn, bc, v, u); 

      /* create linear solver */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(1.0e-14, 500);

			
      int maxIters = 10;
      double tol = 1.0e-12;
      double resid = 1;
      int iter=0;
      while (resid > tol && iter++ < maxIters)
        {
          TSFOut::printf("Picard step# %d", iter);

          /* solve for the improved solution */
          Expr uImproved = prob.solve(solver);
          /* compute the norm of the difference between u0 and uImproved. When this]
           * becomes small enough, we can stop  */
          resid = (uImproved-u0).norm(2);
          TSFOut::printf("resid = %g\n", resid);

          /* replace the old solution with the improved solution */
          u0 = uImproved;

          /* write the current iterate to a matlab file  */
          char fName[100];
          sprintf(fName, "u%d.dat.%d", iter, MPIComm::world().getRank());
					
          FieldWriter writer = new MatlabWriter(fName);
          writer.writeField(u0);

          /* discard the values of the problem's RHS vector. The matrix doesn't need rebuilding,
           * so we keep the matrix values */
          prob.flushVectorValues();
        }

			
			
      /* compare to exact solution */
      Expr exactSoln = 2.0*log(cosh(x/sqrt(2.0)));

      // compute the norm of the error
      double errorNorm = (exactSoln - u0).norm(2);
      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();
}