CoreFlows/Models/src/LinearElasticityModel.cxx

   1 #include "LinearElasticityModel.hxx"
   2 #include "SparseMatrixPetsc.hxx"
   3 #include "Node.hxx"
   4 #include "math.h"
   5 #include <algorithm>
   6 #include <fstream>
   7 #include <sstream>
   8
   9 using namespace std;
  10
  11 LinearElasticityModel::LinearElasticityModel(int dim, bool FECalculation,  double rho, double lambda, double mu,MPI_Comm comm){
  12         /* Initialisation of PETSC */
  13         //check if PETSC is already initialised
  14         PetscBool petscInitialized;
  15         PetscInitialized(&petscInitialized);
  16         if(!petscInitialized)
  17         {//check if MPI is already initialised
  18                 int mpiInitialized;
  19                 MPI_Initialized(&mpiInitialized);
  20                 if(mpiInitialized)
  21                         PETSC_COMM_WORLD = comm;
  22                 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
  23         }
  24         MPI_Comm_rank(PETSC_COMM_WORLD,&_rank);
  25         MPI_Comm_size(PETSC_COMM_WORLD,&_size);
  26         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.
  27         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.
  28         PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
  29
  30     if(lambda < 0.)
  31     {
  32         std::cout<<"First Lamé coefficient="<<lambda<<endl;
  33         throw CdmathException("Error : First Lamé coefficient lambda cannot  be negative");
  34     }
  35     if(2*mu+dim*lambda < 0.)
  36     {
  37         std::cout<<"First Lamé coefficient="<<lambda<<", second Lamé coefficient="<<mu<<", 2*mu+dim*lambda= "<<2*mu+dim*lambda<<endl;
  38         throw CdmathException("Error : 2*mu+dim*lambda cannot  be negative");
  39     }
  40     if(dim<=0)
  41     {
  42         std::cout<<"space dimension="<<dim<<endl;
  43         throw CdmathException("Error : parameter dim cannot  be negative");
  44     }
  45
  46     _FECalculation=FECalculation;
  47     _onlyNeumannBC=false;
  48
  49         _Ndim=dim;
  50         _nVar=dim;
  51         _initializedMemory=false;
  52
  53     //Mesh data
  54     _neibMaxNbCells=0;
  55     _meshSet=false;
  56     _neibMaxNbNodes=0;
  57
  58     //Boundary conditions
  59     _boundaryNodeIds=std::vector< int >(0);
  60     _dirichletNodeIds=std::vector< int >(0);
  61     _NboundaryNodes=0;
  62     _NdirichletNodes=0;
  63     _NunknownNodes=0;
  64     _dirichletValuesSet=false;
  65     _neumannValuesSet=false;
  66
  67     //Linear solver data
  68         _precision=1.e-6;
  69         _precision_Newton=_precision;
  70         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  71         _maxPetscIts=50;
  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 = "LinearElasticityProblem";
  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         _lambda= lambda;
  93         _mu    = mu;
  94         _rho   = rho;
  95         _densityFieldSet=false;
  96 }
  97
  98 void LinearElasticityModel::initialize()
  99 {
 100         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 101
 102         if(!_meshSet)
 103                 throw CdmathException("LinearElasticityModel::initialize() set mesh first");
 104         else
 105     {
 106                 cout<<"!!!! Initialisation of the computation of the elastic deformation of a solid using ";
 107         *_runLogFile<<"!!!!! Initialisation of the computation of the elastic deformation of a solid using ";
 108         if(!_FECalculation)
 109         {
 110             cout<< "Finite volumes method"<<endl<<endl;
 111             *_runLogFile<< "Finite volumes method"<<endl<<endl;
 112         }
 113         else
 114         {
 115             cout<< "Finite elements method"<<endl<<endl;
 116             *_runLogFile<< "Finite elements method"<<endl<<endl;
 117         }
 118     }
 119         /**************** Field creation *********************/
 120
 121         if(!_densityFieldSet){
 122         if(_FECalculation){
 123             _densityField=Field("Density",NODES,_mesh,1);
 124             for(int i =0; i<_Nnodes; i++)
 125                 _densityField(i) = _rho;
 126         }
 127         else{
 128             _densityField=Field("Density",CELLS,_mesh,1);
 129             for(int i =0; i<_Nmailles; i++)
 130                 _densityField(i) = _rho;
 131         }
 132         _densityFieldSet=true;
 133     }
 134
 135     /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 136     if(_FECalculation)
 137     {
 138         if(_NboundaryNodes==_Nnodes)
 139             cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
 140
 141         for(int i=0; i<_NboundaryNodes; i++)
 142         {
 143             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
 144             if( it != _dirichletBoundaryValues.end() )
 145                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 146             else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
 147             {
 148                 cout<<"!!! No boundary value set for boundary node" << _boundaryNodeIds[i]<< endl;
 149                 *_runLogFile<< "!!! No boundary value set for boundary node" << _boundaryNodeIds[i]<<endl;
 150                 _runLogFile->close();
 151                 throw CdmathException("Missing boundary value");
 152             }
 153             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCLinearElasticity)
 154             {
 155                 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
 156                 *_runLogFile<< "!!!No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
 157                 _runLogFile->close();
 158                 throw CdmathException("Missing boundary condition");
 159             }
 160             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletLinearElasticity)
 161                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 162             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannLinearElasticity)
 163             {
 164                 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
 165                 cout<<"!!! Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity << endl;
 166                 *_runLogFile<< "Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
 167                 *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity <<endl;
 168                 _runLogFile->close();
 169                 throw CdmathException("Wrong boundary condition");
 170             }
 171         }
 172         _NdirichletNodes=_dirichletNodeIds.size();
 173         _NunknownNodes=_Nnodes - _NdirichletNodes;
 174         cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 175                 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 176     }
 177
 178         //creation de la matrice
 179     if(!_FECalculation)
 180         MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles*_nVar, _Nmailles*_nVar, (1+_neibMaxNbCells), PETSC_NULL, &_A);
 181     else
 182         MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes*_nVar, _NunknownNodes*_nVar, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
 183
 184         VecCreate(PETSC_COMM_SELF, &_displacements);
 185
 186         VecDuplicate(_displacements, &_b);//RHS of the linear system
 187
 188         //Linear solver
 189         KSPCreate(PETSC_COMM_SELF, &_ksp);
 190         KSPSetType(_ksp, _ksptype);
 191         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 192         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 193         KSPGetPC(_ksp, &_pc);
 194         PCSetType(_pc, _pctype);
 195
 196     //Checking whether all boundary conditions are Neumann boundary condition
 197     if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
 198     {
 199         map<string, LimitFieldLinearElasticity>::iterator it = _limitField.begin();
 200         while(it != _limitField.end() and (it->second).bcType == NeumannLinearElasticity)
 201             it++;
 202         _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
 203     }
 204     else
 205         if(_FECalculation)
 206             _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
 207         else
 208             _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
 209
 210     //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
 211     if(_onlyNeumannBC)
 212     {
 213         std::cout<<"## Warning all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 214         std::cout<<"## As a consequence we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 215         *_runLogFile<<"## Warning all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 216         *_runLogFile<<"## As a consequence we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 217
 218                 //Check that the matrix is symmetric
 219                 PetscBool isSymetric;
 220                 MatIsSymmetric(_A,_precision,&isSymetric);
 221                 if(!isSymetric)
 222                         {
 223                                 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 224                                 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 225                         }
 226                 MatNullSpace nullsp;
 227                 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
 228                 MatSetNullSpace(_A, nullsp);
 229                 MatSetTransposeNullSpace(_A, nullsp);
 230                 MatNullSpaceDestroy(&nullsp);
 231                 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
 232                 //PCFactorSetShiftAmount(_pc,1e-10);
 233     }
 234
 235         _initializedMemory=true;
 236
 237 }
 238
 239 double LinearElasticityModel::computeStiffnessMatrix(bool & stop)
 240 {
 241     double result;
 242
 243     if(_FECalculation)
 244         result=computeStiffnessMatrixFE(stop);
 245     else
 246         result=computeStiffnessMatrixFV(stop);
 247
 248     if(_verbose or _system)
 249         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 250
 251     return  result;
 252 }
 253
 254 double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){
 255         Cell Cj;
 256         string nameOfGroup;
 257         double dn, coeff;
 258         MatZeroEntries(_A);
 259         VecZeroEntries(_b);
 260
 261     Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 262     std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 263     std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 264     std::vector< Node > nodes(_Ndim+1);//cell nodes
 265     int i_int, j_int; //index of nodes j and k considered as unknown nodes
 266     bool dirichletCell_treated;
 267
 268     std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 269     for (int idim=0; idim<_Ndim+1;idim++)
 270         values[idim][idim]=1;
 271
 272     /* parameters for boundary treatment */
 273     vector< Vector > valuesBorder(_Ndim+1);
 274     Vector GradShapeFuncBorder(_Ndim+1);
 275
 276         for (int j=0; j<_Nmailles;j++)
 277     {
 278                 Cj = _mesh.getCell(j);
 279
 280         for (int idim=0; idim<_Ndim+1;idim++){
 281             nodeIds[idim]=Cj.getNodeId(idim);
 282             nodes[idim]=_mesh.getNode(nodeIds[idim]);
 283             for (int jdim=0; jdim<_Ndim;jdim++)
 284                 M(idim,jdim)=nodes[idim].getPoint()[jdim];
 285             M(idim,_Ndim)=1;
 286         }
 287         for (int idim=0; idim<_Ndim+1;idim++)
 288             GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
 289
 290         /* Loop on the edges of the cell */
 291         for (int idim=0; idim<_Ndim+1;idim++)
 292         {
 293             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 294             {//First node of the edge is not Dirichlet node
 295                 i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 296                 dirichletCell_treated=false;
 297                 for (int jdim=0; jdim<_Ndim+1;jdim++)
 298                 {
 299                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 300                     {//Second node of the edge is not Dirichlet node
 301                         j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 302                         MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 303                     }
 304                     else if (!dirichletCell_treated)
 305                     {//Second node of the edge is a Dirichlet node
 306                         dirichletCell_treated=true;
 307                         for (int kdim=0; kdim<_Ndim+1;kdim++)
 308                         {
 309                                                         std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 310                                                         if( it != _dirichletBoundaryValues.end() )
 311                             {
 312                                 if( _dirichletValuesSet )//New way of storing BC
 313                                     valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 314                                 else    //old way of storing BC
 315                                     valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].displacement;
 316                             }
 317                             else
 318                                 valuesBorder[kdim]=Vector(_Ndim);
 319                         }
 320                         GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
 321                         double coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 322                         VecSetValue(_b,i_int,coeff, ADD_VALUES);
 323                     }
 324                 }
 325             }
 326         }
 327         }
 328
 329     //Calcul de la contribution de la condition limite de Neumann au second membre
 330     if( _NdirichletNodes !=_NboundaryNodes)
 331     {
 332         vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
 333         int NboundaryFaces=boundaryFaces.size();
 334         for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
 335         {
 336             Face Fi = _mesh.getFace(boundaryFaces[i]);
 337             for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
 338             {
 339                 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
 340                 {
 341                     j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 342                     if( _neumannValuesSet )
 343                         coeff =Fi.getMeasure()/(_Ndim+1)*_neumannBoundaryValues[Fi.getNodeId(j)];
 344                     else
 345                         coeff =Fi.getMeasure()/(_Ndim+1)*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalForce;
 346                     VecSetValue(_b, j_int, coeff, ADD_VALUES);
 347                 }
 348             }
 349         }
 350     }
 351
 352     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 353         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 354         VecAssemblyBegin(_b);
 355         VecAssemblyEnd(_b);
 356
 357     stop=false ;
 358
 359         return INFINITY;
 360 }
 361
 362 double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
 363         long nbFaces = _mesh.getNumberOfFaces();
 364         Face Fj;
 365         Cell Cell1,Cell2;
 366         string nameOfGroup;
 367         double inv_dxi, inv_dxj;
 368         double barycenterDistance,coeff;
 369         Vector normale(_Ndim);
 370         double dn;
 371         PetscInt idm, idn;
 372         std::vector< int > idCells;
 373         MatZeroEntries(_A);
 374         VecZeroEntries(_b);
 375
 376         for (int j=0; j<nbFaces;j++){
 377                 Fj = _mesh.getFace(j);
 378
 379                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 380                 idCells = Fj.getCellsId();
 381                 Cell1 = _mesh.getCell(idCells[0]);
 382                 idm = idCells[0];
 383         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 384             if (j == Cell1.getFacesId()[l]){
 385                 for (int idim = 0; idim < _Ndim; ++idim)
 386                     normale[idim] = Cell1.getNormalVector(l,idim);
 387                 break;
 388             }
 389         }
 390
 391                 //Compute velocity at the face Fj
 392                 dn=_conductivity*(_DiffusionTensor*normale)*normale;
 393
 394                 // compute 1/dxi = volume of Ci/area of Fj
 395         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 396
 397                 // If Fj is on the boundary
 398                 if (Fj.getNumberOfCells()==1) {
 399                         if(_verbose )
 400                         {
 401                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 402                                 for(int p=0; p<_Ndim; p++)
 403                                         cout << normale[p] << ",";
 404                                 cout << ") "<<endl;
 405                         }
 406
 407             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
 408             if( it != _dirichletBoundaryValues.end() )
 409             {
 410                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 411                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                                     , ADD_VALUES);
 412                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
 413             }
 414             else
 415             {
 416                 nameOfGroup = Fj.getGroupName();
 417
 418                 if (_limitField[nameOfGroup].bcType==NeumannLinearElasticity){
 419                     if( _neumannValuesSet )
 420                         coeff =Fi.getMeasure()*_neumannBoundaryValues[j];
 421                     else
 422                         coeff =Fi.getMeasure()*_limitField[nameOfGroup].normalForce;
 423                     VecSetValue(_b,idm,    coeff, ADD_VALUES);
 424                 }
 425                 else if(_limitField[nameOfGroup].bcType==DirichletLinearElasticity){
 426                     barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 427                     MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 428                     VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].displacement, ADD_VALUES);
 429                 }
 430                 else {
 431                     stop=true ;
 432                     cout<<"!!!!!!!!!!!!!!! Error LinearElasticityModel::computeStiffnessMatrixFV !!!!!!!!!!"<<endl;
 433                     cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 434                     cout<<"Accepted boundary conditions are Neumann "<<NeumannLinearElasticity<< " and Dirichlet "<<DirichletLinearElasticity<<endl;
 435                     *_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 436                     _runLogFile->close();
 437                     throw CdmathException("Boundary condition not accepted");
 438                 }
 439             }
 440                         // if Fj is inside the domain
 441                 } else  if (Fj.getNumberOfCells()==2 ){
 442                         if(_verbose )
 443                         {
 444                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 445                                 cout << " ; vecteur normal=(";
 446                                 for(int p=0; p<_Ndim; p++)
 447                                         cout << normale[p] << ",";
 448                                 cout << ") "<<endl;
 449                         }
 450                         Cell2 = _mesh.getCell(idCells[1]);
 451                         idn = idCells[1];
 452                         if (_Ndim > 1)
 453                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 454                         else
 455                                 inv_dxj = 1/Cell2.getMeasure();
 456
 457                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 458
 459                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 460                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 461                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 462                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 463                 }
 464                 else
 465         {
 466             *_runLogFile<<"LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face"<<endl;
 467                         throw CdmathException("LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face");
 468         }
 469         }
 470
 471         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 472         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 473         VecAssemblyBegin(_b);
 474         VecAssemblyEnd(_b);
 475
 476     stop=false ;
 477
 478     return INFINITY;
 479 }
 480
 481 double LinearElasticityModel::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeStiffnessMatrix)
 482 {
 483         VecAssemblyBegin(_b);
 484
 485     if(!_FECalculation)
 486         for (int i=0; i<_Nmailles;i++)
 487                         for (int j=0; j<_Ndim;j++)
 488                                 VecSetValue(_b,i*nVar+j,_gravity(j)*_densityField(i),ADD_VALUES);
 489     else
 490         {
 491             Cell Ci;
 492             std::vector< int > nodesId;
 493             for (int i=0; i<_Nmailles;i++)
 494             {
 495                 Ci=_mesh.getCell(i);
 496                 nodesId=Ci.getNodesId();
 497                 for (int j=0; j<nodesId.size();j++)
 498                     if(!_mesh.isBorderNode(nodesId[j]))
 499                                                 for (int k=0; k<_Ndim; k++)
 500                                                         VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds)*nVar+k, _gravity(k)*_densityField(j)*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 501             }
 502         }
 503
 504         VecAssemblyEnd(_b);
 505
 506     if(_verbose or _system)
 507         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 508
 509     stop=false ;
 510         return INFINITY;
 511 }
 512
 513 bool LinearElasticityModel::solveLinearSystem()
 514 {
 515         bool resolutionOK;
 516
 517     //Only implicit scheme considered
 518     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 519     MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 520
 521 #if PETSC_VERSION_GREATER_3_5
 522     KSPSetOperators(_ksp, _A, _A);
 523 #else
 524     KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 525 #endif
 526
 527     if(_conditionNumber)
 528         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 529     KSPSolve(_ksp, _b, _displacements);
 530
 531         KSPConvergedReason reason;
 532         KSPGetConvergedReason(_ksp,&reason);
 533     KSPGetIterationNumber(_ksp, &_PetscIts);
 534     double residu;
 535     KSPGetResidualNorm(_ksp,&residu);
 536         if (reason!=2 and reason!=3)
 537     {
 538         cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 539         cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 540         cout<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 541                 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 542         *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 543         *_runLogFile<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 544                 _runLogFile->close();
 545
 546                 resolutionOK = false;
 547     }
 548     else{
 549         if( _MaxIterLinearSolver < _PetscIts)
 550             _MaxIterLinearSolver = _PetscIts;
 551         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;
 552                 *_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;
 553
 554         resolutionOK = true;
 555     }
 556
 557         return resolutionOK;
 558 }
 559
 560 void LinearElasticityModel::setMesh(const Mesh &M)
 561 {
 562         if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
 563         {
 564         cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
 565                 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
 566                 *_runLogFile<<"LinearElasticityModel::setMesh: mesh has incorrect dimension"<<endl;
 567                 _runLogFile->close();
 568                 throw CdmathException("LinearElasticityModel::setMesh: mesh has incorrect  dimension");
 569         }
 570
 571         _mesh=M;
 572         _Nmailles = _mesh.getNumberOfCells();
 573         _Nnodes =   _mesh.getNumberOfNodes();
 574
 575     cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
 576         *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
 577
 578         // find  maximum nb of neibourghs
 579     if(!_FECalculation)
 580     {
 581         _VV=Field ("Displacements", CELLS, _mesh, _Ndim);
 582         _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 583     }
 584     else
 585     {
 586         if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 587                         cout<<"1D Finite element method on segments"<<endl;
 588         else if(_Ndim==2)
 589         {
 590                         if( _mesh.isTriangular() )//Mesh dim=2
 591                                 cout<<"2D Finite element method on triangles"<<endl;
 592                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 593                                 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 594                         else
 595                         {
 596                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 597                                 *_runLogFile<<"LinearElasticityModel::setMesh: mesh has incorrect dimension"<<endl;
 598                                 _runLogFile->close();
 599                                 throw CdmathException("LinearElasticityModel::setMesh: mesh has incorrect cell types");
 600                         }
 601         }
 602         else if(_Ndim==3)
 603         {
 604                         if( _mesh.isTetrahedral() )//Mesh dim=3
 605                                 cout<<"3D Finite element method on tetrahedra"<<endl;
 606                         else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 607                                 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 608                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 609                                 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 610                         else
 611                         {
 612                                 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;
 613                                 *_runLogFile<<"LinearElasticityModel::setMesh: mesh has incorrect dimension"<<endl;
 614                                 _runLogFile->close();
 615                                 throw CdmathException("LinearElasticityModel::setMesh: mesh has incorrect cell types");
 616                         }
 617         }
 618
 619                 _VV=Field ("Temperature", NODES, _mesh, _Ndim);
 620
 621         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 622         _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 623         _NboundaryNodes=_boundaryNodeIds.size();
 624     }
 625
 626         _meshSet=true;
 627 }
 628
 629 void LinearElasticityModel::setLinearSolver(linearSolver kspType, preconditioner pcType)
 630 {
 631         //_maxPetscIts=maxIterationsPetsc;
 632         // set linear solver algorithm
 633         if (kspType==GMRES)
 634                 _ksptype = (char*)&KSPGMRES;
 635         else if (kspType==CG)
 636                 _ksptype = (char*)&KSPCG;
 637         else if (kspType==BCGS)
 638                 _ksptype = (char*)&KSPBCGS;
 639         else {
 640                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 641                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 642                 _runLogFile->close();
 643                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 644         }
 645         // set preconditioner
 646         if (pcType == NONE)
 647                 _pctype = (char*)&PCNONE;
 648         else if (pcType ==LU)
 649                 _pctype = (char*)&PCLU;
 650         else if (pcType == ILU)
 651                 _pctype = (char*)&PCILU;
 652         else if (pcType ==CHOLESKY)
 653                 _pctype = (char*)&PCCHOLESKY;
 654         else if (pcType == ICC)
 655                 _pctype = (char*)&PCICC;
 656         else {
 657                 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 658                 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 659                 _runLogFile->close();
 660                 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 661         }
 662 }
 663
 664 bool LinearElasticityModel::solveStationaryProblem()
 665 {
 666         if(!_initializedMemory)
 667         {
 668                 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
 669                 _runLogFile->close();
 670                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 671         }
 672         bool stop=false; // Does the Problem want to stop (error) ?
 673
 674         cout<< "!!! Running test case "<< _fileName << " using ";
 675         *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
 676
 677     if(!_FECalculation)
 678     {
 679         cout<< "Finite volumes method"<<endl<<endl;
 680                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 681         }
 682     else
 683         {
 684         cout<< "Finite elements method"<<endl<<endl;
 685                 *_runLogFile<< "Finite elements method"<< endl<<endl;
 686         }
 687
 688     computeStiffnessMatrix( stop);
 689     if (stop){
 690         cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 691         *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 692                 _runLogFile->close();
 693        throw CdmathException("Failed computing diffusion matrix");
 694     }
 695     computeRHS(stop);
 696     if (stop){
 697         cout << "Error : failed computing right hand side, stopping calculation"<< endl;
 698         *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
 699         throw CdmathException("Failed computing right hand side");
 700     }
 701     stop = !solveLinearSystem();
 702     if (stop){
 703         cout << "Error : failed solving linear system, stopping calculation"<< endl;
 704         *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
 705                 _runLogFile->close();
 706         throw CdmathException("Failed solving linear system");
 707     }
 708
 709     _computationCompletedSuccessfully=true;
 710     save();
 711
 712     *_runLogFile<< "!!!!!! Computation successful"<< endl;
 713         _runLogFile->close();
 714
 715         return !stop;
 716 }
 717
 718 void LinearElasticityModel::save(){
 719     cout<< "Saving numerical results"<<endl<<endl;
 720     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 721
 722         string resultFile(_path+"/LinearElasticityModel");//Results
 723
 724         resultFile+="_";
 725         resultFile+=_fileName;
 726
 727         // create mesh and component info
 728     string suppress ="rm -rf "+resultFile+"_*";
 729     system(suppress.c_str());//Nettoyage des précédents calculs identiques
 730
 731     if(_verbose or _system)
 732         VecView(_displacements,PETSC_VIEWER_STDOUT_SELF);
 733
 734     double uk;
 735     if(!_FECalculation)
 736         for(int i=0; i<_Nmailles; i++)
 737             {
 738                                 for(int j=0; j<_nVar; j++)
 739                                 {
 740                                         int k=i*_nVar+j;
 741                                         VecGetValues(_displacements, 1, &k, &uk);
 742                                         _VV(i,j)=uk;
 743                                 }
 744             }
 745     else
 746     {
 747         int globalIndex;
 748         for(int i=0; i<_NunknownNodes; i++)
 749         {
 750                         globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
 751                         for(int j=0; j<_nVar; j++)
 752                         {
 753                                 int k=i*_nVar+j;
 754                                 VecGetValues(_displacements, 1, &k, &uk);
 755                                 _VV(globalIndex,j)=uk;
 756                         }
 757         }
 758
 759         Node Ni;
 760         string nameOfGroup;
 761         for(int i=0; i<_NdirichletNodes; i++)
 762         {
 763             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 764             nameOfGroup = Ni.getGroupName();
 765                         for(int j=0; j<_nVar; j++)
 766                                 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].displacement[i];
 767         }
 768     }
 769
 770     switch(_saveFormat)
 771     {
 772         case VTK :
 773             _VV.writeVTK(resultFile);
 774             break;
 775         case MED :
 776             _VV.writeMED(resultFile);
 777             break;
 778         case CSV :
 779             _VV.writeCSV(resultFile);
 780             break;
 781     }
 782 }
 783 Field
 784 LinearElasticityModel::getOutputDisplacementField()
 785 {
 786     if(!_computationCompletedSuccessfully)
 787         throw("Computation not performed yet or failed. No displacement field available");
 788     else
 789         return _VV;
 790 }
 791
 792 void LinearElasticityModel::terminate()
 793 {
 794         VecDestroy(&_displacements);
 795         VecDestroy(&_b);
 796         MatDestroy(&_A);
 797 }
 798 void
 799 LinearElasticityModel::setDirichletValues(map< int, double> dirichletBoundaryValues)
 800 {
 801     _dirichletValuesSet=true;
 802     _dirichletBoundaryValues=dirichletBoundaryValues;
 803 }
 804 void
 805 LinearElasticityModel::setNeumannValues(map< int, double> neumannBoundaryValues)
 806 {
 807     _neumannValuesSet=true;
 808     _neumannBoundaryValues=neumannBoundaryValues;
 809 }
 810
 811 double
 812 LinearElasticityModel::getConditionNumber(bool isSingular, double tol) const
 813 {
 814   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 815   return A.getConditionNumber( isSingular, tol);
 816 }