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                 //Check that the matrix is symmetric
 277                 PetscBool isSymetric;
 278                 MatIsSymmetric(_A,_precision,&isSymetric);
 279                 if(!isSymetric)
 280                         {
 281                                 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 282                                 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 283                         }
 284                 MatNullSpace nullsp;
 285                 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
 286                 MatSetNullSpace(_A, nullsp);
 287                 MatSetTransposeNullSpace(_A, nullsp);
 288                 MatNullSpaceDestroy(&nullsp);
 289                 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
 290                 //PCFactorSetShiftAmount(_pc,1e-10);
 291     }
 292
 293         _initializedMemory=true;
 294
 295 }
 296
 297 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
 298         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 299                 computeDiffusionMatrix(stop);
 300
 301         _dt_src=computeRHS(stop);
 302         stop=false;
 303         return _dt_src;
 304 }
 305
 306 Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
 307     vector< Matrix > matrices(_Ndim);
 308
 309     for (int idim=0; idim<_Ndim;idim++){
 310         matrices[idim]=M.deepCopy();
 311         for (int jdim=0; jdim<_Ndim+1;jdim++)
 312                         matrices[idim](jdim,idim) = values[jdim] ;
 313     }
 314
 315         Vector result(_Ndim);
 316     for (int idim=0; idim<_Ndim;idim++)
 317         result[idim] = matrices[idim].determinant();
 318
 319         return result;
 320 }
 321
 322 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
 323 {
 324     double result;
 325
 326     if(_FECalculation)
 327         result=computeDiffusionMatrixFE(stop);
 328     else
 329         result=computeDiffusionMatrixFV(stop);
 330
 331     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 332     //update value here if variable  heat transfer coefficient
 333     if(_heatTransfertCoeff>_precision)
 334         MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
 335
 336     if(_verbose or _system)
 337         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 338
 339     return  result;
 340 }
 341
 342 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 343         Cell Cj;
 344         string nameOfGroup;
 345         double dn, coeff;
 346         MatZeroEntries(_A);
 347         VecZeroEntries(_b);
 348
 349     Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 350     std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 351     std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 352     std::vector< Node > nodes(_Ndim+1);//cell nodes
 353     int i_int, j_int; //index of nodes j and k considered as unknown nodes
 354     bool dirichletCell_treated;
 355
 356     std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 357     for (int idim=0; idim<_Ndim+1;idim++)
 358         values[idim][idim]=1;
 359
 360     /* parameters for boundary treatment */
 361     vector< double > valuesBorder(_Ndim+1);
 362     Vector GradShapeFuncBorder(_Ndim+1);
 363
 364         for (int j=0; j<_Nmailles;j++)
 365     {
 366                 Cj = _mesh.getCell(j);
 367
 368         for (int idim=0; idim<_Ndim+1;idim++){
 369             nodeIds[idim]=Cj.getNodeId(idim);
 370             nodes[idim]=_mesh.getNode(nodeIds[idim]);
 371             for (int jdim=0; jdim<_Ndim;jdim++)
 372                 M(idim,jdim)=nodes[idim].getPoint()[jdim];
 373             M(idim,_Ndim)=1;
 374         }
 375         for (int idim=0; idim<_Ndim+1;idim++)
 376             GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
 377
 378         /* Loop on the edges of the cell */
 379         for (int idim=0; idim<_Ndim+1;idim++)
 380         {
 381             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 382             {//First node of the edge is not Dirichlet node
 383                 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 384                 dirichletCell_treated=false;
 385                 for (int jdim=0; jdim<_Ndim+1;jdim++)
 386                 {
 387                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 388                     {//Second node of the edge is not Dirichlet node
 389                         j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 390                         MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 391                     }
 392                     else if (!dirichletCell_treated)
 393                     {//Second node of the edge is a Dirichlet node
 394                         dirichletCell_treated=true;
 395                         for (int kdim=0; kdim<_Ndim+1;kdim++)
 396                         {
 397                                                         std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 398                                                         if( it != _dirichletBoundaryValues.end() )
 399                             {
 400                                 if( _dirichletValuesSet )
 401                                     valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 402                                 else
 403                                     valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
 404                             }
 405                             else
 406                                 valuesBorder[kdim]=0;
 407                         }
 408                         GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
 409                         coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 410                         VecSetValue(_b,i_int,coeff, ADD_VALUES);
 411                     }
 412                 }
 413             }
 414         }
 415         }
 416
 417     //Calcul de la contribution de la condition limite de Neumann au second membre
 418     if( _NdirichletNodes !=_NboundaryNodes)
 419     {
 420         vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
 421         int NboundaryFaces=boundaryFaces.size();
 422         for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
 423         {
 424             Face Fi = _mesh.getFace(i);
 425             for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
 426             {
 427                 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
 428                 {
 429                     j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 430                     if( _neumannValuesSet )
 431                         coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
 432                     else
 433                         coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
 434                     VecSetValue(_b, j_int, coeff, ADD_VALUES);
 435                 }
 436             }
 437         }
 438     }
 439     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 440         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 441         VecAssemblyBegin(_b);
 442         VecAssemblyEnd(_b);
 443
 444         _diffusionMatrixSet=true;
 445     stop=false ;
 446
 447         return INFINITY;
 448 }
 449
 450 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 451         long nbFaces = _mesh.getNumberOfFaces();
 452         Face Fj;
 453         Cell Cell1,Cell2;
 454         string nameOfGroup;
 455         double inv_dxi, inv_dxj;
 456         double barycenterDistance;
 457         Vector normale(_Ndim);
 458         double dn;
 459         PetscInt idm, idn;
 460         std::vector< int > idCells;
 461         MatZeroEntries(_A);
 462         VecZeroEntries(_b);
 463
 464         for (int j=0; j<nbFaces;j++){
 465                 Fj = _mesh.getFace(j);
 466
 467                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 468                 idCells = Fj.getCellsId();
 469                 Cell1 = _mesh.getCell(idCells[0]);
 470                 idm = idCells[0];
 471         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 472             if (j == Cell1.getFacesId()[l]){
 473                 for (int idim = 0; idim < _Ndim; ++idim)
 474                     normale[idim] = Cell1.getNormalVector(l,idim);
 475                 break;
 476             }
 477         }
 478
 479                 //Compute velocity at the face Fj
 480                 dn=_conductivity*(_DiffusionTensor*normale)*normale;
 481
 482                 // compute 1/dxi = volume of Ci/area of Fj
 483         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 484
 485                 // If Fj is on the boundary
 486                 if (Fj.getNumberOfCells()==1) {
 487                         if(_verbose )
 488                         {
 489                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 490                                 for(int p=0; p<_Ndim; p++)
 491                                         cout << normale[p] << ",";
 492                                 cout << ") "<<endl;
 493                         }
 494
 495             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
 496             if( it != _dirichletBoundaryValues.end() )
 497             {
 498                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 499                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                                     , ADD_VALUES);
 500                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
 501             }
 502             else
 503             {
 504                 nameOfGroup = Fj.getGroupName();
 505
 506                 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
 507                     VecSetValue(_b,idm,    -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 508                 }
 509                 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
 510                     barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 511                     MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 512                     VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 513                 }
 514                 else {
 515                     stop=true ;
 516                     cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 517                     cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 518                     cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
 519                     *_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 520                     _runLogFile->close();
 521                     throw CdmathException("Boundary condition not accepted");
 522                 }
 523             }
 524                         // if Fj is inside the domain
 525                 } else  if (Fj.getNumberOfCells()==2 ){
 526                         if(_verbose )
 527                         {
 528                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 529                                 cout << " ; vecteur normal=(";
 530                                 for(int p=0; p<_Ndim; p++)
 531                                         cout << normale[p] << ",";
 532                                 cout << ") "<<endl;
 533                         }
 534                         Cell2 = _mesh.getCell(idCells[1]);
 535                         idn = idCells[1];
 536                         if (_Ndim > 1)
 537                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 538                         else
 539                                 inv_dxj = 1/Cell2.getMeasure();
 540
 541                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 542
 543                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 544                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 545                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 546                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 547                 }
 548                 else
 549         {
 550             *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
 551                         throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 552         }
 553         }
 554
 555         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 556         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 557         VecAssemblyBegin(_b);
 558         VecAssemblyEnd(_b);
 559
 560         _diffusionMatrixSet=true;
 561     stop=false ;
 562
 563     return INFINITY;
 564 }
 565
 566 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)
 567 {
 568         VecAssemblyBegin(_b);
 569
 570     if(!_FECalculation)
 571         for (int i=0; i<_Nmailles;i++)
 572             VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
 573     else
 574         {
 575             Cell Ci;
 576             std::vector< int > nodesId;
 577             for (int i=0; i<_Nmailles;i++)
 578             {
 579                 Ci=_mesh.getCell(i);
 580                 nodesId=Ci.getNodesId();
 581                 for (int j=0; j<nodesId.size();j++)
 582                     if(!_mesh.isBorderNode(nodesId[j]))
 583                     {
 584                         double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 585                         VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 586                     }
 587             }
 588         }
 589
 590         VecAssemblyEnd(_b);
 591
 592     if(_verbose or _system)
 593         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 594
 595     stop=false ;
 596         return INFINITY;
 597 }
 598
 599 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
 600 {
 601         bool stop;
 602
 603     //Only implicit scheme considered
 604     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 605     MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 606
 607 #if PETSC_VERSION_GREATER_3_5
 608     KSPSetOperators(_ksp, _A, _A);
 609 #else
 610     KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 611 #endif
 612
 613     if(_conditionNumber)
 614         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 615     KSPSolve(_ksp, _b, _Tk);
 616
 617         KSPConvergedReason reason;
 618         KSPGetConvergedReason(_ksp,&reason);
 619     KSPGetIterationNumber(_ksp, &_PetscIts);
 620     double residu;
 621     KSPGetResidualNorm(_ksp,&residu);
 622         if (reason!=2 and reason!=3)
 623     {
 624         cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 625         cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 626         cout<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 627                 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 628         *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 629         *_runLogFile<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 630                 _runLogFile->close();
 631         converged = false;
 632         stop = true;
 633     }
 634     else{
 635         if( _MaxIterLinearSolver < _PetscIts)
 636             _MaxIterLinearSolver = _PetscIts;
 637         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;
 638                 *_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;
 639         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 640         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 641
 642         _erreur_rel= 0;
 643         double Ti, dTi;
 644
 645         VecAssemblyBegin(_Tk);
 646         VecAssemblyEnd(  _Tk);
 647         VecAssemblyBegin(_deltaT);
 648         VecAssemblyEnd(  _deltaT);
 649
 650         if(!_FECalculation)
 651             for(int i=0; i<_Nmailles; i++)
 652             {
 653                 VecGetValues(_deltaT, 1, &i, &dTi);
 654                 VecGetValues(_Tk, 1, &i, &Ti);
 655                 if(_erreur_rel < fabs(dTi/Ti))
 656                     _erreur_rel = fabs(dTi/Ti);
 657             }
 658         else
 659             for(int i=0; i<_NunknownNodes; i++)
 660             {
 661                 VecGetValues(_deltaT, 1, &i, &dTi);
 662                 VecGetValues(_Tk, 1, &i, &Ti);
 663                 if(_erreur_rel < fabs(dTi/Ti))
 664                     _erreur_rel = fabs(dTi/Ti);
 665             }
 666         stop=false;
 667         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 668     }
 669
 670         VecCopy(_Tk, _Tkm1);
 671
 672         return stop;
 673 }
 674
 675 void StationaryDiffusionEquation::setMesh(const Mesh &M)
 676 {
 677         if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
 678         {
 679         cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
 680                 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
 681                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 682                 _runLogFile->close();
 683                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect  dimension");
 684         }
 685
 686         _mesh=M;
 687         _Nmailles = _mesh.getNumberOfCells();
 688         _Nnodes =   _mesh.getNumberOfNodes();
 689
 690     cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
 691         *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
 692
 693         // find  maximum nb of neibourghs
 694     if(!_FECalculation)
 695     {
 696         _VV=Field ("Temperature", CELLS, _mesh, 1);
 697         _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 698     }
 699     else
 700     {
 701         if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 702                         cout<<"1D Finite element method on segments"<<endl;
 703         else if(_Ndim==2)
 704         {
 705                         if( _mesh.isTriangular() )//Mesh dim=2
 706                                 cout<<"2D Finite element method on triangles"<<endl;
 707                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 708                                 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 709                         else
 710                         {
 711                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 712                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 713                                 _runLogFile->close();
 714                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 715                         }
 716         }
 717         else if(_Ndim==3)
 718         {
 719                         if( _mesh.isTetrahedral() )//Mesh dim=3
 720                                 cout<<"3D Finite element method on tetrahedra"<<endl;
 721                         else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 722                                 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 723                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 724                                 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 725                         else
 726                         {
 727                                 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;
 728                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 729                                 _runLogFile->close();
 730                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 731                         }
 732         }
 733
 734                 _VV=Field ("Temperature", NODES, _mesh, 1);
 735
 736         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 737         _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 738         _NboundaryNodes=_boundaryNodeIds.size();
 739     }
 740
 741         _meshSet=true;
 742 }
 743
 744 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
 745 {
 746         //_maxPetscIts=maxIterationsPetsc;
 747         // set linear solver algorithm
 748         if (kspType==GMRES)
 749                 _ksptype = (char*)&KSPGMRES;
 750         else if (kspType==CG)
 751                 _ksptype = (char*)&KSPCG;
 752         else if (kspType==BCGS)
 753                 _ksptype = (char*)&KSPBCGS;
 754         else {
 755                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 756                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 757                 _runLogFile->close();
 758                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 759         }
 760         // set preconditioner
 761         if (pcType == NONE)
 762                 _pctype = (char*)&PCNONE;
 763         else if (pcType ==LU)
 764                 _pctype = (char*)&PCLU;
 765         else if (pcType == ILU)
 766                 _pctype = (char*)&PCILU;
 767         else if (pcType ==CHOLESKY)
 768                 _pctype = (char*)&PCCHOLESKY;
 769         else if (pcType == ICC)
 770                 _pctype = (char*)&PCICC;
 771         else {
 772                 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 773                 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 774                 _runLogFile->close();
 775                 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 776         }
 777 }
 778
 779 bool StationaryDiffusionEquation::solveStationaryProblem()
 780 {
 781         if(!_initializedMemory)
 782         {
 783                 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
 784                 _runLogFile->close();
 785                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 786         }
 787         bool stop=false; // Does the Problem want to stop (error) ?
 788         bool converged=false; // has the newton scheme converged (end) ?
 789
 790         cout<< "!!! Running test case "<< _fileName << " using ";
 791         *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
 792
 793     if(!_FECalculation)
 794     {
 795         cout<< "Finite volumes method"<<endl<<endl;
 796                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 797         }
 798     else
 799         {
 800         cout<< "Finite elements method"<<endl<<endl;
 801                 *_runLogFile<< "Finite elements method"<< endl<<endl;
 802         }
 803
 804     computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
 805     if (stop){
 806         cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 807         *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 808                 _runLogFile->close();
 809        throw CdmathException("Failed computing diffusion matrix");
 810     }
 811     computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
 812     if (stop){
 813         cout << "Error : failed computing right hand side, stopping calculation"<< endl;
 814         *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
 815         throw CdmathException("Failed computing right hand side");
 816     }
 817     stop = iterateNewtonStep(converged);
 818     if (stop){
 819         cout << "Error : failed solving linear system, stopping calculation"<< endl;
 820         *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
 821                 _runLogFile->close();
 822         throw CdmathException("Failed solving linear system");
 823     }
 824
 825     _computationCompletedSuccessfully=true;
 826     save();
 827
 828         // Newton iteration loop for non linear problems
 829     /*
 830         while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
 831         {
 832         computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
 833         computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
 834         stop = iterateNewtonStep(converged);
 835         _NEWTON_its++;
 836         }
 837     if (stop){
 838         cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 839         *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 840         throw CdmathException("Failed solving a Newton iteration");
 841     }
 842     else if(_NEWTON_its==_maxNewtonIts){
 843         cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 844         *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 845         throw CdmathException("No convergence of Newton scheme");
 846     }
 847     else{
 848         cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 849         *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 850         save();
 851     }
 852     */
 853
 854     *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
 855         _runLogFile->close();
 856
 857         return !stop;
 858 }
 859
 860 void StationaryDiffusionEquation::save(){
 861     cout<< "Saving numerical results"<<endl<<endl;
 862     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 863
 864         string resultFile(_path+"/StationaryDiffusionEquation");//Results
 865
 866         resultFile+="_";
 867         resultFile+=_fileName;
 868
 869         // create mesh and component info
 870     string suppress ="rm -rf "+resultFile+"_*";
 871     system(suppress.c_str());//Nettoyage des précédents calculs identiques
 872
 873     if(_verbose or _system)
 874         VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
 875
 876     double Ti;
 877     if(!_FECalculation)
 878         for(int i=0; i<_Nmailles; i++)
 879             {
 880                 VecGetValues(_Tk, 1, &i, &Ti);
 881                 _VV(i)=Ti;
 882             }
 883     else
 884     {
 885         int globalIndex;
 886         for(int i=0; i<_NunknownNodes; i++)
 887         {
 888             VecGetValues(_Tk, 1, &i, &Ti);
 889             globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 890             _VV(globalIndex)=Ti;
 891         }
 892
 893         Node Ni;
 894         string nameOfGroup;
 895         for(int i=0; i<_NdirichletNodes; i++)
 896         {
 897             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 898             nameOfGroup = Ni.getGroupName();
 899             _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 900         }
 901     }
 902
 903     _VV.setInfoOnComponent(0,"Temperature_(K)");
 904     switch(_saveFormat)
 905     {
 906         case VTK :
 907             _VV.writeVTK(resultFile);
 908             break;
 909         case MED :
 910             _VV.writeMED(resultFile);
 911             break;
 912         case CSV :
 913             _VV.writeCSV(resultFile);
 914             break;
 915     }
 916 }
 917 Field
 918 StationaryDiffusionEquation::getOutputTemperatureField()
 919 {
 920     if(!_computationCompletedSuccessfully)
 921         throw("Computation not performed yet or failed. No temperature field available");
 922     else
 923         return _VV;
 924 }
 925
 926 void StationaryDiffusionEquation::terminate()
 927 {
 928         VecDestroy(&_Tk);
 929         VecDestroy(&_Tkm1);
 930         VecDestroy(&_deltaT);
 931         VecDestroy(&_b);
 932         MatDestroy(&_A);
 933 }
 934 void
 935 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 936 {
 937     _dirichletValuesSet=true;
 938     _dirichletBoundaryValues=dirichletBoundaryValues;
 939 }
 940
 941 void
 942 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
 943 {
 944     _neumannValuesSet=true;
 945     _neumannBoundaryValues=neumannBoundaryValues;
 946 }
 947
 948 double
 949 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
 950 {
 951   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 952   return A.getConditionNumber( isSingular, tol);
 953 }
 954 std::vector< double >
 955 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
 956 {
 957   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 958
 959   if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
 960   {
 961       Vector nodal_volumes(_NunknownNodes);
 962       int j_int;
 963       for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
 964       {
 965         Cell Ci = _mesh.getCell(i);
 966         for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
 967         {
 968             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
 969                         {
 970                 j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 971                 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
 972             }
 973         }
 974        }
 975       for( j_int = 0; j_int< _NunknownNodes ; j_int++)
 976         nodal_volumes[j_int]=1/nodal_volumes[j_int];
 977       A.leftDiagonalScale(nodal_volumes);
 978   }
 979
 980   return A.getEigenvalues( nev, which, tol);
 981 }
 982 std::vector< Vector >
 983 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
 984 {
 985   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 986   return A.getEigenvectors( nev, which, tol);
 987 }
 988 Field
 989 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
 990 {
 991   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 992   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
 993   Field my_eigenfield;
 994
 995   if(_FECalculation)
 996     my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
 997   else
 998     my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
 999
1000   my_eigenfield.setFieldByDataArrayDouble(d);
1001
1002   return my_eigenfield;
1003 }