#include "Sundance.h"


/** \example catenary.cpp
 *
 * catenary.cpp finds the minimum energy configuration of a hanging
 * chain, or equivalently, the minimum area of a surface of
 * revolution. The functional to be minimized is
 * 
 * Integral( u * sqrt(1 + (dx*u)*(dx*u) )
 * 
 * We can derive a differential equation by taking the variation of this
 * functional. The differential equation is nonlinear, so we linearize
 * and solve with Newton's method. 
 *
 * We will solve the problem on the domain [-1, 1] with boundary conditions
 * u(-1) = u(1) = 2.0. With these boundary conditions, the exact solution is 
 * 
 * y = cosh(alpha*x)/alpha
 *
 *  where alpha is the solution of the nonlinear algebraic equation 
 *
 * 2*alpha = cosh(alpha).
 *
 * The auxiliary function catenaryConstant() uses Newton's method to solve
 * this equation for alpha. 
 */

double catenaryConstant(double y0, double x0, double guess)
{
  double a0 = guess;

  for (int i=0; i<10; i++)
    {
      double f = y0*a0 - cosh(a0*x0);
      double df = y0 - x0*sinh(a0*x0);
      if (fabs(f/df) < 1.0e-14) return a0;
      a0 = a0 - f/df;
    }
  TSFError::raise("Newton's method failed to converge in determining "
                  "the exact solution to the catenary problem");
}

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

      int nx = 32;
      double a = 1.0;
      double y0 = 2.0;



      double lintol = 1.0e-14;
      double ftol = 1.0e-13;  // newton solver tolerance on function residual
      double stol = 1.0e-13;  // newton solver tolerance on change in resi
      int maxiterlin = 100;
      int maxiter = 100;

      /* Create a mesh object. */
      MeshGenerator mesher = new PartitionedLineMesher(-a, a, nx);
      Mesh mesh = mesher.getMesh();

      /* compute the exact solution */
      Expr x = new CoordExpr(0);
      double alpha = catenaryConstant(y0, a, 0.75);
      Expr exactSoln = cosh(alpha*x)/alpha;

      /* create boundary cell sets */
      CellSet boundary = new BoundaryCellSet();
      CellSet left = boundary.subset( fabs(x + a) < 1.0e-10 );
      CellSet right = boundary.subset( fabs(x - a) < 1.0e-10 );

      /* Generate a discrete vector space of first-order 
       * Lagrange polynomials */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh,new Lagrange(2));
      Expr initial = 1.0;
      Expr u0 = DiscreteFunction::discretize(discreteSpace, initial );

      /* Define an unknown function and its variation. */

      Expr u = new UnknownFunction(new Lagrange(2), "u");
      Expr w = u.variation();
      
      /* Define differentiation operators */
      Expr dx = new Derivative(0);
      
      /* Write energy functional */
      Expr E = Integral( u*sqrt(1.0 + (dx*u)*(dx*u)) );
      
      /* take variations to obtain a nonlinear PDE */
      Expr nonLinEqn = E.variation(u);
       
      /* linearize */
      Expr linEqn = nonLinEqn.linearization(u,u0);
      Expr du = u.differential();

      /* specify the boundary conditions */
      EssentialBC nonLinBC = EssentialBC(boundary, (u-2.0)*w);

      /* linearize the BCs */
      EssentialBC linBC = nonLinBC.linearization(u,u0);

      /* create a solver */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(prec, lintol, maxiterlin);

      /* create a linearized problem */
      StaticLinearProblem prob(mesh, linEqn, linBC, w, du);
      NewtonLinearization newton(prob, u0, solver);      

      /* solve the nonlinear problem */
      Expr soln = newton.solve(NewtonSolver(solver, maxiter, ftol, stol));
       
      /* write the solution */
      char fName[100];
      sprintf(fName, "catenary.dat.%d", MPIComm::world().getRank());
      FieldWriter writer = new MatlabWriter(fName);
      writer.writeField(soln);

      /* check against the exact solution */
      double errorNorm = (exactSoln - soln).norm(2);  
      double tolerance = 1.0e-6;
      TSFOut::printf("error = %g\n", errorNorm);

      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }
  catch(exception& e)
    {
      TSFOut::println(e.what());
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  Sundance::finalize();
}