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