#include "StationaryDiffusionEquation.hxx"
+#include "DiffusionEquation.hxx"
#include "Node.hxx"
#include "SparseMatrixPetsc.hxx"
#include "math.h"
using namespace std;
-int StationaryDiffusionEquation::fact(int n)
-{
- return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
-}
-int StationaryDiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const
-{//assumes Dirichlet node numbering is strictly increasing
- int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
- int boundarySize=dirichletNodes.size();
- while(j<boundarySize and dirichletNodes[j]<globalIndex)
- j++;
- if(j==boundarySize)
- return globalIndex-boundarySize;
- else if (dirichletNodes[j]>globalIndex)
- return globalIndex-j;
- else
- throw CdmathException("StationaryDiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
-}
-
-int StationaryDiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes) const
-{//assumes Dirichlet boundary node numbering is strictly increasing
- int boundarySize=dirichletNodes.size();
- /* trivial case where all boundary nodes are Neumann BC */
- if(boundarySize==0)
- return unknownNodeIndex;
-
- double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
- int j=0;//indice de parcours des noeuds frontière
- //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
- while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
- {
- unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
- j++;
- }
-
- if(j+1==boundarySize)
- return unknownNodeIndex+boundarySize;
- else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
- return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
-}
-
-StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda){
+StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda,MPI_Comm comm){
+ /* Initialisation of PETSC */
+ //check if PETSC is already initialised
PetscBool petscInitialized;
PetscInitialized(&petscInitialized);
if(!petscInitialized)
- PetscInitialize(NULL,NULL,0,0);
+ {//check if MPI is already initialised
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+ if(mpiInitialized)
+ PETSC_COMM_WORLD = comm;
+ PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
+ }
+ MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
+ MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
+ 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.
+ 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.
+ PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
if(lambda < 0.)
{
std::cout<<"Conductivity="<<lambda<<endl;
- throw CdmathException("Error : conductivity parameter lambda cannot be negative");
+ throw CdmathException("!!!!!!!!Error : conductivity parameter lambda cannot be negative");
}
if(dim<=0)
{
std::cout<<"Space dimension="<<dim<<endl;
- throw CdmathException("Error : parameter dim cannot be negative");
+ throw CdmathException("!!!!!!!!Error : parameter dim cannot be negative");
}
_FECalculation=FECalculation;
_NEWTON_its=0;
int _PetscIts=0;//the number of iterations of the linear solver
_ksptype = (char*)&KSPGMRES;
- _pctype = (char*)&PCLU;
+ _pctype = (char*)&PCILU;
+
_conditionNumber=false;
_erreur_rel= 0;
_heatPowerFieldSet=false;
_heatTransfertCoeff=0;
_heatSource=0;
-}
-
-void StationaryDiffusionEquation::initialize()
-{
- _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
- if(!_meshSet)
- throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
- else
- {
- cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
- *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
- if(!_FECalculation)
- {
- cout<< "Finite volumes method"<<endl<<endl;
- *_runLogFile<< "Finite volumes method"<<endl<<endl;
- }
- else
- {
- cout<< "Finite elements method"<<endl<<endl;
- *_runLogFile<< "Finite elements method"<<endl<<endl;
- }
- }
-
- _DiffusionTensor=Matrix(_Ndim);
+ /* Default diffusion tensor is diagonal */
+ _DiffusionTensor=Matrix(_Ndim);
for(int idim=0;idim<_Ndim;idim++)
- _DiffusionTensor(idim,idim)=1;
- /**************** Field creation *********************/
-
- if(!_heatPowerFieldSet){
- if(_FECalculation){
- _heatPowerField=Field("Heat power",NODES,_mesh,1);
- for(int i =0; i<_Nnodes; i++)
- _heatPowerField(i) = _heatSource;
- }
- else{
- _heatPowerField=Field("Heat power",CELLS,_mesh,1);
- for(int i =0; i<_Nmailles; i++)
- _heatPowerField(i) = _heatSource;
- }
- _heatPowerFieldSet=true;
- }
- if(!_fluidTemperatureFieldSet){
- if(_FECalculation){
- _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
- for(int i =0; i<_Nnodes; i++)
- _fluidTemperatureField(i) = _fluidTemperature;
- }
- else{
- _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
- for(int i =0; i<_Nmailles; i++)
- _fluidTemperatureField(i) = _fluidTemperature;
- }
- _fluidTemperatureFieldSet=true;
- }
+ _DiffusionTensor(idim,idim)=_conductivity;
- /* Détection des noeuds frontière avec une condition limite de Dirichlet */
- if(_FECalculation)
- {
- if(_NboundaryNodes==_Nnodes)
- cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
-
- for(int i=0; i<_NboundaryNodes; i++)
- {
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
- if( it != _dirichletBoundaryValues.end() )
- _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
- {
- cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
- *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
- _runLogFile->close();
- throw CdmathException("Missing boundary group");
- }
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
- {
- cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
- *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
- _runLogFile->close();
- throw CdmathException("Missing boundary condition");
- }
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
- _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
- {
- cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
- *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
- _runLogFile->close();
- throw CdmathException("Wrong boundary condition");
- }
- }
- _NdirichletNodes=_dirichletNodeIds.size();
- _NunknownNodes=_Nnodes - _NdirichletNodes;
- cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
- *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
- }
-
- //creation de la matrice
- if(!_FECalculation)
- MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
+ PetscPrintf(PETSC_COMM_WORLD,"\n Stationary diffusion problem with conductivity %.2f", lambda);
+ if(FECalculation)
+ PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n\n");
else
- MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
-
- VecCreate(PETSC_COMM_SELF, &_Tk);
+ PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n\n");
+}
+void StationaryDiffusionEquation::initialize()
+{
+ if(_mpi_rank==0)
+ {
+ _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
+
+ if(!_meshSet)
+ throw CdmathException("!!!!!!!!StationaryDiffusionEquation::initialize() set mesh first");
+ else
+ {
+ cout<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
+ *_runLogFile<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
+ if(!_FECalculation)
+ {
+ cout<< "Finite volumes method"<<endl<<endl;
+ *_runLogFile<< "Finite volumes method"<<endl<<endl;
+ }
+ else
+ {
+ cout<< "Finite elements method"<<endl<<endl;
+ *_runLogFile<< "Finite elements method"<<endl<<endl;
+ }
+ }
+
+ /**************** Field creation *********************/
+
+ if(!_heatPowerFieldSet){
+ _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _heatPowerField(i) = _heatSource;
+ _heatPowerFieldSet=true;
+ }
+ if(!_fluidTemperatureFieldSet){
+ _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _fluidTemperatureField(i) = _fluidTemperature;
+ _fluidTemperatureFieldSet=true;
+ }
+
+ /* Détection des noeuds frontière avec une condition limite de Dirichlet */
+ if(_FECalculation)
+ {
+ if(_NboundaryNodes==_Nnodes)
+ cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
+
+ for(int i=0; i<_NboundaryNodes; i++)
+ {
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
+ if( it != _dirichletBoundaryValues.end() )
+ _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
+ else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
+ {
+ cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
+ *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
+ _runLogFile->close();
+ throw CdmathException("Missing boundary group");
+ }
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
+ {
+ cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
+ cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
+ *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
+ *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
+ _runLogFile->close();
+ throw CdmathException("Missing boundary condition");
+ }
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
+ _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
+ {
+ cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
+ cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
+ *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
+ *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
+ _runLogFile->close();
+ throw CdmathException("Wrong boundary condition");
+ }
+ }
+ _NdirichletNodes=_dirichletNodeIds.size();
+ _NunknownNodes=_Nnodes - _NdirichletNodes;
+ cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ }
+ }
+
if(!_FECalculation)
- VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
- else
- VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
-
+ _globalNbUnknowns = _Nmailles*_nVar;
+ else{
+ MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ _globalNbUnknowns = _NunknownNodes*_nVar;
+ }
+
+ /* Vectors creations */
+ VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
+ VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
VecSetFromOptions(_Tk);
+ VecGetLocalSize(_Tk, &_localNbUnknowns);
+
VecDuplicate(_Tk, &_Tkm1);
VecDuplicate(_Tk, &_deltaT);
- VecDuplicate(_Tk, &_b);//RHS of the linear system
+ VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
- //Linear solver
- KSPCreate(PETSC_COMM_SELF, &_ksp);
+ /* Matrix creation */
+ MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
+
+ /* Local sequential vector creation */
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0
+ if(_mpi_size==0)
+ _Tk_seq=_Tk;
+ VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
+
+ //PETSc Linear solver
+ KSPCreate(PETSC_COMM_WORLD, &_ksp);
KSPSetType(_ksp, _ksptype);
- // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
KSPGetPC(_ksp, &_pc);
- PCSetType(_pc, _pctype);
+ //PETSc preconditioner
+ if(_mpi_size==1 )
+ PCSetType(_pc, _pctype);
+ else
+ {
+ PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
+ if(_pctype != (char*)&PCILU)//Default pc type is ilu
+ {
+ PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
+ PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
+ //If the above setvalue does not work, try the following
+ /*
+ KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
+ KSP * subKSP;
+ PC subpc;
+ int nlocal;//nb local blocs (should equal 1)
+ PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
+ if(nlocal==1)
+ {
+ KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
+ KSPGetPC(subKSP[0],&subpc);
+ PCSetType(subpc,_pctype);
+ }
+ else
+ throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
+ */
+ }
+ }
//Checking whether all boundary conditions are Neumann boundary condition
//if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
*_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;
*_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;
- //Check that the matrix is symmetric
- PetscBool isSymetric;
- MatIsSymmetric(_A,_precision,&isSymetric);
- if(!isSymetric)
- {
- cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
- throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
- }
MatNullSpace nullsp;
MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
MatSetNullSpace(_A, nullsp);
//PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
//PCFactorSetShiftAmount(_pc,1e-10);
}
-
_initializedMemory=true;
-
}
double StationaryDiffusionEquation::computeTimeStep(bool & stop){
return _dt_src;
}
-Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
- vector< Matrix > matrices(_Ndim);
-
- for (int idim=0; idim<_Ndim;idim++){
- matrices[idim]=M.deepCopy();
- for (int jdim=0; jdim<_Ndim+1;jdim++)
- matrices[idim](jdim,idim) = values[jdim] ;
- }
-
- Vector result(_Ndim);
- for (int idim=0; idim<_Ndim;idim++)
- result[idim] = matrices[idim].determinant();
-
- return result;
-}
-
double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
{
double result;
+ MatZeroEntries(_A);
+ VecZeroEntries(_b);
+
if(_FECalculation)
result=computeDiffusionMatrixFE(stop);
else
result=computeDiffusionMatrixFV(stop);
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+
//Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
//update value here if variable heat transfer coefficient
if(_heatTransfertCoeff>_precision)
MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
if(_verbose or _system)
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
return result;
}
double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
- Cell Cj;
- string nameOfGroup;
- double dn, coeff;
- MatZeroEntries(_A);
- VecZeroEntries(_b);
-
- Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
- std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
- std::vector< int > nodeIds(_Ndim+1);//cell node Ids
- std::vector< Node > nodes(_Ndim+1);//cell nodes
- int i_int, j_int; //index of nodes j and k considered as unknown nodes
- bool dirichletCell_treated;
-
- std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
- for (int idim=0; idim<_Ndim+1;idim++)
- values[idim][idim]=1;
-
- /* parameters for boundary treatment */
- vector< double > valuesBorder(_Ndim+1);
- Vector GradShapeFuncBorder(_Ndim+1);
-
- for (int j=0; j<_Nmailles;j++)
- {
- Cj = _mesh.getCell(j);
-
- for (int idim=0; idim<_Ndim+1;idim++){
- nodeIds[idim]=Cj.getNodeId(idim);
- nodes[idim]=_mesh.getNode(nodeIds[idim]);
- for (int jdim=0; jdim<_Ndim;jdim++)
- M(idim,jdim)=nodes[idim].getPoint()[jdim];
- M(idim,_Ndim)=1;
- }
- for (int idim=0; idim<_Ndim+1;idim++)
- GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
+ if(_mpi_rank == 0)
+ {
+ Cell Cj;
+ string nameOfGroup;
+ double coeff;
+
+ Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
+ std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
+ std::vector< int > nodeIds(_Ndim+1);//cell node Ids
+ std::vector< Node > nodes(_Ndim+1);//cell nodes
+ int i_int, j_int; //index of nodes j and k considered as unknown nodes
+ bool dirichletCell_treated;
+
+ std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
+ for (int idim=0; idim<_Ndim+1;idim++)
+ values[idim][idim]=1;
+
+ /* parameters for boundary treatment */
+ vector< double > valuesBorder(_Ndim+1);
+ Vector GradShapeFuncBorder(_Ndim+1);
+
+ for (int j=0; j<_Nmailles;j++)
+ {
+ Cj = _mesh.getCell(j);
+
+ for (int idim=0; idim<_Ndim+1;idim++){
+ nodeIds[idim]=Cj.getNodeId(idim);
+ nodes[idim]=_mesh.getNode(nodeIds[idim]);
+ for (int jdim=0; jdim<_Ndim;jdim++)
+ M(idim,jdim)=nodes[idim].getPoint()[jdim];
+ M(idim,_Ndim)=1;
+ }
+ for (int idim=0; idim<_Ndim+1;idim++)
+ GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
+
+ /* Loop on the edges of the cell */
+ for (int idim=0; idim<_Ndim+1;idim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
+ {//First node of the edge is not Dirichlet node
+ i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ dirichletCell_treated=false;
+ for (int jdim=0; jdim<_Ndim+1;jdim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
+ {//Second node of the edge is not Dirichlet node
+ j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
+ }
+ else if (!dirichletCell_treated)
+ {//Second node of the edge is a Dirichlet node
+ dirichletCell_treated=true;
+ for (int kdim=0; kdim<_Ndim+1;kdim++)
+ {
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
+ if( it != _dirichletBoundaryValues.end() )
+ {
+ if( _dirichletValuesSet )
+ valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
+ else
+ valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
+ }
+ else
+ valuesBorder[kdim]=0;
+ }
+ GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
+ coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
+ VecSetValue(_b,i_int,coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+ }
+
+ //Calcul de la contribution de la condition limite de Neumann au second membre
+ if( _NdirichletNodes !=_NboundaryNodes)
+ {
+ vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
+ int NboundaryFaces=boundaryFaces.size();
+ for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
+ {
+ Face Fi = _mesh.getFace(boundaryFaces[i]);
+ for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
+ {
+ j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ if( _neumannValuesSet )
+ coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
+ else
+ coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
+ VecSetValue(_b, j_int, coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+ }
- /* Loop on the edges of the cell */
- for (int idim=0; idim<_Ndim+1;idim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
- {//First node of the edge is not Dirichlet node
- i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- dirichletCell_treated=false;
- for (int jdim=0; jdim<_Ndim+1;jdim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
- {//Second node of the edge is not Dirichlet node
- j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
- }
- else if (!dirichletCell_treated)
- {//Second node of the edge is a Dirichlet node
- dirichletCell_treated=true;
- for (int kdim=0; kdim<_Ndim+1;kdim++)
- {
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
- if( it != _dirichletBoundaryValues.end() )
- {
- if( _dirichletValuesSet )
- valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
- else
- valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
- }
- else
- valuesBorder[kdim]=0;
- }
- GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
- coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
- VecSetValue(_b,i_int,coeff, ADD_VALUES);
- }
- }
- }
- }
+ if(_onlyNeumannBC) //Check that the matrix is symmetric
+ {
+ PetscBool isSymetric;
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+ MatIsSymmetric(_A,_precision,&isSymetric);
+ if(!isSymetric)
+ {
+ cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
+ throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
+ }
}
-
- //Calcul de la contribution de la condition limite de Neumann au second membre
- if( _NdirichletNodes !=_NboundaryNodes)
- {
- vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
- int NboundaryFaces=boundaryFaces.size();
- for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
- {
- Face Fi = _mesh.getFace(i);
- for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
- {
- j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
- if( _neumannValuesSet )
- coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
- else
- coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
- VecSetValue(_b, j_int, coeff, ADD_VALUES);
- }
- }
- }
- }
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
_diffusionMatrixSet=true;
stop=false ;
}
double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
- long nbFaces = _mesh.getNumberOfFaces();
- Face Fj;
- Cell Cell1,Cell2;
- string nameOfGroup;
- double inv_dxi, inv_dxj;
- double barycenterDistance;
- Vector normale(_Ndim);
- double dn;
- PetscInt idm, idn;
- std::vector< int > idCells;
- MatZeroEntries(_A);
- VecZeroEntries(_b);
-
- for (int j=0; j<nbFaces;j++){
- Fj = _mesh.getFace(j);
-
- // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
- idCells = Fj.getCellsId();
- Cell1 = _mesh.getCell(idCells[0]);
- idm = idCells[0];
- for(int l=0; l<Cell1.getNumberOfFaces(); l++){
- if (j == Cell1.getFacesId()[l]){
- for (int idim = 0; idim < _Ndim; ++idim)
- normale[idim] = Cell1.getNormalVector(l,idim);
- break;
- }
- }
-
- //Compute velocity at the face Fj
- dn=_conductivity*(_DiffusionTensor*normale)*normale;
-
- // compute 1/dxi = volume of Ci/area of Fj
- inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
-
- // If Fj is on the boundary
- if (Fj.getNumberOfCells()==1) {
- if(_verbose )
- {
- cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
- }
-
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
- if( it != _dirichletBoundaryValues.end() )
- {
- barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
- MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
- }
- else
- {
- nameOfGroup = Fj.getGroupName();
-
- if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
- VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
- }
- else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
- barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
- MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
- }
- else {
- stop=true ;
- cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
- cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
- *_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- _runLogFile->close();
- throw CdmathException("Boundary condition not accepted");
- }
- }
- // if Fj is inside the domain
- } else if (Fj.getNumberOfCells()==2 ){
- if(_verbose )
- {
- cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
- cout << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
+ if(_mpi_rank == 0)
+ {
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Cell1,Cell2;
+ string nameOfGroup;
+ double inv_dxi, inv_dxj;
+ double barycenterDistance;
+ Vector normale(_Ndim);
+ double dn;
+ PetscInt idm, idn;
+ std::vector< int > idCells;
+
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+
+ // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
+ idCells = Fj.getCellsId();
+ Cell1 = _mesh.getCell(idCells[0]);
+ idm = idCells[0];
+ for(int l=0; l<Cell1.getNumberOfFaces(); l++){
+ if (j == Cell1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ normale[idim] = Cell1.getNormalVector(l,idim);
+ break;
+ }
+ }
+
+ //Compute velocity at the face Fj
+ dn=(_DiffusionTensor*normale)*normale;
+
+ // compute 1/dxi = volume of Ci/area of Fj
+ inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
+
+ // If Fj is on the boundary
+ if (Fj.getNumberOfCells()==1) {
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
+ if( it != _dirichletBoundaryValues.end() )
+ {
+ barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
+ MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
+ }
+ else
+ {
+ nameOfGroup = Fj.getGroupName();
+
+ if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
+ VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
+ }
+ else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
+ barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
+ MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
+ }
+ else {
+ stop=true ;
+ cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
+ cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
+ *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ _runLogFile->close();
+ throw CdmathException("Boundary condition not set");
+ }
+ }
+ // if Fj is inside the domain
+ } else if (Fj.getNumberOfCells()==2 ){
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+ Cell2 = _mesh.getCell(idCells[1]);
+ idn = idCells[1];
+ if (_Ndim > 1)
+ inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
+ else
+ inv_dxj = 1/Cell2.getMeasure();
+
+ barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
+
+ MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
}
- Cell2 = _mesh.getCell(idCells[1]);
- idn = idCells[1];
- if (_Ndim > 1)
- inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
else
- inv_dxj = 1/Cell2.getMeasure();
-
- barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
-
- MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ {
+ *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
+ throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
+ }
}
- else
- {
- *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
- throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
- }
}
-
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
+ if(_onlyNeumannBC) //Check that the matrix is symmetric
+ {
+ PetscBool isSymetric;
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+ MatIsSymmetric(_A,_precision,&isSymetric);
+ if(!isSymetric)
+ {
+ cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
+ throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
+ }
+ }
+
_diffusionMatrixSet=true;
stop=false ;
return INFINITY;
}
-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)
+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
{
- VecAssemblyBegin(_b);
- if(!_FECalculation)
- for (int i=0; i<_Nmailles;i++)
- VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
- else
- {
- Cell Ci;
- std::vector< int > nodesId;
- for (int i=0; i<_Nmailles;i++)
- {
- Ci=_mesh.getCell(i);
- nodesId=Ci.getNodesId();
- for (int j=0; j<nodesId.size();j++)
- if(!_mesh.isBorderNode(nodesId[j]))
- {
- double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
- VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
- }
- }
- }
-
+ if(_mpi_rank == 0)
+ {
+ if(!_FECalculation)
+ for (int i=0; i<_Nmailles;i++)
+ VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
+ else
+ {
+ Cell Ci;
+ std::vector< int > nodesId;
+ for (int i=0; i<_Nmailles;i++)
+ {
+ Ci=_mesh.getCell(i);
+ nodesId=Ci.getNodesId();
+ for (int j=0; j<nodesId.size();j++)
+ if(!_mesh.isBorderNode(nodesId[j]))
+ {
+ double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
+ VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
+ }
+ }
+ }
+ }
+ VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
if(_verbose or _system)
- VecView(_b,PETSC_VIEWER_STDOUT_SELF);
+ VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
stop=false ;
return INFINITY;
bool stop;
//Only implicit scheme considered
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
#if PETSC_VERSION_GREATER_3_5
KSPSetOperators(_ksp, _A, _A);
KSPGetResidualNorm(_ksp,&residu);
if (reason!=2 and reason!=3)
{
- cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
- cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
- cout<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc.");
+ PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Itérations maximales %d atteintes, résidu = %1.2e, précision demandée= %1.2e",_maxPetscIts,residu,_precision);
+ PetscPrintf(PETSC_COMM_WORLD,"Solver used %s, preconditioner %s, Final number of iteration = %d",_ksptype,_pctype,_PetscIts);
*_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
*_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
*_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
else{
if( _MaxIterLinearSolver < _PetscIts)
_MaxIterLinearSolver = _PetscIts;
- 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;
+ 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);
*_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;
VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
- _erreur_rel= 0;
- double Ti, dTi;
-
VecAssemblyBegin(_Tk);
VecAssemblyEnd( _Tk);
VecAssemblyBegin(_deltaT);
VecAssemblyEnd( _deltaT);
- if(!_FECalculation)
- for(int i=0; i<_Nmailles; i++)
- {
- VecGetValues(_deltaT, 1, &i, &dTi);
- VecGetValues(_Tk, 1, &i, &Ti);
- if(_erreur_rel < fabs(dTi/Ti))
- _erreur_rel = fabs(dTi/Ti);
- }
- else
- for(int i=0; i<_NunknownNodes; i++)
- {
- VecGetValues(_deltaT, 1, &i, &dTi);
- VecGetValues(_Tk, 1, &i, &Ti);
- if(_erreur_rel < fabs(dTi/Ti))
- _erreur_rel = fabs(dTi/Ti);
- }
+ if(_verbose)
+ PetscPrintf(PETSC_COMM_WORLD,"Début calcul de l'erreur maximale");
+
+ VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
+
+ if(_verbose)
+ PetscPrintf(PETSC_COMM_WORLD,"Fin calcul de la variation relative, erreur maximale : %1.2e", _erreur_rel );
stop=false;
converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
void StationaryDiffusionEquation::setMesh(const Mesh &M)
{
- if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
+ if(_mpi_rank==0)
{
- cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
- *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
+ if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
+ {
+ cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
+ *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
+ }
+
+ _mesh=M;
+ _Nmailles = _mesh.getNumberOfCells();
+ _Nnodes = _mesh.getNumberOfNodes();
+
+ cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
+ *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
+
+ // find maximum nb of neibourghs
+ if(!_FECalculation)
+ {
+ _VV=Field ("Temperature", CELLS, _mesh, 1);
+ _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
+ }
+ else
+ {
+ if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
+ cout<<"1D Finite element method on segments"<<endl;
+ else if(_Ndim==2)
+ {
+ if( _mesh.isTriangular() )//Mesh dim=2
+ cout<<"2D Finite element method on triangles"<<endl;
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else
+ {
+ cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
+ }
+ }
+ else if(_Ndim==3)
+ {
+ if( _mesh.isTetrahedral() )//Mesh dim=3
+ cout<<"3D Finite element method on tetrahedra"<<endl;
+ else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
+ cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else
+ {
+ 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;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
+ }
+ }
+
+ _VV=Field ("Temperature", NODES, _mesh, 1);
+
+ _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
+ _boundaryNodeIds = _mesh.getBoundaryNodeIds();
+ _NboundaryNodes=_boundaryNodeIds.size();
+ }
}
-
- _mesh=M;
- _Nmailles = _mesh.getNumberOfCells();
- _Nnodes = _mesh.getNumberOfNodes();
-
- cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
- *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
-
- // find maximum nb of neibourghs
- if(!_FECalculation)
- {
- _VV=Field ("Temperature", CELLS, _mesh, 1);
- _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
- }
- else
- {
- if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
- cout<<"1D Finite element method on segments"<<endl;
- else if(_Ndim==2)
- {
- if( _mesh.isTriangular() )//Mesh dim=2
- cout<<"2D Finite element method on triangles"<<endl;
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else
- {
- cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
- }
- }
- else if(_Ndim==3)
- {
- if( _mesh.isTetrahedral() )//Mesh dim=3
- cout<<"3D Finite element method on tetrahedra"<<endl;
- else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
- cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else
- {
- 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;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
- }
- }
-
- _VV=Field ("Temperature", NODES, _mesh, 1);
-
- _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
- _boundaryNodeIds = _mesh.getBoundaryNodeIds();
- _NboundaryNodes=_boundaryNodeIds.size();
- }
+
+ /* MPI distribution parameters */
+ int nbVoisinsMax;//Mettre en attribut ?
+ if(!_FECalculation){
+ MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbCells;
+ }
+ else{
+ MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbNodes;
+ }
+ _d_nnz = (nbVoisinsMax+1)*_nVar;
+ _o_nnz = nbVoisinsMax *_nVar;
_meshSet=true;
}
throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
}
// set preconditioner
- if (pcType == NONE)
+ if (pcType == NOPC)
_pctype = (char*)&PCNONE;
else if (pcType ==LU)
_pctype = (char*)&PCLU;
else if (pcType == ICC)
_pctype = (char*)&PCICC;
else {
- cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
- *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
+ cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
+ *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
_runLogFile->close();
- throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
+ throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
}
}
bool stop=false; // Does the Problem want to stop (error) ?
bool converged=false; // has the newton scheme converged (end) ?
- cout<< "!!! Running test case "<< _fileName << " using ";
+ PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName);
*_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
if(!_FECalculation)
{
- cout<< "Finite volumes method"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
*_runLogFile<< "Finite volumes method"<<endl<<endl;
}
else
{
- cout<< "Finite elements method"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
*_runLogFile<< "Finite elements method"<< endl<<endl;
}
computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
if (stop){
- cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing diffusion matrix, stopping calculation\n");
*_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
_runLogFile->close();
throw CdmathException("Failed computing diffusion matrix");
}
computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
if (stop){
- cout << "Error : failed computing right hand side, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n");
*_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
throw CdmathException("Failed computing right hand side");
}
stop = iterateNewtonStep(converged);
if (stop){
- cout << "Error : failed solving linear system, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n");
*_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
_runLogFile->close();
throw CdmathException("Failed solving linear system");
}
void StationaryDiffusionEquation::save(){
- cout<< "Saving numerical results"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
*_runLogFile<< "Saving numerical results"<< endl<<endl;
string resultFile(_path+"/StationaryDiffusionEquation");//Results
string suppress ="rm -rf "+resultFile+"_*";
system(suppress.c_str());//Nettoyage des précédents calculs identiques
- if(_verbose or _system)
- VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
-
- double Ti;
- if(!_FECalculation)
- for(int i=0; i<_Nmailles; i++)
- {
- VecGetValues(_Tk, 1, &i, &Ti);
- _VV(i)=Ti;
- }
- else
- {
- int globalIndex;
- for(int i=0; i<_NunknownNodes; i++)
- {
- VecGetValues(_Tk, 1, &i, &Ti);
- globalIndex = globalNodeIndex(i, _dirichletNodeIds);
- _VV(globalIndex)=Ti;
- }
-
- Node Ni;
- string nameOfGroup;
- for(int i=0; i<_NdirichletNodes; i++)
- {
- Ni=_mesh.getNode(_dirichletNodeIds[i]);
- nameOfGroup = Ni.getGroupName();
- _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
- }
- }
+ if(_mpi_size>1){
+ VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
+ VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
+ }
- _VV.setInfoOnComponent(0,"Temperature_(K)");
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile);
- break;
- case MED :
- _VV.writeMED(resultFile);
- break;
- case CSV :
- _VV.writeCSV(resultFile);
- break;
- }
-}
-Field
-StationaryDiffusionEquation::getOutputTemperatureField()
-{
- if(!_computationCompletedSuccessfully)
- throw("Computation not performed yet or failed. No temperature field available");
- else
- return _VV;
+ if(_verbose or _system)
+ VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD);
+
+ if(_mpi_rank==0){
+ double Ti;
+ if(!_FECalculation)
+ for(int i=0; i<_Nmailles; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tk_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tk , 1, &i, &Ti);
+ _VV(i)=Ti;
+ }
+ else
+ {
+ int globalIndex;
+ for(int i=0; i<_NunknownNodes; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tk_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tk , 1, &i, &Ti);
+ globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
+ _VV(globalIndex)=Ti;
+ }
+
+ Node Ni;
+ string nameOfGroup;
+ for(int i=0; i<_NdirichletNodes; i++)
+ {
+ Ni=_mesh.getNode(_dirichletNodeIds[i]);
+ nameOfGroup = Ni.getGroupName();
+ _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
+ }
+ }
+
+ _VV.setInfoOnComponent(0,"Temperature_(K)");
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile);
+ break;
+ case MED :
+ _VV.writeMED(resultFile);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile);
+ break;
+ }
+ }
}
void StationaryDiffusionEquation::terminate()
VecDestroy(&_deltaT);
VecDestroy(&_b);
MatDestroy(&_A);
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecDestroy(&_Tk_seq);
+
+ PetscBool petscInitialized;
+ PetscInitialized(&petscInitialized);
+ if(petscInitialized)
+ PetscFinalize();
+
+ delete _runLogFile;
}
void
StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
{
if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
{
- j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ j_int=DiffusionEquation::unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
}
}
MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
Field my_eigenfield;
- if(_FECalculation)
- my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
- else
- my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
+ my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
my_eigenfield.setFieldByDataArrayDouble(d);
return my_eigenfield;
}
+
+Field&
+StationaryDiffusionEquation::getOutputTemperatureField()
+{
+ if(!_computationCompletedSuccessfully)
+ throw("Computation not performed yet or failed. No temperature field available");
+ else
+ return _VV;
+}
+
+Field&
+StationaryDiffusionEquation::getRodTemperatureField()
+{
+ return getOutputTemperatureField();
+}
+
+vector<string>
+StationaryDiffusionEquation::getInputFieldsNames()
+{
+ vector<string> result(2);
+
+ result[0]="FluidTemperature";
+ result[1]="HeatPower";
+
+ return result;
+}
+vector<string>
+StationaryDiffusionEquation::getOutputFieldsNames()
+{
+ vector<string> result(1);
+
+ result[0]="RodTemperature";
+
+ return result;
+}
+
+Field&
+StationaryDiffusionEquation::getOutputField(const string& nameField )
+{
+ if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
+ return getRodTemperatureField();
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
+ throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
+ }
+}
+
+void
+StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField )
+{
+ if(!_meshSet)
+ throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setInputField set the mesh first");
+
+ if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
+ return setFluidTemperatureField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name");
+ }
+}
+
+void
+StationaryDiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
+ if(!_meshSet)
+ throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setFluidTemperatureField set initial field first");
+
+ coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
+ _fluidTemperatureField=coupledTemperatureField;
+ _fluidTemperatureFieldSet=true;
+};
+
+void
+StationaryDiffusionEquation::setHeatPowerField(Field heatPower){
+ if(!_meshSet)
+ throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
+
+ heatPower.getMesh().checkFastEquivalWith(_mesh);
+ _heatPowerField=heatPower;
+ _heatPowerFieldSet=true;
+}
+
+void
+StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){
+ if(!_meshSet)
+ throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
+
+ _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel);
+ _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
+ _heatPowerFieldSet=true;
+}
+
+void
+StationaryDiffusionEquation::setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error StationaryDiffusionEquation::setDirichletBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
+
+void
+StationaryDiffusionEquation::setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error StationaryDiffusionEquation::setNeumannBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
