#include "Sundance.h"

void writeToFile(const Expr& f, const char *name, int i);

int main()
{

  try
    {
      //ExprBase::verbosePrintOn();
      int nx = 20;
      double tLeft = 1.0;
      double tRight = 0.5;
      Mesh mesh = lineMesh(0.0, 1.0, nx);

      CellSet left(0, "left");
      CellSet right(0, "right");

      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), "blah");
      writeToFile(TStep0, "blah", 0);
      TStep0 = new DiscreteFunction(mesh, tInit, Lagrange(1), "TStep0");

      Expr varT = new TestFunction(Lagrange(1), "varT");
      Expr dT = new UnknownFunction(Lagrange(1), "dT");
			
      Expr dx = new Derivative(0,1);
			
      Expr newton = (dx*varT)*(dx*(pow(T0, 4) + 4.0*pow(T0, 3)*dT));

      cerr << "newton eqn = " << newton  << endl << endl;

      Expr eqn = Integral(newton);
			
      EssentialBC bc = 
        EssentialBC(left, (T0+dT-tLeft)*varT) 
        && EssentialBC(right, (T0+dT-tRight)*varT) ;
			
      LinearSolver solver = new BICGSTABSolver(1.e-12, 5000);
      StaticLinearProblem prob(mesh, eqn, bc, varT, dT, solver);

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

      cerr << "resid eqn = " << residEqn << endl << endl;

      EssentialBC residBC = 
        EssentialBC(left, varT*(residExpr + (tLeft - TStep0)))
        && EssentialBC(right, varT*(residExpr + ( tRight - TStep0)));
			
      StaticLinearProblem residProb(mesh, residEqn, residBC, varT, 
                                    residExpr, 
                                    new BICGSTABSolver(1.e-12, 5000));

      Expr residSoln;			
      residProb.solve(residSoln);

      writeToFile(residSoln, "r", 0);
      writeToFile(TStep0, "tsInit", 0);
      double prevResid = residSoln.norm();
      cerr << "initial residual is " << prevResid << endl;

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

      while (tResid > tol && iter < maxIters)
        {
          cerr << "newton step# " << iter << endl;

          cerr << "probBC " << bc << endl;
          Expr dTSoln;
          prob.solve(dTSoln);

          double step = 1.0;
          bool decrease = false;
          int b = 0;
					

					
          while (!decrease)
            {
              cerr << "newton step length = " << step << endl;
              char tmpName[100];
              sprintf(tmpName, "Ts%d", b);
              TStep0 = new DiscreteFunction(mesh, T0, 
                                            Lagrange(1), "blah");
              residProb.solve(residSoln);
              resid = residSoln.norm();
              writeToFile(residSoln, "s", b);
              writeToFile(TStep0, "tp", b);
							
              TStep0 = new DiscreteFunction(mesh, T0+step*dTSoln, 
                                            Lagrange(1), tmpName);
              writeToFile(TStep0, "ts", b);
              b++;
              cerr << "residBC " << residBC << endl;
              if (b > 20) TSFError::raise("too many backtracks");
              residProb.solve(residSoln);
              resid = residSoln.norm();
              writeToFile(residSoln, "r", b);
              cerr << " resid=" << resid << " prevResid=" << prevResid << endl;
							
              if (resid < prevResid) 
                {
                  decrease = true;
                  prevResid = resid;
                }
              else
                {
                  decrease = false;
                  step = 0.5*step;
                  cerr << "backtracking to step = " << step << endl;
                }
            }

          iter++;

          T0 = TStep0;
					
          tResid = step*step*dTSoln.norm();
          cerr << "norm(step*dtSoln) = " << tResid << endl;
					
          writeToFile(T0, "T0", iter);
        }
			
      cerr << "backtracking newton converged" << endl;
			
    }
  catch(exception& e)
    {
      e.print();
      exit(1);
    }
				

}


void writeToFile(const Expr& f, const char *name, int i)
{
  char fName[100];
  sprintf(fName, "%s%d.dat", name, i);
  ofstream of(fName);
  f.matlabDump(of);
}