#include "Sundance.h"


/**
 *\example coupled1D.cpp
 * Solves two problems coupled at an interface
 */


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

      /*
       * Create a mesh object. 
       */
      int n = 2;
      MeshGenerator mesher = new LineMesher(-1.0, 1.0, n);
      Mesh mesh = mesher.getMesh().getSubmesh();

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

      /* get the left and right boundary cell sets */
      CellSet boundary = new BoundaryCellSet();
      CellSet leftEnd = boundary.subset( fabs(x + 1.0) < 1.0e-10 );
      CellSet rightEnd = boundary.subset( fabs(x - 1.0) < 1.0e-10 );

      /* get the interface cell set */
      CellSet zeroCells = new DimensionalCellSet(0);
      CellSet interface = zeroCells.subset(x==0.0);

      /* get the left and right zone cell sets */
      CellSet maximal = CellSet::maximalCells();
      CellSet leftZone = maximal.subset( x <= 0.0 ) ;
      CellSet rightZone = maximal.subset( x >= 0.0 ) ;
			



      /* define field variables. 
       * u1 is the field on the left side
       * u2 is the field on the right side
       * lambda is the lagrange multiplier used to impose the
       * interface continuity condition 
       */
      Expr u1 = new UnknownFunction(new Lagrange(1));
      Expr u2 = new UnknownFunction(new Lagrange(1));
      Expr lambda = new UnknownFunction(new Lagrange(1));
      Expr v1 = u1.variation();
      Expr v2 = u2.variation();
      Expr mu = lambda.variation();

      /* define the derivative */
      Expr dx = new Derivative(0);
			
      /* write the variational equation for both domains 
       * plus the interface condition */
      Expr eqn = Integral(leftZone, (dx*u1)*(dx*v1))
        + Integral(rightZone,  (dx*u2)*(dx*v2))
        + Integral(interface, v1*lambda - v2*lambda)
        + Integral(interface, mu*(u1-u2));

      /* boundary conditions */
      EssentialBC bc = 
        EssentialBC(leftEnd, v1*u1) && EssentialBC(rightEnd, v2*(u2-1.0));

      /*
        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(v1,v2,mu), 
                               List(u1,u2,lambda));
      cerr << prob.getOperator() << endl;
      cerr << prob.getRHS() << endl;

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

      TSFVector u1Vec;
      TSFVector u2Vec;
      soln[0].getVector(u1Vec);
      soln[1].getVector(u2Vec);
			
      cerr << u1Vec << endl;
      cerr << u2Vec << endl;

      cerr << "u1 norm = " << soln[0].norm(2, leftZone) << endl;
      cerr << "u2 norm = " << soln[1].norm(2, rightZone) << endl;

      FieldWriter w1 = new MatlabWriter("u1.dat", leftZone);
      FieldWriter w2 = new MatlabWriter("u2.dat", rightZone);

      w1.writeField(u1);
      w2.writeField(soln[1]);
    }
  catch(exception& e)
    {
      TSFOut::println(e.what());
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  Sundance::finalize();
}