CoreFlows/Models/src/StationaryDiffusionEquation.cxx

   1 #include "StationaryDiffusionEquation.hxx"
   2 #include "Node.hxx"
   3 #include "SparseMatrixPetsc.hxx"
   4 #include "math.h"
   5 #include <algorithm>
   6 #include <fstream>
   7 #include <sstream>
   8
   9 using namespace std;
  10
  11 int StationaryDiffusionEquation::fact(int n)
  12 {
  13   return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
  14 }
  15 int StationaryDiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const
  16 {//assumes Dirichlet node numbering is strictly increasing
  17     int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
  18     int boundarySize=dirichletNodes.size();
  19     while(j<boundarySize and dirichletNodes[j]<globalIndex)
  20         j++;
  21     if(j==boundarySize)
  22         return globalIndex-boundarySize;
  23     else if (dirichletNodes[j]>globalIndex)
  24         return globalIndex-j;
  25     else
  26         throw CdmathException("StationaryDiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
  27 }
  28
  29 int StationaryDiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes) const
  30 {//assumes Dirichlet boundary node numbering is strictly increasing
  31     int boundarySize=dirichletNodes.size();
  32     /* trivial case where all boundary nodes are Neumann BC */
  33     if(boundarySize==0)
  34         return unknownNodeIndex;
  35
  36     double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
  37     int j=0;//indice de parcours des noeuds frontière
  38     //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
  39     while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
  40     {
  41         unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
  42         j++;
  43     }
  44
  45     if(j+1==boundarySize)
  46         return unknownNodeIndex+boundarySize;
  47     else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
  48         return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
  49 }
  50
  51 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda){
  52         PetscBool petscInitialized;
  53         PetscInitialized(&petscInitialized);
  54         if(!petscInitialized)
  55                 PetscInitialize(NULL,NULL,0,0);
  56
  57     if(lambda < 0.)
  58     {
  59         std::cout<<"Conductivity="<<lambda<<endl;
  60         throw CdmathException("Error : conductivity parameter lambda cannot  be negative");
  61     }
  62     if(dim<=0)
  63     {
  64         std::cout<<"Space dimension="<<dim<<endl;
  65         throw CdmathException("Error : parameter dim cannot  be negative");
  66     }
  67
  68     _FECalculation=FECalculation;
  69     _onlyNeumannBC=false;
  70
  71         _Ndim=dim;
  72         _nVar=1;//scalar prolem
  73         _dt_src=0;
  74         _diffusionMatrixSet=false;
  75         _initializedMemory=false;
  76
  77     //Mesh data
  78     _neibMaxNbCells=0;
  79     _meshSet=false;
  80     _neibMaxNbNodes=0;
  81
  82     //Boundary conditions
  83     _boundaryNodeIds=std::vector< int >(0);
  84     _dirichletNodeIds=std::vector< int >(0);
  85     _NboundaryNodes=0;
  86     _NdirichletNodes=0;
  87     _NunknownNodes=0;
  88     _dirichletValuesSet=false;
  89     _neumannValuesSet=false;
  90
  91     //Linear solver data
  92         _precision=1.e-6;
  93         _precision_Newton=_precision;
  94         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  95         _maxPetscIts=50;
  96         _maxNewtonIts=50;
  97         _NEWTON_its=0;
  98         int _PetscIts=0;//the number of iterations of the linear solver
  99         _ksptype = (char*)&KSPGMRES;
 100         _pctype = (char*)&PCLU;
 101         _conditionNumber=false;
 102         _erreur_rel= 0;
 103
 104     //parameters for monitoring simulation
 105         _verbose = false;
 106         _system = false;
 107         _runLogFile=new ofstream;
 108
 109         //result save parameters
 110         _fileName = "StationaryDiffusionProblem";
 111         char result[ PATH_MAX ];//extracting current directory
 112         getcwd(result, PATH_MAX );
 113         _path=string( result );
 114         _saveFormat=VTK;
 115     _computationCompletedSuccessfully=false;
 116
 117     //heat transfer parameters
 118         _conductivity=lambda;
 119         _fluidTemperatureFieldSet=false;
 120         _fluidTemperature=0;
 121         _heatPowerFieldSet=false;
 122         _heatTransfertCoeff=0;
 123         _heatSource=0;
 124 }
 125
 126 void StationaryDiffusionEquation::initialize()
 127 {
 128         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 129
 130         if(!_meshSet)
 131                 throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
 132         else
 133     {
 134                 cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
 135         *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
 136         if(!_FECalculation)
 137         {
 138             cout<< "Finite volumes method"<<endl<<endl;
 139             *_runLogFile<< "Finite volumes method"<<endl<<endl;
 140         }
 141         else
 142         {
 143             cout<< "Finite elements method"<<endl<<endl;
 144             *_runLogFile<< "Finite elements method"<<endl<<endl;
 145         }
 146     }
 147
 148         _DiffusionTensor=Matrix(_Ndim);
 149         for(int idim=0;idim<_Ndim;idim++)
 150                 _DiffusionTensor(idim,idim)=1;
 151         /**************** Field creation *********************/
 152
 153         if(!_heatPowerFieldSet){
 154         if(_FECalculation){
 155             _heatPowerField=Field("Heat power",NODES,_mesh,1);
 156             for(int i =0; i<_Nnodes; i++)
 157                 _heatPowerField(i) = _heatSource;
 158         }
 159         else{
 160             _heatPowerField=Field("Heat power",CELLS,_mesh,1);
 161             for(int i =0; i<_Nmailles; i++)
 162                 _heatPowerField(i) = _heatSource;
 163         }
 164         _heatPowerFieldSet=true;
 165     }
 166         if(!_fluidTemperatureFieldSet){
 167         if(_FECalculation){
 168             _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
 169             for(int i =0; i<_Nnodes; i++)
 170                 _fluidTemperatureField(i) = _fluidTemperature;
 171         }
 172         else{
 173             _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
 174             for(int i =0; i<_Nmailles; i++)
 175                 _fluidTemperatureField(i) = _fluidTemperature;
 176         }
 177         _fluidTemperatureFieldSet=true;
 178         }
 179
 180     /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 181     if(_FECalculation)
 182     {
 183         if(_NboundaryNodes==_Nnodes)
 184             cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
 185
 186         for(int i=0; i<_NboundaryNodes; i++)
 187         {
 188             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
 189             if( it != _dirichletBoundaryValues.end() )
 190                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 191             else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
 192             {
 193                 cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
 194                 *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
 195                 _runLogFile->close();
 196                 throw CdmathException("Missing boundary group");
 197             }
 198             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
 199             {
 200                 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
 201                 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
 202                 *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
 203                 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
 204                 _runLogFile->close();
 205                 throw CdmathException("Missing boundary condition");
 206             }
 207             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
 208                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 209             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
 210             {
 211                 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
 212                 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
 213                 *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
 214                 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
 215                 _runLogFile->close();
 216                 throw CdmathException("Wrong boundary condition");
 217             }
 218         }
 219         _NdirichletNodes=_dirichletNodeIds.size();
 220         _NunknownNodes=_Nnodes - _NdirichletNodes;
 221         cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 222                 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 223     }
 224
 225         //creation de la matrice
 226     if(!_FECalculation)
 227         MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
 228     else
 229         MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
 230
 231         VecCreate(PETSC_COMM_SELF, &_Tk);
 232
 233     if(!_FECalculation)
 234         VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
 235     else
 236         VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
 237
 238         VecSetFromOptions(_Tk);
 239         VecDuplicate(_Tk, &_Tkm1);
 240         VecDuplicate(_Tk, &_deltaT);
 241         VecDuplicate(_Tk, &_b);//RHS of the linear system
 242
 243         //Linear solver
 244         KSPCreate(PETSC_COMM_SELF, &_ksp);
 245         KSPSetType(_ksp, _ksptype);
 246         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 247         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 248         KSPGetPC(_ksp, &_pc);
 249         PCSetType(_pc, _pctype);
 250
 251     //Checking whether all boundary conditions are Neumann boundary condition
 252     //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
 253     if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
 254     {
 255         map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
 256         while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
 257             it++;
 258         _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
 259     }
 260     else
 261         if(_FECalculation)
 262             _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
 263         else
 264             _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
 265
 266     //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
 267     if(_onlyNeumannBC)
 268     {
 269         std::cout<<"### Warning : all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 270         std::cout<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
 271         std::cout<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 272         *_runLogFile<<"### Warning : all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 273         *_runLogFile<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 274         *_runLogFile<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
 275
 276                 MatNullSpace nullsp;
 277                 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
 278                 MatSetNullSpace(_A, nullsp);
 279                 MatSetTransposeNullSpace(_A, nullsp);
 280                 MatNullSpaceDestroy(&nullsp);
 281                 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
 282                 //PCFactorSetShiftAmount(_pc,1e-10);
 283     }
 284         _initializedMemory=true;
 285 }
 286
 287 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
 288         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 289                 computeDiffusionMatrix(stop);
 290
 291         _dt_src=computeRHS(stop);
 292         stop=false;
 293         return _dt_src;
 294 }
 295
 296 Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
 297     vector< Matrix > matrices(_Ndim);
 298
 299     for (int idim=0; idim<_Ndim;idim++){
 300         matrices[idim]=M.deepCopy();
 301         for (int jdim=0; jdim<_Ndim+1;jdim++)
 302                         matrices[idim](jdim,idim) = values[jdim] ;
 303     }
 304
 305         Vector result(_Ndim);
 306     for (int idim=0; idim<_Ndim;idim++)
 307         result[idim] = matrices[idim].determinant();
 308
 309         return result;
 310 }
 311
 312 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
 313 {
 314     double result;
 315
 316     if(_FECalculation)
 317         result=computeDiffusionMatrixFE(stop);
 318     else
 319         result=computeDiffusionMatrixFV(stop);
 320
 321     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 322     //update value here if variable  heat transfer coefficient
 323     if(_heatTransfertCoeff>_precision)
 324         MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
 325
 326     if(_verbose or _system)
 327         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 328
 329     return  result;
 330 }
 331
 332 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 333         Cell Cj;
 334         string nameOfGroup;
 335         double dn, coeff;
 336         MatZeroEntries(_A);
 337         VecZeroEntries(_b);
 338
 339     Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 340     std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 341     std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 342     std::vector< Node > nodes(_Ndim+1);//cell nodes
 343     int i_int, j_int; //index of nodes j and k considered as unknown nodes
 344     bool dirichletCell_treated;
 345
 346     std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 347     for (int idim=0; idim<_Ndim+1;idim++)
 348         values[idim][idim]=1;
 349
 350     /* parameters for boundary treatment */
 351     vector< double > valuesBorder(_Ndim+1);
 352     Vector GradShapeFuncBorder(_Ndim+1);
 353
 354         for (int j=0; j<_Nmailles;j++)
 355     {
 356                 Cj = _mesh.getCell(j);
 357
 358         for (int idim=0; idim<_Ndim+1;idim++){
 359             nodeIds[idim]=Cj.getNodeId(idim);
 360             nodes[idim]=_mesh.getNode(nodeIds[idim]);
 361             for (int jdim=0; jdim<_Ndim;jdim++)
 362                 M(idim,jdim)=nodes[idim].getPoint()[jdim];
 363             M(idim,_Ndim)=1;
 364         }
 365         for (int idim=0; idim<_Ndim+1;idim++)
 366             GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
 367
 368         /* Loop on the edges of the cell */
 369         for (int idim=0; idim<_Ndim+1;idim++)
 370         {
 371             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 372             {//First node of the edge is not Dirichlet node
 373                 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 374                 dirichletCell_treated=false;
 375                 for (int jdim=0; jdim<_Ndim+1;jdim++)
 376                 {
 377                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 378                     {//Second node of the edge is not Dirichlet node
 379                         j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 380                         MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 381                     }
 382                     else if (!dirichletCell_treated)
 383                     {//Second node of the edge is a Dirichlet node
 384                         dirichletCell_treated=true;
 385                         for (int kdim=0; kdim<_Ndim+1;kdim++)
 386                         {
 387                                                         std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 388                                                         if( it != _dirichletBoundaryValues.end() )
 389                             {
 390                                 if( _dirichletValuesSet )
 391                                     valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 392                                 else
 393                                     valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
 394                             }
 395                             else
 396                                 valuesBorder[kdim]=0;
 397                         }
 398                         GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
 399                         coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 400                         VecSetValue(_b,i_int,coeff, ADD_VALUES);
 401                     }
 402                 }
 403             }
 404         }
 405         }
 406
 407     //Calcul de la contribution de la condition limite de Neumann au second membre
 408     if( _NdirichletNodes !=_NboundaryNodes)
 409     {
 410         vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
 411         int NboundaryFaces=boundaryFaces.size();
 412         for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
 413         {
 414             Face Fi = _mesh.getFace(i);
 415             for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
 416             {
 417                 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
 418                 {
 419                     j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 420                     if( _neumannValuesSet )
 421                         coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
 422                     else
 423                         coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
 424                     VecSetValue(_b, j_int, coeff, ADD_VALUES);
 425                 }
 426             }
 427         }
 428     }
 429     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 430         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 431         VecAssemblyBegin(_b);
 432         VecAssemblyEnd(_b);
 433
 434         if(_onlyNeumannBC)      //Check that the matrix is symmetric
 435         {
 436                 PetscBool isSymetric;
 437                 MatIsSymmetric(_A,_precision,&isSymetric);
 438                 if(!isSymetric)
 439                 {
 440                         cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 441                         throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 442                 }
 443         }
 444
 445         _diffusionMatrixSet=true;
 446     stop=false ;
 447
 448         return INFINITY;
 449 }
 450
 451 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 452         long nbFaces = _mesh.getNumberOfFaces();
 453         Face Fj;
 454         Cell Cell1,Cell2;
 455         string nameOfGroup;
 456         double inv_dxi, inv_dxj;
 457         double barycenterDistance;
 458         Vector normale(_Ndim);
 459         double dn;
 460         PetscInt idm, idn;
 461         std::vector< int > idCells;
 462         MatZeroEntries(_A);
 463         VecZeroEntries(_b);
 464
 465         for (int j=0; j<nbFaces;j++){
 466                 Fj = _mesh.getFace(j);
 467
 468                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 469                 idCells = Fj.getCellsId();
 470                 Cell1 = _mesh.getCell(idCells[0]);
 471                 idm = idCells[0];
 472         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 473             if (j == Cell1.getFacesId()[l]){
 474                 for (int idim = 0; idim < _Ndim; ++idim)
 475                     normale[idim] = Cell1.getNormalVector(l,idim);
 476                 break;
 477             }
 478         }
 479
 480                 //Compute velocity at the face Fj
 481                 dn=_conductivity*(_DiffusionTensor*normale)*normale;
 482
 483                 // compute 1/dxi = volume of Ci/area of Fj
 484         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 485
 486                 // If Fj is on the boundary
 487                 if (Fj.getNumberOfCells()==1) {
 488                         if(_verbose )
 489                         {
 490                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 491                                 for(int p=0; p<_Ndim; p++)
 492                                         cout << normale[p] << ",";
 493                                 cout << ") "<<endl;
 494                         }
 495
 496             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
 497             if( it != _dirichletBoundaryValues.end() )
 498             {
 499                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 500                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                                     , ADD_VALUES);
 501                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
 502             }
 503             else
 504             {
 505                 nameOfGroup = Fj.getGroupName();
 506
 507                 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
 508                     VecSetValue(_b,idm,    -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 509                 }
 510                 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
 511                     barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 512                     MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 513                     VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 514                 }
 515                 else {//Problem for 3D tetrahera
 516                                         //cout<< "Fj.getGroupNames().size()= "<<Fj.getGroupNames().size()<<", Fj.x()= "<<Fj.x()<<", Fj.y()= "<<Fj.y()<<", Fj.z()= "<<Fj.z()<<endl;
 517                     stop=true ;
 518                     cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 519                     cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 520                     cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
 521                     *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 522                     _runLogFile->close();
 523                     throw CdmathException("Boundary condition not set");
 524                 }
 525             }
 526                         // if Fj is inside the domain
 527                 } else  if (Fj.getNumberOfCells()==2 ){
 528                         if(_verbose )
 529                         {
 530                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 531                                 cout << " ; vecteur normal=(";
 532                                 for(int p=0; p<_Ndim; p++)
 533                                         cout << normale[p] << ",";
 534                                 cout << ") "<<endl;
 535                         }
 536                         Cell2 = _mesh.getCell(idCells[1]);
 537                         idn = idCells[1];
 538                         if (_Ndim > 1)
 539                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 540                         else
 541                                 inv_dxj = 1/Cell2.getMeasure();
 542
 543                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 544
 545                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 546                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 547                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 548                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 549                 }
 550                 else
 551         {
 552             *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
 553                         throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 554         }
 555         }
 556
 557         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 558         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 559         VecAssemblyBegin(_b);
 560         VecAssemblyEnd(_b);
 561
 562         if(_onlyNeumannBC)      //Check that the matrix is symmetric
 563         {
 564                 PetscBool isSymetric;
 565                 MatIsSymmetric(_A,_precision,&isSymetric);
 566                 if(!isSymetric)
 567                 {
 568                         cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 569                         throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 570                 }
 571         }
 572
 573         _diffusionMatrixSet=true;
 574     stop=false ;
 575
 576     return INFINITY;
 577 }
 578
 579 double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix)
 580 {
 581         VecAssemblyBegin(_b);
 582
 583     if(!_FECalculation)
 584         for (int i=0; i<_Nmailles;i++)
 585             VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
 586     else
 587         {
 588             Cell Ci;
 589             std::vector< int > nodesId;
 590             for (int i=0; i<_Nmailles;i++)
 591             {
 592                 Ci=_mesh.getCell(i);
 593                 nodesId=Ci.getNodesId();
 594                 for (int j=0; j<nodesId.size();j++)
 595                     if(!_mesh.isBorderNode(nodesId[j]))
 596                     {
 597                         double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 598                         VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 599                     }
 600             }
 601         }
 602
 603         VecAssemblyEnd(_b);
 604
 605     if(_verbose or _system)
 606         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 607
 608     stop=false ;
 609         return INFINITY;
 610 }
 611
 612 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
 613 {
 614         bool stop;
 615
 616     //Only implicit scheme considered
 617     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 618     MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 619
 620 #if PETSC_VERSION_GREATER_3_5
 621     KSPSetOperators(_ksp, _A, _A);
 622 #else
 623     KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 624 #endif
 625
 626     if(_conditionNumber)
 627         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 628     KSPSolve(_ksp, _b, _Tk);
 629
 630         KSPConvergedReason reason;
 631         KSPGetConvergedReason(_ksp,&reason);
 632     KSPGetIterationNumber(_ksp, &_PetscIts);
 633     double residu;
 634     KSPGetResidualNorm(_ksp,&residu);
 635         if (reason!=2 and reason!=3)
 636     {
 637         cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 638         cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 639         cout<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 640                 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 641         *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 642         *_runLogFile<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 643                 _runLogFile->close();
 644         converged = false;
 645         stop = true;
 646     }
 647     else{
 648         if( _MaxIterLinearSolver < _PetscIts)
 649             _MaxIterLinearSolver = _PetscIts;
 650         cout<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<<  _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
 651                 *_runLogFile<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<<  _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
 652         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 653         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 654
 655         _erreur_rel= 0;
 656         double Ti, dTi;
 657
 658         VecAssemblyBegin(_Tk);
 659         VecAssemblyEnd(  _Tk);
 660         VecAssemblyBegin(_deltaT);
 661         VecAssemblyEnd(  _deltaT);
 662
 663         if(!_FECalculation)
 664             for(int i=0; i<_Nmailles; i++)
 665             {
 666                 VecGetValues(_deltaT, 1, &i, &dTi);
 667                 VecGetValues(_Tk, 1, &i, &Ti);
 668                 if(_erreur_rel < fabs(dTi/Ti))
 669                     _erreur_rel = fabs(dTi/Ti);
 670             }
 671         else
 672             for(int i=0; i<_NunknownNodes; i++)
 673             {
 674                 VecGetValues(_deltaT, 1, &i, &dTi);
 675                 VecGetValues(_Tk, 1, &i, &Ti);
 676                 if(_erreur_rel < fabs(dTi/Ti))
 677                     _erreur_rel = fabs(dTi/Ti);
 678             }
 679         stop=false;
 680         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 681     }
 682
 683         VecCopy(_Tk, _Tkm1);
 684
 685         return stop;
 686 }
 687
 688 void StationaryDiffusionEquation::setMesh(const Mesh &M)
 689 {
 690         if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
 691         {
 692         cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
 693                 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
 694                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 695                 _runLogFile->close();
 696                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect  dimension");
 697         }
 698
 699         _mesh=M;
 700         _Nmailles = _mesh.getNumberOfCells();
 701         _Nnodes =   _mesh.getNumberOfNodes();
 702
 703     cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
 704         *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
 705
 706         // find  maximum nb of neibourghs
 707     if(!_FECalculation)
 708     {
 709         _VV=Field ("Temperature", CELLS, _mesh, 1);
 710         _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 711     }
 712     else
 713     {
 714         if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 715                         cout<<"1D Finite element method on segments"<<endl;
 716         else if(_Ndim==2)
 717         {
 718                         if( _mesh.isTriangular() )//Mesh dim=2
 719                                 cout<<"2D Finite element method on triangles"<<endl;
 720                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 721                                 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 722                         else
 723                         {
 724                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 725                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 726                                 _runLogFile->close();
 727                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 728                         }
 729         }
 730         else if(_Ndim==3)
 731         {
 732                         if( _mesh.isTetrahedral() )//Mesh dim=3
 733                                 cout<<"3D Finite element method on tetrahedra"<<endl;
 734                         else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 735                                 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 736                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 737                                 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 738                         else
 739                         {
 740                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either tetrahedral, either a triangularised surface or 1D network"<<endl;
 741                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 742                                 _runLogFile->close();
 743                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 744                         }
 745         }
 746
 747                 _VV=Field ("Temperature", NODES, _mesh, 1);
 748
 749         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 750         _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 751         _NboundaryNodes=_boundaryNodeIds.size();
 752     }
 753
 754         _meshSet=true;
 755 }
 756
 757 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
 758 {
 759         //_maxPetscIts=maxIterationsPetsc;
 760         // set linear solver algorithm
 761         if (kspType==GMRES)
 762                 _ksptype = (char*)&KSPGMRES;
 763         else if (kspType==CG)
 764                 _ksptype = (char*)&KSPCG;
 765         else if (kspType==BCGS)
 766                 _ksptype = (char*)&KSPBCGS;
 767         else {
 768                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 769                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 770                 _runLogFile->close();
 771                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 772         }
 773         // set preconditioner
 774         if (pcType == NONE)
 775                 _pctype = (char*)&PCNONE;
 776         else if (pcType ==LU)
 777                 _pctype = (char*)&PCLU;
 778         else if (pcType == ILU)
 779                 _pctype = (char*)&PCILU;
 780         else if (pcType ==CHOLESKY)
 781                 _pctype = (char*)&PCCHOLESKY;
 782         else if (pcType == ICC)
 783                 _pctype = (char*)&PCICC;
 784         else {
 785                 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 786                 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 787                 _runLogFile->close();
 788                 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 789         }
 790 }
 791
 792 bool StationaryDiffusionEquation::solveStationaryProblem()
 793 {
 794         if(!_initializedMemory)
 795         {
 796                 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
 797                 _runLogFile->close();
 798                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 799         }
 800         bool stop=false; // Does the Problem want to stop (error) ?
 801         bool converged=false; // has the newton scheme converged (end) ?
 802
 803         cout<< "!!! Running test case "<< _fileName << " using ";
 804         *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
 805
 806     if(!_FECalculation)
 807     {
 808         cout<< "Finite volumes method"<<endl<<endl;
 809                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 810         }
 811     else
 812         {
 813         cout<< "Finite elements method"<<endl<<endl;
 814                 *_runLogFile<< "Finite elements method"<< endl<<endl;
 815         }
 816
 817     computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
 818     if (stop){
 819         cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 820         *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 821                 _runLogFile->close();
 822        throw CdmathException("Failed computing diffusion matrix");
 823     }
 824     computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
 825     if (stop){
 826         cout << "Error : failed computing right hand side, stopping calculation"<< endl;
 827         *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
 828         throw CdmathException("Failed computing right hand side");
 829     }
 830     stop = iterateNewtonStep(converged);
 831     if (stop){
 832         cout << "Error : failed solving linear system, stopping calculation"<< endl;
 833         *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
 834                 _runLogFile->close();
 835         throw CdmathException("Failed solving linear system");
 836     }
 837
 838     _computationCompletedSuccessfully=true;
 839     save();
 840
 841         // Newton iteration loop for non linear problems
 842     /*
 843         while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
 844         {
 845         computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
 846         computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
 847         stop = iterateNewtonStep(converged);
 848         _NEWTON_its++;
 849         }
 850     if (stop){
 851         cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 852         *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 853         throw CdmathException("Failed solving a Newton iteration");
 854     }
 855     else if(_NEWTON_its==_maxNewtonIts){
 856         cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 857         *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 858         throw CdmathException("No convergence of Newton scheme");
 859     }
 860     else{
 861         cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 862         *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 863         save();
 864     }
 865     */
 866
 867     *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
 868         _runLogFile->close();
 869
 870         return !stop;
 871 }
 872
 873 void StationaryDiffusionEquation::save(){
 874     cout<< "Saving numerical results"<<endl<<endl;
 875     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 876
 877         string resultFile(_path+"/StationaryDiffusionEquation");//Results
 878
 879         resultFile+="_";
 880         resultFile+=_fileName;
 881
 882         // create mesh and component info
 883     string suppress ="rm -rf "+resultFile+"_*";
 884     system(suppress.c_str());//Nettoyage des précédents calculs identiques
 885
 886     if(_verbose or _system)
 887         VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
 888
 889     double Ti;
 890     if(!_FECalculation)
 891         for(int i=0; i<_Nmailles; i++)
 892             {
 893                 VecGetValues(_Tk, 1, &i, &Ti);
 894                 _VV(i)=Ti;
 895             }
 896     else
 897     {
 898         int globalIndex;
 899         for(int i=0; i<_NunknownNodes; i++)
 900         {
 901             VecGetValues(_Tk, 1, &i, &Ti);
 902             globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 903             _VV(globalIndex)=Ti;
 904         }
 905
 906         Node Ni;
 907         string nameOfGroup;
 908         for(int i=0; i<_NdirichletNodes; i++)
 909         {
 910             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 911             nameOfGroup = Ni.getGroupName();
 912             _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 913         }
 914     }
 915
 916     _VV.setInfoOnComponent(0,"Temperature_(K)");
 917     switch(_saveFormat)
 918     {
 919         case VTK :
 920             _VV.writeVTK(resultFile);
 921             break;
 922         case MED :
 923             _VV.writeMED(resultFile);
 924             break;
 925         case CSV :
 926             _VV.writeCSV(resultFile);
 927             break;
 928     }
 929 }
 930 Field
 931 StationaryDiffusionEquation::getOutputTemperatureField()
 932 {
 933     if(!_computationCompletedSuccessfully)
 934         throw("Computation not performed yet or failed. No temperature field available");
 935     else
 936         return _VV;
 937 }
 938
 939 void StationaryDiffusionEquation::terminate()
 940 {
 941         VecDestroy(&_Tk);
 942         VecDestroy(&_Tkm1);
 943         VecDestroy(&_deltaT);
 944         VecDestroy(&_b);
 945         MatDestroy(&_A);
 946 }
 947 void
 948 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 949 {
 950     _dirichletValuesSet=true;
 951     _dirichletBoundaryValues=dirichletBoundaryValues;
 952 }
 953
 954 void
 955 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
 956 {
 957     _neumannValuesSet=true;
 958     _neumannBoundaryValues=neumannBoundaryValues;
 959 }
 960
 961 double
 962 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
 963 {
 964   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 965   return A.getConditionNumber( isSingular, tol);
 966 }
 967 std::vector< double >
 968 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
 969 {
 970   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 971
 972   if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
 973   {
 974       Vector nodal_volumes(_NunknownNodes);
 975       int j_int;
 976       for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
 977       {
 978         Cell Ci = _mesh.getCell(i);
 979         for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
 980         {
 981             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
 982                         {
 983                 j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 984                 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
 985             }
 986         }
 987        }
 988       for( j_int = 0; j_int< _NunknownNodes ; j_int++)
 989         nodal_volumes[j_int]=1/nodal_volumes[j_int];
 990       A.leftDiagonalScale(nodal_volumes);
 991   }
 992
 993   return A.getEigenvalues( nev, which, tol);
 994 }
 995 std::vector< Vector >
 996 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
 997 {
 998   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 999   return A.getEigenvectors( nev, which, tol);
1000 }
1001 Field
1002 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
1003 {
1004   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
1005   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
1006   Field my_eigenfield;
1007
1008   if(_FECalculation)
1009     my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
1010   else
1011     my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
1012
1013   my_eigenfield.setFieldByDataArrayDouble(d);
1014
1015   return my_eigenfield;
1016 }