#include "Sundance.h"


/**
\example optDistribPoisson1D.cpp



 */


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

      /* mesh the line segment [0, pi] */
      int n = 10;
      const double pi = 4.0*atan(1.0);
      MeshGenerator mesher = new LineMesher(0.0, pi, 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 - pi) < 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 v = u.variation();

      Expr lambda = new UnknownFunction(new Lagrange(2));
      Expr mu = lambda.variation();

      Expr alpha = new UnknownFunction(new Lagrange(2));
      Expr beta = alpha.variation();
			
      /* Define the differentiation operator of order 1 in direction 0 */
      Expr dx = new Derivative(0);

      /* define a target function */
      Expr target = sin(x);			

      /* control penalty weight */
      double R = 0.1;

      /* write the least-squares objective function plus a control penalty */
			
      Expr objectiveFunction = 0.5*Integral(pow(u-target, 2.0))
        + 0.5*R*Integral(alpha*alpha);

      /* write the lagrangian */
      Expr lagrangian = objectiveFunction + Integral(-(dx*u)*(dx*lambda))
        + Integral(-lambda*alpha);

      /* take variations to derive the first-order necessary condition */
      Expr eqn = lagrangian.variation(List(u, lambda, alpha));

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

      EssentialBC bc = EssentialBC(left, u*mu) && EssentialBC(right, u*mu) 
        && EssentialBC(left, lambda*v) && EssentialBC(right, lambda*v);

			


      /*
        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, List(mu, v, beta),
                               List(u, lambda, alpha));

      /*
        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();
      cerr << "u" << endl;
      writer.writeField(soln[0]);
      cerr << "lambda" << endl;
      writer.writeField(soln[1]);
      cerr << "alpha" << endl;
      writer.writeField(soln[2]);

			
    }
  catch(exception& e)
    {
      TSFOut::println(e.what());
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  Sundance::finalize();
}