NitschePoisson2D.cpp

Go to the documentation of this file.

#include "Sundance.hpp"
#include "SundanceProblemTesting.hpp"


class NitschePoisson2DTest : public LPRectTestBase
{
public:
  /** */
  NitschePoisson2DTest(const Array<int>& n, int k, double C)
    : LPRectTestBase(n), k_(k), C_(C) {}

  /** */
  std::string name() const {return "NitschePoisson2D";}

  /** */
  Expr exactSoln() const 
    {
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);
      const double pi = 4.0*atan(1.0);
      return cos(pi*x/2.0)*sin(pi*y);    
    }

  /** */
  Array<int> pExpected() const {return tuple<int>(k_+1);}

  /** */
  LinearProblem prob(const Mesh& mesh) const
    {
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);
      Expr grad = List( dx , dy );

      CellFilter allBdry = domain().north() + domain().south() 
        + domain().east() + domain().west();

      BasisFamily L = new Lagrange( k_ );

      Expr u = new UnknownFunction( L , "u" );
      Expr v = new TestFunction( L , "v" );
      QuadratureFamily quad = new GaussianQuadrature( 2 * k_ );
      
      const double pi = 4.0*atan(1.0);
      Expr force = 1.25*pi*pi*cos(pi*x/2.0)*sin(pi*y);
      Expr eqn = Integral( interior(), 
        (grad*v) * (grad*u) - force * v , quad )
        + NitschePoissonDirichletBC(2, allBdry, quad, 
          1.0, v, u, exactSoln(), C_);
      

      Expr bc;

      return LinearProblem(mesh, eqn, bc, v, u, vecType());
    }

  /**
   * Specify which linear solvers are to be used.
   */
  Array<LPTestSpec> specs() const
    {
      double tol = 0.05;
      return tuple(
        LPTestSpec("amesos.xml", tol, makeSet<int>(1))/*,
                                                        LPTestSpec("aztec-ifpack.xml", tol), 
                                                        LPTestSpec("aztec-ml.xml", tol),
                                          LPTestSpec("belos-ml.xml", tol),
                                          LPTestSpec("bicgstab.xml", tol) */
        );
    }

  void postRunCallback(int meshID, const Mesh& mesh,
    const string& solverFile,
    const Expr& soln) const
    {
      string solverName = StrUtils::before(solverFile, ".");
      string file = "NitschePoisson2D-mesh-" + Teuchos::toString(meshID)
        + "-" + solverName + "-p-" + Teuchos::toString(k_)
        + "-C-" + Teuchos::toString(C_);

      DiscreteSpace discSpace(mesh, new Lagrange(1), vecType());
      Expr uEx = exactSoln();
      L2Projector proj(discSpace, soln - uEx);
      Expr err = proj.project();

      FieldWriter w = new VTKWriter(file);
      w.addMesh(mesh);
      w.addField("soln", new ExprFieldWrapper(soln));
      w.addField("error", new ExprFieldWrapper(err));
      w.write();
    }
private:
  int k_;
  double C_;
  
};


int main( int argc , char **argv )
{

  try
  {
    Sundance::init(&argc, &argv);
    Tabs::showDepth() = false;
    LinearSolveDriver::solveFailureIsFatal() = false;
    int np = MPIComm::world().getNProc();

    LPTestSuite tests;

    Array<int> nx = tuple(4, 8, 16, 32, 64);
    double C = 100.0;

    for (int p=1; p<=3; p++)
    {
      tests.registerTest(rcp(new NitschePoisson2DTest(nx, p, C)));
    }

    bool pass = tests.run();

    Out::root() << "total test status: " << pass << std::endl;

    Sundance::passFailTest(pass);
  }
  catch(std::exception& e)
  {
    Sundance::handleException(e);
  }
  Sundance::finalize(); 
  return Sundance::testStatus(); 
}