#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 */
      BasisFamily lagr1 = new Lagrange(1);
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, tuple(lagr1, lagr1),
                                  new DenseSerialVectorType());

      /* initial guess is u=1.0 */

      Expr u0 = DiscreteFunction::discretize(discreteSpace, List(Expr(1.0), Expr(1.0)));
      Expr v0 = u0[0];
      Expr w0 = u0[1];

      /*
        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 dv = new UnknownFunction(new Lagrange(1));
      Expr varV = new TestFunction(new Lagrange(1));
      Expr dw = new UnknownFunction(new Lagrange(1));
      Expr varW = new TestFunction(new Lagrange(1));
      Expr du = List(dv, dw);
      Expr varU = List(varV, varW);

      Expr v = v0 + dv;
      Expr w = w0 + dw;

      /*
       * 
       */
      Expr eqn = Integral(varV*(-1.0-1.0/v0 + v + dv/v0/v0)) 
        + Integral(varW*(w0*w0 - 2.0 + 2.0*w0*dw));

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


      /* Linearized problem */
      StaticLinearProblem prob(mesh, eqn, varU, du,
                               new DenseSerialVectorType());

      NewtonLinearization newton(prob, u0, solver);
      Expr soln = newton.solve(NewtonSolver(solver, 50, 1.0e-14, 1.0e-14));

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

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

      /*
        compute the norm of the error
      */
      double errorNorm = (soln[0] - exactSoln[0]).norm(2) 
        + (soln[1] - exactSoln[1]).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();
}