#include "Sundance.h"
#include "NewtonSolver.h"



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


      /* create a simple mesh on the unit line */
      double L=1.0;
      int nx = 2;
      int ny = 2;
      int npx = 1;
      int npy = 1;
      MeshGenerator mesher 
        = new PartitionedRectangleMesher(-L, L, nx, npx, -L, L, ny, npy);
      Mesh mesh = mesher.getMesh();

      /* define an expression representing the x-coordinate function */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      /* create a cell set on the middle */
      CellSet lines = new DimensionalCellSet(1);
      CellSet interface = lines.subset( x == 0.0 );
      CellSet maximal = CellSet::maximalCells();
      CellSet left = maximal.subset( x <= 0.0 );
      CellSet right = maximal.subset( x >= 0.0 );

      /* create a discrete space on the mesh */
      TSFVectorSpace leftSpace = new SundanceVectorSpace(mesh, new Lagrange(2), left);
      TSFVectorSpace rightSpace = new SundanceVectorSpace(mesh, new Lagrange(2), right);

      Expr f0 = new DiscreteFunction(leftSpace, y);
      Expr g0 = new DiscreteFunction(rightSpace, y);

      Expr f = new UnknownFunction(new Lagrange(2));
      Expr g = new UnknownFunction(new Lagrange(2));
      Expr df = f.variation();
      Expr dg = g.variation();

      Expr eqn = Integral(interface, df*(f0 - f) + dg*(g0 - g));

      WorkSet::verboseOn();

      StaticLinearProblem prob(mesh, eqn, List(df, dg), List(f, g));

      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(prec, 1.0e-14, 300);
      Expr soln = prob.solve(solver);

      double errorNorm = /*(soln[0]-f0).norm(2, interface) + */(soln[1]-g0).norm(2, interface);

      
      double tolerance = 1.0e-10;
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();

}