#include "Sundance.h"


void writeToFile(const Expr& f, const string& name);

int main(int argc, void** argv)
{

  try
    {
      PMachine::init(&argc, &argv);

      int nx = 30;
      int ny = 6;
      double tLeft = 1.0;
      double tRight = 0.5;
      Mesh mesh = rectMesh(0.0, 1.0, nx, 0.0, 1.0, ny, true);
      mesh.getSubmesh();
      string filename = "test.msh." + toString(PMachine::getRank());
      mesh.textWrite(filename);

      CellSet left(1, "left");
      CellSet right(1, "right");
      CellSet top(1, "top");
      CellSet bottom(1, "bottom");
			

      Expr x = new CoordExpr(0, "x");
      Expr tInit = tLeft + (tRight-tLeft)*x;

      Expr T0 = new DiscreteFunction(mesh, tInit, Lagrange(1), "T0");

      Expr TStep0 = new DiscreteFunction(mesh, T0, Lagrange(1), "TStep0");

      Expr varT = new TestFunction(Lagrange(1), "varT");
      Expr dT = new UnknownFunction(Lagrange(1), "dT");
			
      Expr dx = new Derivative(0,1);
      Expr dy = new Derivative(1,1);
      Expr grad = List(dx, dy);

      Expr newton = (grad*varT)*(grad*(pow(T0, 4) + 4.0*pow(T0, 3)*dT));

      Expr eqn = Integral(newton, GaussLegendre(8));
			
      EssentialBC bc = 
        EssentialBC(left, (T0+dT-tLeft)*varT) 
        && EssentialBC(right, (T0+dT-tRight)*varT) ;
			
      AztecOptions azOpt;
      azOpt.setParameter(AZ_tol, 1.0e-10);
      azOpt.setOption(AZ_max_iter, 1500);
      azOpt.setOption(AZ_solver, AZ_tfqmr);
      azOpt.setOption(AZ_output, AZ_last);
      LinearSolver solver = new AztecSolver(azOpt);
			
      StaticLinearProblem prob(mesh, eqn, bc, varT, dT, solver);

      Expr residExpr = new UnknownFunction(Lagrange(1), "resid");
      Expr residEqn 
        = varT*residExpr + (grad*varT)*(grad*(pow(TStep0,4)));



      EssentialBC residBC = 
        EssentialBC(left, varT*(residExpr + (tLeft - TStep0)))
        && EssentialBC(right, varT*(residExpr + ( tRight - TStep0)));




			
      StaticLinearProblem residProb(mesh, residEqn, residBC, varT, 
                                    residExpr, 
                                    new AztecSolver(azOpt));

      Expr residSoln;			
      residProb.solve(residSoln);

      double prevResid = residSoln.norm();

      TSFOut::printf("initial residual = %g", prevResid);

      int maxIters = 100;
      double tol = 1.0e-8;
      double resid;
      double tResid = 1.0;
      int iter=0;

      while (tResid > tol && iter < maxIters)
        {
          TSFOut::printf("newton step# %d", iter);
          Expr dTSoln;
          prob.solve(dTSoln);

          double step = 1.0;
          bool decrease = false;
          int b = 0;
					
          while (!decrease)
            {
              TSFOut::printf("\tfraction of full step = %g", step);
              TStep0 = new DiscreteFunction(mesh, T0+step*dTSoln, 
                                            Lagrange(1), "T" 
                                            + toString(iter) + "." + toString(b));

              residProb.solve(residSoln);
              resid = residSoln.norm();
              writeToFile(residSoln, "r" + toString(iter) + "." + toString(b)
                          + ".dat");
              writeToFile(TStep0, "tp" +toString(iter) + "." + toString(b)
                          + ".dat" );
              b++;
              if (b > 10) TSFError::raise("too many backtracks");
              TSFOut::printf("\tresid = %g, prevResid = %g", 
                             resid, prevResid);
							
              if (resid < prevResid) 
                {
                  decrease = true;
                  prevResid = resid;
                }
              else
                {
                  decrease = false;
                  step = 0.5*step;
                }
            }

          iter++;

          T0 = TStep0;
					
          tResid = step*dTSoln.norm();
          TSFOut::printf("norm(step*dtSoln) = %g", tResid);
					
          writeToFile(T0, "T0" + toString(iter) + ".dat");
        }
			

      if (iter >= maxIters)
        {
          TSFError::raise("backtracking newton failed to converge");
        }

      TSFOut::printf("backtracking newton converged after %d iterations",
                     iter);

      Expr exactSoln = pow(tLeft,4.0) + (pow(tRight,4.0)-pow(tLeft,4.0))*x;
			
      double errorNorm 
        = definiteIntegral(mesh, pow(exactSoln-pow(T0, 4.0), 2.0),
                           GaussLegendre(8));
			
      // compute the norm of the error
      double tolerance = 1.0e-6;
			
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  catch(exception& e)
    {
      e.print();
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  PMachine::finalize();
}


void writeToFile(const Expr& f, const string& name)
{
  ofstream of(name.c_str());
  f.matlabDump(of);
}