#include "Sundance.h"


/**
\example heat1D.cpp
Solve the heat equation with a constant source term on the unit line.
This example will illustrate the creation of a simple mesh, the definition
of a symbolic representation of a finite-element problem, and the
solution of that problem.

 */


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 nProc = MPIComm::world().getNProc();
      int n = 10*nProc;
      MeshGenerator mesher = new PartitionedLineMesher(0.0, 1.0, n);
      Mesh mesh = mesher.getMesh();


      /* 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));



      /*
        Define the differentiation operator of order 1 in direction 0.
      */
      Expr dx = new Derivative(0);
      
      

      /*
        Create a discrete function with constant value 1.0, discretized
        on our mesh object with first-order Lagrange basis. We will use
        this discrete function as the source term. (We could simply use
        a constant expression 1.0 instead, but we wanted to illustrate
        how to create a discrete function). 
      */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(1));

      Expr f = new DiscreteFunction(discreteSpace, 1.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 = Integral(-(dx*U)*(dx*varU) + f*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 solver = new BICGSTABSolver(prec, 1.0e-14, 300);



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

      /*
        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);

      /*
        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 - exactSoln).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();
}