CoreFlows/examples/C/StationaryDiffusionEquation_2DFV_StructuredTriangles_Neumann.cxx

   1 #include "StationaryDiffusionEquation.hxx"
   2 #include "Node.hxx"
   3 #include "math.h"
   4 #include <assert.h>
   5
   6 double pi = M_PI;
   7 using namespace std;
   8
   9 int main(int argc, char** argv)
  10 {
  11         int spaceDim = 2;
  12
  13         /* Mesh data */
  14         double xinf=0.0;
  15         double xsup=1.0;
  16         double yinf=0.0;
  17         double ysup=1.0;
  18         int nx=20;
  19         int ny=20;
  20
  21     /* Mesh construction */
  22         Mesh M(xinf,xsup,nx,yinf,ysup,ny,0); //Regular triangular mesh
  23
  24         /* set the limit field for each boundary */
  25         double eps=1e-6;
  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");
  30
  31     /* set the boundary values for each boundary */
  32         double T1=0;
  33         double T2=0;
  34         double T3=0;
  35         double T4=0;
  36
  37         cout<< "Building of a regular triangular 2D mesh from a square mesh with "<< nx<<"x" <<ny<< " cells"<<endl;
  38
  39     /* Create the problem */
  40     bool FEComputation=false;
  41         StationaryDiffusionEquation myProblem(spaceDim,FEComputation);
  42         myProblem.setMesh(M);
  43
  44     /* set the boundary conditions */
  45         myProblem.setNeumannBoundaryCondition("Bord1");
  46         myProblem.setNeumannBoundaryCondition("Bord2");
  47         myProblem.setNeumannBoundaryCondition("Bord3");
  48         myProblem.setNeumannBoundaryCondition("Bord4");
  49
  50         /* Set the right hand side function*/
  51         Field my_RHSfield("RHS_field", CELLS, M, 1);
  52     Cell Ci;
  53     double x, y;
  54         for(int i=0; i< M.getNumberOfCells(); i++)
  55     {
  56                 Ci= M.getCell(i);
  57                 x = Ci.x();
  58                 y = Ci.y();
  59
  60                 my_RHSfield[i]=2*pi*pi*cos(pi*x)*cos(pi*y);//mettre la fonction definie au second membre de l'edp
  61         }
  62         myProblem.setHeatPowerField(my_RHSfield);
  63         myProblem.setLinearSolver(GMRES,ILU);
  64
  65     /* name the result file */
  66         string fileName = "StationnaryDiffusion_2DFV_RegularTriangles_Neumann";
  67         myProblem.setFileName(fileName);
  68
  69         /* Run the computation */
  70         myProblem.initialize();
  71         bool ok = myProblem.solveStationaryProblem();
  72         if (!ok)
  73                 cout << "Simulation of "<<fileName<<" failed !" << endl;
  74         else
  75         {
  76                 /********************** Postprocessing and measure od numerical error******************************/
  77                 Field my_ResultField = myProblem.getOutputTemperatureField();
  78                 /* The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi */
  79                 double max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi);
  80                 double max_sol_num=my_ResultField.max();
  81                 double min_sol_num=my_ResultField.min();
  82                 double erreur_abs=0;
  83                 for(int i=0; i< M.getNumberOfCells() ; i++)
  84                         if( erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) )
  85                                 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]);
  86
  87                 cout<<"Absolute error = max(| exact solution - numerical solution |) = "<<erreur_abs <<endl;
  88                 cout<<"Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = "<<erreur_abs/max_abs_sol_exacte<<endl;
  89                 cout<<"Maximum numerical solution = "<< max_sol_num<< " Minimum numerical solution = "<< min_sol_num << endl;
  90
  91                 assert( erreur_abs/max_abs_sol_exacte <1.);
  92
  93         cout << "Simulation of "<<fileName<<" is successful ! " << endl;
  94     }
  95
  96         cout << "------------ End of calculation !!! -----------" << endl;
  97         myProblem.terminate();
  98
  99         return EXIT_SUCCESS;
 100 }