1 #include "StationaryDiffusionEquation.hxx"
9 int main(int argc, char** argv)
24 /* Mesh construction : splitting policy to 0 yield all nodes considered boundary nodes */
25 Mesh M(xinf,xsup,nx,yinf,ysup,ny,zinf,zsup,nz,1); //Regular tetrahadral mesh
27 /* set the limit field for each boundary */
29 M.setGroupAtPlan(xsup,0,eps,"Bord1");
30 M.setGroupAtPlan(xinf,0,eps,"Bord2");
31 M.setGroupAtPlan(ysup,1,eps,"Bord3");
32 M.setGroupAtPlan(yinf,1,eps,"Bord4");
33 M.setGroupAtPlan(zsup,2,eps,"Bord5");
34 M.setGroupAtPlan(zinf,2,eps,"Bord6");
36 /* set the boundary values for each boundary */
44 cout<< "Built a regular tetrahedral 3D mesh from a cube mesh with "<< nx<<"x" <<ny<<"x" <<nz<< " cells"<<endl;
46 /* Create the problem */
47 bool FEComputation=false;
48 StationaryDiffusionEquation myProblem(spaceDim,FEComputation);
51 /* set the boundary conditions */
52 myProblem.setDirichletBoundaryCondition("Bord1",T1);
53 myProblem.setDirichletBoundaryCondition("Bord2",T2);
54 myProblem.setDirichletBoundaryCondition("Bord3",T3);
55 myProblem.setDirichletBoundaryCondition("Bord4",T4);
56 myProblem.setDirichletBoundaryCondition("Bord5",T5);
57 myProblem.setDirichletBoundaryCondition("Bord6",T6);
59 /* Set the right hand side function*/
60 Field my_RHSfield("RHS_field", CELLS, M, 1);
63 for(int i=0; i< M.getNumberOfCells(); i++)
70 my_RHSfield[i]=2*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z);//mettre la fonction definie au second membre de l'edp
72 myProblem.setHeatPowerField(my_RHSfield);
73 myProblem.setLinearSolver(GMRES,ILU);
75 /* name the result file */
76 string fileName = "StationaryDiffusion_3DFV_StructuredTetrahedra";
77 myProblem.setFileName(fileName);
79 /* Run the computation */
80 myProblem.initialize();
81 bool ok = myProblem.solveStationaryProblem();
83 cout << "Simulation of "<<fileName<<" failed !" << endl;
86 /********************** Postprocessing and measure od numerical error******************************/
87 Field my_ResultField = myProblem.getOutputTemperatureField();
88 /* The following formulas use the fact that the exact solution is equal the right hand side divided by 3*pi*pi */
89 double max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(3*pi*pi);
90 double max_sol_num=my_ResultField.max();
91 double min_sol_num=my_ResultField.min();
93 for(int i=0; i< M.getNumberOfNodes() ; i++)
94 if( erreur_abs < abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i]) )
95 erreur_abs = abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i]);
97 cout<<"Absolute error = max(| exact solution - numerical solution |) = "<<erreur_abs <<endl;
98 cout<<"Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = "<<erreur_abs/max_abs_sol_exacte<<endl;
99 cout<<"Maximum numerical solution = "<< max_sol_num<< " Minimum numerical solution = "<< min_sol_num << endl;
101 assert( erreur_abs/max_abs_sol_exacte <1.);
103 cout << "Simulation of "<<fileName<<" is successful ! " << endl;
105 cout << "------------ End of calculation !!! -----------" << endl;
106 myProblem.terminate();