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 #include "Sundance.hpp"
00026
00027
00028 int main(int argc, char** argv)
00029 {
00030 try
00031 {
00032 const double pi = 4.0*atan(1.0);
00033 double lambda = 1.25*pi*pi;
00034
00035 int nx = 32;
00036 int nt = 10;
00037 double tFinal = 1.0/lambda;
00038
00039 Sundance::setOption("nx", nx, "Number of elements");
00040 Sundance::setOption("nt", nt, "Number of timesteps");
00041 Sundance::setOption("tFinal", tFinal, "Final time");
00042
00043 Sundance::init(&argc, &argv);
00044
00045
00046 VectorType<double> vecType = new EpetraVectorType();
00047
00048
00049 MeshType meshType = new BasicSimplicialMeshType();
00050
00051 MeshSource meshSrc = new PartitionedRectangleMesher(
00052 0.0, 1.0, nx,
00053 0.0, 1.0, nx,
00054 meshType);
00055 Mesh mesh = meshSrc.getMesh();
00056
00057
00058
00059
00060 CellFilter interior = new MaximalCellFilter();
00061 CellFilter edges = new DimensionalCellFilter(1);
00062 CellFilter west = edges.coordSubset(0, 0.0);
00063 CellFilter east = edges.coordSubset(0, 1.0);
00064 CellFilter south = edges.coordSubset(1, 0.0);
00065 CellFilter north = edges.coordSubset(1, 1.0);
00066
00067
00068 BasisFamily basis = new Lagrange(1);
00069 Expr u = new UnknownFunction(basis, "u");
00070 Expr v = new TestFunction(basis, "v");
00071
00072
00073 Expr grad = gradient(2);
00074
00075 Expr x = new CoordExpr(0);
00076 Expr y = new CoordExpr(1);
00077
00078 Expr t = new Sundance::Parameter(0.0);
00079 Expr tPrev = new Sundance::Parameter(0.0);
00080
00081
00082 DiscreteSpace discSpace(mesh, basis, vecType);
00083 Expr uExact = cos(0.5*pi*y)*sin(pi*x)*exp(-lambda*t);
00084 L2Projector proj(discSpace, uExact);
00085 Expr uPrev = proj.project();
00086
00087
00088
00089
00090
00091 QuadratureFamily quad = new GaussianQuadrature(2);
00092
00093 double deltaT = tFinal/nt;
00094
00095 Expr gWest = -pi*exp(-lambda*t)*cos(0.5*pi*y);
00096 Expr gWestPrev = -pi*exp(-lambda*tPrev)*cos(0.5*pi*y);
00097
00098
00099 Expr eqn = Integral(interior, v*(u-uPrev)/deltaT
00100 + 0.5*(grad*v)*(grad*u + grad*uPrev), quad)
00101 + Integral(west, -0.5*v*(gWest+gWestPrev), quad);
00102
00103 Expr bc = EssentialBC(east + north, v*u, quad);
00104
00105
00106 LinearProblem prob(mesh, eqn, bc, v, u, vecType);
00107
00108
00109 LinearSolver<double> solver
00110 = LinearSolverBuilder::createSolver("amesos.xml");
00111
00112 FieldWriter w0 = new VTKWriter("TransientHeat2D-0");
00113 w0.addMesh(mesh);
00114 w0.addField("T", new ExprFieldWrapper(uPrev[0]));
00115 w0.write();
00116
00117 for (int i=0; i<nt; i++)
00118 {
00119 t.setParameterValue((i+1)*deltaT);
00120 tPrev.setParameterValue(i*deltaT);
00121 Out::root() << "t=" << (i+1)*deltaT << endl;
00122 Expr uNext = prob.solve(solver);
00123
00124 ostringstream oss;
00125 oss << "TransientHeat2D-" << i+1;
00126 FieldWriter w = new VTKWriter(oss.str());
00127 w.addMesh(mesh);
00128 w.addField("T", new ExprFieldWrapper(uNext[0]));
00129 w.write();
00130
00131 updateDiscreteFunction(uNext, uPrev);
00132 }
00133
00134
00135
00136 double err = L2Norm(mesh, interior, uExact-uPrev, quad);
00137 Out::root() << "error norm=" << err << endl;
00138
00139 double h = 1.0/(nx-1.0);
00140 double tol = 0.1*(pow(h,2.0) + pow(lambda*deltaT, 2.0));
00141 Out::root() << "tol=" << tol << endl;
00142
00143
00144 Sundance::passFailTest(err, tol);
00145 }
00146 catch(std::exception& e)
00147 {
00148 Sundance::handleException(e);
00149 }
00150 Sundance::finalize();
00151 return Sundance::testStatus();
00152 }
00153