1 #include "StationaryDiffusionEquation.hxx"
2 #include "DiffusionEquation.hxx"
4 #include "SparseMatrixPetsc.hxx"
12 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda,MPI_Comm comm){
13 /* Initialisation of PETSC */
14 //check if PETSC is already initialised
15 PetscBool petscInitialized;
16 PetscInitialized(&petscInitialized);
18 {//check if MPI is already initialised
20 MPI_Initialized(&mpiInitialized);
22 PETSC_COMM_WORLD = comm;
23 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
25 MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
26 MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
27 PetscPrintf(PETSC_COMM_WORLD,"\n Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
28 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
29 PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
33 std::cout<<"Conductivity="<<lambda<<endl;
34 throw CdmathException("!!!!!!!!Error : conductivity parameter lambda cannot be negative");
38 std::cout<<"Space dimension="<<dim<<endl;
39 throw CdmathException("!!!!!!!!Error : parameter dim cannot be negative");
42 _FECalculation=FECalculation;
46 _nVar=1;//scalar prolem
48 _diffusionMatrixSet=false;
49 _initializedMemory=false;
57 _boundaryNodeIds=std::vector< int >(0);
58 _dirichletNodeIds=std::vector< int >(0);
62 _dirichletValuesSet=false;
63 _neumannValuesSet=false;
67 _precision_Newton=_precision;
68 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
72 int _PetscIts=0;//the number of iterations of the linear solver
73 _ksptype = (char*)&KSPGMRES;
74 _pctype = (char*)&PCILU;
76 _conditionNumber=false;
79 //parameters for monitoring simulation
82 _runLogFile=new ofstream;
84 //result save parameters
85 _fileName = "StationaryDiffusionProblem";
86 char result[ PATH_MAX ];//extracting current directory
87 getcwd(result, PATH_MAX );
88 _path=string( result );
90 _computationCompletedSuccessfully=false;
92 //heat transfer parameters
94 _fluidTemperatureFieldSet=false;
96 _heatPowerFieldSet=false;
97 _heatTransfertCoeff=0;
100 /* Default diffusion tensor is diagonal */
101 _DiffusionTensor=Matrix(_Ndim);
102 for(int idim=0;idim<_Ndim;idim++)
103 _DiffusionTensor(idim,idim)=_conductivity;
105 PetscPrintf(PETSC_COMM_WORLD,"\n Stationary diffusion problem with conductivity %.2f", lambda);
107 PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n\n");
109 PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n\n");
112 void StationaryDiffusionEquation::initialize()
116 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
119 throw CdmathException("!!!!!!!!StationaryDiffusionEquation::initialize() set mesh first");
122 cout<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
123 *_runLogFile<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
126 cout<< "Finite volumes method"<<endl<<endl;
127 *_runLogFile<< "Finite volumes method"<<endl<<endl;
131 cout<< "Finite elements method"<<endl<<endl;
132 *_runLogFile<< "Finite elements method"<<endl<<endl;
136 /**************** Field creation *********************/
138 if(!_heatPowerFieldSet){
139 _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
140 for(int i =0; i<_VV.getNumberOfElements(); i++)
141 _heatPowerField(i) = _heatSource;
142 _heatPowerFieldSet=true;
144 if(!_fluidTemperatureFieldSet){
145 _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
146 for(int i =0; i<_VV.getNumberOfElements(); i++)
147 _fluidTemperatureField(i) = _fluidTemperature;
148 _fluidTemperatureFieldSet=true;
151 /* Détection des noeuds frontière avec une condition limite de Dirichlet */
154 if(_NboundaryNodes==_Nnodes)
155 cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
157 for(int i=0; i<_NboundaryNodes; i++)
159 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
160 if( it != _dirichletBoundaryValues.end() )
161 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
162 else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
164 cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
165 *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
166 _runLogFile->close();
167 throw CdmathException("Missing boundary group");
169 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
171 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
172 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
173 *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
174 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
175 _runLogFile->close();
176 throw CdmathException("Missing boundary condition");
178 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
179 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
180 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
182 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
183 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
184 *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
185 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
186 _runLogFile->close();
187 throw CdmathException("Wrong boundary condition");
190 _NdirichletNodes=_dirichletNodeIds.size();
191 _NunknownNodes=_Nnodes - _NdirichletNodes;
192 cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
193 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
198 _globalNbUnknowns = _Nmailles*_nVar;
200 MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
201 _globalNbUnknowns = _NunknownNodes*_nVar;
204 /* Vectors creations */
205 VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
206 VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
207 VecSetFromOptions(_Tk);
208 VecGetLocalSize(_Tk, &_localNbUnknowns);
210 VecDuplicate(_Tk, &_Tkm1);
211 VecDuplicate(_Tk, &_deltaT);
212 VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
214 /* Matrix creation */
215 MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
217 /* Local sequential vector creation */
218 if(_mpi_size>1 && _mpi_rank == 0)
219 VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0
222 VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
224 //PETSc Linear solver
225 KSPCreate(PETSC_COMM_WORLD, &_ksp);
226 KSPSetType(_ksp, _ksptype);
227 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
228 KSPGetPC(_ksp, &_pc);
229 //PETSc preconditioner
231 PCSetType(_pc, _pctype);
234 PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
235 if(_pctype != (char*)&PCILU)//Default pc type is ilu
237 PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
238 PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
239 //If the above setvalue does not work, try the following
241 KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
244 int nlocal;//nb local blocs (should equal 1)
245 PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
248 KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
249 KSPGetPC(subKSP[0],&subpc);
250 PCSetType(subpc,_pctype);
253 throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
258 //Checking whether all boundary conditions are Neumann boundary condition
259 //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
260 if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
262 map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
263 while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
265 _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
269 _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
271 _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
273 //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
276 std::cout<<"### Warning : all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
277 std::cout<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
278 std::cout<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
279 *_runLogFile<<"### Warning : all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
280 *_runLogFile<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
281 *_runLogFile<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
284 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
285 MatSetNullSpace(_A, nullsp);
286 MatSetTransposeNullSpace(_A, nullsp);
287 MatNullSpaceDestroy(&nullsp);
288 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
289 //PCFactorSetShiftAmount(_pc,1e-10);
291 _initializedMemory=true;
294 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
295 if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
296 computeDiffusionMatrix(stop);
298 _dt_src=computeRHS(stop);
303 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
311 result=computeDiffusionMatrixFE(stop);
313 result=computeDiffusionMatrixFV(stop);
315 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
316 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
318 //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
319 //update value here if variable heat transfer coefficient
320 if(_heatTransfertCoeff>_precision)
321 MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
323 if(_verbose or _system)
324 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
329 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
336 Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
337 std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
338 std::vector< int > nodeIds(_Ndim+1);//cell node Ids
339 std::vector< Node > nodes(_Ndim+1);//cell nodes
340 int i_int, j_int; //index of nodes j and k considered as unknown nodes
341 bool dirichletCell_treated;
343 std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
344 for (int idim=0; idim<_Ndim+1;idim++)
345 values[idim][idim]=1;
347 /* parameters for boundary treatment */
348 vector< double > valuesBorder(_Ndim+1);
349 Vector GradShapeFuncBorder(_Ndim+1);
351 for (int j=0; j<_Nmailles;j++)
353 Cj = _mesh.getCell(j);
355 for (int idim=0; idim<_Ndim+1;idim++){
356 nodeIds[idim]=Cj.getNodeId(idim);
357 nodes[idim]=_mesh.getNode(nodeIds[idim]);
358 for (int jdim=0; jdim<_Ndim;jdim++)
359 M(idim,jdim)=nodes[idim].getPoint()[jdim];
362 for (int idim=0; idim<_Ndim+1;idim++)
363 GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
365 /* Loop on the edges of the cell */
366 for (int idim=0; idim<_Ndim+1;idim++)
368 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
369 {//First node of the edge is not Dirichlet node
370 i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
371 dirichletCell_treated=false;
372 for (int jdim=0; jdim<_Ndim+1;jdim++)
374 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
375 {//Second node of the edge is not Dirichlet node
376 j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
377 MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
379 else if (!dirichletCell_treated)
380 {//Second node of the edge is a Dirichlet node
381 dirichletCell_treated=true;
382 for (int kdim=0; kdim<_Ndim+1;kdim++)
384 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
385 if( it != _dirichletBoundaryValues.end() )
387 if( _dirichletValuesSet )
388 valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
390 valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
393 valuesBorder[kdim]=0;
395 GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
396 coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
397 VecSetValue(_b,i_int,coeff, ADD_VALUES);
404 //Calcul de la contribution de la condition limite de Neumann au second membre
405 if( _NdirichletNodes !=_NboundaryNodes)
407 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
408 int NboundaryFaces=boundaryFaces.size();
409 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
411 Face Fi = _mesh.getFace(boundaryFaces[i]);
412 for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
414 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
416 j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
417 if( _neumannValuesSet )
418 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
420 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
421 VecSetValue(_b, j_int, coeff, ADD_VALUES);
428 if(_onlyNeumannBC) //Check that the matrix is symmetric
430 PetscBool isSymetric;
431 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
432 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
433 MatIsSymmetric(_A,_precision,&isSymetric);
436 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
437 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
441 _diffusionMatrixSet=true;
447 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
450 long nbFaces = _mesh.getNumberOfFaces();
454 double inv_dxi, inv_dxj;
455 double barycenterDistance;
456 Vector normale(_Ndim);
459 std::vector< int > idCells;
461 for (int j=0; j<nbFaces;j++){
462 Fj = _mesh.getFace(j);
464 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
465 idCells = Fj.getCellsId();
466 Cell1 = _mesh.getCell(idCells[0]);
468 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
469 if (j == Cell1.getFacesId()[l]){
470 for (int idim = 0; idim < _Ndim; ++idim)
471 normale[idim] = Cell1.getNormalVector(l,idim);
476 //Compute velocity at the face Fj
477 dn=(_DiffusionTensor*normale)*normale;
479 // compute 1/dxi = volume of Ci/area of Fj
480 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
482 // If Fj is on the boundary
483 if (Fj.getNumberOfCells()==1) {
486 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
487 for(int p=0; p<_Ndim; p++)
488 cout << normale[p] << ",";
492 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
493 if( it != _dirichletBoundaryValues.end() )
495 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
496 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
497 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
501 nameOfGroup = Fj.getGroupName();
503 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
504 VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
506 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
507 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
508 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
509 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
513 cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
514 cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
515 cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
516 *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
517 _runLogFile->close();
518 throw CdmathException("Boundary condition not set");
521 // if Fj is inside the domain
522 } else if (Fj.getNumberOfCells()==2 ){
525 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
526 cout << " ; vecteur normal=(";
527 for(int p=0; p<_Ndim; p++)
528 cout << normale[p] << ",";
531 Cell2 = _mesh.getCell(idCells[1]);
534 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
536 inv_dxj = 1/Cell2.getMeasure();
538 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
540 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
541 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
542 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
543 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
547 *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
548 throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
553 if(_onlyNeumannBC) //Check that the matrix is symmetric
555 PetscBool isSymetric;
556 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
557 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
558 MatIsSymmetric(_A,_precision,&isSymetric);
561 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
562 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
566 _diffusionMatrixSet=true;
572 double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix
578 for (int i=0; i<_Nmailles;i++)
579 VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
583 std::vector< int > nodesId;
584 for (int i=0; i<_Nmailles;i++)
587 nodesId=Ci.getNodesId();
588 for (int j=0; j<nodesId.size();j++)
589 if(!_mesh.isBorderNode(nodesId[j]))
591 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
592 VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
597 VecAssemblyBegin(_b);
600 if(_verbose or _system)
601 VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
607 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
611 //Only implicit scheme considered
613 #if PETSC_VERSION_GREATER_3_5
614 KSPSetOperators(_ksp, _A, _A);
616 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
620 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
621 KSPSolve(_ksp, _b, _Tk);
623 KSPConvergedReason reason;
624 KSPGetConvergedReason(_ksp,&reason);
625 KSPGetIterationNumber(_ksp, &_PetscIts);
627 KSPGetResidualNorm(_ksp,&residu);
628 if (reason!=2 and reason!=3)
630 PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc.");
631 PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Itérations maximales %d atteintes, résidu = %1.2e, précision demandée= %1.2e",_maxPetscIts,residu,_precision);
632 PetscPrintf(PETSC_COMM_WORLD,"Solver used %s, preconditioner %s, Final number of iteration = %d",_ksptype,_pctype,_PetscIts);
633 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
634 *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
635 *_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
636 _runLogFile->close();
641 if( _MaxIterLinearSolver < _PetscIts)
642 _MaxIterLinearSolver = _PetscIts;
643 PetscPrintf(PETSC_COMM_WORLD,"## Système linéaire résolu en %d itérations par le solveur %s et le preconditioneur %s, précision demandée = %1.2e",_PetscIts,_ksptype,_pctype,_precision);
644 *_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;
645 VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
646 VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
648 VecAssemblyBegin(_Tk);
649 VecAssemblyEnd( _Tk);
650 VecAssemblyBegin(_deltaT);
651 VecAssemblyEnd( _deltaT);
654 PetscPrintf(PETSC_COMM_WORLD,"Début calcul de l'erreur maximale");
656 VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
659 PetscPrintf(PETSC_COMM_WORLD,"Fin calcul de la variation relative, erreur maximale : %1.2e", _erreur_rel );
661 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
669 void StationaryDiffusionEquation::setMesh(const Mesh &M)
673 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
675 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
676 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
677 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
678 _runLogFile->close();
679 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
683 _Nmailles = _mesh.getNumberOfCells();
684 _Nnodes = _mesh.getNumberOfNodes();
686 cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
687 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
689 // find maximum nb of neibourghs
692 _VV=Field ("Temperature", CELLS, _mesh, 1);
693 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
697 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
698 cout<<"1D Finite element method on segments"<<endl;
701 if( _mesh.isTriangular() )//Mesh dim=2
702 cout<<"2D Finite element method on triangles"<<endl;
703 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
704 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
707 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
708 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
709 _runLogFile->close();
710 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
715 if( _mesh.isTetrahedral() )//Mesh dim=3
716 cout<<"3D Finite element method on tetrahedra"<<endl;
717 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
718 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
719 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
720 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
723 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;
724 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
725 _runLogFile->close();
726 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
730 _VV=Field ("Temperature", NODES, _mesh, 1);
732 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
733 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
734 _NboundaryNodes=_boundaryNodeIds.size();
738 /* MPI distribution parameters */
739 int nbVoisinsMax;//Mettre en attribut ?
741 MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
742 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
743 nbVoisinsMax = _neibMaxNbCells;
746 MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
747 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
748 nbVoisinsMax = _neibMaxNbNodes;
750 _d_nnz = (nbVoisinsMax+1)*_nVar;
751 _o_nnz = nbVoisinsMax *_nVar;
756 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
758 //_maxPetscIts=maxIterationsPetsc;
759 // set linear solver algorithm
761 _ksptype = (char*)&KSPGMRES;
762 else if (kspType==CG)
763 _ksptype = (char*)&KSPCG;
764 else if (kspType==BCGS)
765 _ksptype = (char*)&KSPBCGS;
767 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
768 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
769 _runLogFile->close();
770 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
772 // set preconditioner
774 _pctype = (char*)&PCNONE;
775 else if (pcType ==LU)
776 _pctype = (char*)&PCLU;
777 else if (pcType == ILU)
778 _pctype = (char*)&PCILU;
779 else if (pcType ==CHOLESKY)
780 _pctype = (char*)&PCCHOLESKY;
781 else if (pcType == ICC)
782 _pctype = (char*)&PCICC;
784 cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
785 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
786 _runLogFile->close();
787 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
791 bool StationaryDiffusionEquation::solveStationaryProblem()
793 if(!_initializedMemory)
795 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
796 _runLogFile->close();
797 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
799 bool stop=false; // Does the Problem want to stop (error) ?
800 bool converged=false; // has the newton scheme converged (end) ?
802 PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName);
803 *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
807 PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
808 *_runLogFile<< "Finite volumes method"<<endl<<endl;
812 PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
813 *_runLogFile<< "Finite elements method"<< endl<<endl;
816 computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
818 PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing diffusion matrix, stopping calculation\n");
819 *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
820 _runLogFile->close();
821 throw CdmathException("Failed computing diffusion matrix");
823 computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
825 PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n");
826 *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
827 throw CdmathException("Failed computing right hand side");
829 stop = iterateNewtonStep(converged);
831 PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n");
832 *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
833 _runLogFile->close();
834 throw CdmathException("Failed solving linear system");
837 _computationCompletedSuccessfully=true;
840 // Newton iteration loop for non linear problems
842 while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
844 computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
845 computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
846 stop = iterateNewtonStep(converged);
850 cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
851 *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
852 throw CdmathException("Failed solving a Newton iteration");
854 else if(_NEWTON_its==_maxNewtonIts){
855 cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
856 *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
857 throw CdmathException("No convergence of Newton scheme");
860 cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
861 *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
866 *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
867 _runLogFile->close();
872 void StationaryDiffusionEquation::save(){
873 PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
874 *_runLogFile<< "Saving numerical results"<< endl<<endl;
876 string resultFile(_path+"/StationaryDiffusionEquation");//Results
879 resultFile+=_fileName;
881 // create mesh and component info
882 string suppress ="rm -rf "+resultFile+"_*";
883 system(suppress.c_str());//Nettoyage des précédents calculs identiques
886 VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
887 VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
890 if(_verbose or _system)
891 VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD);
896 for(int i=0; i<_Nmailles; i++)
899 VecGetValues(_Tk_seq, 1, &i, &Ti);
901 VecGetValues(_Tk , 1, &i, &Ti);
907 for(int i=0; i<_NunknownNodes; i++)
910 VecGetValues(_Tk_seq, 1, &i, &Ti);
912 VecGetValues(_Tk , 1, &i, &Ti);
913 globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
919 for(int i=0; i<_NdirichletNodes; i++)
921 Ni=_mesh.getNode(_dirichletNodeIds[i]);
922 nameOfGroup = Ni.getGroupName();
923 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
927 _VV.setInfoOnComponent(0,"Temperature_(K)");
931 _VV.writeVTK(resultFile);
934 _VV.writeMED(resultFile);
937 _VV.writeCSV(resultFile);
943 void StationaryDiffusionEquation::terminate()
947 VecDestroy(&_deltaT);
950 if(_mpi_size>1 && _mpi_rank == 0)
951 VecDestroy(&_Tk_seq);
953 PetscBool petscInitialized;
954 PetscInitialized(&petscInitialized);
961 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
963 _dirichletValuesSet=true;
964 _dirichletBoundaryValues=dirichletBoundaryValues;
968 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
970 _neumannValuesSet=true;
971 _neumannBoundaryValues=neumannBoundaryValues;
975 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
977 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
978 return A.getConditionNumber( isSingular, tol);
980 std::vector< double >
981 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
983 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
985 if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
987 Vector nodal_volumes(_NunknownNodes);
989 for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
991 Cell Ci = _mesh.getCell(i);
992 for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
994 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
996 j_int=DiffusionEquation::unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
997 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
1001 for( j_int = 0; j_int< _NunknownNodes ; j_int++)
1002 nodal_volumes[j_int]=1/nodal_volumes[j_int];
1003 A.leftDiagonalScale(nodal_volumes);
1006 return A.getEigenvalues( nev, which, tol);
1008 std::vector< Vector >
1009 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
1011 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
1012 return A.getEigenvectors( nev, which, tol);
1015 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
1017 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
1018 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
1019 Field my_eigenfield;
1021 my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
1023 my_eigenfield.setFieldByDataArrayDouble(d);
1025 return my_eigenfield;
1029 StationaryDiffusionEquation::getOutputTemperatureField()
1031 if(!_computationCompletedSuccessfully)
1032 throw("Computation not performed yet or failed. No temperature field available");
1038 StationaryDiffusionEquation::getRodTemperatureField()
1040 return getOutputTemperatureField();
1044 StationaryDiffusionEquation::getInputFieldsNames()
1046 vector<string> result(2);
1048 result[0]="FluidTemperature";
1049 result[1]="HeatPower";
1054 StationaryDiffusionEquation::getOutputFieldsNames()
1056 vector<string> result(1);
1058 result[0]="RodTemperature";
1064 StationaryDiffusionEquation::getOutputField(const string& nameField )
1066 if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
1067 return getRodTemperatureField();
1070 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
1071 throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
1076 StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField )
1079 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setInputField set the mesh first");
1081 if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
1082 return setFluidTemperatureField( inputField) ;
1083 else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
1084 return setHeatPowerField( inputField );
1087 cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
1088 throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name");
1093 StationaryDiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
1095 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setFluidTemperatureField set initial field first");
1097 coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
1098 _fluidTemperatureField=coupledTemperatureField;
1099 _fluidTemperatureFieldSet=true;
1103 StationaryDiffusionEquation::setHeatPowerField(Field heatPower){
1105 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1107 heatPower.getMesh().checkFastEquivalWith(_mesh);
1108 _heatPowerField=heatPower;
1109 _heatPowerFieldSet=true;
1113 StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){
1115 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1117 _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel);
1118 _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
1119 _heatPowerFieldSet=true;
1123 StationaryDiffusionEquation::setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
1124 if(_FECalculation && field_support_type != NODES)
1125 cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
1126 else if(!_FECalculation && field_support_type == NODES)
1127 cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
1131 switch(field_support_type)
1134 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
1137 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
1140 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
1143 std::ostringstream message;
1144 message << "Error StationaryDiffusionEquation::setDirichletBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
1145 throw CdmathException(message.str().c_str());
1147 /* For the moment the boundary value is taken constant equal to zero */
1148 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
1152 StationaryDiffusionEquation::setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
1153 if(_FECalculation && field_support_type != NODES)
1154 cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
1155 else if(!_FECalculation && field_support_type == NODES)
1156 cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
1160 switch(field_support_type)
1163 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
1166 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
1169 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
1172 std::ostringstream message;
1173 message << "Error StationaryDiffusionEquation::setNeumannBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
1174 throw CdmathException(message.str().c_str());
1176 /* For the moment the boundary value is taken constant equal to zero */
1177 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab