v0.14.0
ADOLCPlasticityMaterialModels.hpp
Go to the documentation of this file.
1 /** \file ADOLCPlasticityMaterialModels.hpp
2  * \ingroup adoc_plasticity
3  * \brief Matetial models for plasticity
4  *
5  * \example ADOLCPlasticityMaterialModels.hpp
6  *
7  */
8 
9 #ifndef __ADOLC_MATERIAL_MODELS_HPP
10 #define __ADOLC_MATERIAL_MODELS_HPP
11 
12 #ifndef WITH_ADOL_C
13 #error "MoFEM need to be compiled with ADOL-C"
14 #endif
15 
16 namespace ADOLCPlasticity {
17 
18 /** \brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
19  * \ingroup adoc_plasticity
20  */
21 
22 /**
23  * @brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
24  *
25  * @tparam DIM model dimension (2 - is plane stress)
26  */
27 template <int DIM> struct J2Plasticity;
28 
29 /**
30  * @brief J2 (Von Misses) plasticity
31  * @ingroup adoc_plasticity
32  *
33  * Free energy:
34  *
35  * \f[
36  * \psi(\varepsilon^e, \alpha, \beta) =
37  * \frac{1}{2} \lambda \textrm{tr}[\varepsilon^e]^2
38  * + \mu \varepsilon^e : \varepsilon^e
39  * +
40  * \sigma_y\alpha
41  * +
42  * \frac{1}{2} \phi H \alpha^2
43  * +
44  * \frac{1}{2} (1-\phi)K \beta^2
45  * \f]
46  *
47  * Isotropic hardening variable \f$\alpha\f$ is in first index of
48  * internal variables vector. Kinematic hardening variable is in
49  * the remaining indices of internal variables vector.
50  *
51  * Yield function:
52  *
53  * Deviator of actual stress:
54  * \f[
55  * \eta=\textrm{dev}[\sigma]-\overline{\beta}
56  * \f]
57  * where \f$\overline{\beta}\f$ is back stress.
58  *
59  * Square of the deviator norm
60  * \f[
61  * f = \frac{3}{2} \eta:\eta
62  * \f]
63  *
64  * Yield function:
65  * \f[
66  * y = \sqrt{f} - \overline{\alpha}
67  * \f]
68  * where \f$f\f$ is defined in \ref eva
69  *
70  * This is associated potential model, so flow potential is equal to yield
71  *
72  */
73 template <> struct J2Plasticity<3> : public ClosestPointProjection {
74 
76  nbInternalVariables = 7;
77  activeVariablesW.resize(12 + nbInternalVariables, false);
78  }
79 
83 
91 
93 
94  double mu; //< Lamé parameter
95  double lambda; //< Lamé parameters
96  double H; //< Isotropic hardening
97  double K; //< Kinematic hardening
98  double phi; //< Combined isotropic/kinematic hardening
99  double sigmaY; //< Yield stress
100 
101  enum ModelParams { LAMBDA, MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM };
102 
103  using ParamsArray = std::array<double, LAST_PARAM>;
104  boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
105  boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
106  boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
107 
110 
111  double young_modulus = 1;
112  double poisson_ratio = 0.25;
113  sigmaY = 1;
114  H = 0.1;
115  K = 0;
116  phi = 1;
117 
118  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-young_modulus", &young_modulus,
119  PETSC_NULL);
120  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-poisson_ratio", &poisson_ratio,
121  0);
122  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_y", &sigmaY, PETSC_NULL);
123  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-H", &H, PETSC_NULL);
124  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-K", &K, PETSC_NULL);
125  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-phi", &phi, PETSC_NULL);
126 
127  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
128  << "Young modulus: " << young_modulus;
129  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
130  << "Poisson ratio: " << poisson_ratio;
131 
134 
135  defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
136  (*defaultParamsArrayPtr)[LAMBDA] = lambda;
137  (*defaultParamsArrayPtr)[MU] = mu;
138  (*defaultParamsArrayPtr)[H] = H;
139  (*defaultParamsArrayPtr)[K] = K;
140  (*defaultParamsArrayPtr)[PHI] = phi;
141  (*defaultParamsArrayPtr)[SIGMA_Y] = sigmaY;
142 
143  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
144  << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
145  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
146  << "MU: " << (*defaultParamsArrayPtr)[MU];
147  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
148  << "H: " << (*defaultParamsArrayPtr)[ISOH];
149  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
150  << "K: " << (*defaultParamsArrayPtr)[KINK];
151  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
152  << "PHI: " << (*defaultParamsArrayPtr)[PHI];
153  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
154  << "SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
155 
156  paramsArrayPtr = boost::make_shared<ParamsArray>();
157  std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
158  paramsArrayPtr->begin());
159  oldParamsArrayPtr = boost::make_shared<ParamsArray>();
160  std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
161 
163  };
164 
165  /**
166  * @brief Get material parameters from block
167  *
168  * @param m_field
169  * @param pip
170  * @param block_name
171  * @param sev
172  * @return MoFEMErrorCode
173  */
176  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
177  std::string block_name, Sev sev) {
179 
180  struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
181 
183 
184  OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
185  boost::shared_ptr<ParamsArray> def_params_array_otr,
186  MoFEM::Interface &m_field, Sev sev,
187  std::vector<const CubitMeshSets *> meshset_vec_ptr)
188  : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
189  defParamsArray(def_params_array_otr) {
190  CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
191  "Can not get data from block");
192  }
193 
194  MoFEMErrorCode doWork(int side, EntityType type,
197 
198  for (auto &b : blockData) {
199 
200  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
201  std::copy(b.paramsArray.begin(), b.paramsArray.end(),
202  paramsArrayPtr->begin());
204  }
205  }
206 
207  std::copy(defParamsArray->begin(), defParamsArray->end(),
208  paramsArrayPtr->begin());
209 
211  }
212 
213  private:
214  boost::shared_ptr<ParamsArray> paramsArrayPtr;
215  boost::shared_ptr<ParamsArray> defParamsArray;
216 
217  struct BlockData {
218  ParamsArray paramsArray;
219  Range blockEnts;
220  };
221 
222  std::vector<BlockData> blockData;
223 
225  extractBlockData(MoFEM::Interface &m_field,
226  std::vector<const CubitMeshSets *> meshset_vec_ptr,
227  Sev sev) {
229 
230  for (auto m : meshset_vec_ptr) {
231  MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "JP2 parameters") << *m;
232  std::vector<double> block_data;
233  CHKERR m->getAttributes(block_data);
234  if (block_data.size() != 6) {
235  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
236  "Expected that block has two attribute");
237  }
238  auto get_block_ents = [&]() {
239  Range ents;
240  CHKERR
241  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
242  return ents;
243  };
244 
245  blockData.push_back({*defParamsArray, get_block_ents()});
246  std::copy(block_data.begin(), block_data.end(),
247  blockData.back().paramsArray.begin());
248  const auto young_modulus = blockData.back().paramsArray[LAMBDA];
249  const auto poisson_ratio = blockData.back().paramsArray[MU];
250 
251  // Is assumed that uset provide young_modulus and poisson_ratio in
252  // fist two argiumens of the block
253 
254  MOFEM_LOG("ADOLCPlasticityWold", sev)
255  << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
256  MOFEM_LOG("ADOLCPlasticityWold", sev)
257  << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
258 
259  blockData.back().paramsArray[LAMBDA] =
261  blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
262 
263  MOFEM_LOG("ADOLCPlasticityWold", sev)
264  << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
265  MOFEM_LOG("ADOLCPlasticityWold", sev)
266  << "MU: " << (blockData.back().paramsArray)[MU];
267  MOFEM_LOG("ADOLCPlasticityWold", sev)
268  << "H: " << (blockData.back().paramsArray)[ISOH];
269  MOFEM_LOG("ADOLCPlasticityWold", sev)
270  << "K: " << (blockData.back().paramsArray)[KINK];
271  MOFEM_LOG("ADOLCPlasticityWold", sev)
272  << "PHI: " << (blockData.back().paramsArray)[PHI];
273  MOFEM_LOG("ADOLCPlasticityWold", sev)
274  << "SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
275  }
277  }
278  };
279 
280  auto cubit_meshsets_vec_ptr =
281  m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
282  std::regex((boost::format("%s(.*)") % block_name).str()));
283 
284  pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
285  m_field, sev, cubit_meshsets_vec_ptr));
286 
288  }
289 
290  /**
291  * @brief Set material parameters at integration point
292  *
293  * @param tag
294  * @param recalculate_elastic_tangent
295  * @return MoFEMErrorCode
296  */
297  MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
299  switch (tag) {
300  case TypesTags::W: {
301  set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
302  paramsArrayPtr->data());
303  for (auto p = 0; p != LAST_PARAM; p++) {
304  if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
305  recalculate_elastic_tangent = true;
306  std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
307  oldParamsArrayPtr->begin());
308  break;
309  }
310  }
311 
312  } break;
313  case TypesTags::Y:
314  case TypesTags::H:
315  break;
316  default:
317  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
318  "Unknown tag for setParams");
319  }
321  }
322 
323  //! [Free energy J2]
326 
327  // ADOL-C register parameters (those can varied for integration points)
328  auto p_lambda = mkparam(lambda); // 0
329  auto p_mu = mkparam(mu); // 1
330  auto p_H = mkparam(H); // 2
331  auto p_K = mkparam(K); // 3
332  auto p_phi = mkparam(phi); // 4
333  auto p_sigma_y = mkparam(sigmaY); // 5
334 
335  // Operator converting Voight and tensor notation
336  auto t_voight_strain_op = voight_to_strain_op();
337 
338  // Variable in voight notation
339  auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
340  auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
341  auto t_voight_internal_tensor =
342  getFTensor1FromPtr<6>(&a_internalVariables[1]);
343 
344  // Convert variables to tensor notation
345  tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
346  tPlasticStrain(i, j) =
347  t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
348  tInternalTensor(i, j) =
349  t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
350 
351  // Evaluate elastic strain
352  tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
353  tR = tElasticStrain(i, i);
354 
355  // Strain energy
356  a_w = 0.5 * p_lambda * tR * tR +
357  p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
358  // Isotropic strain energy
359  a_w +=
360  a_internalVariables[0] * p_sigma_y +
361  0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
362  // Kinematic strain energy
363  a_w += (0.5 * (1 - p_phi) * p_K) *
364  (tInternalTensor(i, j) * tInternalTensor(i, j));
365 
367  }
368  //! [Free energy J2]
369 
371 
372  //! [Yield function J2]
375  auto t_voight_stress_op = voight_to_stress_op();
376  auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
377  auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
378  tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
379  tBackStress(i, j) = t_voight_stress_op(i, j, Z) * t_back_stress(Z);
380 
381  constexpr double third = boost::math::constants::third<double>();
383  tR = tStress(i, i) / 3.;
384  tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j) - tBackStress(i, j);
385 
386  // Fix the constant to mach uniaxial test
387  constexpr double c = 3. / 2.;
388  f = c * tDeviator(i, j) * tDeviator(i, j);
389 
391  }
392 
395  CHKERR evalF();
396  a_y = sqrt(f) - a_internalFluxes[0];
398  }
399  //! [Yield function J2]
400 
403  CHKERR evalF();
404  a_h = sqrt(f) - a_internalFluxes[0];
406  }
407 };
408 
409 //! [J2 2D]
410 template <> struct J2Plasticity<2> : public J2Plasticity<3> {
412 
415 
416  auto p_lambda = mkparam(lambda); // 0
417  auto p_mu = mkparam(mu); // 1
418  auto p_H = mkparam(H); // 2
419  auto p_K = mkparam(K); // 3
420  auto p_phi = mkparam(phi); // 4
421  auto p_sigma_y = mkparam(sigmaY); // 5
422 
423  auto t_voight_strain_op = voight_to_strain_op();
424  auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
425  auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
426  auto t_voight_internal_tensor =
427  getFTensor1FromPtr<6>(&a_internalVariables[1]);
428 
429  tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
430  tPlasticStrain(i, j) =
431  t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
432  tInternalTensor(i, j) =
433  t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
434 
435  tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
436 
437  //! [Plane stress]
438  // Calculate st rain at ezz enforcing plane strain case
439  tElasticStrain(2, 2) =
440  (-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
441  (p_lambda + 2 * p_mu);
442  //! [Plane stress]
443 
444  tR = tElasticStrain(i, i);
445 
446  a_w = 0.5 * p_lambda * tR * tR +
447  p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
448 
449  // plastic part
450  // isotropic
451  a_w +=
452  a_internalVariables[0] * p_sigma_y +
453  0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
454  // kinematic
455  a_w += (0.5 * (1 - p_phi) * p_K) *
456  (tInternalTensor(i, j) * tInternalTensor(i, j));
457 
459  }
460 };
461 //! [J2 2D]
462 
463 /** \brief Paraboloidal yield criterion
464  *
465  * See paper for reference \cite ullah2017three and \cite stier2015comparing
466  *
467  * Linear hardening
468  * \f[
469  * \psi =
470  * \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
471  * \varepsilon
472  * +
473  * \sigma_{yt}\alpha_0
474  * +
475  * \frac{1}{2} H_t \alpha_0^2
476  * +
477  * \sigma_{yc}\alpha_1
478  * +
479  * \frac{1}{2} H_c \alpha_1^2
480  * \f]
481  *
482  * Exponential hardening
483  * \f[
484  * \psi =
485  * \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
486  * \varepsilon\\
487  * +
488  * (\sigma_{yt}+H_t)\alpha_0
489  * +
490  * \frac{H_t}{n_t} \exp{(-n_t \alpha_0)}\\
491  * +
492  * (\sigma_{yc}+H_c)\alpha_0
493  * +
494  * \frac{H_c}{n_c} \exp{(-n_c \alpha_1)}
495  * \f]
496  *
497  * \f[
498  * I_1 = \textrm{tr} (\boldsymbol{\sigma})
499  * \f]
500  *
501  * \f[
502  * \eta=\textrm{dev}[\boldsymbol{\sigma}]
503  * \f]
504  *
505  * \f[
506  * J_2 = \frac{1}{2} \eta:\eta
507  * \f]
508  *
509  * \f[
510  * y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) -
511  * 2\overline{\alpha_0} \,\overline{\alpha_1} \f]
512  * where
513  * \f[
514  * \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t
515  * \alpha_0
516  * \f]
517  *
518  * \f[
519  * \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c
520  * \alpha_1
521  * \f]
522  *
523  * \f[
524  * \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right)
525  * - 2\overline{\alpha_0} \,\overline{\alpha_1} \f] \f[ \alpha=
526  * \frac{1-2\nu_p}{1+\nu_p}
527  * \f]
528  *
529  * \ingroup adoc_plasticity
530  */
532 
535  activeVariablesW.resize(12 + nbInternalVariables, false);
536  }
537 
541 
547 
549 
550  double lambda; //< Lamé parameters
551  double mu; //< Lamé parameter
552  double nup; //< Plastic Poisson’s ratio
553  double Ht, Hc; //< Isotropic hardening for tension and compression
554  double sigmaYt, sigmaYc; //< Yield stress in tension and compression
555  double nt, nc; //< Hardening parameters for tension and compression
556 
557  enum ModelParams {
559  MU,
563  HT,
564  HC,
565  NT,
566  NC,
568  };
569 
570  using ParamsArray = std::array<double, LAST_PARAM>;
571  boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
572 
575 
576  double young_modulus = 1;
577  double poisson_ratio = 0.25;
578 
579  sigmaYt = 1;
580  sigmaYc = 1;
581  Ht = 0.1;
582  Hc = 0.1;
583  nc = 1.;
584  nt = 1.;
585  nup = 0.;
586 
587  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-young_modulus", &young_modulus,
588  PETSC_NULL);
589  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-poisson_ratio", &poisson_ratio,
590  0);
591  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_yt", &sigmaYt, PETSC_NULL);
592  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_yc", &sigmaYc, PETSC_NULL);
593  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-Ht", &Ht, PETSC_NULL);
594  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-Hc", &Hc, PETSC_NULL);
595  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nt", &nt, PETSC_NULL);
596  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nc", &nc, PETSC_NULL);
597  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nup", &nup, PETSC_NULL);
598 
599  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
600  << "Young modulus: " << young_modulus;
601  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
602  << "Poisson ratio: " << poisson_ratio;
603 
606 
607  defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
608  (*defaultParamsArrayPtr)[LAMBDA] = lambda;
609  (*defaultParamsArrayPtr)[MU] = mu;
610  (*defaultParamsArrayPtr)[NUP] = nup;
611  (*defaultParamsArrayPtr)[SIGMA_YT] = sigmaYt;
612  (*defaultParamsArrayPtr)[SIGMA_YC] = sigmaYc;
613  (*defaultParamsArrayPtr)[HT] = Ht;
614  (*defaultParamsArrayPtr)[HC] = Hc;
615  (*defaultParamsArrayPtr)[NT] = nt;
616  (*defaultParamsArrayPtr)[NC] = nc;
617 
618  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
619  << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
620  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
621  << "MU: " << (*defaultParamsArrayPtr)[MU];
622  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
623  << "NUP: " << (*defaultParamsArrayPtr)[NUP];
624  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
625  << "SIGMA_YT: " << (*defaultParamsArrayPtr)[SIGMA_YT];
626  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
627  << "SIGMA_YC: " << (*defaultParamsArrayPtr)[SIGMA_YC];
628  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
629  << "HT: " << (*defaultParamsArrayPtr)[HT];
630  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
631  << "HC: " << (*defaultParamsArrayPtr)[HC];
632  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
633  << "NT: " << (*defaultParamsArrayPtr)[NT];
634  MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
635  << "NC: " << (*defaultParamsArrayPtr)[NC];
636 
638  };
639 
642  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
643  std::string block_name, Sev sev) {
644  return 0;
645  }
646 
647  MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
648  recalculate_elastic_tangent = true;
649  return 0;
650  }
651 
652  //! [Paraboloidal free energy]
655 
656  auto t_voight_strain_op = voight_to_strain_op();
657  auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
658  auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
659 
660  tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
661  tPlasticStrain(i, j) =
662  t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
664  tR = tElasticStrain(i, i);
665 
666  a_w = 0.5 * lambda * tR * tR +
667  mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
668 
673 
675  }
676  //! [Paraboloidal free energy]
677 
679 
680 //! [Paraboloidal yield function]
683  auto t_voight_stress_op = voight_to_stress_op();
684  auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
685  tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
686  tR = tStress(i, i);
687  I1 = tR;
689  tR /= 3.;
690  tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j);
691  J2 = 0.5 * (tDeviator(i, j) * tDeviator(i, j));
693  }
694 
697  CHKERR evalF();
698  a_y = 6 * J2 + 2 * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
701  }
702 //! [Paraboloidal yield function]
703 
704 //! [Paraboloidal flow potential]
707  CHKERR evalF();
708  double alpha =
709  (1 - 2 * nup) /
710  (1 + nup); // relation between alpha and plastic poisson's ratio
711  a_h = 6 * J2 +
712  2 * alpha * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
715  }
716 //! [Paraboloidal flow potential]
717 
718 };
719 
720 template <typename MODEL>
721 inline auto createMaterial(
722  std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
723  auto cp_ptr = boost::make_shared<MODEL>();
724  cp_ptr->getDefaultMaterialParameters();
725  cp_ptr->createMatAVecR();
726  cp_ptr->snesCreate();
727  cp_ptr->recordTapes();
728  cp_ptr->tapesTags = tapes_tags;
729  return cp_ptr;
730 }
731 
732 } // namespace ADOLCPlasticity
733 
734 #endif //__ADOLC_MATERIAL_MODELS_HPP
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
ADOLCPlasticity::ClosestPointProjection::a_h
adouble a_h
Definition: ADOLCPlasticity.hpp:330
ADOLCPlasticity::ParaboloidalPlasticity::defaultParamsArrayPtr
boost::shared_ptr< ParamsArray > defaultParamsArrayPtr
Definition: ADOLCPlasticityMaterialModels.hpp:571
ADOLCPlasticity::J2Plasticity< 3 >::tStress
FTensor::Tensor2_symmetric< adouble, 3 > tStress
Definition: ADOLCPlasticityMaterialModels.hpp:88
ADOLCPlasticity::J2Plasticity< 3 >::f
adouble f
[Free energy J2]
Definition: ADOLCPlasticityMaterialModels.hpp:370
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
ADOLCPlasticity::ParaboloidalPlasticity::tElasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tElasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:544
ADOLCPlasticity::J2Plasticity< 3 >::lambda
double lambda
Definition: ADOLCPlasticityMaterialModels.hpp:95
ADOLCPlasticity::ParaboloidalPlasticity::SIGMA_YC
@ SIGMA_YC
Definition: ADOLCPlasticityMaterialModels.hpp:562
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:123
ADOLCPlasticity::J2Plasticity
J2 plasticity (Kinematic Isotropic (Linear) Hardening)
Definition: ADOLCPlasticityMaterialModels.hpp:27
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
ADOLCPlasticity::ClosestPointProjection::a_plasticStrain
ublas::vector< adouble > a_plasticStrain
Definition: ADOLCPlasticity.hpp:325
ADOLCPlasticity::ParaboloidalPlasticity::HT
@ HT
Definition: ADOLCPlasticityMaterialModels.hpp:563
ADOLCPlasticity::ParaboloidalPlasticity::nt
double nt
Definition: ADOLCPlasticityMaterialModels.hpp:555
ADOLCPlasticity::ParaboloidalPlasticity::tR
adouble tR
Definition: ADOLCPlasticityMaterialModels.hpp:548
ADOLCPlasticity::ParaboloidalPlasticity::j
FTensor::Index< 'j', 3 > j
Definition: ADOLCPlasticityMaterialModels.hpp:539
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ADOLCPlasticity::ParaboloidalPlasticity::Ht
double Ht
Definition: ADOLCPlasticityMaterialModels.hpp:553
ADOLCPlasticity::ParaboloidalPlasticity::sigmaYt
double sigmaYt
Definition: ADOLCPlasticityMaterialModels.hpp:554
ADOLCPlasticity::ParaboloidalPlasticity::Hc
double Hc
Definition: ADOLCPlasticityMaterialModels.hpp:553
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
ADOLCPlasticity::ClosestPointProjection::a_sTrain
ublas::vector< adouble > a_sTrain
Definition: ADOLCPlasticity.hpp:324
ADOLCPlasticity::J2Plasticity< 3 >::i
FTensor::Index< 'i', 3 > i
Definition: ADOLCPlasticityMaterialModels.hpp:80
ADOLCPlasticity::ParaboloidalPlasticity::NUP
@ NUP
Definition: ADOLCPlasticityMaterialModels.hpp:560
W
double W
Definition: free_surface.cpp:166
FTensor::Tensor2_symmetric< adouble, 3 >
ADOLCPlasticity::J2Plasticity< 3 >::ParamsArray
std::array< double, LAST_PARAM > ParamsArray
Definition: ADOLCPlasticityMaterialModels.hpp:103
ADOLCPlasticity::ParaboloidalPlasticity::flowPotential
MoFEMErrorCode flowPotential()
[Paraboloidal yield function]
Definition: ADOLCPlasticityMaterialModels.hpp:705
ADOLCPlasticity::ParaboloidalPlasticity::mu
double mu
Definition: ADOLCPlasticityMaterialModels.hpp:551
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ADOLCPlasticity::ParaboloidalPlasticity::NT
@ NT
Definition: ADOLCPlasticityMaterialModels.hpp:565
ADOLCPlasticity::ParaboloidalPlasticity
[J2 2D]
Definition: ADOLCPlasticityMaterialModels.hpp:531
ADOLCPlasticity::ParaboloidalPlasticity::ModelParams
ModelParams
Definition: ADOLCPlasticityMaterialModels.hpp:557
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
ADOLCPlasticity::ParaboloidalPlasticity::nc
double nc
Definition: ADOLCPlasticityMaterialModels.hpp:555
ADOLCPlasticity::J2Plasticity< 3 >::tR
adouble tR
Definition: ADOLCPlasticityMaterialModels.hpp:92
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ADOLCPlasticity::ClosestPointProjection
Closest point projection algorithm.
Definition: ADOLCPlasticity.hpp:139
ADOLCPlasticity::J2Plasticity< 3 >::freeHemholtzEnergy
MoFEMErrorCode freeHemholtzEnergy()
[Free energy J2]
Definition: ADOLCPlasticityMaterialModels.hpp:324
ADOLCPlasticity::J2Plasticity< 3 >::J2Plasticity
J2Plasticity()
Definition: ADOLCPlasticityMaterialModels.hpp:75
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
ADOLCPlasticity::ParaboloidalPlasticity::getDefaultMaterialParameters
MoFEMErrorCode getDefaultMaterialParameters()
Definition: ADOLCPlasticityMaterialModels.hpp:573
BlockData
Definition: gel_analysis.cpp:74
ADOLCPlasticity::J2Plasticity< 3 >::Z
FTensor::Index< 'Z', 6 > Z
Definition: ADOLCPlasticityMaterialModels.hpp:82
ADOLCPlasticity::ParaboloidalPlasticity::evalF
MoFEMErrorCode evalF()
[Paraboloidal yield function]
Definition: ADOLCPlasticityMaterialModels.hpp:681
ADOLCPlasticity::ParaboloidalPlasticity::I1
adouble I1
[Paraboloidal free energy]
Definition: ADOLCPlasticityMaterialModels.hpp:678
convert.type
type
Definition: convert.py:64
ADOLCPlasticity::J2Plasticity< 3 >::j
FTensor::Index< 'j', 3 > j
Definition: ADOLCPlasticityMaterialModels.hpp:81
ADOLCPlasticity::ParaboloidalPlasticity::i
FTensor::Index< 'i', 3 > i
Definition: ADOLCPlasticityMaterialModels.hpp:538
ADOLCPlasticity::J2Plasticity< 2 >::freeHemholtzEnergy
MoFEMErrorCode freeHemholtzEnergy()
[Free energy J2]
Definition: ADOLCPlasticityMaterialModels.hpp:413
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
ADOLCPlasticity::J2Plasticity< 3 >::tDeviator
FTensor::Tensor2_symmetric< adouble, 3 > tDeviator
Definition: ADOLCPlasticityMaterialModels.hpp:90
ADOLCPlasticity::createMaterial
auto createMaterial(std::array< int, ClosestPointProjection::LAST_TAPE > tapes_tags={0, 1, 2})
Definition: ADOLCPlasticityMaterialModels.hpp:721
ADOLCPlasticity::voight_to_strain_op
FTensor::Dg< double, 3, 6 > voight_to_strain_op()
Op convert Vight strain vector to strain tensor.
Definition: ADOLCPlasticity.hpp:29
ADOLCPlasticity::J2Plasticity< 3 >::K
double K
Definition: ADOLCPlasticityMaterialModels.hpp:97
ADOLCPlasticity::ParaboloidalPlasticity::Z
FTensor::Index< 'Z', 6 > Z
Definition: ADOLCPlasticityMaterialModels.hpp:540
ADOLCPlasticity::ParaboloidalPlasticity::ParaboloidalPlasticity
ParaboloidalPlasticity()
Definition: ADOLCPlasticityMaterialModels.hpp:533
ADOLCPlasticity::J2Plasticity< 3 >::yieldFunction
MoFEMErrorCode yieldFunction()
Set yield function.
Definition: ADOLCPlasticityMaterialModels.hpp:393
ADOLCPlasticity::ParaboloidalPlasticity::tPlasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tPlasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:543
ADOLCPlasticity::ParaboloidalPlasticity::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get model parameters from blocks.
Definition: ADOLCPlasticityMaterialModels.hpp:641
ADOLCPlasticity::ClosestPointProjection::a_w
adouble a_w
Definition: ADOLCPlasticity.hpp:328
ADOLCPlasticity
Definition: ADOLCPlasticity.hpp:24
ADOLCPlasticity::ParaboloidalPlasticity::nup
double nup
Definition: ADOLCPlasticityMaterialModels.hpp:552
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ADOLCPlasticity::ClosestPointProjection::a_internalFluxes
ublas::vector< adouble > a_internalFluxes
Definition: ADOLCPlasticity.hpp:327
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
ADOLCPlasticity::ClosestPointProjection::nbInternalVariables
int nbInternalVariables
Definition: ADOLCPlasticity.hpp:141
ADOLCPlasticity::ClosestPointProjection::activeVariablesW
VectorAdaptor activeVariablesW
Definition: ADOLCPlasticity.hpp:173
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
ADOLCPlasticity::J2Plasticity< 3 >::tElasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tElasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:86
ADOLCPlasticity::ParaboloidalPlasticity::tDeviator
FTensor::Tensor2_symmetric< adouble, 3 > tDeviator
Definition: ADOLCPlasticityMaterialModels.hpp:546
ADOLCPlasticity::J2Plasticity< 3 >::SIGMA_Y
@ SIGMA_Y
Definition: ADOLCPlasticityMaterialModels.hpp:101
ADOLCPlasticity::ParaboloidalPlasticity::LAST_PARAM
@ LAST_PARAM
Definition: ADOLCPlasticityMaterialModels.hpp:567
Range
ADOLCPlasticity::J2Plasticity< 3 >::mu
double mu
Definition: ADOLCPlasticityMaterialModels.hpp:94
adouble
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
ADOLCPlasticity::J2Plasticity< 3 >::tInternalTensor
FTensor::Tensor2_symmetric< adouble, 3 > tInternalTensor
Definition: ADOLCPlasticityMaterialModels.hpp:87
ADOLCPlasticity::J2Plasticity< 3 >::tBackStress
FTensor::Tensor2_symmetric< adouble, 3 > tBackStress
Definition: ADOLCPlasticityMaterialModels.hpp:89
ADOLCPlasticity::ParaboloidalPlasticity::tTotalStrain
FTensor::Tensor2_symmetric< adouble, 3 > tTotalStrain
Definition: ADOLCPlasticityMaterialModels.hpp:542
ADOLCPlasticity::ClosestPointProjection::a_internalVariables
ublas::vector< adouble > a_internalVariables
Definition: ADOLCPlasticity.hpp:326
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
ADOLCPlasticity::ClosestPointProjection::a_y
adouble a_y
Definition: ADOLCPlasticity.hpp:329
ADOLCPlasticity::J2Plasticity< 3 >::evalF
MoFEMErrorCode evalF()
[Yield function J2]
Definition: ADOLCPlasticityMaterialModels.hpp:373
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
ADOLCPlasticity::ParaboloidalPlasticity::LAMBDA
@ LAMBDA
Definition: ADOLCPlasticityMaterialModels.hpp:558
ADOLCPlasticity::J2Plasticity< 3 >::phi
double phi
Definition: ADOLCPlasticityMaterialModels.hpp:98
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ADOLCPlasticity::voight_to_stress_op
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
Definition: ADOLCPlasticity.hpp:49
ADOLCPlasticity::ClosestPointProjection::a_sTress
ublas::vector< adouble > a_sTress
Definition: ADOLCPlasticity.hpp:327
ADOLCPlasticity::ParaboloidalPlasticity::ParamsArray
std::array< double, LAST_PARAM > ParamsArray
Definition: ADOLCPlasticityMaterialModels.hpp:570
ADOLCPlasticity::ParaboloidalPlasticity::MU
@ MU
Definition: ADOLCPlasticityMaterialModels.hpp:559
ADOLCPlasticity::ParaboloidalPlasticity::sigmaYc
double sigmaYc
Definition: ADOLCPlasticityMaterialModels.hpp:554
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
ADOLCPlasticity::ParaboloidalPlasticity::setParams
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set parameters for ADOL-C of constitutive relations.
Definition: ADOLCPlasticityMaterialModels.hpp:647
ADOLCPlasticity::ParaboloidalPlasticity::yieldFunction
MoFEMErrorCode yieldFunction()
Set yield function.
Definition: ADOLCPlasticityMaterialModels.hpp:695
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
ADOLCPlasticity::ParaboloidalPlasticity::SIGMA_YT
@ SIGMA_YT
Definition: ADOLCPlasticityMaterialModels.hpp:561
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
ADOLCPlasticity::ParaboloidalPlasticity::NC
@ NC
Definition: ADOLCPlasticityMaterialModels.hpp:566
ADOLCPlasticity::J2Plasticity< 3 >::tPlasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tPlasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:85
ADOLCPlasticity::ParaboloidalPlasticity::HC
@ HC
Definition: ADOLCPlasticityMaterialModels.hpp:564
ADOLCPlasticity::J2Plasticity< 3 >::getDefaultMaterialParameters
MoFEMErrorCode getDefaultMaterialParameters()
Definition: ADOLCPlasticityMaterialModels.hpp:108
ADOLCPlasticity::J2Plasticity< 3 >::ModelParams
ModelParams
Definition: ADOLCPlasticityMaterialModels.hpp:101
ADOLCPlasticity::ParaboloidalPlasticity::tStress
FTensor::Tensor2_symmetric< adouble, 3 > tStress
Definition: ADOLCPlasticityMaterialModels.hpp:545
ADOLCPlasticity::J2Plasticity< 3 >::setParams
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set material parameters at integration point.
Definition: ADOLCPlasticityMaterialModels.hpp:297
ADOLCPlasticity::ParaboloidalPlasticity::freeHemholtzEnergy
MoFEMErrorCode freeHemholtzEnergy()
[Paraboloidal free energy]
Definition: ADOLCPlasticityMaterialModels.hpp:653
OP
ADOLCPlasticity::J2Plasticity< 3 >::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get material parameters from block.
Definition: ADOLCPlasticityMaterialModels.hpp:175
ADOLCPlasticity::ParaboloidalPlasticity::lambda
double lambda
Definition: ADOLCPlasticityMaterialModels.hpp:550
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
ADOLCPlasticity::J2Plasticity< 3 >::flowPotential
MoFEMErrorCode flowPotential()
[Yield function J2]
Definition: ADOLCPlasticityMaterialModels.hpp:401
H
double H
Hardening.
Definition: plastic.cpp:124
ADOLCPlasticity::J2Plasticity< 3 >::tTotalStrain
FTensor::Tensor2_symmetric< adouble, 3 > tTotalStrain
Definition: ADOLCPlasticityMaterialModels.hpp:84
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
ADOLCPlasticity::J2Plasticity< 3 >::sigmaY
double sigmaY
Definition: ADOLCPlasticityMaterialModels.hpp:99
ADOLCPlasticity::J2Plasticity< 3 >::H
double H
Definition: ADOLCPlasticityMaterialModels.hpp:96
MU
#define MU(E, NU)
Definition: fem_tools.h:23
ADOLCPlasticity::ParaboloidalPlasticity::J2
adouble J2
Definition: ADOLCPlasticityMaterialModels.hpp:678