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