#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 n = 10;
      MeshGenerator mesher = new LineMesher(-1.0, 1.0, n);
      Mesh mesh = mesher.getMesh();

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

      mesh.labelMaximalCells("left", x <= 0.0);
      mesh.labelMaximalCells("right", x >= 0.0);
			
      Expr f = new RegionalExpr("f");
      f.setRegionalExpr("left", x);
      f.setRegionalExpr("right", x*x);
			
      /*
        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 v = new TestFunction(new Lagrange(2));

      Expr eqn = Integral(v*(u-f));

      /*
        Combine the geometry, the variational form, the BCs, and the solver
        to form a complete problem.
      */
      StaticLinearProblem prob(mesh, eqn, v, u,
                               new DenseSerialVectorType());

      /*
        solve the problem, obtaining the solution as a (discrete) Expr object
      */
      TSFLinearSolver solver = new DirectSolver();
      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 = f;

      /*
        compute the norm of the error
      */
      double errorNorm = (soln - f).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();
}