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 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; }
66 void setMesh(const Mesh &M);
67 void setFileName(string fileName){ _fileName = fileName; }
68 bool solveStationaryProblem();
69 Field getOutputDisplacementField();
71 //Linear system and spectrum
72 void setLinearSolver(linearSolver kspType, preconditioner pcType);
73 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
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
82 /* Boundary conditions */
83 void setBoundaryFields(map<string, LimitFieldLinearElasticity> boundaryFields){
84 _limitField = boundaryFields;
86 /** \fn setDirichletBoundaryCondition
87 * \brief adds a new boundary condition of type Dirichlet
89 * \param [in] string : the name of the boundary
90 * \param [in] double : the value of the displacement at the boundary
93 void setDirichletBoundaryCondition(string groupName,double displacement){
94 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,displacement,-1);
97 /** \fn setNeumannBoundaryCondition
98 * \brief adds a new boundary condition of type Neumann
100 * \param [in] string : the name of the boundary
101 * \param [in] double : outward normal force
104 void setNeumannBoundaryCondition(string groupName, double normalForce=0){
105 _limitField[groupName]=LimitFieldLinearElasticity(DirichletLinearElasticity,-1, normalForce);
108 void setDirichletValues(map< int, double> dirichletBoundaryValues);
109 void setNeumannValues (map< int, double> neumannBoundaryValues);
116 int _Ndim;//space dimension
117 int _nVar;//Number of equations to solve=1
122 bool _initializedMemory;
123 int _Nmailles;//number of cells for FV calculation
124 int _neibMaxNbCells;//maximum number of cells around a cell
127 double _precision_Newton;
128 double _erreur_rel;//norme(Uk+1-Uk)
129 bool _computationCompletedSuccessfully;
131 //Linear solver and petsc
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
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
148 Vec _displacements;//unknown of the linear system
150 //Physical parameterss
151 double _lambda, _mu;//Lamé coefficients
152 double _rho;//constantDensity
153 Field _densityField;//For non constant density field
154 bool _densityFieldSet;
158 bool _verbose, _system;
159 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
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
165 double computeRHS(bool & stop);
166 double computeStiffnessMatrixFV(bool & stop);
168 /************ Data for FE calculation *************/
170 int _Nnodes;/* number of nodes for FE calculation */
171 int _neibMaxNbNodes;/* maximum number of nodes around a node */
172 int _NunknownNodes;/* number of unknown nodes for FE calculation */
173 int _NboundaryNodes;/* total number of boundary nodes */
174 int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
175 std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
176 std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
178 /*********** Functions for finite element method ***********/
179 Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
180 double computeStiffnessMatrixFE(bool & stop);
182 int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
183 int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
185 /********* Possibility to set a boundary field as Dirichlet boundary condition *********/
186 bool _dirichletValuesSet;
187 bool _neumannValuesSet;
188 std::map< int, double> _dirichletBoundaryValues;
189 std::map< int, double> _neumannBoundaryValues;
192 #endif /* LinearElasticityModel_HXX_ */