1 #include "StationaryDiffusionEquation.hxx"
9 int main(int argc, char** argv)
21 /* Mesh construction */
22 Mesh M(xinf,xsup,nx,yinf,ysup,ny,0); //Regular triangular mesh
24 /* set the limit field for each boundary */
26 M.setGroupAtPlan(xsup,0,eps,"Bord1");
27 M.setGroupAtPlan(xinf,0,eps,"Bord2");
28 M.setGroupAtPlan(ysup,1,eps,"Bord3");
29 M.setGroupAtPlan(yinf,1,eps,"Bord4");
31 cout<< "Built a regular triangular 2D mesh from a square mesh with "<< nx<<"x" <<ny<< " cells"<<endl;
33 /* Create the problem */
34 bool FEComputation=true;
35 StationaryDiffusionEquation myProblem(spaceDim,FEComputation);
38 /* set the boundary conditions */
39 myProblem.setNeumannBoundaryCondition("Bord1");
40 myProblem.setNeumannBoundaryCondition("Bord2");
41 myProblem.setNeumannBoundaryCondition("Bord3");
42 myProblem.setNeumannBoundaryCondition("Bord4");
44 /* Set the right hand side function*/
45 Field my_RHSfield("RHS_field", NODES, M, 1);
48 for(int i=0; i< M.getNumberOfNodes(); i++)
54 my_RHSfield[i]=2*pi*pi*cos(pi*x)*cos(pi*y);//mettre la fonction definie au second membre de l'edp
56 myProblem.setHeatPowerField(my_RHSfield);
57 myProblem.setLinearSolver(GMRES,ILU);
59 /* name the result file */
60 string fileName = "StationnaryDiffusion_2DFE_StructuredTriangles_Neumann";
61 myProblem.setFileName(fileName);
63 /* Run the computation */
64 myProblem.initialize();
65 bool ok = myProblem.solveStationaryProblem();
67 cout << "Simulation of "<<fileName<<" failed !" << endl;
70 /********************** Postprocessing and measure od numerical error******************************/
71 Field my_ResultField = myProblem.getOutputTemperatureField();
72 /* The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi */
73 double max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi);
74 double max_sol_num=my_ResultField.max();
75 double min_sol_num=my_ResultField.min();
77 for(int i=0; i< M.getNumberOfNodes() ; i++)
78 if( erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) )
79 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]);
81 cout<<"Absolute error = max(| exact solution - numerical solution |) = "<<erreur_abs <<endl;
82 cout<<"Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = "<<erreur_abs/max_abs_sol_exacte<<endl;
83 cout<<"Maximum numerical solution = "<< max_sol_num<< " Minimum numerical solution = "<< min_sol_num << endl;
85 assert( erreur_abs/max_abs_sol_exacte <1.);
87 cout << "Simulation of "<<fileName<<" is successful ! " << endl;
89 cout << "------------ End of calculation !!! -----------" << endl;
90 myProblem.terminate();