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