1 #include "LinearElasticityModel.hxx"
2 #include "SparseMatrixPetsc.hxx"
11 int LinearElasticityModel::fact(int n)
13 return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
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)
22 return globalIndex-boundarySize;
23 else if (dirichletNodes[j]>globalIndex)
26 throw CdmathException("LinearElasticityModel::unknownNodeIndex : Error : node is a Dirichlet boundary node");
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 */
34 return unknownNodeIndex;
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)
41 unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
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;
51 LinearElasticityModel::LinearElasticityModel(int dim, bool FECalculation, double rho, double lambda, double mu){
52 PetscBool petscInitialized;
53 PetscInitialized(&petscInitialized);
55 PetscInitialize(NULL,NULL,0,0);
59 std::cout<<"First Lamé coefficient="<<lambda<<endl;
60 throw CdmathException("Error : First Lamé coefficient lambda cannot be negative");
62 if(2*mu+dim*lambda < 0.)
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");
69 std::cout<<"space dimension="<<dim<<endl;
70 throw CdmathException("Error : parameter dim cannot be negative");
73 _FECalculation=FECalculation;
78 _initializedMemory=false;
86 _boundaryNodeIds=std::vector< int >(0);
87 _dirichletNodeIds=std::vector< int >(0);
91 _dirichletValuesSet=false;
92 _neumannValuesSet=false;
96 _precision_Newton=_precision;
97 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
99 int _PetscIts=0;//the number of iterations of the linear solver
100 _ksptype = (char*)&KSPGMRES;
101 _pctype = (char*)&PCLU;
102 _conditionNumber=false;
105 //parameters for monitoring simulation
108 _runLogFile=new ofstream;
110 //result save parameters
111 _fileName = "LinearElasticityProblem";
112 char result[ PATH_MAX ];//extracting current directory
113 getcwd(result, PATH_MAX );
114 _path=string( result );
116 _computationCompletedSuccessfully=false;
118 //heat transfer parameters
122 _densityFieldSet=false;
125 void LinearElasticityModel::initialize()
127 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
130 throw CdmathException("LinearElasticityModel::initialize() set mesh first");
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 ";
137 cout<< "Finite volumes method"<<endl<<endl;
138 *_runLogFile<< "Finite volumes method"<<endl<<endl;
142 cout<< "Finite elements method"<<endl<<endl;
143 *_runLogFile<< "Finite elements method"<<endl<<endl;
146 /**************** Field creation *********************/
148 if(!_densityFieldSet){
150 _densityField=Field("Density",NODES,_mesh,1);
151 for(int i =0; i<_Nnodes; i++)
152 _densityField(i) = _rho;
155 _densityField=Field("Density",CELLS,_mesh,1);
156 for(int i =0; i<_Nmailles; i++)
157 _densityField(i) = _rho;
159 _densityFieldSet=true;
162 /* Détection des noeuds frontière avec une condition limite de Dirichlet */
165 if(_NboundaryNodes==_Nnodes)
166 cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
168 for(int i=0; i<_NboundaryNodes; i++)
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 )
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");
180 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCLinearElasticity)
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");
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)
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");
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;
205 //creation de la matrice
207 MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles*_nVar, _Nmailles*_nVar, (1+_neibMaxNbCells), PETSC_NULL, &_A);
209 MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes*_nVar, _NunknownNodes*_nVar, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
211 VecCreate(PETSC_COMM_SELF, &_displacements);
213 VecDuplicate(_displacements, &_b);//RHS of the linear system
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);
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
227 map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
228 while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
230 _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
234 _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
236 _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
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)
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
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;
251 //Check that the matrix is symmetric
252 PetscBool isSymetric;
253 MatIsSymmetric(_A,_precision,&isSymetric);
256 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
257 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
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);
268 _initializedMemory=true;
272 Vector LinearElasticityModel::gradientNodal(Matrix M, vector< double > values){
273 vector< Matrix > matrices(_Ndim);
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] ;
281 Vector result(_Ndim);
282 for (int idim=0; idim<_Ndim;idim++)
283 result[idim] = matrices[idim].determinant();
288 double LinearElasticityModel::computeStiffnessMatrix(bool & stop)
293 result=computeStiffnessMatrixFE(stop);
295 result=computeStiffnessMatrixFV(stop);
297 if(_verbose or _system)
298 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
303 double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){
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;
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;
321 /* parameters for boundary treatment */
322 vector< Vector > valuesBorder(_Ndim+1);
323 Vector GradShapeFuncBorder(_Ndim+1);
325 for (int j=0; j<_Nmailles;j++)
327 Cj = _mesh.getCell(j);
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];
336 for (int idim=0; idim<_Ndim+1;idim++)
337 GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
339 /* Loop on the edges of the cell */
340 for (int idim=0; idim<_Ndim+1;idim++)
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++)
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);
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++)
358 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
359 if( it != _dirichletBoundaryValues.end() )
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;
367 valuesBorder[kdim]=Vector(_Ndim);
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);
378 //Calcul de la contribution de la condition limite de Neumann au second membre
379 if( _NdirichletNodes !=_NboundaryNodes)
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
385 Face Fi = _mesh.getFace(i);
386 for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
388 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
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)];
394 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalForce;
395 VecSetValue(_b, j_int, coeff, ADD_VALUES);
401 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
402 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
403 VecAssemblyBegin(_b);
411 double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
412 long nbFaces = _mesh.getNumberOfFaces();
416 double inv_dxi, inv_dxj;
417 double barycenterDistance;
418 Vector normale(_Ndim);
421 std::vector< int > idCells;
425 for (int j=0; j<nbFaces;j++){
426 Fj = _mesh.getFace(j);
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]);
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);
440 //Compute velocity at the face Fj
441 dn=_conductivity*(_DiffusionTensor*normale)*normale;
443 // compute 1/dxi = volume of Ci/area of Fj
444 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
446 // If Fj is on the boundary
447 if (Fj.getNumberOfCells()==1) {
450 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
451 for(int p=0; p<_Ndim; p++)
452 cout << normale[p] << ",";
456 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
457 if( it != _dirichletBoundaryValues.end() )
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);
465 nameOfGroup = Fj.getGroupName();
467 if (_limitField[nameOfGroup].bcType==NeumannLinearElasticity){//Nothing to do
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);
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");
484 // if Fj is inside the domain
485 } else if (Fj.getNumberOfCells()==2 ){
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] << ",";
494 Cell2 = _mesh.getCell(idCells[1]);
497 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
499 inv_dxj = 1/Cell2.getMeasure();
501 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
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);
510 *_runLogFile<<"LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face"<<endl;
511 throw CdmathException("LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face");
515 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
516 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
517 VecAssemblyBegin(_b);
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)
527 VecAssemblyBegin(_b);
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);
536 std::vector< int > nodesId;
537 for (int i=0; i<_Nmailles;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);
550 if(_verbose or _system)
551 VecView(_b,PETSC_VIEWER_STDOUT_SELF);
557 bool LinearElasticityModel::solveLinearSystem()
561 //Only implicit scheme considered
562 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
563 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
565 #if PETSC_VERSION_GREATER_3_5
566 KSPSetOperators(_ksp, _A, _A);
568 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
572 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
573 KSPSolve(_ksp, _b, _displacements);
575 KSPConvergedReason reason;
576 KSPGetConvergedReason(_ksp,&reason);
577 KSPGetIterationNumber(_ksp, &_PetscIts);
579 KSPGetResidualNorm(_ksp,&residu);
580 if (reason!=2 and reason!=3)
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();
590 resolutionOK = false;
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;
604 void LinearElasticityModel::setMesh(const Mesh &M)
606 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
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");
616 _Nmailles = _mesh.getNumberOfCells();
617 _Nnodes = _mesh.getNumberOfNodes();
619 cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
620 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
622 // find maximum nb of neibourghs
625 _VV=Field ("Displacements", CELLS, _mesh, _Ndim);
626 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
630 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
631 cout<<"1D Finite element method on segments"<<endl;
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;
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");
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;
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");
663 _VV=Field ("Temperature", NODES, _mesh, _Ndim);
665 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
666 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
667 _NboundaryNodes=_boundaryNodeIds.size();
673 void LinearElasticityModel::setLinearSolver(linearSolver kspType, preconditioner pcType)
675 //_maxPetscIts=maxIterationsPetsc;
676 // set linear solver algorithm
678 _ksptype = (char*)&KSPGMRES;
679 else if (kspType==CG)
680 _ksptype = (char*)&KSPCG;
681 else if (kspType==BCGS)
682 _ksptype = (char*)&KSPBCGS;
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 !!!");
689 // set preconditioner
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;
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 !!!" );
708 bool LinearElasticityModel::solveStationaryProblem()
710 if(!_initializedMemory)
712 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
713 _runLogFile->close();
714 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
716 bool stop=false; // Does the Problem want to stop (error) ?
718 cout<< "!!! Running test case "<< _fileName << " using ";
719 *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
723 cout<< "Finite volumes method"<<endl<<endl;
724 *_runLogFile<< "Finite volumes method"<<endl<<endl;
728 cout<< "Finite elements method"<<endl<<endl;
729 *_runLogFile<< "Finite elements method"<< endl<<endl;
732 computeStiffnessMatrix( 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");
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");
745 stop = !solveLinearSystem();
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");
753 _computationCompletedSuccessfully=true;
756 *_runLogFile<< "!!!!!! Computation successful"<< endl;
757 _runLogFile->close();
762 void LinearElasticityModel::save(){
763 cout<< "Saving numerical results"<<endl<<endl;
764 *_runLogFile<< "Saving numerical results"<< endl<<endl;
766 string resultFile(_path+"/LinearElasticityModel");//Results
769 resultFile+=_fileName;
771 // create mesh and component info
772 string suppress ="rm -rf "+resultFile+"_*";
773 system(suppress.c_str());//Nettoyage des précédents calculs identiques
775 if(_verbose or _system)
776 VecView(_displacements,PETSC_VIEWER_STDOUT_SELF);
780 for(int i=0; i<_Nmailles; i++)
782 for(int j=0; j<_nVar; j++)
785 VecGetValues(_displacements, 1, &k, &uk);
792 for(int i=0; i<_NunknownNodes; i++)
794 globalIndex = globalNodeIndex(i, _dirichletNodeIds);
795 for(int j=0; j<_nVar; j++)
798 VecGetValues(_displacements, 1, &k, &uk);
799 _VV(globalIndex,j)=uk;
805 for(int i=0; i<_NdirichletNodes; i++)
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];
817 _VV.writeVTK(resultFile);
820 _VV.writeMED(resultFile);
823 _VV.writeCSV(resultFile);
828 LinearElasticityModel::getOutputDisplacementField()
830 if(!_computationCompletedSuccessfully)
831 throw("Computation not performed yet or failed. No displacement field available");
836 void LinearElasticityModel::terminate()
838 VecDestroy(&_displacements);
843 LinearElasticityModel::setDirichletValues(map< int, double> dirichletBoundaryValues)
845 _dirichletValuesSet=true;
846 _dirichletBoundaryValues=dirichletBoundaryValues;
849 LinearElasticityModel::setNeumannValues(map< int, double> neumannBoundaryValues)
851 _neumannValuesSet=true;
852 _neumannBoundaryValues=neumannBoundaryValues;
856 LinearElasticityModel::getConditionNumber(bool isSingular, double tol) const
858 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
859 return A.getConditionNumber( isSingular, tol);