CoreFlows/examples/C/StationaryDiffusionEquation_2DEF_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=3;
  19         int ny=3;
  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         cout<< "Built a regular triangular 2D mesh from a square mesh with "<< nx<<"x" <<ny<< " cells"<<endl;
  32
  33     /* Create the problem */
  34     bool FEComputation=true;
  35         StationaryDiffusionEquation myProblem(spaceDim,FEComputation);
  36         myProblem.setMesh(M);
  37
  38     /* set the boundary conditions */
  39         myProblem.setNeumannBoundaryCondition("Bord1");
  40         myProblem.setNeumannBoundaryCondition("Bord2");
  41         myProblem.setNeumannBoundaryCondition("Bord3");
  42         myProblem.setNeumannBoundaryCondition("Bord4");
  43
  44         /* Set the right hand side function*/
  45         Field my_RHSfield("RHS_field", NODES, M, 1);
  46     Node Ni;
  47     double x, y;
  48         for(int i=0; i< M.getNumberOfNodes(); i++)
  49     {
  50                 Ni= M.getNode(i);
  51                 x = Ni.x();
  52                 y = Ni.y();
  53
  54                 my_RHSfield[i]=2*pi*pi*cos(pi*x)*cos(pi*y);//mettre la fonction definie au second membre de l'edp
  55         }
  56         myProblem.setHeatPowerField(my_RHSfield);
  57         myProblem.setLinearSolver(GMRES,ILU);
  58
  59     /* name the result file */
  60         string fileName = "StationnaryDiffusion_2DFE_StructuredTriangles_Neumann";
  61         myProblem.setFileName(fileName);
  62
  63         /* Run the computation */
  64         myProblem.initialize();
  65         bool ok = myProblem.solveStationaryProblem();
  66         if (!ok)
  67                 cout << "Simulation of "<<fileName<<" failed !" << endl;
  68         else
  69         {
  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();
  76                 double erreur_abs=0;
  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]);
  80
  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;
  84
  85                 assert( erreur_abs/max_abs_sol_exacte <1.);
  86
  87         cout << "Simulation of "<<fileName<<" is successful ! " << endl;
  88     }
  89         cout << "------------ End of calculation !!! -----------" << endl;
  90         myProblem.terminate();
  91
  92         return EXIT_SUCCESS;
  93 }