00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "Sundance.hpp"
00024
00025
00026 #if defined(HAVE_SUNDANCE_EXODUS) && defined(Trilinos_DATA_DIR)
00027
00028 using namespace Sundance;
00029
00030
00031
00032
00033
00034
00035 int main(int argc, char** argv)
00036 {
00037 try
00038 {
00039
00040 std::string meshFile="plateWithHole3D-1";
00041 std::string solverFile = "aztec-ml.xml";
00042 Sundance::setOption("meshFile", meshFile, "mesh file");
00043 Sundance::setOption("solver", solverFile,
00044 "name of XML file for solver");
00045
00046
00047 Sundance::init(&argc, &argv);
00048
00049
00050 VectorType<double> vecType = new EpetraVectorType();
00051
00052
00053 MeshType meshType = new BasicSimplicialMeshType();
00054 MeshSource meshSrc
00055 = new ExodusMeshReader(meshFile, meshType);
00056 Mesh mesh = meshSrc.getMesh();
00057
00058
00059
00060
00061 CellFilter interior = new MaximalCellFilter();
00062
00063
00064 CellFilter edges = new DimensionalCellFilter(2);
00065
00066 CellFilter south = edges.labeledSubset(1);
00067 CellFilter east = edges.labeledSubset(2);
00068 CellFilter north = edges.labeledSubset(3);
00069 CellFilter west = edges.labeledSubset(4);
00070 CellFilter hole = edges.labeledSubset(5);
00071 CellFilter down = edges.labeledSubset(6);
00072 CellFilter up = edges.labeledSubset(7);
00073
00074
00075
00076
00077 BasisFamily basis = new Lagrange(1);
00078 Expr u = new UnknownFunction(basis, "u");
00079 Expr v = new TestFunction(basis, "v");
00080
00081
00082 Expr dx = new Derivative(0);
00083 Expr dy = new Derivative(1);
00084 Expr dz = new Derivative(2);
00085 Expr grad = List(dx, dy, dz);
00086
00087
00088 QuadratureFamily quad1 = new GaussianQuadrature(1);
00089 QuadratureFamily quad2 = new GaussianQuadrature(2);
00090
00091
00092
00093 Expr eqn
00094 = Integral(interior, (grad*u)*(grad*v), quad1)
00095 + Integral(east, v, quad1);
00096
00097
00098 Expr h = new CellDiameterExpr();
00099 Expr bc = EssentialBC(west, v*u/h, quad2);
00100
00101
00102 LinearProblem prob(mesh, eqn, bc, v, u, vecType);
00103
00104
00105
00106
00107
00108 LinearSolver<double> solver
00109 = LinearSolverBuilder::createSolver(solverFile);
00110
00111
00112
00113 Expr soln = prob.solve(solver);
00114
00115
00116
00117
00118 DiscreteSpace discSpace(mesh, List(basis, basis, basis), vecType);
00119 L2Projector proj(discSpace, grad*soln);
00120 Expr gradU = proj.project();
00121
00122
00123 FieldWriter w = new VTKWriter("LaplaceDemo3D");
00124 w.addMesh(mesh);
00125 w.addField("soln", new ExprFieldWrapper(soln[0]));
00126 w.addField("du_dx", new ExprFieldWrapper(gradU[0]));
00127 w.addField("du_dy", new ExprFieldWrapper(gradU[1]));
00128 w.addField("du_dz", new ExprFieldWrapper(gradU[2]));
00129 w.write();
00130
00131
00132 Expr n = CellNormalExpr(3, "n");
00133 CellFilter wholeBdry = east+west+north+south+up+down+hole;
00134 Expr fluxExpr
00135 = Integral(wholeBdry, (n*grad)*soln, quad1);
00136 double flux = evaluateIntegral(mesh, fluxExpr);
00137 Out::root() << "numerical flux = " << flux << std::endl;
00138
00139
00140
00141
00142
00143 Expr x = new CoordExpr(0);
00144 Expr y = new CoordExpr(1);
00145 Expr z = new CoordExpr(2);
00146
00147 Expr xCMExpr = Integral(interior, x, quad1);
00148 Expr yCMExpr = Integral(interior, y, quad1);
00149 Expr zCMExpr = Integral(interior, z, quad1);
00150 Expr volExpr = Integral(interior, 1.0, quad1);
00151
00152 double vol = evaluateIntegral(mesh, volExpr);
00153 double xCM = evaluateIntegral(mesh, xCMExpr)/vol;
00154 double yCM = evaluateIntegral(mesh, yCMExpr)/vol;
00155 double zCM = evaluateIntegral(mesh, zCMExpr)/vol;
00156 Out::root() << "centroid = (" << xCM << ", " << yCM
00157 << ", " << zCM << ")" << std::endl;
00158
00159
00160
00161 Expr r = sqrt(x*x + y*y);
00162 Expr sinPhi = y/r;
00163
00164
00165 QuadratureFamily quad4 = new GaussianQuadrature(4);
00166
00167 Expr fourierSin1Expr = Integral(hole, sinPhi*soln, quad4);
00168 Expr fourierDenomExpr = Integral(hole, sinPhi*sinPhi, quad2);
00169 double fourierSin1 = evaluateIntegral(mesh, fourierSin1Expr);
00170 double fourierDenom = evaluateIntegral(mesh, fourierDenomExpr);
00171 Out::root() << "fourier sin m=1 = " << fourierSin1/fourierDenom << std::endl;
00172
00173
00174 Expr L2NormExpr = Integral(interior, soln*soln, quad2);
00175 double l2Norm_method1 = sqrt(evaluateIntegral(mesh, L2NormExpr));
00176 Out::os() << "method #1: ||soln|| = " << l2Norm_method1 << endl;
00177
00178
00179 double l2Norm_method2 = L2Norm(mesh, interior, soln, quad2);
00180 Out::os() << "method #2: ||soln|| = " << l2Norm_method2 << endl;
00181
00182
00183
00184
00185
00186
00187
00188 Sundance::passFailTest(fabs(flux), 1.0e-2);
00189 }
00190 catch(std::exception& e)
00191 {
00192 Sundance::handleException(e);
00193 }
00194 Sundance::finalize();
00195
00196 return Sundance::testStatus();
00197 }
00198
00199
00200
00201
00202
00203 #else
00204
00205
00206 int main(int argc, char** argv)
00207 {
00208 Sundance::init(&argc, &argv);
00209 std::cout << "dummy Laplace3D PASSED. Enable exodus to run the actual test" << std::endl;
00210 Sundance::finalize();
00211 return 0;
00212 }
00213
00214
00215 #endif
00216