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
00033 using Sundance::List;
00034
00035 using namespace Sundance;
00036 using namespace Playa;
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 int main(int argc, char** argv)
00050 {
00051
00052 try
00053 {
00054 Sundance::init(&argc, &argv);
00055
00056
00057 VectorType<double> vecType = new EpetraVectorType();
00058
00059
00060
00061 MeshType meshType = new BasicSimplicialMeshType();
00062 int nx = 32;
00063
00064 MeshSource mesher = new PartitionedRectangleMesher(
00065 0.0, 1.0, nx,
00066 0.0, 1.0, nx,
00067 meshType);
00068 Mesh mesh = mesher.getMesh();
00069
00070
00071
00072 CellFilter interior = new MaximalCellFilter();
00073
00074
00075 CellFilter edges = new DimensionalCellFilter(1);
00076 CellFilter west = edges.coordSubset(0, 0.0);
00077 CellFilter east = edges.coordSubset(0, 1.0);
00078
00079
00080 Expr u = new UnknownFunction(new Lagrange(1), "u");
00081 Expr v = new TestFunction(new Lagrange(1), "v");
00082
00083
00084 Expr x = new CoordExpr(0);
00085 Expr grad = gradient(1);
00086
00087
00088 QuadratureFamily quad = new GaussianQuadrature(4);
00089
00090
00091 Expr xi = new Sundance::Parameter(0.0);
00092
00093
00094
00095 Expr uEx = x*(1.0-x)*(1.0+xi*exp(x));
00096 Expr f = -(-20 - exp(x)*xi*(1 + 32*x + 10*x*x +
00097 exp(x)*(-1 + 2*x*(2 + x))*xi))/10.0;
00098
00099
00100
00101 Expr eqn = Integral(interior,
00102 (1.0 + 0.1*xi*exp(x))*(grad*v)*(grad*u) - f*v, quad);
00103
00104
00105 Expr h = new CellDiameterExpr();
00106 Expr bc = EssentialBC(east + west, v*u/h, quad);
00107
00108
00109
00110 LinearProblem prob(mesh, eqn, bc, v, u, vecType);
00111
00112
00113
00114 DiscreteSpace ds(mesh, new Lagrange(1), vecType);
00115 L2Projector proj(ds, uEx);
00116
00117
00118 LinearSolver<double> solver = LinearSolverBuilder::createSolver("aztec-ml.xml");
00119 Expr soln;
00120 SolverState<double> state;
00121
00122
00123 int nSteps = 10;
00124 double xiMax = 2.0;
00125
00126
00127 Array<double> err(nSteps);
00128
00129
00130 for (int n=0; n<nSteps; n++)
00131 {
00132
00133 double xiVal = xiMax*n/(nSteps - 1.0);
00134 xi.setParameterValue(xiVal);
00135 Out::root() << "step n=" << n << " of " << nSteps << " xi=" << xiVal;
00136
00137
00138 state = prob.solve(solver, soln);
00139
00140 TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
00141 std::runtime_error,
00142 "solve failed!");
00143
00144
00145
00146 Expr uEx0 = proj.project();
00147
00148
00149 FieldWriter w = new VTKWriter("ParameterSweep-" + Teuchos::toString(n));
00150 w.addMesh(mesh);
00151 w.addField("u", new ExprFieldWrapper(soln[0]));
00152 w.addField("uEx", new ExprFieldWrapper(uEx0[0]));
00153 w.write();
00154
00155
00156 err[n] = L2Norm(mesh, interior, soln-uEx, quad);
00157 Out::root() << " L2 error = " << err[n] << endl;
00158 }
00159
00160
00161 double hVal = 1.0/(nx-1.0);
00162 double fudge = 2.0;
00163 double tol = fudge*hVal*hVal;
00164
00165
00166 double maxErr = *std::max_element(err.begin(), err.end());
00167
00168
00169 Sundance::passFailTest(maxErr, tol);
00170 }
00171 catch(std::exception& e)
00172 {
00173 Sundance::handleException(e);
00174 }
00175 Sundance::finalize();
00176 return Sundance::testStatus();
00177 }