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 #include "Teuchos_XMLParameterListWriter.hpp"
00034
00035 using Sundance::List;
00036
00037
00038
00039
00040
00041 CELL_PREDICATE(LeftPointTest, {return fabs(x[0]) < 1.0e-10;})
00042 CELL_PREDICATE(BottomPointTest, {return fabs(x[1]) < 1.0e-10;})
00043 CELL_PREDICATE(RightPointTest, {return fabs(x[0]-1.0) < 1.0e-10;})
00044 CELL_PREDICATE(TopPointTest, {return fabs(x[1]-1.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 = 32;
00061 int ny = 32;
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 edges = new DimensionalCellFilter(1);
00072
00073 CellFilter left = edges.coordSubset(0,0.0);
00074 CellFilter right = edges.coordSubset(0,1.0);
00075 CellFilter top = edges.subset(new TopPointTest());
00076 CellFilter bottom = edges.subset(new BottomPointTest());
00077
00078
00079
00080
00081 Expr psi = new UnknownFunction(new Lagrange(1), "psi");
00082 Expr vPsi = new TestFunction(new Lagrange(1), "vPsi");
00083 Expr omega = new UnknownFunction(new Lagrange(1), "omega");
00084 Expr vOmega = new TestFunction(new Lagrange(1), "vOmega");
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 Expr eqn = Integral(interior, (grad*vPsi)*(grad*psi)
00099 + (grad*vOmega)*(grad*omega) + vPsi*omega, quad2)
00100 + Integral(top, -1.0*vPsi, quad4);
00101
00102 Expr bc = EssentialBC(bottom, vOmega*psi, quad2)
00103 + EssentialBC(top, vOmega*psi, quad2)
00104 + EssentialBC(left, vOmega*psi, quad2)
00105 + EssentialBC(right, vOmega*psi, quad2);
00106
00107
00108
00109
00110 LinearProblem prob(mesh, eqn, bc, List(vPsi, vOmega),
00111 List(psi, omega), vecType);
00112
00113
00114 LinearSolver<double> solver
00115 = LinearSolverBuilder::createSolver("bicgstab.xml");
00116
00117
00118 std::cerr << "starting solve..." << std::endl;
00119
00120 Expr soln = prob.solve(solver);
00121
00122
00123 FieldWriter w = new VTKWriter("VorticityStokes2D");
00124 w.addMesh(mesh);
00125 w.addField("psi", new ExprFieldWrapper(soln[0]));
00126 w.addField("omega", new ExprFieldWrapper(soln[1]));
00127 w.write();
00128
00129
00130
00131
00132 Expr totalVorticityExpr = Integral(interior, soln[1], quad2);
00133 double totalVorticity = evaluateIntegral(mesh, totalVorticityExpr);
00134 std::cerr << "total vorticity = " << totalVorticity << std::endl;
00135
00136 double tol = 1.0e-4;
00137 Sundance::passFailTest(fabs(totalVorticity-1.0), tol);
00138 }
00139 catch(std::exception& e)
00140 {
00141 std::cerr << e.what() << std::endl;
00142 }
00143 Sundance::finalize(); return Sundance::testStatus();
00144 }