#include "Sundance.h"

/** 
 * Solve Laplace eqn on the Sandia thunderbird
 * logo. The thunderbird has been meshed with Jonathan Shewchuk's
 * Triangle Delaunay mesher; the mesh is read into Sundance using a
 * ShewchukMeshReader object.  
 *
 * We apply Dirichlet boundary conditions on the left and right sides
 * of the thunderbird's tail. The left side of the tail is the line
 * segment between (-3,-10) and (-2,0), and the right side is the line
 * segment between (3,-10) and (2,0). These can be identified as the
 * boundary subsets on the lines y = 20 +/- 10 x.  We set u=-1 on the
 * left and u=1 on the right.
 */

int main(int argc, void** argv)
{
  MPIComm::init(&argc, &argv);
  
  try
    {
      Timer::start();
      
      /* read the mesh */
      MeshReader reader = new ShewchukMeshReader("tBird.1");  
      Mesh mesh = reader.getMesh();

      /* define expressions representing coordinate functions. These will
      * be used in the identification of boundary subsets. */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      CellSet boundary = new BoundaryCellSet();
      CellSet leftInlet = boundary.subset( x < 0.0 && y == -7.0 );
      CellSet rightInlet = boundary.subset( x > 0.0 && y == -7.0 );
      CellSet topOutlet = boundary.subset( x == -4.0 );
      CellSet bottomOutlet = boundary.subset( y == -10.0 );
      
      CellSet inflow = leftInlet + rightInlet;
      CellSet outflow = topOutlet + bottomOutlet;

      CellSet walls = boundary - inflow - outflow;

      CellSet bottomLeftWall = walls.subset( x >= -6.5 && y < 1.0 && x < 0.0 );
      CellSet bottomRightWall = walls.subset( x <= 6.5 && y < 1.0 && x > 0.0 );
      CellSet topLeftWall = walls.subset( x == -7.0 
                                          || (x < 0.0 && y > 1.0 && y < 6.0) );
      CellSet topRightWall = walls.subset( y == 7.0 || x==7.0 || x==2.0 
                                           || (x > 0.0 && y==3.0));

      /* create test and unknown functions */
      Expr v = new TestFunction(new Lagrange(1));
      Expr u = new UnknownFunction(new Lagrange(1));

      /* create gradient operator */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);
      
      /* write weak equation */
      Integral eqn = Integral(-(grad*v)*(grad*u));

      /* set dirichlet boundary conditions */
      EssentialBC bc = EssentialBC(bottomLeftWall, v*u)
        && EssentialBC(topLeftWall, v*(u-1.0))
        && EssentialBC(topRightWall, v*(u-2.0))
        && EssentialBC(bottomRightWall, v*(u-3.0));
      
      /* create a preconditioner and solver */
      TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(2);
      TSFLinearSolver solver = new BICGSTABSolver(precond, 1.e-14, 500);

      /* stuff everything into a problem object */
      StaticLinearProblem prob(mesh, eqn, bc, v, u,
                               new PetraVectorType()); 

      /* get the solution */
      Expr solnU = prob.solve(solver);

      /* write the solution in matlab format */
      ofstream of("tBirdHeat.dat");
      solnU.matlabDump(of);

      Timer::report();
    }
  catch(exception& e)
    {
      TSFOut::println(e.what());
    }
  MPIComm::finalize();

}
tBirdHeat.cpp
