1 //============================================================================
3 * \file LinearElasticityModel.hxx
4 * \author Michael NDJINGA
7 * \brief Stationary linear elasticity model
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.
14 //============================================================================
16 /*! \class LinearElasticityModel LinearElasticityModel.hxx "LinearElasticityModel.hxx"
17 * \brief Linear Elasticity Model solved with either finite element or finite volume method.
19 * \sigma = 2 \mu e(u) + \lambda Tr(e(u)) I_d
22 #ifndef LinearElasticityModel_HXX_
23 #define LinearElasticityModel_HXX_
25 #include "ProblemCoreFlows.hxx"
29 /*! Boundary condition type */
30 enum BoundaryTypeLinearElasticity { NeumannLinearElasticity, DirichletLinearElasticity, NoneBCLinearElasticity};
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;
40 BoundaryTypeLinearElasticity bcType;
41 double displacement; //for Dirichlet
42 double normalForce; //for Neumann
45 class LinearElasticityModel
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
58 LinearElasticityModel( int dim, bool FECalculation=true, double rho, double lambda, double mu);
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; }
66 void setMesh(const Mesh &M);
67 void setFileName(string fileName){
70 bool solveStationaryProblem();
71 Field getOutputDisplacementField();
73 //Linear system and spectrum
74 void setLinearSolver(linearSolver kspType, preconditioner pcType);
75 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
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
84 /* Boundary conditions */
85 void setBoundaryFields(map<string, LimitFieldLinearElasticity> boundaryFields){
86 _limitField = boundaryFields;
88 /** \fn setDirichletBoundaryCondition
89 * \brief adds a new boundary condition of type Dirichlet
91 * \param [in] string : the name of the boundary
92 * \param [in] double : the value of the displacement at the boundary
95 void setDirichletBoundaryCondition(string groupName,double displacement){
96 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,displacement,-1);
99 /** \fn setNeumannBoundaryCondition
100 * \brief adds a new boundary condition of type Neumann
102 * \param [in] string : the name of the boundary
103 * \param [in] double : outward normal force
106 void setNeumannBoundaryCondition(string groupName, double normalForce=0){
107 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,-1, normalForce);
110 void setDirichletValues(map< int, double> dirichletBoundaryValues);
111 void setNeumannValues (map< int, double> neumannBoundaryValues);
118 int _Ndim;//space dimension
119 int _nVar;//Number of equations to solve=1
124 bool _initializedMemory;
125 int _Nmailles;//number of cells for FV calculation
126 int _neibMaxNbCells;//maximum number of cells around a cell
129 double _precision_Newton;
130 double _erreur_rel;//norme(Uk+1-Uk)
131 bool _computationCompletedSuccessfully;
133 //Linear solver and petsc
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
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
150 Vec _displacements;//unknown of the linear system
152 //Physical parameterss
153 double _lambda, _mu;//Lamé coefficients
154 double _rho;//constantDensity
155 Field _densityField;//For non constant density field
156 bool _densityFieldSet;
160 bool _verbose, _system;
161 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
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
167 double computeRHS(bool & stop);
168 double computeStiffnessMatrixFV(bool & stop);
170 /************ Data for FE calculation *************/
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 */
180 /*********** Functions for finite element method ***********/
181 Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
182 double computeStiffnessMatrixFE(bool & stop);
184 int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
185 int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
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;
194 #endif /* LinearElasticityModel_HXX_ */