v0.10.0
PlDataAtGaussPts.hpp
Go to the documentation of this file.
1 
2 #ifndef __PLASTICITYCOMMONDATA_HPP__
3 #define __PLASTICITYCOMMONDATA_HPP__
4 
5 
6 struct CommonData {
7 
8 
9  int cOunter = 0;
11  int maxDiv;
12  double adaptPowCoef;
13  double stepRed;
14  double hardMod; ///< Isotropic Hardening modulus \f$ H^{i} \f$
15  double kinMod; ///< Kinematic Hardening modulus \f$ H^{k} \f$
16  double kAppa; ///< isostatic elasticity modulus \f$ \kappa =
17  ///< \frac{E\nu}{(1-2\nu)(1 + \nu)} + \frac{2}{3}\mu \f$
18  double yOung; ///< Young's modulus \f$ E \f$
19  double pOisson; ///< Poisson's ratio \f$ \nu \f$
20  double Mu; ///< shear modulus \f$ \mu = \frac{E}{2(1+\nu)} \f$
21  double cOefficient;
22  double
23  yieldStress; ///< Prescribed value of the yield stress \f$ \sigma_{Y} \f$
24  double yieldCondition; ///< Yield criterion for Von Mises plasticity
25  double plasticMultiplier; ///< Plastic multiplier \f$ \dot{\gamma} \f$
26  double tOlerance;
27  double
28  lAmbda; ///< lame parameter \f$ \lambda = \frac{E \nu}{(1+\nu)(1-2\nu)}
29  ///< \f$
30 
32  tD; ///< Elastic stiffness tensor (4th rank
33  ///< tensor with minor and major
34  ///< symmetry) \f$ \mathbb{D} =
35  ///< \mu(\mathbb{I} + \bar{\mathbb{I}}) -
36  ///< \lambda(\mathbf{I} \otimes
37  ///< \mathbf{I}) \f$
38 
39  boost::shared_ptr<MatrixDouble> gradDispPtr; ///< Total strain tensor
40  ///< \f$
41  ///< \boldsymbol{\varepsilon}
42  ///< \f$
43 
44  MatrixDouble consistentMatrix; ///< \f$ \mathbb{C} = \mathbb{D}
45  ///< - \frac{\mathbb{D} :
46  ///< \mathbf{n} \otimes
47  ///< \mathbb{D} :
48  ///< \mathbf{n}}{\mathbf{n} :
49  ///< \mathbb{D} : \mathbf{n} +
50  ///< \mathbf{K}^{2}} \f$
51 
52  boost::shared_ptr<MatrixDouble> stressPtr; ///< Stress tensor \f$
53  ///< \boldsymbol{\sigma} \f$
54 
55  boost::shared_ptr<MatrixDouble>
56  plasticStrainPtr; ///< Plastic Strain \f$
57  ///< \boldsymbol{\varepsilon}^{p} =
58  ///< \boldsymbol{\varepsilon} -
59  ///< \boldsymbol{\varepsilon}^{e} \f$
60 
61  boost::shared_ptr<VectorDouble> hardParamPtr; ///< Vector containing the
62  ///< isotrtopic hardening
63  ///< parameter at each Gauss
64  ///< point
65 
66  boost::shared_ptr<VectorDouble> kinParamPtr; ///< Vector containing the
67  ///< kinematic hardening
68  ///< parameter at each Gauss
69  ///< point
70 
71  boost::shared_ptr<VectorDouble>
72  kinParamPrevPtr; ///< container that keeps the
73  ///< previous value of the
74  ///< kinetic hardening
75  ///< parameter
76 
77  boost::shared_ptr<MatrixDouble>
78  backStressPtr; ///< Back stress tensor \f$ \boldsymbol{\beta} \f$
79 
80  PetscBool calcReactions;
81 
82  PetscInt testNb;
83 
84  PetscBool isAtomTest;
85 
86  CommonData(MoFEM::Interface &m_field) : mField(m_field) {
87  // default base parameters
88 
89  kinMod = 0.1;
90  hardMod = 0.1;
91  plasticMultiplier = 0.0;
92  yieldStress = 1000;
93  pOisson = 0.1;
94  yOung = 50000;
95  tOlerance = 1e-9;
96  maxDiv = 100;
97  stepRed = 0.2;
98  adaptPowCoef = 1. / 4.;
99 
100  gradDispPtr = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
101  stressPtr = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
102  plasticStrainPtr = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
103  hardParamPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
104  kinParamPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
105  kinParamPrevPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
106  backStressPtr = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
107 
108  ierr = createTags(thPlasticStrain, "PLASTIC_STRAIN", 9);
109  ierr = createTags(thHardParam, "HARDENING_PARAMETER", 1);
110  ierr = createTags(thBackStress, "BACK_STRESS", 9);
111  ierr = createTags(thKinParam, "KIN_HARDENING_PARAMETER", 1);
112 
113  CHKERRABORT(PETSC_COMM_WORLD, ierr);
114  }
115 
117  const int nb_gauss_pts,
118  boost::shared_ptr<MatrixDouble> tag_ptr,
119  Tag tag_ref) {
121  double *tag_data;
122  int tag_size;
123  CHKERR mField.get_moab().tag_get_by_ptr(tag_ref, &fe_ent, 1,
124  (const void **)&tag_data,
125  &tag_size);
126 
127  if (tag_size == 1) {
128  tag_ptr->resize(9, nb_gauss_pts , false);
129  tag_ptr->clear();
130  void const *tag_data[] = {&*tag_ptr->data().begin()};
131  const int tag_size = tag_ptr->data().size();
132  CHKERR mField.get_moab().tag_set_by_ptr(tag_ref, &fe_ent, 1,
133  tag_data, &tag_size);
134  } else if (tag_size != nb_gauss_pts * 9) {
135  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
136  "Wrong size of the tag, data inconsistency");
137  } else {
138  MatrixAdaptor tag_vec = MatrixAdaptor(
139  9, nb_gauss_pts,
140  ublas::shallow_array_adaptor<double>(tag_size, tag_data));
141 
142  *tag_ptr = tag_vec;
143  }
144 
145 
147  }
148 
150  const int nb_gauss_pts,
151  boost::shared_ptr<VectorDouble> tag_ptr,
152  Tag tag_ref) {
154  double *tag_data;
155  int tag_size;
156  CHKERR mField.get_moab().tag_get_by_ptr(
157  tag_ref, &fe_ent, 1, (const void **)&tag_data, &tag_size);
158 
159  if (tag_size == 1) {
160  tag_ptr->resize(nb_gauss_pts);
161  tag_ptr->clear();
162  void const *tag_data[] = {&*tag_ptr->begin()};
163  const int tag_size = tag_ptr->size();
164  CHKERR mField.get_moab().tag_set_by_ptr(tag_ref, &fe_ent, 1,
165  tag_data, &tag_size);
166  } else if (tag_size != nb_gauss_pts) {
167  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
168  "Wrong size of the tag, data inconsistency");
169  } else {
170  VectorAdaptor tag_vec = VectorAdaptor(
171  tag_size, ublas::shallow_array_adaptor<double>(tag_size, tag_data));
172  *tag_ptr = tag_vec;
173  }
175  }
176  template<class T>
177  MoFEMErrorCode setTags(const EntityHandle fe_ent, const int nb_gauss_pts,
178  T tag_ptr, Tag tag_ref){
180  void const *tag_data[] = {&*tag_ptr->data().begin()};
181  const int tag_size = tag_ptr->data().size();
182  CHKERR mField.get_moab().tag_set_by_ptr(tag_ref, &fe_ent, 1,
183  tag_data, &tag_size);
185  }
186 
187  /**
188  * @brief Get the Parameters object
189  *
190  * @return MoFEMErrorCode
191  *
192  * create the function to pass some command lines to the main file.
193  */
196 
197  reactionsSidesetId = 0;
198  calcReactions = PETSC_TRUE;
199  testNb = 1;
200  isAtomTest = PETSC_FALSE;
201 
202  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
203  "none");
204 
205  CHKERR PetscOptionsScalar("-young", "Young modulus", "", yOung, &yOung,
206  PETSC_NULL);
207 
208  CHKERR PetscOptionsScalar("-poisson", "Poisson ratio", "", pOisson,
209  &pOisson, PETSC_NULL);
210 
211  CHKERR PetscOptionsScalar("-tolerance", "tolerance", "", tOlerance,
212  &tOlerance, PETSC_NULL);
213 
214  CHKERR PetscOptionsScalar("-yield", "yield stress", "", yieldStress,
215  &yieldStress, PETSC_NULL);
216 
217  CHKERR PetscOptionsScalar("-HardMod", "Hardening Modulus", "", hardMod,
218  &hardMod, PETSC_NULL);
219 
220  CHKERR PetscOptionsScalar("-KinMod", "Kinematic Hardening Modulus", "",
221  kinMod, &kinMod, PETSC_NULL);
222 
223  CHKERR PetscOptionsScalar("-step_red", "step reductcion while diverge", "",
224  stepRed, &stepRed, PETSC_NULL);
225 
226  CHKERR PetscOptionsScalar("-adapt_coef",
227  "code will be smarter depending on the value", "",
228  adaptPowCoef, &adaptPowCoef, PETSC_NULL);
229 
230  CHKERR PetscOptionsInt(
231  "-calc_reactions", "calculate reactions for blocksets", "",
233 
234  CHKERR PetscOptionsInt("-is_atom_test", "Give number of the atom test", "",
235  testNb, &testNb, &isAtomTest);
236 
237  CHKERR PetscOptionsInt("-steps_max_div", "number of steps", "", maxDiv,
238  &maxDiv, PETSC_NULL);
239 
240  cOefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
241  Mu = yOung / (2 * (1 + pOisson));
242  kAppa = cOefficient * pOisson + (2. / 3.) * Mu;
248  tD(i, j, k, l) = 0.;
249 
250  tD(0, 0, 0, 0) = lAmbda + 2 * Mu;
251  tD(1, 1, 1, 1) = lAmbda + 2 * Mu;
252  tD(2, 2, 2, 2) = lAmbda + 2 * Mu;
253  tD(0, 0, 1, 1) = lAmbda;
254  tD(0, 0, 2, 2) = lAmbda;
255  tD(1, 1, 2, 2) = lAmbda;
256  tD(2, 2, 1, 1) = lAmbda;
257  tD(1, 1, 2, 2) = lAmbda;
258  tD(2, 2, 0, 0) = lAmbda;
259  tD(1, 1, 0, 0) = lAmbda;
260  tD(0, 1, 0, 1) = Mu;
261  tD(0, 2, 0, 2) = Mu;
262  tD(1, 2, 1, 2) = Mu;
263  tD(0, 1, 1, 0) = Mu;
264  tD(0, 2, 2, 0) = Mu;
265  tD(1, 2, 2, 1) = Mu;
266  tD(2, 0, 0, 2) = Mu;
267  tD(1, 0, 0, 1) = Mu;
268  tD(2, 1, 1, 2) = Mu;
269  tD(2, 0, 2, 0) = Mu;
270  tD(1, 0, 1, 0) = Mu;
271  tD(2, 1, 2, 1) = Mu;
272 
273 
274  ierr = PetscOptionsEnd();
275  CHKERRQ(ierr);
276 
278  }
279 
284 
285  // MoFEMErrorCode
286  // calcReactionForces(VectorDouble &reactions, DM &dm,
287  // boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhs,
288  // int &reactionsSidesetId, MoFEM::Interface &mField);
289 
290 private:
292 
294  const string tag_string, const int &dim){
295 
297  std::vector<double> def_val(dim, 0);
298  CHKERR mField.get_moab().tag_get_handle(
299  tag_string.c_str(), 1, MB_TYPE_DOUBLE, tag_ref,
300  MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, &def_val[0]);
302  }
303 
304  public:
306  calcReactionForces(VectorDouble &reactions, DM &dm,
307  boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhs,
308  int &reactionsSidesetId);
309 
310  MoFEMErrorCode atomTests(const int step_number,
311  const VectorDouble &reactions);
312 
313 };
314 
315 #endif //__PLASTICITYCOMMONDATA_HPP__
PetscBool calcReactions
MoFEMErrorCode atomTests(const int step_number, const VectorDouble &reactions)
boost::shared_ptr< MatrixDouble > plasticStrainPtr
Deprecated interface functions.
boost::shared_ptr< MatrixDouble > gradDispPtr
virtual moab::Interface & get_moab()=0
double yieldCondition
Yield criterion for Von Mises plasticity.
PetscBool isAtomTest
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
double Mu
shear modulus
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MatrixDouble consistentMatrix
MoFEMErrorCode calcReactionForces(VectorDouble &reactions, DM &dm, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, int &reactionsSidesetId)
static Index< 'l', 3 > l
boost::shared_ptr< VectorDouble > hardParamPtr
MoFEMErrorCode getTags(const EntityHandle fe_ent, const int nb_gauss_pts, boost::shared_ptr< MatrixDouble > tag_ptr, Tag tag_ref)
const int dim
MoFEM::Interface & mField
double yOung
Young's modulus .
boost::shared_ptr< VectorDouble > kinParamPtr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
double yieldStress
Prescribed value of the yield stress .
static Index< 'i', 3 > i
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:126
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
static Index< 'j', 3 > j
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:109
double plasticMultiplier
Plastic multiplier .
#define CHKERR
Inline error check.
Definition: definitions.h:604
double pOisson
Poisson's ratio .
boost::shared_ptr< VectorDouble > kinParamPrevPtr
double hardMod
Isotropic Hardening modulus .
static Index< 'k', 3 > k
CommonData(MoFEM::Interface &m_field)
MoFEMErrorCode getTags(const EntityHandle fe_ent, const int nb_gauss_pts, boost::shared_ptr< VectorDouble > tag_ptr, Tag tag_ref)
double kinMod
Kinematic Hardening modulus .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MoFEMErrorCode createTags(Tag &tag_ref, const string tag_string, const int &dim)
MoFEMErrorCode setTags(const EntityHandle fe_ent, const int nb_gauss_pts, T tag_ptr, Tag tag_ref)
boost::shared_ptr< MatrixDouble > stressPtr
FTensor::Tensor4< double, 3, 3, 3, 3 > tD
boost::shared_ptr< MatrixDouble > backStressPtr
Back stress tensor .
MoFEMErrorCode getParameters()
Get the Parameters object.