export PETSC_LIBRARIES=@PETSC_DIR@/@PETSC_ARCH@/lib
export MEDFILE_ROOT_DIR=@MEDFILE_ROOT_DIR@
export MEDFILE_INCLUDE_DIRS=@MEDFILE_INCLUDE_DIRS@
-export MEDFILE_LIBRARIES=@MEDFILE_LIBRARIES@
+export MEDFILE_LIBRARIES=@MEDFILE_LIBRARIES_INSTALL@
export MEDCOUPLING_ROOT_DIR=@MEDCOUPLING_ROOT_DIR@
export MEDCOUPLING_INCLUDE_DIR=@MEDCOUPLING_INCLUDE_DIR@
export MEDCOUPLING_LIBRARIES=@MEDCOUPLING_LIBRARIES@
TransportEquation.hxx
DiffusionEquation.hxx
StationaryDiffusionEquation.hxx
+ LinearElasticityModel.hxx
Fluide.h
ProblemFluid.hxx
utilitaire_algebre.h
* \brief Stationary linear elasticity model
* -div \sigma = f
* with the stress \sigma given by the Hooke's law
- * \sigma=2\mu e(u)+\lambda Tr(e(u)) I_d
- * solved with either finite elements or finite volume method
- * Dirichlet (fixed boundary) or Neumann (free boundary) boundary conditions
+ * \sigma = 2 \mu e(u) + \lambda Tr(e(u)) I_d
+ * solved with either finite element or finite volume method
+ * Dirichlet (fixed boundary) or Neumann (free boundary) boundary conditions.
* */
//============================================================================
/*! \class LinearElasticityModel LinearElasticityModel.hxx "LinearElasticityModel.hxx"
- * \brief Linear Elasticity Model solved with either finite elements or finite volume method.
+ * \brief Linear Elasticity Model solved with either finite element or finite volume method.
* -div \sigma = f
- * \sigma=2\mu e(u)+\lambda Tr(e(u)) I_d
+ * \sigma = 2 \mu e(u) + \lambda Tr(e(u)) I_d
*/
+
#ifndef LinearElasticityModel_HXX_
#define LinearElasticityModel_HXX_
using namespace std;
+/*! Boundary condition type */
+enum BoundaryTypeLinearElasticity { NeumannLinearElasticity, DirichletLinearElasticity, NoneBCLinearElasticity};
+
+/** \struct LimitField
+ * \brief value of some fields on the boundary */
+struct LimitFieldLinearElasticity{
+ LimitFieldLinearElasticity(){bcType=NoneBCLinearElasticity; displacement=0; normalForce=0;}
+ LimitFieldLinearElasticity(BoundaryTypeLinearElasticity _bcType, double _displacement, double _normalForce){
+ bcType=_bcType; displacement=_displacement; normalForce=_normalForce;
+ }
+
+ BoundaryTypeLinearElasticity bcType;
+ double displacement; //for Dirichlet
+ double normalForce; //for Neumann
+};
+
class LinearElasticityModel
{
/** \fn LinearElasticityModel
* \brief Constructor for the linear elasticity in a solid
* \param [in] int : space dimension
- * \param [in] double : numerical method
+ * \param [in] bool : numerical method
* \param [in] double : solid density
* \param [in] double : first Lamé coefficient
* \param [in] double : second Lamé coefficient
void setDensityField(Field densityField) { _densityField=densityField; _densityFieldSet=true;}
void setLameCoefficient(double lambda, double mu) { _lambda = lambda; _mu = mu;}
void setYoungAndPoissonModuli(double E, double nu) { _lambda = E*nu/(1+nu)/(1-2*nu); _mu = E/2/(1+nu);}
- void setGravity(Vector gravite ) { _gravite=gravite; }
+ void setGravity(Vector gravite ) { _gravity=gravite; }
void setMesh(const Mesh &M);
void setFileName(string fileName){
void save();
/* Boundary conditions */
- void setBoundaryFields(map<string, LimitField> boundaryFields){
+ void setBoundaryFields(map<string, LimitFieldLinearElasticity> boundaryFields){
_limitField = boundaryFields;
};
/** \fn setDirichletBoundaryCondition
* \brief adds a new boundary condition of type Dirichlet
* \details
* \param [in] string : the name of the boundary
- * \param [in] double : the value of the temperature at the boundary
+ * \param [in] double : the value of the displacement at the boundary
* \param [out] void
* */
- void setDirichletBoundaryCondition(string groupName,double Temperature){
- _limitField[groupName]=LimitField(Dirichlet,-1, vector<double>(_Ndim,0),vector<double>(_Ndim,0),
- vector<double>(_Ndim,0),Temperature,-1,-1,-1);
+ void setDirichletBoundaryCondition(string groupName,double displacement){
+ _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,displacement,-1);
};
/** \fn setNeumannBoundaryCondition
* \brief adds a new boundary condition of type Neumann
* \details
* \param [in] string : the name of the boundary
+ * \param [in] double : outward normal force
* \param [out] void
* */
- void setNeumannBoundaryCondition(string groupName){
- _limitField[groupName]=LimitField(Neumann,-1, vector<double>(0),vector<double>(0),
- vector<double>(0),-1,-1,-1,-1);
+ void setNeumannBoundaryCondition(string groupName, double normalForce=0){
+ _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,-1, normalForce);
};
void setDirichletValues(map< int, double> dirichletBoundaryValues);
+ void setNeumannValues (map< int, double> neumannBoundaryValues);
protected :
double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
bool _conditionNumber;//computes an estimate of the condition number
- map<string, LimitField> _limitField;
+ map<string, LimitFieldLinearElasticity> _limitField;
bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
Vector _normale;
/********* Possibility to set a boundary field as Dirichlet boundary condition *********/
bool _dirichletValuesSet;
+ bool _neumannValuesSet;
std::map< int, double> _dirichletBoundaryValues;
+ std::map< int, double> _neumannBoundaryValues;
};
#endif /* LinearElasticityModel_HXX_ */
* \date June 2019
* \brief Stationary heat diffusion equation solved with either finite elements or finite volume method.
* -\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T)
- * Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions
+ * Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions.
* */
//============================================================================
* \details see \ref StationaryDiffusionEqPage for more details
* -\lambda\Delta T=\Phi(T) + \lambda_{sf} (T_{fluid}-T)
*/
+
#ifndef StationaryDiffusionEquation_HXX_
#define StationaryDiffusionEquation_HXX_
#include "ProblemCoreFlows.hxx"
-/* for the laplacian spectrum */
+/* For the laplacian spectrum */
#include <slepceps.h>
#include <slepcsvd.h>
using namespace std;
-//! enumeration BoundaryType
/*! Boundary condition type */
enum BoundaryTypeStationaryDiffusion { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
/** \fn StationaryDiffusionEquation
* \brief Constructor for the temperature diffusion in a solid
* \param [in] int : space dimension
+ * \param [in] bool : numerical method
* \param [in] double : solid conductivity
* */
* \brief adds a new boundary condition of type Neumann
* \details
* \param [in] string : the name of the boundary
+ * \param [in] double : outward normal flux
* \param [out] void
* */
void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
TransportEquation.cxx
DiffusionEquation.cxx
StationaryDiffusionEquation.cxx
+# LinearElasticityModel.cxx
Fluide.cxx
ProblemFluid.cxx
utilitaire_algebre.cxx
_NdirichletNodes=0;
_NunknownNodes=0;
_dirichletValuesSet=false;
+ _neumannValuesSet=false;
//Linear solver data
_precision=1.e-6;
_runLogFile->close();
throw CdmathException("Missing boundary value");
}
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoTypeSpecified)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCLinearElasticity)
{
cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
*_runLogFile<< "!!!No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
_runLogFile->close();
throw CdmathException("Missing boundary condition");
}
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==Dirichlet)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletLinearElasticity)
_dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=Neumann)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannLinearElasticity)
{
cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann << endl;
+ cout<<"!!! Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity << endl;
*_runLogFile<< "Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann <<endl;
+ *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity <<endl;
_runLogFile->close();
throw CdmathException("Wrong boundary condition");
}
else
MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes*_nVar, _NunknownNodes*_nVar, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
- VecCreate(PETSC_COMM_SELF, &_displacments);
+ VecCreate(PETSC_COMM_SELF, &_displacements);
VecDuplicate(_displacements, &_b);//RHS of the linear system
KSPGetPC(_ksp, &_pc);
PCSetType(_pc, _pctype);
+ //Checking whether all boundary conditions are Neumann boundary condition
+ //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
+ if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
+ {
+ map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
+ while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
+ it++;
+ _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
+ }
+ else
+ if(_FECalculation)
+ _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
+ else
+ _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
+
//Checking whether all boundaries are Neumann boundaries
- map<string, LimitField>::iterator it = _limitField.begin();
- while(it != _limitField.end() and (it->second).bcType == Neumann)
+ map<string, LimitFieldLinearElasticity>::iterator it = _limitField.begin();
+ while(it != _limitField.end() and (it->second).bcType == NeumannLinearElasticity)
it++;
_onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);
//If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
return result;
}
-double LinearElasticityModel::computeDiffusionMatrix(bool & stop)
+double LinearElasticityModel::computeStiffnessMatrix(bool & stop)
{
double result;
if(_FECalculation)
- result=computeDiffusionMatrixFE(stop);
+ result=computeStiffnessMatrixFE(stop);
else
- result=computeDiffusionMatrixFV(stop);
+ result=computeStiffnessMatrixFV(stop);
if(_verbose or _system)
MatView(_A,PETSC_VIEWER_STDOUT_SELF);
return result;
}
-double LinearElasticityModel::computeDiffusionMatrixFE(bool & stop){
+double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){
Cell Cj;
string nameOfGroup;
- double dn;
+ double dn, coeff;
MatZeroEntries(_A);
VecZeroEntries(_b);
if( _dirichletValuesSet )//New way of storing BC
valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
else //old way of storing BC
- valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].Displacements;
+ valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].displacement;
}
else
valuesBorder[kdim]=Vector(_Ndim);
}
}
+ //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()].normalForce;
+ VecSetValue(_b, j_int, coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+
MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
VecAssemblyBegin(_b);
return INFINITY;
}
-double LinearElasticityModel::computeDiffusionMatrixFV(bool & stop){
+double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
long nbFaces = _mesh.getNumberOfFaces();
Face Fj;
Cell Cell1,Cell2;
{
nameOfGroup = Fj.getGroupName();
- if (_limitField[nameOfGroup].bcType==Neumann){//Nothing to do
+ if (_limitField[nameOfGroup].bcType==NeumannLinearElasticity){//Nothing to do
}
- else if(_limitField[nameOfGroup].bcType==Dirichlet){
+ else if(_limitField[nameOfGroup].bcType==DirichletLinearElasticity){
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);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].displacement, ADD_VALUES);
}
else {
stop=true ;
- cout<<"!!!!!!!!!!!!!!! Error LinearElasticityModel::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
+ cout<<"!!!!!!!!!!!!!!! Error LinearElasticityModel::computeStiffnessMatrixFV !!!!!!!!!!"<<endl;
cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- cout<<"Accepted boundary conditions are Neumann "<<Neumann<< " and Dirichlet "<<Dirichlet<<endl;
+ cout<<"Accepted boundary conditions are Neumann "<<NeumannLinearElasticity<< " and Dirichlet "<<DirichletLinearElasticity<<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");
}
else
{
- *_runLogFile<<"LinearElasticityModel::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
- throw CdmathException("LinearElasticityModel::computeDiffusionMatrixFV(): incompatible number of cells around a face");
+ *_runLogFile<<"LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face"<<endl;
+ throw CdmathException("LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face");
}
}
if(_conditionNumber)
KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
- KSPSolve(_ksp, _b, displacements);
+ KSPSolve(_ksp, _b, _displacements);
KSPConvergedReason reason;
KSPGetConvergedReason(_ksp,&reason);
*_runLogFile<< "Finite elements method"<< endl<<endl;
}
- computeDiffusionMatrix( stop);
+ computeStiffnessMatrix( stop);
if (stop){
cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
*_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
_dirichletValuesSet=true;
_dirichletBoundaryValues=dirichletBoundaryValues;
}
+void
+LinearElasticityModel::setNeumannValues(map< int, double> neumannBoundaryValues)
+{
+ _neumannValuesSet=true;
+ _neumannBoundaryValues=neumannBoundaryValues;
+}
double
LinearElasticityModel::getConditionNumber(bool isSingular, double tol) const
#!/bin/bash
-export CDMATH_DIR=@CDMATH_INSTALL@
-source ${CMAKE_INSTALL_PREFIX}/env_SOLVERLAB.sh
+export CDMATH_INSTALL=@CMAKE_INSTALL_PREFIX@
+source @CMAKE_INSTALL_PREFIX@/env_SOLVERLAB.sh
export CoreFlows_INSTALL=@CMAKE_INSTALL_PREFIX@
export PETSC_DIR=@PETSC_DIR@
#-------------------------------------------------------------------------------------------------------------------
export CoreFlows=$CoreFlows_INSTALL/bin/CoreFlowsMainExe
-export LD_LIBRARY_PATH=$CoreFlows_INSTALL/lib:$CDMATH_DIR/lib:${PETSC_DIR}/${PETSC_ARCH}/lib:${MEDCOUPLING_LIBRARIES}:${MEDFILE_C_LIBRARIES}:${LD_LIBRARY_PATH}
-export PYTHONPATH=$CoreFlows_INSTALL/lib:$CoreFlows_INSTALL/lib/coreflows:$CoreFlows_INSTALL/bin/coreflows:$CoreFlows_INSTALL/lib/python2.7/site-packages/salome:$CDMATH_DIR/lib/cdmath:$CDMATH_DIR/bin/cdmath:$CDMATH_DIR/bin/cdmath/postprocessing:${PETSC_DIR}/${PETSC_ARCH}/lib:${MEDCOUPLING_LIBRARIES}:${MEDFILE_C_LIBRARIES}:${PYTHONPATH}
+export LD_LIBRARY_PATH=$CoreFlows_INSTALL/lib:$CDMATH_INSTALL/lib:${PETSC_DIR}/${PETSC_ARCH}/lib:${MEDCOUPLING_LIBRARIES}:${MEDFILE_C_LIBRARIES}:${LD_LIBRARY_PATH}
+export PYTHONPATH=$CoreFlows_INSTALL/lib:$CoreFlows_INSTALL/lib/coreflows:$CoreFlows_INSTALL/bin/coreflows:$CoreFlows_INSTALL/lib/python2.7/site-packages/salome:$CDMATH_INSTALL/lib/cdmath:$CDMATH_INSTALL/bin/cdmath:$CDMATH_INSTALL/bin/cdmath/postprocessing:${PETSC_DIR}/${PETSC_ARCH}/lib:${MEDCOUPLING_LIBRARIES}:${MEDFILE_C_LIBRARIES}:${PYTHONPATH}
export CoreFlowsGUI=$CoreFlows_INSTALL/bin/salome/CoreFlows_Standalone.py
- New preconditioners for implicit methods for two phase flows
- The coupling of fluid models or multiphysics coupling (eg thermal hydraulics and neutronics or thermal hydraulics and solid thermics)
-The library is currently developed for linux distributions and is maintained on Ubuntu 16.04 LTS and 18.04 LTS, as well as on Fedora 24, 26, 28, 29 and 30.
+The library is currently developed for linux distributions and is maintained on Ubuntu 16.04 LTS and 18.04 LTS, as well as on Fedora 24, 26, 28, 30 and 32.
Examples of use
---------------
> This will download and build the following dependencies
> - PETSc from http://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.13.5.tar.gz
> - SLEPc from https://slepc.upv.es/download/distrib/slepc-3.13.4.tar.gz
+> - F2CBLASLAPACK from http://ftp.mcs.anl.gov/pub/petsc/externalpackages/f2cblaslapack-3.4.2.q4.tar.gz
> - HDF5 https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.3/src/hdf5-1.10.3.tar.gz
> - MEDFILE from http://files.salome-platform.org/Salome/other/med-4.1.0.tar.gz
> - MEDCOUPLING from http://files.salome-platform.org/Salome/other/medCoupling-9.4.0.tar.gz
