CoreFlows/Models/src/StationaryDiffusionEquation.cxx

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