From fc396aa59f8a7e038312e65c12916e9ba2d9e7f0 Mon Sep 17 00:00:00 2001 From: michael Date: Sat, 19 Sep 2020 11:05:00 +0200 Subject: [PATCH] Updated environment files --- CDMATH/env_CDMATH.sh | 2 +- CoreFlows/Models/inc/CMakeLists.txt | 1 + .../Models/inc/LinearElasticityModel.hxx | 51 ++++++---- .../inc/StationaryDiffusionEquation.hxx | 8 +- CoreFlows/Models/src/CMakeLists.txt | 1 + .../Models/src/LinearElasticityModel.cxx | 93 ++++++++++++++----- CoreFlows/env_CoreFlows.sh | 8 +- README.md | 3 +- 8 files changed, 118 insertions(+), 49 deletions(-) diff --git a/CDMATH/env_CDMATH.sh b/CDMATH/env_CDMATH.sh index b358940..2b2a1fb 100755 --- a/CDMATH/env_CDMATH.sh +++ b/CDMATH/env_CDMATH.sh @@ -7,7 +7,7 @@ export PETSC_INCLUDES=@PETSC_INCLUDES_INSTALL@ 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@ diff --git a/CoreFlows/Models/inc/CMakeLists.txt b/CoreFlows/Models/inc/CMakeLists.txt index 71cd0e9..a165746 100755 --- a/CoreFlows/Models/inc/CMakeLists.txt +++ b/CoreFlows/Models/inc/CMakeLists.txt @@ -7,6 +7,7 @@ SET(include_models_HXX TransportEquation.hxx DiffusionEquation.hxx StationaryDiffusionEquation.hxx + LinearElasticityModel.hxx Fluide.h ProblemFluid.hxx utilitaire_algebre.h diff --git a/CoreFlows/Models/inc/LinearElasticityModel.hxx b/CoreFlows/Models/inc/LinearElasticityModel.hxx index 4a431b7..741a306 100755 --- a/CoreFlows/Models/inc/LinearElasticityModel.hxx +++ b/CoreFlows/Models/inc/LinearElasticityModel.hxx @@ -7,17 +7,18 @@ * \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_ @@ -25,6 +26,22 @@ 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 { @@ -32,7 +49,7 @@ public : /** \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 @@ -44,7 +61,7 @@ public : 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){ @@ -65,33 +82,33 @@ public : void save(); /* Boundary conditions */ - void setBoundaryFields(map boundaryFields){ + void setBoundaryFields(map 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(_Ndim,0),vector(_Ndim,0), - vector(_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(0),vector(0), - vector(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 : @@ -126,7 +143,7 @@ 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 _limitField; + map _limitField; bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector Vector _normale; @@ -169,7 +186,9 @@ protected : /********* 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_ */ diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx index 7b8ced5..6fd4220 100755 --- a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx +++ b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx @@ -6,7 +6,7 @@ * \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. * */ //============================================================================ @@ -15,18 +15,18 @@ * \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 #include using namespace std; -//! enumeration BoundaryType /*! Boundary condition type */ enum BoundaryTypeStationaryDiffusion { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion}; @@ -50,6 +50,7 @@ public : /** \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 * */ @@ -96,6 +97,7 @@ public : * \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){ diff --git a/CoreFlows/Models/src/CMakeLists.txt b/CoreFlows/Models/src/CMakeLists.txt index b6084d4..f07e759 100755 --- a/CoreFlows/Models/src/CMakeLists.txt +++ b/CoreFlows/Models/src/CMakeLists.txt @@ -14,6 +14,7 @@ SET(src_models_CXX TransportEquation.cxx DiffusionEquation.cxx StationaryDiffusionEquation.cxx +# LinearElasticityModel.cxx Fluide.cxx ProblemFluid.cxx utilitaire_algebre.cxx diff --git a/CoreFlows/Models/src/LinearElasticityModel.cxx b/CoreFlows/Models/src/LinearElasticityModel.cxx index 6bb90d3..a81eebf 100755 --- a/CoreFlows/Models/src/LinearElasticityModel.cxx +++ b/CoreFlows/Models/src/LinearElasticityModel.cxx @@ -89,6 +89,7 @@ LinearElasticityModel::LinearElasticityModel(int dim, bool FECalculation, doubl _NdirichletNodes=0; _NunknownNodes=0; _dirichletValuesSet=false; + _neumannValuesSet=false; //Linear solver data _precision=1.e-6; @@ -176,21 +177,21 @@ void LinearElasticityModel::initialize() _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]<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]<close(); throw CdmathException("Wrong boundary condition"); } @@ -207,7 +208,7 @@ void LinearElasticityModel::initialize() 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 @@ -219,9 +220,24 @@ void LinearElasticityModel::initialize() 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::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::iterator it = _limitField.begin(); - while(it != _limitField.end() and (it->second).bcType == Neumann) + map::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 @@ -269,14 +285,14 @@ Vector LinearElasticityModel::gradientNodal(Matrix M, vector< double > values){ 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); @@ -284,10 +300,10 @@ double LinearElasticityModel::computeDiffusionMatrix(bool & stop) 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); @@ -345,7 +361,7 @@ double LinearElasticityModel::computeDiffusionMatrixFE(bool & stop){ 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); @@ -359,6 +375,29 @@ double LinearElasticityModel::computeDiffusionMatrixFE(bool & stop){ } } + //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); @@ -369,7 +408,7 @@ double LinearElasticityModel::computeDiffusionMatrixFE(bool & stop){ return INFINITY; } -double LinearElasticityModel::computeDiffusionMatrixFV(bool & stop){ +double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){ long nbFaces = _mesh.getNumberOfFaces(); Face Fj; Cell Cell1,Cell2; @@ -425,18 +464,18 @@ double LinearElasticityModel::computeDiffusionMatrixFV(bool & stop){ { 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 !!!!!!!!!!"< dirichletBoundaryVal _dirichletValuesSet=true; _dirichletBoundaryValues=dirichletBoundaryValues; } +void +LinearElasticityModel::setNeumannValues(map< int, double> neumannBoundaryValues) +{ + _neumannValuesSet=true; + _neumannBoundaryValues=neumannBoundaryValues; +} double LinearElasticityModel::getConditionNumber(bool isSingular, double tol) const diff --git a/CoreFlows/env_CoreFlows.sh b/CoreFlows/env_CoreFlows.sh index 19cb5d5..312107c 100755 --- a/CoreFlows/env_CoreFlows.sh +++ b/CoreFlows/env_CoreFlows.sh @@ -1,7 +1,7 @@ #!/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@ @@ -10,6 +10,6 @@ export PETSC_INCLUDES=@PETSC_INCLUDES_INSTALL@ #------------------------------------------------------------------------------------------------------------------- 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 diff --git a/README.md b/README.md index 3adbb6d..b37e485 100755 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ The main research objectives of SOLVERLAB are the study of - 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 --------------- @@ -72,6 +72,7 @@ Simpler build for a minimum version: > 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 -- 2.39.2