#include "Sundance.h"


/**


 */


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

      /*
        Create a mesh object. In this example, we will use a built-in method
        to create a uniform mesh on the unit line. In more realistic problems
        we would use a mesher to create a mesh, and then read the mesh using
        a MeshReader object. 
      */
      int n = 1;
      MeshGenerator mesher = new LineMesher(0.0, 1.0, n);
      Mesh mesh = mesher.getMesh();

      Defaults::setVectorType(new DenseSerialVectorType());


      /* create a discrete space on the mesh */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(1));

      /* initial guess is u=1.0 */

      Expr U0 = new DiscreteFunction(discreteSpace, 1.0, "u0");

      /*
        Define an unknown function and its variation. The constructor
        argument is the basis family with which the function will be
        represented, in this case second-order Lagrange (nodal) polynomials.
      */
      Expr U = new UnknownFunction(new Lagrange(1));
      Expr varU = new TestFunction(new Lagrange(1));

      /*
       * 
       */
      Expr eqn = Integral(varU*(U0 - (1.0+1.0/U0)) 
                          + varU*(1.0+1.0/U0/U0));

			
      /*
       * Create a solver object: LAPACK direct method
       */
      TSFLinearSolver solver = new DirectSolver();


      /* Linearized problem */
      StaticLinearProblem prob(mesh, eqn, varU, U);

      PicardLinearization newton(prob, U0, solver);
      Expr soln = newton.solve(NewtonSolver(50, 1.0e-12));

      /*
        write the solution in a form readable by matlab
      */
      FieldWriter writer = new MatlabWriter();
      writer.writeField(soln);

      /*
        compute the error and represent as a discrete function
      */
      Expr exactSoln = 0.5*(1.0 + sqrt(5.0));

      /*
        compute the norm of the error
      */
      double errorNorm = (soln - exactSoln).norm(2);
      double tolerance = 1.0e-10;

      /*
        decide if the error is within tolerance
      */
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);

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