v0.15.0
Loading...
Searching...
No Matches
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
16namespace 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 */
27template <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 */
73template <> 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_NULLPTR, "-young_modulus", &young_modulus,
119 PETSC_NULLPTR);
120 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio", &poisson_ratio,
121 0);
122 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_y", &sigmaY, PETSC_NULLPTR);
123 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-H", &H, PETSC_NULLPTR);
124 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-K", &K, PETSC_NULLPTR);
125 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-phi", &phi, PETSC_NULLPTR);
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)[ISOH] = H;
139 (*defaultParamsArrayPtr)[KINK] = 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 */
174 MoFEMErrorCode
176 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
177 std::string block_name, Sev sev) {
179
180 struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
181
182 using OP = ForcesAndSourcesCore::UserDataOperator;
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,
195 EntitiesFieldData::EntData &data) {
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
224 MoFEMErrorCode
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 // It is assumed that user provide young_modulus and poisson_ratio in
252 // first two arguments 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]
324 MoFEMErrorCode freeHemholtzEnergy() {
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]
373 MoFEMErrorCode evalF() {
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
393 MoFEMErrorCode yieldFunction() {
395 CHKERR evalF();
396 a_y = sqrt(f) - a_internalFluxes[0];
398 }
399 //! [Yield function J2]
400
401 MoFEMErrorCode flowPotential() {
403 CHKERR evalF();
404 a_h = sqrt(f) - a_internalFluxes[0];
406 }
407
408 MoFEMErrorCode codedHessianW(vector<double *> hessian)
409 {
411 lambda = (*paramsArrayPtr)[0];
412 mu = (*paramsArrayPtr)[1];
413 H = (*paramsArrayPtr)[5];
414
415 double lambda_plus_2mu = lambda + 2 * mu;
416 // first 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon^p^2}
417 // lower symmetric part
418 // Diagonal elements
419 hessian[0][0] = lambda_plus_2mu;
420 hessian[1][1] = lambda_plus_2mu;
421 hessian[2][2] = lambda_plus_2mu;
422
423 // Off-diagonal elements
424 hessian[1][0] = lambda;
425 hessian[2][0] = lambda;
426 hessian[2][1] = lambda;
427
428 // Shear components
429 for (int i = 3; i <= 5; ++i) {
430 hessian[i][i] = mu;
431 }
432
433 //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon^p}
434 double neg_lambda_plus_2mu = -(lambda + 2 * mu);
435 double neg_lambda = -lambda;
436 double neg_mu = -mu;
437
438 for (int i = 6; i <= 8; ++i) {
439 for (int j = 0; j <= 2; ++j) {
440 hessian[i][j] = (i - 6 == j) ? neg_lambda_plus_2mu : neg_lambda;
441 }
442 }
443
444 hessian[9][3] = neg_mu;
445 hessian[10][4] = neg_mu;
446 hessian[11][5] = neg_mu;
447
448 //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon}
449 hessian[6][6] = lambda_plus_2mu;
450 hessian[7][6] = lambda;
451 hessian[7][7] = lambda_plus_2mu;
452 hessian[8][6] = lambda;
453 hessian[8][7] = lambda;
454 hessian[8][8] = lambda_plus_2mu;
455 hessian[9][9] = mu;
456 hessian[10][10] = mu;
457 hessian[11][11] = mu;
458
459 // internalVariables is \frac{\partial^2 \psi}{\partial \alpha^2}
460 hessian[12][12] = H;
461
463 }
464};
465
466//! [J2 2D]
467template <> struct J2Plasticity<2> : public J2Plasticity<3> {
468 using J2Plasticity<3>::J2Plasticity;
469
470 MoFEMErrorCode freeHemholtzEnergy() {
472
473 auto p_lambda = mkparam(lambda); // 0
474 auto p_mu = mkparam(mu); // 1
475 auto p_H = mkparam(H); // 2
476 auto p_K = mkparam(K); // 3
477 auto p_phi = mkparam(phi); // 4
478 auto p_sigma_y = mkparam(sigmaY); // 5
479
480 auto t_voight_strain_op = voight_to_strain_op();
481 auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
482 auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
483 auto t_voight_internal_tensor =
484 getFTensor1FromPtr<6>(&a_internalVariables[1]);
485
486 tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
487 tPlasticStrain(i, j) =
488 t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
489 tInternalTensor(i, j) =
490 t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
491
492 tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
493
494 //! [Plane stress]
495 // Calculate st rain at ezz enforcing plane strain case
496 tElasticStrain(2, 2) =
497 (-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
498 (p_lambda + 2 * p_mu);
499 //! [Plane stress]
500
501 tR = tElasticStrain(i, i);
502
503 a_w = 0.5 * p_lambda * tR * tR +
504 p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
505
506 // plastic part
507 // isotropic
508 a_w +=
509 a_internalVariables[0] * p_sigma_y +
510 0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
511 // kinematic
512 a_w += (0.5 * (1 - p_phi) * p_K) *
513 (tInternalTensor(i, j) * tInternalTensor(i, j));
514
516 }
517
518 MoFEMErrorCode codedHessianW(vector<double *> hessian) {
520 // Not implemented yet
521 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
522 "2D Energy Hessian not yet implemented");
524 }
525};
526//! [J2 2D]
527
528/** \brief Paraboloidal yield criterion
529 *
530 * See paper for reference \cite ullah2017three and \cite stier2015comparing
531 *
532 * Linear hardening
533 * \f[
534 * \psi =
535 * \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
536 * \varepsilon
537 * +
538 * \sigma_{yt}\alpha_0
539 * +
540 * \frac{1}{2} H_t \alpha_0^2
541 * +
542 * \sigma_{yc}\alpha_1
543 * +
544 * \frac{1}{2} H_c \alpha_1^2
545 * \f]
546 *
547 * Exponential hardening
548 * \f[
549 * \psi =
550 * \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
551 * \varepsilon\\
552 * +
553 * (\sigma_{yt}+H_t)\alpha_0
554 * +
555 * \frac{H_t}{n_t} \exp{(-n_t \alpha_0)}\\
556 * +
557 * (\sigma_{yc}+H_c)\alpha_0
558 * +
559 * \frac{H_c}{n_c} \exp{(-n_c \alpha_1)}
560 * \f]
561 *
562 * \f[
563 * I_1 = \textrm{tr} (\boldsymbol{\sigma})
564 * \f]
565 *
566 * \f[
567 * \eta=\textrm{dev}[\boldsymbol{\sigma}]
568 * \f]
569 *
570 * \f[
571 * J_2 = \frac{1}{2} \eta:\eta
572 * \f]
573 *
574 * \f[
575 * y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) -
576 * 2\overline{\alpha_0} \,\overline{\alpha_1} \f]
577 * where
578 * \f[
579 * \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t
580 * \alpha_0
581 * \f]
582 *
583 * \f[
584 * \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c
585 * \alpha_1
586 * \f]
587 *
588 * \f[
589 * \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right)
590 * - 2\overline{\alpha_0} \,\overline{\alpha_1} \f] \f[ \alpha=
591 * \frac{1-2\nu_p}{1+\nu_p}
592 * \f]
593 *
594 * \ingroup adoc_plasticity
595 */
597
602
606
612
614
615 double lambda; //< Lamé parameters
616 double mu; //< Lamé parameter
617 double nup; //< Plastic Poisson’s ratio
618 double Ht, Hc; //< Isotropic hardening for tension and compression
619 double sigmaYt, sigmaYc; //< Yield stress in tension and compression
620 double nt, nc; //< Hardening parameters for tension and compression
621
634
635 using ParamsArray = std::array<double, LAST_PARAM>;
636 boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
637 boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
638 boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
639
642
643 double young_modulus = 1;
644 double poisson_ratio = 0.25;
645
646 sigmaYt = 1;
647 sigmaYc = 1;
648 Ht = 0.1;
649 Hc = 0.1;
650 nc = 1.;
651 nt = 1.;
652 nup = 0.;
653
654 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-young_modulus", &young_modulus,
655 PETSC_NULLPTR);
656 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio", &poisson_ratio,
657 0);
658 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_yt", &sigmaYt, PETSC_NULLPTR);
659 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_yc", &sigmaYc, PETSC_NULLPTR);
660 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-Ht", &Ht, PETSC_NULLPTR);
661 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-Hc", &Hc, PETSC_NULLPTR);
662 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nt", &nt, PETSC_NULLPTR);
663 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nc", &nc, PETSC_NULLPTR);
664 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nup", &nup, PETSC_NULLPTR);
665
666 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
667 << "Young modulus: " << young_modulus;
668 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
669 << "Poisson ratio: " << poisson_ratio;
670
673
674 defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
675 (*defaultParamsArrayPtr)[LAMBDA] = lambda;
676 (*defaultParamsArrayPtr)[MU] = mu;
677 (*defaultParamsArrayPtr)[NUP] = nup;
678 (*defaultParamsArrayPtr)[SIGMA_YT] = sigmaYt;
679 (*defaultParamsArrayPtr)[SIGMA_YC] = sigmaYc;
680 (*defaultParamsArrayPtr)[HT] = Ht;
681 (*defaultParamsArrayPtr)[HC] = Hc;
682 (*defaultParamsArrayPtr)[NT] = nt;
683 (*defaultParamsArrayPtr)[NC] = nc;
684
685 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
686 << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
687 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
688 << "MU: " << (*defaultParamsArrayPtr)[MU];
689 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
690 << "NUP: " << (*defaultParamsArrayPtr)[NUP];
691 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
692 << "SIGMA_YT: " << (*defaultParamsArrayPtr)[SIGMA_YT];
693 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
694 << "SIGMA_YC: " << (*defaultParamsArrayPtr)[SIGMA_YC];
695 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
696 << "HT: " << (*defaultParamsArrayPtr)[HT];
697 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
698 << "HC: " << (*defaultParamsArrayPtr)[HC];
699 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
700 << "NT: " << (*defaultParamsArrayPtr)[NT];
701 MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
702 << "NC: " << (*defaultParamsArrayPtr)[NC];
703
704 paramsArrayPtr = boost::make_shared<ParamsArray>();
705 std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
706 paramsArrayPtr->begin());
707 oldParamsArrayPtr = boost::make_shared<ParamsArray>();
708 std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
709
711 };
712
713 MoFEMErrorCode
715 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
716 std::string block_name, Sev sev) {
718
719 struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
720
721 using OP = ForcesAndSourcesCore::UserDataOperator;
722
723 OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
724 boost::shared_ptr<ParamsArray> def_params_array_ptr,
725 MoFEM::Interface &m_field, Sev sev,
726 std::vector<const CubitMeshSets *> meshset_vec_ptr)
727 : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
728 defParamsArray(def_params_array_ptr) {
729 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
730 "Can not get data from block");
731 }
732
733 MoFEMErrorCode doWork(int side, EntityType type,
734 EntitiesFieldData::EntData &data) {
736
737 for (auto &b : blockData) {
738
739 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
740 std::copy(b.paramsArray.begin(), b.paramsArray.end(),
741 paramsArrayPtr->begin());
743 }
744 }
745
746 std::copy(defParamsArray->begin(), defParamsArray->end(),
747 paramsArrayPtr->begin());
748
750 }
751
752 private:
753 boost::shared_ptr<ParamsArray> paramsArrayPtr;
754 boost::shared_ptr<ParamsArray> defParamsArray;
755
756 struct BlockData {
757 ParamsArray paramsArray;
758 Range blockEnts;
759 };
760
761 std::vector<BlockData> blockData;
762
763 MoFEMErrorCode
764 extractBlockData(MoFEM::Interface &m_field,
765 std::vector<const CubitMeshSets *> meshset_vec_ptr,
766 Sev sev) {
768
769 for (auto m : meshset_vec_ptr) {
770 MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "Paraboloidal") << *m;
771 std::vector<double> block_data;
772 CHKERR m->getAttributes(block_data);
773 if (block_data.size() != 9) {
774 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
775 "Expected that block has nine attributes");
776 }
777 auto get_block_ents = [&]() {
778 Range ents;
779 CHKERR
780 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
781 return ents;
782 };
783
784 blockData.push_back({*defParamsArray, get_block_ents()});
785 std::copy(block_data.begin(), block_data.end(),
786 blockData.back().paramsArray.begin());
787
788 const auto young_modulus = blockData.back().paramsArray[LAMBDA];
789 const auto poisson_ratio = blockData.back().paramsArray[MU];
790
791 // It is assumed that user provide young_modulus and poisson_ratio in
792 // first two arguments of the block
793
794 MOFEM_LOG("ADOLCPlasticityWold", sev)
795 << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
796 MOFEM_LOG("ADOLCPlasticityWold", sev)
797 << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
798
799 blockData.back().paramsArray[LAMBDA] =
801 blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
802
803 MOFEM_LOG("ADOLCPlasticityWold", sev)
804 << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
805 MOFEM_LOG("ADOLCPlasticityWold", sev)
806 << "MU: " << (blockData.back().paramsArray)[MU];
807 MOFEM_LOG("ADOLCPlasticityWold", sev)
808 << "NUP: " << (blockData.back().paramsArray)[NUP];
809 MOFEM_LOG("ADOLCPlasticityWold", sev)
810 << "SIGMA_YT: " << (blockData.back().paramsArray)[SIGMA_YT];
811 MOFEM_LOG("ADOLCPlasticityWold", sev)
812 << "SIGMA_YC: " << (blockData.back().paramsArray)[SIGMA_YC];
813 MOFEM_LOG("ADOLCPlasticityWold", sev)
814 << "HT: " << (blockData.back().paramsArray)[HT];
815 MOFEM_LOG("ADOLCPlasticityWold", sev)
816 << "HC: " << (blockData.back().paramsArray)[HC];
817 MOFEM_LOG("ADOLCPlasticityWold", sev)
818 << "NT: " << (blockData.back().paramsArray)[NT];
819 MOFEM_LOG("ADOLCPlasticityWold", sev)
820 << "NC: " << (blockData.back().paramsArray)[NC];
821 }
823 }
824 };
825
826 auto cubit_meshsets_vec_ptr =
827 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
828 std::regex((boost::format("%s(.*)") % block_name).str()));
829
830 pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
831 m_field, sev, cubit_meshsets_vec_ptr));
832
834 }
835
836 MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
838 switch (tag) {
839 case TypesTags::W: {
840 set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
841 paramsArrayPtr->data());
842 for (auto p = 0; p != LAST_PARAM; p++) {
843 if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
844 recalculate_elastic_tangent = true;
845 std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
846 oldParamsArrayPtr->begin());
847 break;
848 }
849 }
850
851 } break;
852 case TypesTags::Y:
853 case TypesTags::H:
854 break;
855 default:
856 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
857 "Unknown tag for setParams");
858 }
860 return 0;
861 }
862
863 //! [Paraboloidal free energy]
864 MoFEMErrorCode freeHemholtzEnergy() {
866
867 auto p_lambda = mkparam(lambda); // 0
868 auto p_mu = mkparam(mu); // 1
869 auto p_nup = mkparam(nup); // 2
870 auto p_sigma_yt = mkparam(sigmaYt); // 3
871 auto p_sigma_yc = mkparam(sigmaYc); // 4
872 auto p_Ht = mkparam(Ht); // 5
873 auto p_Hc = mkparam(Hc); // 6
874 auto p_nt = mkparam(nt); // 7
875 auto p_nc = mkparam(nc); // 8
876
877 auto t_voight_strain_op = voight_to_strain_op();
878 auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
879 auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
880
881 tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
882 tPlasticStrain(i, j) =
883 t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
885 tR = tElasticStrain(i, i);
886
887 a_w = 0.5 * p_lambda * tR * tR +
888 p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
889
890 a_w += p_sigma_yt * a_internalVariables[0] +
891 0.5 * p_Ht * a_internalVariables[0] * a_internalVariables[0];
892 a_w += p_sigma_yc * a_internalVariables[1] +
893 0.5 * p_Hc * a_internalVariables[1] * a_internalVariables[1];
894
896 }
897 //! [Paraboloidal free energy]
898
899
900 MoFEMErrorCode codedHessianW(vector<double *> hessian)
901 {
903 // first 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon^p^2}
904 // lower symmetric part
905
906 lambda = (*paramsArrayPtr)[0];
907 mu = (*paramsArrayPtr)[1];
908
909 Ht = (*paramsArrayPtr)[5];
910 Hc = (*paramsArrayPtr)[6];
911
912 double lambda_plus_2mu = lambda + 2 * mu;
913
914 // Diagonal elements
915 hessian[0][0] = lambda_plus_2mu;
916 hessian[1][1] = lambda_plus_2mu;
917 hessian[2][2] = lambda_plus_2mu;
918
919 // Off-diagonal elements
920 hessian[1][0] = lambda;
921 hessian[2][0] = lambda;
922 hessian[2][1] = lambda;
923
924 // Shear components
925 for (int i = 3; i <= 5; ++i) {
926 hessian[i][i] = mu;
927 }
928
929 //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon^p}
930 double neg_lambda_plus_2mu = -(lambda + 2 * mu);
931 double neg_lambda = -lambda;
932 double neg_mu = -mu;
933
934 for (int i = 6; i <= 8; ++i) {
935 for (int j = 0; j <= 2; ++j) {
936 hessian[i][j] = (i - 6 == j) ? neg_lambda_plus_2mu : neg_lambda;
937 }
938 }
939
940 hessian[9][3] = neg_mu;
941 hessian[10][4] = neg_mu;
942 hessian[11][5] = neg_mu;
943
944 //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon}
945 hessian[6][6] = lambda_plus_2mu;
946 hessian[7][6] = lambda;
947 hessian[7][7] = lambda_plus_2mu;
948 hessian[8][6] = lambda;
949 hessian[8][7] = lambda;
950 hessian[8][8] = lambda_plus_2mu;
951 hessian[9][9] = mu;
952 hessian[10][10] = mu;
953 hessian[11][11] = mu;
954
955 // internalVariables is \frac{\partial^2 \psi}{\partial \alpha^2}
956 hessian[12][12] = Ht;
957 hessian[13][13] = Hc;
958
960 }
961
963
964//! [Paraboloidal yield function]
965 MoFEMErrorCode evalF() {
967 auto t_voight_stress_op = voight_to_stress_op();
968 auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
969 tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
970 tR = tStress(i, i);
971 I1 = tR;
973 tR /= 3.;
974 tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j);
975 J2 = 0.5 * (tDeviator(i, j) * tDeviator(i, j));
977 }
978
979 MoFEMErrorCode yieldFunction() {
981 CHKERR evalF();
982 a_y = 6 * J2 + 2 * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
985 }
986//! [Paraboloidal yield function]
987
988//! [Paraboloidal flow potential]
989 MoFEMErrorCode flowPotential() {
991 CHKERR evalF();
992 double alpha =
993 (1 - 2 * nup) /
994 (1 + nup); // relation between alpha and plastic poisson's ratio
995 a_h = 6 * J2 +
996 2 * alpha * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
999 }
1000//! [Paraboloidal flow potential]
1001};
1002
1004
1006
1007 MoFEMErrorCode
1009 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1010 std::string block_name, Sev sev) override {
1012
1013 struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
1014
1015 using OP = ForcesAndSourcesCore::UserDataOperator;
1016
1017 double poissonRatio = 0.25;
1018 Tag thYoungModulus;
1019 Tag thCompressiveYieldStress;
1020 Tag thTensionYieldStress;
1021 double youngsModulus;
1022 double yieldStrengthC;
1023 double yieldStrengthT;
1024
1025 OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
1026 boost::shared_ptr<ParamsArray> def_params_array_ptr,
1027 MoFEM::Interface &m_field, Sev sev,
1028 std::vector<const CubitMeshSets *> meshset_vec_ptr)
1029 : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
1030 defParamsArray(def_params_array_ptr), mField(m_field) {
1031 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio",
1032 &poissonRatio, 0);
1033 rval = mField.get_moab().tag_get_handle(
1034 "YOUNGS_MODULUS", 1, MB_TYPE_DOUBLE, thYoungModulus, MB_TAG_SPARSE);
1035 if (rval != MB_SUCCESS) {
1036 MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
1037 << "YOUNGS_MODULUS tag does not exist using default value";
1038 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-young_modulus",
1039 &youngsModulus, 0);
1040 }
1041
1042 rval = mField.get_moab().tag_get_handle(
1043 "Yield_Strength_C", 1, MB_TYPE_DOUBLE, thCompressiveYieldStress,
1044 MB_TAG_SPARSE);
1045
1046 if (rval != MB_SUCCESS) {
1047 MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
1048 << "Yield_Strength_C tag does not exist using default value";
1049 yieldStrengthC = defParamsArray->at(SIGMA_YC);
1050 }
1051
1052 rval = mField.get_moab().tag_get_handle(
1053 "Yield_Strength_T", 1, MB_TYPE_DOUBLE, thTensionYieldStress,
1054 MB_TAG_SPARSE);
1055
1056 if (rval != MB_SUCCESS) {
1057 MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
1058 << "Yield_Strength_T tag does not exist using default value";
1059 yieldStrengthT = defParamsArray->at(SIGMA_YT);
1060 }
1061
1062 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
1063 "Can not get data from block");
1064 }
1065
1066 MoFEMErrorCode doWork(int side, EntityType type,
1067 EntitiesFieldData::EntData &data) {
1069
1070 for (auto &b : blockData) {
1071
1072 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
1073 EntityHandle fe_ent = getFEEntityHandle();
1074
1075 if (thYoungModulus != 0)
1076 CHKERR mField.get_moab().tag_get_data(thYoungModulus, &fe_ent, 1,
1077 &youngsModulus);
1078
1079 if (thCompressiveYieldStress != 0)
1080 CHKERR mField.get_moab().tag_get_data(
1081 thCompressiveYieldStress, &fe_ent, 1, &yieldStrengthC);
1082
1083 if (thTensionYieldStress != 0)
1084 CHKERR mField.get_moab().tag_get_data(
1085 thTensionYieldStress, &fe_ent, 1, &yieldStrengthT);
1086
1087 std::copy(b.paramsArray.begin(), b.paramsArray.end(),
1088 paramsArrayPtr->begin());
1089
1090 (*paramsArrayPtr)[LAMBDA] = LAMBDA(youngsModulus, poissonRatio);
1091 (*paramsArrayPtr)[MU] = MU(youngsModulus, poissonRatio);
1092 (*paramsArrayPtr)[SIGMA_YC] = yieldStrengthC;
1093 (*paramsArrayPtr)[SIGMA_YT] = yieldStrengthT;
1094
1096 }
1097 }
1098
1099 std::copy(defParamsArray->begin(), defParamsArray->end(),
1100 paramsArrayPtr->begin());
1101
1103 }
1104
1105 private:
1106 boost::shared_ptr<ParamsArray> paramsArrayPtr;
1107 boost::shared_ptr<ParamsArray> defParamsArray;
1108 MoFEM::Interface &mField;
1109
1110 struct BlockData {
1111 ParamsArray paramsArray;
1112 Range blockEnts;
1113 };
1114
1115 std::vector<BlockData> blockData;
1116
1117 MoFEMErrorCode
1118 extractBlockData(MoFEM::Interface &m_field,
1119 std::vector<const CubitMeshSets *> meshset_vec_ptr,
1120 Sev sev) {
1122
1123 for (auto m : meshset_vec_ptr) {
1124 MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev,
1125 "Heterogeneous Paraboloidal")
1126 << *m;
1127 std::vector<double> block_data;
1128 CHKERR m->getAttributes(block_data);
1129 if (block_data.size() != 0) {
1130 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1131 "Expected that block has zero attributes");
1132 }
1133 auto get_block_ents = [&]() {
1134 Range ents;
1135 CHKERR
1136 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
1137 return ents;
1138 };
1139
1140 blockData.push_back({*defParamsArray, get_block_ents()});
1141 std::copy(block_data.begin(), block_data.end(),
1142 blockData.back().paramsArray.begin());
1143 const auto young_modulus = blockData.back().paramsArray[LAMBDA];
1144 const auto poisson_ratio = blockData.back().paramsArray[MU];
1145
1146 // It is assumed that user provide young_modulus and poisson_ratio in
1147 // first two arguments of the block
1148
1149 MOFEM_LOG("ADOLCPlasticityWold", sev)
1150 << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
1151 MOFEM_LOG("ADOLCPlasticityWold", sev)
1152 << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
1153
1154 blockData.back().paramsArray[LAMBDA] =
1156 blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
1157
1158 MOFEM_LOG("ADOLCPlasticityWold", sev)
1159 << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
1160 MOFEM_LOG("ADOLCPlasticityWold", sev)
1161 << "MU: " << (blockData.back().paramsArray)[MU];
1162 MOFEM_LOG("ADOLCPlasticityWold", sev)
1163 << "NUP: " << (blockData.back().paramsArray)[NUP];
1164 MOFEM_LOG("ADOLCPlasticityWold", sev)
1165 << "SIGMA_YT: " << (blockData.back().paramsArray)[SIGMA_YT];
1166 MOFEM_LOG("ADOLCPlasticityWold", sev)
1167 << "SIGMA_YC: " << (blockData.back().paramsArray)[SIGMA_YC];
1168 MOFEM_LOG("ADOLCPlasticityWold", sev)
1169 << "HT: " << (blockData.back().paramsArray)[HT];
1170 MOFEM_LOG("ADOLCPlasticityWold", sev)
1171 << "HC: " << (blockData.back().paramsArray)[HC];
1172 MOFEM_LOG("ADOLCPlasticityWold", sev)
1173 << "NT: " << (blockData.back().paramsArray)[NT];
1174 MOFEM_LOG("ADOLCPlasticityWold", sev)
1175 << "NC: " << (blockData.back().paramsArray)[NC];
1176 }
1178 }
1179 };
1180
1181 auto cubit_meshsets_vec_ptr =
1182 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1183 std::regex((boost::format("%s(.*)") % block_name).str()));
1184
1185 pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
1186 m_field, sev, cubit_meshsets_vec_ptr));
1187
1189 }
1190};
1191
1192template <typename MODEL>
1193inline auto createMaterial(
1194 std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
1195 auto cp_ptr = boost::make_shared<MODEL>();
1196 cp_ptr->getDefaultMaterialParameters();
1197 cp_ptr->createMatAVecR();
1198 cp_ptr->snesCreate();
1199 cp_ptr->recordTapes();
1200 cp_ptr->tapesTags = tapes_tags;
1201 return cp_ptr;
1202}
1203
1204} // namespace ADOLCPlasticity
1205
1206#endif //__ADOLC_MATERIAL_MODELS_HPP
constexpr double third
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MU(E, NU)
Definition fem_tools.h:23
#define LAMBDA(E, NU)
Definition fem_tools.h:22
constexpr auto t_kd
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
auto createMaterial(std::array< int, ClosestPointProjection::LAST_TAPE > tapes_tags={0, 1, 2})
FTensor::Dg< double, 3, 6 > voight_to_strain_op()
Op convert Vight strain vector to strain tensor.
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
double H
Hardening.
Definition plastic.cpp:128
double sigmaY
Yield stress.
Definition plastic.cpp:127
static double phi
FTensor::Index< 'm', 3 > m
Closest point projection algorithm.
std::array< int, LAST_TAPE > tapesTags
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev) override
Get model parameters from blocks.
MoFEMErrorCode codedHessianW(vector< double * > hessian)
MoFEMErrorCode freeHemholtzEnergy()
Set Hemholtz energy.
MoFEMErrorCode flowPotential()
[Yield function J2]
MoFEMErrorCode codedHessianW(vector< double * > hessian)
MoFEMErrorCode evalF()
[Yield function J2]
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get material parameters from block.
FTensor::Tensor2_symmetric< adouble, 3 > tInternalTensor
MoFEMErrorCode yieldFunction()
Set yield function.
FTensor::Tensor2_symmetric< adouble, 3 > tPlasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tTotalStrain
FTensor::Tensor2_symmetric< adouble, 3 > tDeviator
FTensor::Tensor2_symmetric< adouble, 3 > tElasticStrain
MoFEMErrorCode freeHemholtzEnergy()
[Free energy J2]
FTensor::Tensor2_symmetric< adouble, 3 > tStress
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set material parameters at integration point.
FTensor::Tensor2_symmetric< adouble, 3 > tBackStress
J2 plasticity (Kinematic Isotropic (Linear) Hardening)
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get model parameters from blocks.
MoFEMErrorCode codedHessianW(vector< double * > hessian)
[Paraboloidal free energy]
FTensor::Tensor2_symmetric< adouble, 3 > tDeviator
MoFEMErrorCode evalF()
[Paraboloidal yield function]
FTensor::Tensor2_symmetric< adouble, 3 > tPlasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tTotalStrain
FTensor::Tensor2_symmetric< adouble, 3 > tElasticStrain
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set parameters for ADOL-C of constitutive relations.
MoFEMErrorCode freeHemholtzEnergy()
[Paraboloidal free energy]
FTensor::Tensor2_symmetric< adouble, 3 > tStress
MoFEMErrorCode flowPotential()
[Paraboloidal yield function]
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.