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