#include "Sundance.h"

/**
 * 
 */

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

      /* create a simple mesh on the rectangle */

      int nx = 160;
      int ny = 160;

      MeshGenerator mesher 
        = new RectangleMesher(0.0, 1.0, nx, 
                              0.0, 1.0, ny);

      Mesh mesh = mesher.getMesh();
      CellSet boundary = new BoundaryCellSet();

      /* define coordinate functions for x and y coordinates */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      const double pi = 4.0*atan(1.0);


      /* f, g, and h are related by:
       * Laplacian(f) = g
       * div(h) = g
       */
      Expr f = cos(pi*x)*cos(pi*y);
      Expr g = -2.0*pi*pi*f;
      Expr h = -pi * List(sin(pi*x)*cos(pi*y), cos(pi*x)*sin(pi*y));
			
			

      /* Create a vector space factory, used to 
       * specify the low-level linear algebra representation */
      TSFVectorType petra = new PetraVectorType();

      TSFVectorSpace discreteSpace = new SundanceVectorSpace(mesh, 
                                                             new Lagrange(2),
                                                             petra);
      Expr f0 = new DiscreteFunction(discreteSpace, f);
      Expr g0 = new DiscreteFunction(discreteSpace, g);
      Expr hx0 = new DiscreteFunction(discreteSpace, h[0]);
      Expr hy0 = new DiscreteFunction(discreteSpace, h[1]);
      Expr h0 = List(hx0, hy0);
			
      /* create symbolic objects for test and unknown functions */
      Expr v = new TestFunction(new Lagrange(2));
      Expr r = new UnknownFunction(new Lagrange(2));


      /* create symbolic differential operators */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);

      ExprArray eqns;
      eqns.append(Integral(v*r + (grad*f)*(grad*v) + v*g));
      eqns.append(Integral(v*r + v*g + (grad*f0)*(grad*v)));
      eqns.append(Integral(v*r + (grad*f0)*(grad*v) + v*g));
      eqns.append(Integral(v*r + (grad*f)*(grad*v) + v*(grad*h0)));
      eqns.append(Integral(v*r + (grad*h0)*v + (grad*f)*(grad*v)));
      eqns.append(Integral(v*r + (grad*f)*(grad*v) + v*g0));
      eqns.append(Integral(v*r + v*g0 + (grad*f)*(grad*v)));

      double err = 0.0;
      for (int i=0; i<eqns.length(); i++)
        {
          TSFOut::printf("-------------- Eqn %d -----------------\n",
                         i);
          StaticLinearProblem prob(mesh, eqns[i], v, r);
          double resid = prob.getRHS().norm2();
          TSFOut::printf("resid = %g\n", resid);
          err += resid;

        }
      double tolerance = 1.0e-6;
      TSFOut::printf("error = %g\n", err);

      Testing::passFailCheck(__FILE__, err, tolerance);


    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}