#include "Sundance.h"


int main(int argc, void** argv)
{
  try
    {
      Sundance::init(&argc, &argv);
			
      double Lx = 1.0;
      double Ly = 1.0;
      double Lz = 1.0;
			
      int nx = 4;
      int ny = 6;
      int nz = 8;

      MeshGenerator mesher = new BrickMesher(0.0, Lx, nx,
                                             0.0, Ly, ny,
                                             0.0, Lz, nz);
      Mesh mesh = mesher.getMesh();

			
			
      /* define coordinate functions for x and y coordinates */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);
      Expr z = new CoordExpr(2);
      Expr u0 = x + sqrt(2.0)*y + atan(1.0)*z;

      /* define cells sets for each of the six sides of the box */
      CellSet boundary = new BoundaryCellSet();
      CellSet left = boundary.subset( fabs(x) < 1.0e-9 );
      CellSet right = boundary.subset( fabs(x-Lx) < 1.0e-9 );
      CellSet front = boundary.subset( fabs(y) < 1.0e-9 );
      CellSet back = boundary.subset( fabs(y-Ly) < 1.0e-9 );
      CellSet bottom = boundary.subset( fabs(z) < 1.0e-9 );
      CellSet top = boundary.subset( fabs(z-Lz) < 1.0e-9 );
			
      /* Create a vector space factory, used to 
       * specify the low-level linear algebra representation */
      TSFVectorType petra = new PetraVectorType();
			

      Expr varU = new TestFunction(new Lagrange(1));
      Expr U = new UnknownFunction(new Lagrange(1));

      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr dz = new Derivative(2);

      Expr grad = List(dx, dy, dz);     

      Expr e = (grad*varU)*(grad*U);      

      Expr eqn = Integral(e, new GaussianQuadrature(2));
      
      EssentialBC bc = EssentialBC(boundary, varU*(U-u0));
      //      EssentialBC bc = EssentialBC(left, varU*U) && EssentialBC(right, varU*(U-1));

      /* Assemble everything into a problem object, with a specification that
       * Petra be used as the low-level linear algebra representation */
      StaticLinearProblem prob(mesh, eqn, bc, varU, U, petra);
			
      /* create a preconditioner and solver */
      TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(precond, 1.e-14, 500);

      /* solve the problem and return the solution as a symbolic object */
      Expr soln = prob.solve(solver);

      /* write to matlab */
      FieldWriter writer = new MatlabWriter("heat3D.dat");
      writer.writeField(soln);

      /* compare to known solution */
      Expr exactSoln = u0;

      // compute the norm of the error
      double errorNorm = (soln-exactSoln).norm(2, new GaussianQuadrature(2));
      double tolerance = 1.0e-9;
				
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}