CoreFlows/examples/C/StationaryDiffusionEquation_3DFV_StructuredTetrahedra.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 = 3;
  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         double zinf=0.0;
  19         double zsup=1.0;
  20         int nx=20;
  21         int ny=20;
  22         int nz=20;
  23
  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
  26
  27         /* set the limit field for each boundary */
  28         double eps=1e-6;
  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");
  35
  36     /* set the boundary values for each boundary */
  37         double T1=0;
  38         double T2=0;
  39         double T3=0;
  40         double T4=0;
  41         double T5=0;
  42         double T6=0;
  43
  44         cout<< "Built a regular tetrahedral 3D mesh from a cube mesh with "<< nx<<"x" <<ny<<"x" <<nz<< " cells"<<endl;
  45
  46     /* Create the problem */
  47     bool FEComputation=false;
  48         StationaryDiffusionEquation myProblem(spaceDim,FEComputation);
  49         myProblem.setMesh(M);
  50
  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);
  58
  59         /* Set the right hand side function*/
  60         Field my_RHSfield("RHS_field", CELLS, M, 1);
  61     Cell Ci;
  62     double x, y, z;
  63         for(int i=0; i< M.getNumberOfCells(); i++)
  64     {
  65                 Ci= M.getCell(i);
  66                 x = Ci.x();
  67                 y = Ci.y();
  68                 z = Ci.z();
  69
  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
  71         }
  72         myProblem.setHeatPowerField(my_RHSfield);
  73         myProblem.setLinearSolver(GMRES,ILU);
  74
  75     /* name the result file */
  76         string fileName = "StationaryDiffusion_3DFV_StructuredTetrahedra";
  77         myProblem.setFileName(fileName);
  78
  79         /* Run the computation */
  80         myProblem.initialize();
  81         bool ok = myProblem.solveStationaryProblem();
  82         if (!ok)
  83                 cout << "Simulation of "<<fileName<<" failed !" << endl;
  84         else
  85         {
  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();
  92                 double erreur_abs=0;
  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]);
  96
  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;
 100
 101                 assert( erreur_abs/max_abs_sol_exacte <1.);
 102
 103         cout << "Simulation of "<<fileName<<" is successful ! " << endl;
 104     }
 105         cout << "------------ End of calculation !!! -----------" << endl;
 106         myProblem.terminate();
 107
 108         return EXIT_SUCCESS;
 109 }