00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "Sundance.hpp"
00032 #include "SundanceEvaluator.hpp"
00033
00034 using Sundance::List;
00035
00036
00037
00038
00039 CELL_PREDICATE(LeftPointTest, {return fabs(x[0]) < 1.0e-10;})
00040 CELL_PREDICATE(BottomPointTest, {return fabs(x[1]) < 1.0e-10;})
00041 CELL_PREDICATE(RightPointTest, {return fabs(x[0]-1.0) < 1.0e-10;})
00042 CELL_PREDICATE(TopPointTest, {return fabs(x[1]-1.0) < 1.0e-10;})
00043
00044 CELL_PREDICATE(CornerPointTest, {return fabs(x[1]) < 1.0e-10 && fabs(x[0])<1.0e-10;})
00045
00046
00047 int main(int argc, char** argv)
00048 {
00049
00050 try
00051 {
00052 Sundance::init(&argc, &argv);
00053 int np = MPIComm::world().getNProc();
00054
00055
00056 VectorType<double> vecType = new EpetraVectorType();
00057
00058
00059
00060 int nx = 8;
00061 int ny = 8;
00062 MeshType meshType = new BasicSimplicialMeshType();
00063 MeshSource mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, np,
00064 0.0, 1.0, ny, 1,
00065 meshType);
00066 Mesh mesh = mesher.getMesh();
00067
00068
00069
00070 CellFilter interior = new MaximalCellFilter();
00071 CellFilter bdry = new BoundaryCellFilter();
00072 CellFilter nodes = new DimensionalCellFilter(0);
00073 CellFilter corner = nodes.subset(new CornerPointTest());
00074
00075
00076
00077 Expr ux = new UnknownFunction(new Lagrange(2), "u_x");
00078 Expr vx = new TestFunction(new Lagrange(2), "v_x");
00079 Expr uy = new UnknownFunction(new Lagrange(2), "u_y");
00080 Expr vy = new TestFunction(new Lagrange(2), "v_y");
00081 Expr p = new UnknownFunction(new Lagrange(1), "p");
00082 Expr q = new TestFunction(new Lagrange(1), "q");
00083 Expr u = List(ux, uy);
00084 Expr v = List(vx, vy);
00085
00086
00087 Expr dx = new Derivative(0);
00088 Expr dy = new Derivative(1);
00089 Expr grad = List(dx, dy);
00090 Expr x = new CoordExpr(0);
00091 Expr y = new CoordExpr(1);
00092
00093
00094 QuadratureFamily quad2 = new GaussianQuadrature(2);
00095 QuadratureFamily quad4 = new GaussianQuadrature(4);
00096
00097
00098
00099 double pi = 4.0*atan(1.0);
00100 Expr sx = sin(pi*x);
00101 Expr cx = cos(pi*x);
00102 Expr sy = sin(pi*y);
00103 Expr cy = cos(pi*y);
00104 Expr psiExact = pow(pi, -3.0) * sx*sy;
00105 Expr uExact = pow(pi, -2.0)*List(-sx*cy, cx*sy);
00106 Expr fy = 4.0*cx*sy;
00107
00108 Expr eqn = Integral(interior, (grad*vx)*(grad*ux)
00109 + (grad*vy)*(grad*uy) - p*(dx*vx+dy*vy)
00110 + q*(dx*ux+dy*uy) - vy*fy,
00111 quad2);
00112
00113
00114 Expr bc = EssentialBC(bdry, v*(u-uExact), quad4)
00115 + EssentialBC(corner, q*p, quad4);
00116
00117
00118 LinearProblem prob(mesh, eqn, bc, List(vx, vy, q),
00119 List(ux, uy, p), vecType);
00120
00121
00122 #ifdef HAVE_CONFIG_H
00123 ParameterXMLFileReader reader(searchForFile("SolverParameters/aztec-native.xml"));
00124 #else
00125 ParameterXMLFileReader reader("aztec-native.xml");
00126 #endif
00127 ParameterList solverParams = reader.getParameters();
00128 LinearSolver<double> solver
00129 = LinearSolverBuilder::createSolver(solverParams);
00130
00131 Expr soln = prob.solve(solver);
00132
00133
00134
00135 FieldWriter w = new VTKWriter("TaylorHoodStokes2d");
00136 w.addMesh(mesh);
00137 w.addField("ux", new ExprFieldWrapper(soln[0]));
00138 w.addField("uy", new ExprFieldWrapper(soln[1]));
00139 w.addField("p", new ExprFieldWrapper(soln[2]));
00140 w.write();
00141
00142 Expr err = List(soln[0], soln[1]) - uExact;
00143 Expr errExpr = Integral(interior,
00144 err*err,
00145 quad4);
00146
00147 FunctionalEvaluator errInt(mesh, errExpr);
00148
00149 double errorSq = errInt.evaluate();
00150 std::cerr << "velocity error norm = " << sqrt(errorSq) << std::endl << std::endl;
00151
00152 Sundance::passFailTest(sqrt(errorSq), 1.0e-3);
00153
00154 }
00155 catch(std::exception& e)
00156 {
00157 Sundance::handleException(e);
00158 }
00159 Sundance::finalize(); return Sundance::testStatus();
00160 }