00001 #include "Sundance.hpp"
00002 #include "SundanceProblemTesting.hpp"
00003
00004
00005 class NitschePoisson2DTest : public LPRectTestBase
00006 {
00007 public:
00008
00009 NitschePoisson2DTest(const Array<int>& n, int k, double C)
00010 : LPRectTestBase(n), k_(k), C_(C) {}
00011
00012
00013 std::string name() const {return "NitschePoisson2D";}
00014
00015
00016 Expr exactSoln() const
00017 {
00018 Expr x = new CoordExpr(0);
00019 Expr y = new CoordExpr(1);
00020 const double pi = 4.0*atan(1.0);
00021 return cos(pi*x/2.0)*sin(pi*y);
00022 }
00023
00024
00025 Array<int> pExpected() const {return tuple<int>(k_+1);}
00026
00027
00028 LinearProblem prob(const Mesh& mesh) const
00029 {
00030 Expr dx = new Derivative(0);
00031 Expr dy = new Derivative(1);
00032 Expr x = new CoordExpr(0);
00033 Expr y = new CoordExpr(1);
00034 Expr grad = List( dx , dy );
00035
00036 CellFilter allBdry = domain().north() + domain().south()
00037 + domain().east() + domain().west();
00038
00039 BasisFamily L = new Lagrange( k_ );
00040
00041 Expr u = new UnknownFunction( L , "u" );
00042 Expr v = new TestFunction( L , "v" );
00043 QuadratureFamily quad = new GaussianQuadrature( 2 * k_ );
00044
00045 const double pi = 4.0*atan(1.0);
00046 Expr force = 1.25*pi*pi*cos(pi*x/2.0)*sin(pi*y);
00047 Expr eqn = Integral( interior(),
00048 (grad*v) * (grad*u) - force * v , quad )
00049 + NitschePoissonDirichletBC(2, allBdry, quad,
00050 1.0, v, u, exactSoln(), C_);
00051
00052
00053 Expr bc;
00054
00055 return LinearProblem(mesh, eqn, bc, v, u, vecType());
00056 }
00057
00058
00059
00060
00061 Array<LPTestSpec> specs() const
00062 {
00063 double tol = 0.05;
00064 return tuple(
00065 LPTestSpec("amesos.xml", tol, makeSet<int>(1))
00066
00067
00068
00069
00070 );
00071 }
00072
00073 void postRunCallback(int meshID, const Mesh& mesh,
00074 const string& solverFile,
00075 const Expr& soln) const
00076 {
00077 string solverName = StrUtils::before(solverFile, ".");
00078 string file = "NitschePoisson2D-mesh-" + Teuchos::toString(meshID)
00079 + "-" + solverName + "-p-" + Teuchos::toString(k_)
00080 + "-C-" + Teuchos::toString(C_);
00081
00082 DiscreteSpace discSpace(mesh, new Lagrange(1), vecType());
00083 Expr uEx = exactSoln();
00084 L2Projector proj(discSpace, soln - uEx);
00085 Expr err = proj.project();
00086
00087 FieldWriter w = new VTKWriter(file);
00088 w.addMesh(mesh);
00089 w.addField("soln", new ExprFieldWrapper(soln));
00090 w.addField("error", new ExprFieldWrapper(err));
00091 w.write();
00092 }
00093 private:
00094 int k_;
00095 double C_;
00096
00097 };
00098
00099
00100 int main( int argc , char **argv )
00101 {
00102
00103 try
00104 {
00105 Sundance::init(&argc, &argv);
00106 Tabs::showDepth() = false;
00107 LinearSolveDriver::solveFailureIsFatal() = false;
00108 int np = MPIComm::world().getNProc();
00109
00110 LPTestSuite tests;
00111
00112 Array<int> nx = tuple(4, 8, 16, 32, 64);
00113 double C = 100.0;
00114
00115 for (int p=1; p<=3; p++)
00116 {
00117 tests.registerTest(rcp(new NitschePoisson2DTest(nx, p, C)));
00118 }
00119
00120 bool pass = tests.run();
00121
00122 Out::root() << "total test status: " << pass << std::endl;
00123
00124 Sundance::passFailTest(pass);
00125 }
00126 catch(std::exception& e)
00127 {
00128 Sundance::handleException(e);
00129 }
00130 Sundance::finalize();
00131 return Sundance::testStatus();
00132 }