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,MPI_Comm comm = MPI_COMM_WORLD);
  59
  60         void setConstantDensity(double rho) { _rho=rho; }
  61         void setDensityField(Field densityField) { _densityField=densityField; _densityFieldSet=true;}
  62         void setLameCoefficients(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){  _fileName = fileName;    }
  68     bool solveStationaryProblem();
  69     Field getOutputDisplacementField();
  70
  71     //Linear system and spectrum
  72     void setLinearSolver(linearSolver kspType, preconditioner pcType);
  73     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
  74
  75         //Gestion du calcul
  76         void initialize();
  77         void terminate();//vide la mémoire et enregistre le résultat final
  78         double computeStiffnessMatrix(bool & stop);
  79         bool solveLinearSystem();//return true if resolution successfull
  80         void save();
  81
  82     /* Boundary conditions */
  83         void setBoundaryFields(map<string, LimitFieldLinearElasticity> boundaryFields){
  84                 _limitField = boundaryFields;
  85     };
  86         /** \fn setDirichletBoundaryCondition
  87                          * \brief adds a new boundary condition of type Dirichlet
  88                          * \details
  89                          * \param [in] string : the name of the boundary
  90                          * \param [in] double : the value of the displacement at the boundary
  91                          * \param [out] void
  92                          *  */
  93         void setDirichletBoundaryCondition(string groupName,double displacement){
  94                 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,displacement,-1);
  95         };
  96
  97         /** \fn setNeumannBoundaryCondition
  98                          * \brief adds a new boundary condition of type Neumann
  99                          * \details
 100                          * \param [in] string : the name of the boundary
 101                          * \param [in] double : outward normal force
 102                          * \param [out] void
 103                          *  */
 104         void setNeumannBoundaryCondition(string groupName, double normalForce=0){
 105                 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,-1, normalForce);
 106         };
 107
 108         void setDirichletValues(map< int, double> dirichletBoundaryValues);
 109         void setNeumannValues  (map< int, double> neumannBoundaryValues);
 110
 111
 112 protected :
 113         //Main unknown field
 114         Field _VV;
 115
 116         int _Ndim;//space dimension
 117         int _nVar;//Number of equations to solve=1
 118
 119     //Mesh data
 120         Mesh _mesh;
 121     bool _meshSet;
 122         bool _initializedMemory;
 123         int _Nmailles;//number of cells for FV calculation
 124         int _neibMaxNbCells;//maximum number of cells around a cell
 125
 126         double _precision;
 127         double _precision_Newton;
 128         double _erreur_rel;//norme(Uk+1-Uk)
 129     bool _computationCompletedSuccessfully;
 130
 131         //Linear solver and petsc
 132         KSP _ksp;
 133         KSPType _ksptype;
 134         PC _pc;
 135         PCType _pctype;
 136         string _pc_hypre;
 137         int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
 138         int _PetscIts;//the number of iterations of the linear solver
 139         Mat  _A;//Linear system matrix
 140         Vec _b;//Linear system right hand side
 141         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
 142         bool _conditionNumber;//computes an estimate of the condition number
 143
 144         map<string, LimitFieldLinearElasticity> _limitField;
 145     bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
 146
 147         Vector _normale;
 148     Vec _displacements;//unknown of the linear system
 149
 150         //Physical parameterss
 151         double _lambda, _mu;//Lamé coefficients
 152         double _rho;//constantDensity
 153         Field _densityField;//For non constant density field
 154         bool _densityFieldSet;
 155         Vector _gravity;
 156
 157         //Display variables
 158         bool _verbose, _system;
 159         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 160     //saving parameters
 161         string _fileName;//name of the calculation
 162         string _path;//path to execution directory used for saving results
 163         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
 164
 165         double computeRHS(bool & stop);
 166         double computeStiffnessMatrixFV(bool & stop);
 167         double computeStiffnessMatrixFE(bool & stop);
 168
 169     /************ Data for FE calculation *************/
 170     bool _FECalculation;
 171         int _Nnodes;/* number of nodes for FE calculation */
 172         int _neibMaxNbNodes;/* maximum number of nodes around a node */
 173         int _NunknownNodes;/* number of unknown nodes for FE calculation */
 174         int _NboundaryNodes;/* total number of boundary nodes */
 175         int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
 176     std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
 177     std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
 178
 179     /********* Possibility to set a boundary field as Dirichlet boundary condition *********/
 180     bool _dirichletValuesSet;
 181     bool _neumannValuesSet;
 182     std::map< int, double> _dirichletBoundaryValues;
 183     std::map< int, double> _neumannBoundaryValues;
 184
 185         //MPI variables
 186         PetscMPIInt    _size;        /* size of communicator */
 187         PetscMPIInt    _rank;        /* processor rank */
 188 };
 189
 190 #endif /* LinearElasticityModel_HXX_ */