#include "Sundance.h"
#include "UnknownParameter.h"
#include "TestParameter.h"
#include "BlockBacksolver.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 = 10;
      MeshGenerator mesher = new LineMesher(0.0, 1.0, n);
      Mesh mesh = mesher.getMesh().getSubmesh();

      /* Define a symbolic object to represent the x coordinate function. */
      Expr x = new CoordExpr(0);

      /*
       * Define a cell set that contains all boundary cells 
       */
      CellSet boundary = new BoundaryCellSet();
      /*
       *	Define a cell set that includes all cells at position x=0.
       */
      CellSet left = boundary.subset( fabs(x - 0.0) < 1.0e-10 );

      /*
       *	Define a cell set that includes all cells at position x=1.
       */
      CellSet right = boundary.subset( fabs(x - 1.0) < 1.0e-10 );



      /*
        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(2));
      Expr varU = new TestFunction(new Lagrange(2));

      Expr alpha = new UnknownParameter();
      Expr beta = new UnknownParameter();
      Expr varAlpha = new TestParameter();
      Expr varBeta = new TestParameter();

      TSFVectorType denseVT = new DenseSerialVectorType();
      TSFVectorType vt = new PetraVectorType();

      TSFArray<Block> unkBlocks = tuple(Block(U, vt), 
                                        Block(List(alpha,beta), denseVT));
      TSFArray<Block> varBlocks = tuple(Block(varU, vt), 
                                        Block(List(varAlpha,varBeta), denseVT));

      /*
        Define the differentiation operator of order 1 in direction 0.
      */
      Expr dx = new Derivative(0);
			
			
      /*
        Having constructed the unknown, variation, source, 
        and differentiation operator,
        we can now define the variational form of the problem on the interior
        of the domain. 
      */
      Expr eqn = varAlpha*(2.0*alpha + beta - 17.0) 
        + varBeta*(alpha - beta + 2.0) 
        + Integral(-(dx*U)*(dx*varU) + (3.0*alpha -2.0*beta)*varU);

			

      /*
        Now specify the boundary conditions on the left and right CellSets.
      */

      EssentialBC bc = 
        EssentialBC(left, varU*U) && EssentialBC(right, varU*U);


      /*
        Create a solver object: stablized biconjugate gradient solver
      */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver0 = new BICGSTABSolver(1.0e-14, 300);
      TSFLinearSolver solver1 = new DirectSolver();
      TSFLinearSolver solver = new BlockBacksolver(solver0, solver1);

      /*
        Combine the geometry, the variational form, the BCs, and the solver
        to form a complete problem.
      */
      StaticLinearProblem prob(mesh, eqn, bc, varBlocks, unkBlocks);

      /*
        solve the problem, obtaining the solution as a (discrete) Expr object
      */
      Expr soln = prob.solve(solver);



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

      cerr << "params = " << soln[1][0] << " " << soln[1][1] << endl;
      /*
        compute the error and represent as a discrete function
      */
      Expr exactSoln = 0.5*x*(1.0-x);

      /*
        compute the norm of the error
      */
      double errorNorm = (soln[0] - exactSoln).norm(2)
        + sqrt(pow(soln[1][0].value()-5.0, 2.0)) 
        + sqrt(pow(soln[1][1].value()-7.0, 2.0)); 

      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();
}