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

/** \example radDiffusion2D.cpp
 * Solve the radiation diffusion equation \f$\nabla^2 T^4 = 0\f$ on the 
 * unit rectangle.  
 */

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

      /* create mesh */
      double L=1.0;
      int nx = 32;
      int ny = 32;
      MeshGenerator mesher = new RectangleMesher(0.0, L, nx, 0.0, L, ny);
      Mesh mesh = mesher.getMesh();

      /* define coordinate expressions */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      /* exact solution */
      Expr exactSoln = pow(1.0 + cos(x)*exp(y), 0.25);

      /* create a cell set representing the right boundary */
      CellSet boundary = new BoundaryCellSet();

      /* 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. Try T0(x,y)=1 + x + y.
       */
      Expr T0 = new DiscreteFunction(discreteSpace, 1.0 + x + y);
      Expr bcFunc = exactSoln;

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

      /* 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 equations and boundary conditions */
      Expr eqn = Integral(-(grad*pow(T,4.0))*(grad*v));
      EssentialBC bc = EssentialBC(boundary, v*(T-bcFunc));

      /* linearize the equations and BC */
      Expr linearizedEqn = eqn.linearization(T, T0);
      EssentialBC linearizedBC = bc.linearization(T, T0);
      Expr dT = T.differential();

      /* linear problem for the step dT */
      StaticLinearProblem prob(mesh, linearizedEqn, linearizedBC, v, dT); 

      /* create linear solver */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(prec, 1.0e-14, 1000);
			
      NewtonLinearization newton(prob, T0, solver);
      Expr soln = newton.solve(NewtonSolver(solver, 10, 1.0e-12, 1.0e-12));
			
      Expr error = new DiscreteFunction(discreteSpace, exactSoln-soln);
      cerr << "exact soln = " << exactSoln << endl;
      /* write to matlab */
      FieldWriter writer = new MatlabWriter("radDiffusion2D.dat");
      FieldWriter writer2 = new MatlabWriter("error.dat");
      writer.writeField(soln);
      writer2.writeField(error);

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

}