#include "Sundance.h"
#include <cstdio>

/**
 * \example fisher1D.cpp
 */

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

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

      /* create unknown and variational functions */
      Expr U = new UnknownFunction(new Lagrange(2),"u");
      Expr delU = U.variation();

      /* create a differentiation operator */
      Expr dx = new Derivative(0);

      /* the initial conditions will be u0(x,t=0) = sin(pi*x).
       * create a coordinate expression to represent x, then
       * create sin(pi*x), and then project it onto a discrete function. */

      const double pi = 4.0*atan(1.0);
      Expr x = new CoordExpr(0);

      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(2));

      Expr one = new ParameterExpr(1.0);
      Expr zero = new ParameterExpr(0.0);
			
      Expr U0 = new DiscreteFunction(discreteSpace, one - 0.9*sin(pi*x/2.0));
      Expr UStart = new DiscreteFunction(discreteSpace, one - 0.9*sin(pi*x/2.0));

      /* 
         set up backward Euler timestepping 
      */

      Expr deltaT = new ParameterExpr(0.1, "dt");

      Expr dU = U.differential();

      Expr uAvg = (U + UStart)/2.0;

      /* Define BCs to be one at both ends */
      CellSet boundary = new BoundaryCellSet();
      CellSet left = boundary.subset( fabs(x - 0.0) < 1.0e-10 );
      CellSet right = boundary.subset( fabs(x - 1.0) < 1.0e-10 );

      Expr nonlinEqn = Integral(delU*(U - UStart)/deltaT 
                                + (dx*delU)*(dx*uAvg) 
                                - delU*(uAvg*(one - uAvg)));
      //				- Integral(right, delU*(dx*UStart)/2.0);

      Expr linEqn = nonlinEqn.linearization(U, U0);


      EssentialBC nonLinBC = EssentialBC(left, (U - one)*delU);
      //				&& EssentialBC(right, (U-one)*delU);
      EssentialBC linBC = nonLinBC.linearization(U,U0);


      /* create a solver object */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(prec, 1.0e-12, 300);
      solver.setVerbosityLevel(0);

      /* 
         put the time-discretized eqn into a StaticLinearProblem object
         which will do the spatial discretization.
      */
      StaticLinearProblem prob(mesh, linEqn, linBC, delU, dU);
      //NewtonLinearization newton(prob, UNew, solver); 
      NewtonLinearization newton(prob, U0, solver); 
      /* 
         Now, loop over timesteps, solving the elliptic problem for u at each
         step. At the end of each step, assign the solution solnU into u0.
         Because Exprs are stored by reference, the updating of u0 propagates
         to the copies of u0 in the equation set and in the 
         StaticLinearProblem. The same StaticLinearProblem can be reused
         at all timesteps.
      */
      int nSteps = 150;
      for (int i=0; i<nSteps; i++)
        {
          cerr << "-------- timestep " << i << " ---------" << endl;
          /* solve the problem */
          Expr soln = newton.solve(NewtonSolver(solver, 100, 1.0e-6, 1.0e-6));
          TSFVector solnVec;
          soln.getVector(solnVec);
          UStart.setVector(solnVec.copy()); 
          U0.setVector(solnVec); 
          /* write the solution at step i to a file */
          char fName[20];
          sprintf(fName, "timeStep%d.dat", i);
          ofstream of(fName);
          FieldWriter writer = new MatlabWriter(fName);
          writer.writeField(U0);
          /* flush the matrix and RHS values */
          prob.flushMatrixValues();
        }

      /* compute the exact solution and the error */
      Expr tFinal = nSteps * deltaT;
      Expr exactSoln = one;
			
      /*
        compute the norm of the error
      */
      double errorNorm = (exactSoln-U0).norm(2);
      double tolerance = 1.0e-4;

      /*
        decide if the error is within tolerance
      */
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}