CoreFlows/Models/inc/LinearElasticityModel.hxx

   1 //============================================================================
   2 /**
   3  * \file LinearElasticityModel.hxx
   4  * \author Michael NDJINGA
   5  * \version 1.0
   6  * \date August 2020
   7  * \brief Stationary linear elasticity model
   8  * -div \sigma = f
   9  * with the stress \sigma given by the Hooke's law
  10  * \sigma=2\mu e(u)+\lambda Tr(e(u)) I_d
  11  * solved with either finite elements or finite volume method
  12  * Dirichlet (fixed boundary) or Neumann (free boundary) boundary conditions
  13  * */
  14 //============================================================================
  15
  16 /*! \class LinearElasticityModel LinearElasticityModel.hxx "LinearElasticityModel.hxx"
  17  *  \brief Linear Elasticity Model solved with either finite elements or finite volume method.
  18  * -div \sigma = f
  19  * \sigma=2\mu e(u)+\lambda Tr(e(u)) I_d
  20  */
  21 #ifndef LinearElasticityModel_HXX_
  22 #define LinearElasticityModel_HXX_
  23
  24 #include "ProblemCoreFlows.hxx"
  25
  26 using namespace std;
  27
  28 class LinearElasticityModel
  29 {
  30
  31 public :
  32         /** \fn LinearElasticityModel
  33                          * \brief Constructor for the linear elasticity in a solid
  34                          * \param [in] int : space dimension
  35                          * \param [in] double : numerical method
  36                          * \param [in] double : solid density
  37                          * \param [in] double : first Lamé coefficient
  38                          * \param [in] double : second  Lamé coefficient
  39                          *  */
  40
  41         LinearElasticityModel( int dim, bool FECalculation=true, double rho, double lambda, double mu);
  42
  43         void setConstantDensity(double rho) { _rho=rho; }
  44         void setDensityField(Field densityField) { _densityField=densityField; _densityFieldSet=true;}
  45         void setLameCoefficient(double lambda, double mu) { _lambda = lambda; _mu = mu;}
  46         void setYoungAndPoissonModuli(double E, double nu) { _lambda = E*nu/(1+nu)/(1-2*nu); _mu = E/2/(1+nu);}
  47         void setGravity(Vector gravite ) { _gravite=gravite; }
  48
  49     void setMesh(const Mesh &M);
  50     void setFileName(string fileName){
  51         _fileName = fileName;
  52     }
  53     bool solveStationaryProblem();
  54     Field getOutputDisplacementField();
  55
  56     //Linear system and spectrum
  57     void setLinearSolver(linearSolver kspType, preconditioner pcType);
  58     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
  59
  60         //Gestion du calcul
  61         void initialize();
  62         void terminate();//vide la mémoire et enregistre le résultat final
  63         double computeStiffnessMatrix(bool & stop);
  64         bool solveLinearSystem();//return true if resolution successfull
  65         void save();
  66
  67     /* Boundary conditions */
  68         void setBoundaryFields(map<string, LimitField> boundaryFields){
  69                 _limitField = boundaryFields;
  70     };
  71         /** \fn setDirichletBoundaryCondition
  72                          * \brief adds a new boundary condition of type Dirichlet
  73                          * \details
  74                          * \param [in] string : the name of the boundary
  75                          * \param [in] double : the value of the temperature at the boundary
  76                          * \param [out] void
  77                          *  */
  78         void setDirichletBoundaryCondition(string groupName,double Temperature){
  79                 _limitField[groupName]=LimitField(Dirichlet,-1, vector<double>(_Ndim,0),vector<double>(_Ndim,0),
  80                                                         vector<double>(_Ndim,0),Temperature,-1,-1,-1);
  81         };
  82
  83         /** \fn setNeumannBoundaryCondition
  84                          * \brief adds a new boundary condition of type Neumann
  85                          * \details
  86                          * \param [in] string : the name of the boundary
  87                          * \param [out] void
  88                          *  */
  89         void setNeumannBoundaryCondition(string groupName){
  90                 _limitField[groupName]=LimitField(Neumann,-1, vector<double>(0),vector<double>(0),
  91                                                       vector<double>(0),-1,-1,-1,-1);
  92         };
  93
  94         void setDirichletValues(map< int, double> dirichletBoundaryValues);
  95
  96
  97 protected :
  98         //Main unknown field
  99         Field _VV;
 100
 101         int _Ndim;//space dimension
 102         int _nVar;//Number of equations to solve=1
 103
 104     //Mesh data
 105         Mesh _mesh;
 106     bool _meshSet;
 107         bool _initializedMemory;
 108         int _Nmailles;//number of cells for FV calculation
 109         int _neibMaxNbCells;//maximum number of cells around a cell
 110
 111         double _precision;
 112         double _precision_Newton;
 113         double _erreur_rel;//norme(Uk+1-Uk)
 114     bool _computationCompletedSuccessfully;
 115
 116         //Linear solver and petsc
 117         KSP _ksp;
 118         KSPType _ksptype;
 119         PC _pc;
 120         PCType _pctype;
 121         string _pc_hypre;
 122         int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
 123         int _PetscIts;//the number of iterations of the linear solver
 124         Mat  _A;//Linear system matrix
 125         Vec _b;//Linear system right hand side
 126         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
 127         bool _conditionNumber;//computes an estimate of the condition number
 128
 129         map<string, LimitField> _limitField;
 130     bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
 131
 132         Vector _normale;
 133     Vec _displacements;//unknown of the linear system
 134
 135         //Physical parameterss
 136         double _lambda, _mu;//Lamé coefficients
 137         double _rho;//constantDensity
 138         Field _densityField;//For non constant density field
 139         bool _densityFieldSet;
 140         Vector _gravity;
 141
 142         //Display variables
 143         bool _verbose, _system;
 144         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 145     //saving parameters
 146         string _fileName;//name of the calculation
 147         string _path;//path to execution directory used for saving results
 148         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
 149
 150         double computeRHS(bool & stop);
 151         double computeStiffnessMatrixFV(bool & stop);
 152
 153     /************ Data for FE calculation *************/
 154     bool _FECalculation;
 155         int _Nnodes;/* number of nodes for FE calculation */
 156         int _neibMaxNbNodes;/* maximum number of nodes around a node */
 157         int _NunknownNodes;/* number of unknown nodes for FE calculation */
 158         int _NboundaryNodes;/* total number of boundary nodes */
 159         int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
 160     std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
 161     std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
 162
 163     /*********** Functions for finite element method ***********/
 164     Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
 165         double computeStiffnessMatrixFE(bool & stop);
 166     int fact(int n);
 167     int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
 168     int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
 169
 170     /********* Possibility to set a boundary field as Dirichlet boundary condition *********/
 171     bool _dirichletValuesSet;
 172     std::map< int, double> _dirichletBoundaryValues;
 173 };
 174
 175 #endif /* LinearElasticityModel_HXX_ */