v0.14.0
Loading...
Searching...
No Matches
PlasticOpsGeneric.hpp
Go to the documentation of this file.
1
2
3/** \file PlasticOpsGeneric.hpp
4 * \example PlasticOpsGeneric.hpp
5 */
6
7namespace PlasticOps {
8
9//! [Lambda functions]
10template <int DIM> inline auto diff_tensor(FTensor::Number<DIM>) {
17 t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
18 return t_diff;
19};
20
21template <int DIM> inline auto symm_L_tensor(FTensor::Number<DIM>) {
24 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
27 t_L(i, j, L) = 0;
28 if constexpr (DIM == 2) {
29 t_L(0, 0, 0) = 1;
30 t_L(1, 0, 1) = 1;
31 t_L(1, 1, 2) = 1;
32 } else {
33 t_L(0, 0, 0) = 1;
34 t_L(1, 0, 1) = 1;
35 t_L(2, 0, 2) = 1;
36 t_L(1, 1, 3) = 1;
37 t_L(2, 1, 4) = 1;
38 t_L(2, 2, 5) = 1;
39 }
40 return t_L;
41}
42
43template <int DIM> inline auto diff_symmetrize(FTensor::Number<DIM>) {
44
49
51
52 t_diff(i, j, k, l) = 0;
53 t_diff(0, 0, 0, 0) = 1;
54 t_diff(1, 1, 1, 1) = 1;
55
56 t_diff(1, 0, 1, 0) = 0.5;
57 t_diff(1, 0, 0, 1) = 0.5;
58
59 t_diff(0, 1, 0, 1) = 0.5;
60 t_diff(0, 1, 1, 0) = 0.5;
61
62 if constexpr (DIM == 3) {
63 t_diff(2, 2, 2, 2) = 1;
64
65 t_diff(2, 0, 2, 0) = 0.5;
66 t_diff(2, 0, 0, 2) = 0.5;
67 t_diff(0, 2, 0, 2) = 0.5;
68 t_diff(0, 2, 2, 0) = 0.5;
69
70 t_diff(2, 1, 2, 1) = 0.5;
71 t_diff(2, 1, 1, 2) = 0.5;
72 t_diff(1, 2, 1, 2) = 0.5;
73 t_diff(1, 2, 2, 1) = 0.5;
74 }
75
76 return t_diff;
77};
78
79template <typename T>
80inline double trace(FTensor::Tensor2_symmetric<T, 2> &t_stress) {
81 constexpr double third = boost::math::constants::third<double>();
82 return (t_stress(0, 0) + t_stress(1, 1)) * third;
83};
84
85template <typename T>
86inline double trace(FTensor::Tensor2_symmetric<T, 3> &t_stress) {
87 constexpr double third = boost::math::constants::third<double>();
88 return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) * third;
89};
90
91template <typename T, int DIM>
92inline auto deviator(
93
96
97) {
99 t_dev(I, J) = 0;
100 for (int ii = 0; ii != DIM; ++ii)
101 for (int jj = ii; jj != DIM; ++jj)
102 t_dev(ii, jj) = t_stress(ii, jj);
103 t_dev(0, 0) -= trace;
104 t_dev(1, 1) -= trace;
105 t_dev(2, 2) -= trace;
106 for (int ii = 0; ii != DIM; ++ii)
107 for (int jj = ii; jj != DIM; ++jj)
108 t_dev(ii, jj) -= t_alpha(ii, jj);
109 return t_dev;
110};
111
112template <typename T>
113inline auto deviator(FTensor::Tensor2_symmetric<T, 2> &t_stress, double trace,
115 return deviator(t_stress, trace, t_alpha, FTensor::Number<2>());
116};
117
118template <typename T>
119inline auto deviator(FTensor::Tensor2_symmetric<T, 3> &t_stress, double trace,
121 return deviator(t_stress, trace, t_alpha, FTensor::Number<3>());
122};
123
124template <int DIM>
129 FTensor::Ddg<double, 3, DIM> t_diff_deviator;
130 t_diff_deviator(I, J, k, l) = 0;
131 for (int ii = 0; ii != DIM; ++ii)
132 for (int jj = ii; jj != DIM; ++jj)
133 for (int kk = 0; kk != DIM; ++kk)
134 for (int ll = kk; ll != DIM; ++ll)
135 t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
136
137 constexpr double third = boost::math::constants::third<double>();
138
139 t_diff_deviator(0, 0, 0, 0) -= third;
140 t_diff_deviator(0, 0, 1, 1) -= third;
141
142 t_diff_deviator(1, 1, 0, 0) -= third;
143 t_diff_deviator(1, 1, 1, 1) -= third;
144
145 t_diff_deviator(2, 2, 0, 0) -= third;
146 t_diff_deviator(2, 2, 1, 1) -= third;
147
148 if constexpr (DIM == 3) {
149 t_diff_deviator(0, 0, 2, 2) -= third;
150 t_diff_deviator(1, 1, 2, 2) -= third;
151 t_diff_deviator(2, 2, 2, 2) -= third;
152 }
153
154 return t_diff_deviator;
155};
156
157inline auto diff_deviator(FTensor::Ddg<double, 2, 2> &&t_diff_stress) {
158 return diff_deviator(std::move(t_diff_stress), FTensor::Number<2>());
159}
160
161inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
162 return diff_deviator(std::move(t_diff_stress), FTensor::Number<3>());
163}
164
165/**
166 *
167
168\f[
169\begin{split}
170f&=\sqrt{s_{ij}s_{ij}}\\
171A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
172\frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
173\frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
174\sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
175s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
176-A_{mn}A_{ij}
177\right)\\
178\frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
179\\
180\frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
181\frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
182\sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
183\sigma_{mn}} D_{mnkl}
184\end{split}
185\f]
186
187 */
188inline double
190 return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
191 std::numeric_limits<double>::epsilon();
192};
193
194template <int DIM>
195inline auto plastic_flow(long double f,
197 FTensor::Ddg<double, 3, DIM> &&t_diff_deviator) {
201 t_diff_f(k, l) =
202 (1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l))) / f;
203 return t_diff_f;
204};
205
206template <typename T, int DIM>
207inline auto
210 FTensor::Ddg<double, 3, DIM> &&t_diff_deviator) {
216 t_diff_flow(i, j, k, l) =
217 (1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
218 (2. / 3.) * t_flow(i, j) * t_flow(k, l))) /
219 f;
220 return t_diff_flow;
221};
222
223template <typename T, int DIM>
226 FTensor::Ddg<double, DIM, DIM> &&t_diff_plastic_flow_dstress) {
234 t_diff_flow(i, j, k, l) =
235 t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
236 return t_diff_flow;
237};
238
239// inline double constrain_diff_sign(double x) {
240// const auto y = x / zeta;
241// if (y > std::numeric_limits<float>::max_exponent10 ||
242// y < std::numeric_limits<float>::min_exponent10) {
243// return 0;
244// } else {
245// const auto e = std::exp(y);
246// const auto ep1 = e + 1;
247// return (2 / zeta) * (e / (ep1 * ep1));
248// }
249// };
250
251inline double constrian_sign(double x) {
252 const auto y = x / zeta;
253 if (y > std::numeric_limits<float>::max_exponent10 ||
254 y < std::numeric_limits<float>::min_exponent10) {
255 if (x > 0)
256 return 1.;
257 else
258 return -1.;
259 } else {
260 const auto e = std::exp(y);
261 return (e - 1) / (1 + e);
262 }
263};
264
265inline double constrain_abs(double x) {
266 const auto y = -x / zeta;
267 if (y > std::numeric_limits<float>::max_exponent10 ||
268 y < std::numeric_limits<float>::min_exponent10) {
269 return std::abs(x);
270 } else {
271 const double e = std::exp(y);
272 return x + 2 * zeta * std::log1p(e);
273 }
274};
275
276inline double w(double eqiv, double dot_tau, double f, double sigma_y,
277 double sigma_Y) {
278 return (f - sigma_y) / sigma_Y + cn1 * (eqiv);
279};
280
281/**
282
283\f[
284\dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
285\| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0 \\
286c_n \sigma_y \dot{\tau} - \frac{1}{2}\left\{c_n\sigma_y \dot{\tau} +
287(f(\pmb\sigma) - \sigma_y) +
288\| c_n \sigma_y \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
289\f]
290
291 */
292inline double constraint(double eqiv, double dot_tau, double f, double sigma_y,
293 double abs_w, double vis_H, double sigma_Y) {
294 return vis_H * dot_tau +
295 (sigma_Y / 2) * ((cn0 * (dot_tau - eqiv) + cn1 * (eqiv) -
296 (f - sigma_y) / sigma_Y) -
297 abs_w);
298};
299
300inline double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau,
301 double vis_H, double sigma_Y) {
302 return vis_H + (sigma_Y / 2) * (cn0);
303};
304
305inline double diff_constrain_deqiv(double sign, double eqiv, double dot_tau,
306 double sigma_Y) {
307 return (sigma_Y / 2) * (-cn0 + cn1 * (1 - sign));
308};
309
310inline auto diff_constrain_df(double sign) { return (-1 - sign) / 2; };
311
312inline auto diff_constrain_dsigma_y(double sign) { return (1 + sign) / 2; }
313
314template <typename T, int DIM>
315inline auto
317 FTensor::Tensor2_symmetric<T, DIM> &t_plastic_flow) {
320 FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstress;
321 t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
322 return t_diff_constrain_dstress;
323};
324
325template <typename T1, typename T2, int DIM>
326inline auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress) {
331 FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstrain;
332 t_diff_constrain_dstrain(k, l) =
333 t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
334 return t_diff_constrain_dstrain;
335};
336
337template <typename T, int DIM>
339 FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain_dot) {
342 constexpr double A = 2. / 3;
343 return std::sqrt(A * t_plastic_strain_dot(i, j) *
344 t_plastic_strain_dot(i, j)) +
345 std::numeric_limits<double>::epsilon();
346};
347
348template <typename T1, typename T2, typename T3, int DIM>
349inline auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot,
350 T3 &t_diff_plastic_strain,
356 constexpr double A = 2. / 3;
358 t_diff_eqiv(i, j) = A * (t_plastic_strain_dot(k, l) / eqiv) *
359 t_diff_plastic_strain(k, l, i, j);
360 return t_diff_eqiv;
361};
362
363//! [Lambda functions]
364
365//! [Auxiliary functions functions
366static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
369 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
370 &mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
371}
372
373static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
376 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
377 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
378 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
379 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
380 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
381 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
382}
383
384//! [Auxiliary functions functions
385
386template <int DIM, typename DomainEleOp>
388 : public DomainEleOp {
390 boost::shared_ptr<CommonData> common_data_ptr);
391 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
392
393protected:
394 boost::shared_ptr<CommonData> commonDataPtr;
395};
396
397template <int DIM, typename DomainEleOp>
400 boost::shared_ptr<CommonData> common_data_ptr)
402 commonDataPtr(common_data_ptr) {
403 // Operator is only executed for vertices
404 std::fill(&DomainEleOp::doEntities[MBEDGE],
405 &DomainEleOp::doEntities[MBMAXTYPE], false);
406}
407
408template <int DIM, typename DomainEleOp>
410 int side, EntityType type, EntData &data) {
412
415
416 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
417 auto t_stress =
418 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
419 auto t_plastic_strain =
420 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
421
422 commonDataPtr->plasticSurface.resize(nb_gauss_pts, false);
423 commonDataPtr->plasticFlow.resize(size_symm, nb_gauss_pts, false);
424 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
425
426 auto &params = commonDataPtr->blockParams;
427
428 for (auto &f : commonDataPtr->plasticSurface) {
429
430 f = platsic_surface(
431
432 deviator(
433 t_stress, trace(t_stress),
434 kinematic_hardening(t_plastic_strain, params[CommonData::C1_k]))
435
436 );
437
438 auto t_flow_tmp =
439 plastic_flow(f,
440
441 deviator(t_stress, trace(t_stress),
442 kinematic_hardening(t_plastic_strain,
443 params[CommonData::C1_k])),
444
446 t_flow(i, j) = t_flow_tmp(i, j);
447
448 ++t_plastic_strain;
449 ++t_flow;
450 ++t_stress;
451 }
452
454}
455
456template <int DIM, typename DomainEleOp>
458 OpCalculatePlasticityImpl(const std::string field_name,
459 boost::shared_ptr<CommonData> common_data_ptr,
460 boost::shared_ptr<MatrixDouble> m_D_ptr);
461 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
462
463protected:
464 boost::shared_ptr<CommonData> commonDataPtr;
465 boost::shared_ptr<MatrixDouble> mDPtr;
466};
467
468template <int DIM, typename DomainEleOp>
470 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
471 boost::shared_ptr<MatrixDouble> m_D_ptr)
473 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
474 // Opetor is only executed for vertices
475 std::fill(&DomainEleOp::doEntities[MBEDGE],
476 &DomainEleOp::doEntities[MBMAXTYPE], false);
477}
478
479template <int DIM, typename DomainEleOp>
481 int side, EntityType type, EntData &data) {
483
490
491 auto &params = commonDataPtr->blockParams; ///< material parameters
492
493 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
495 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
496 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
497 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
498 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
499 auto t_plastic_strain =
500 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
501 auto t_plastic_strain_dot =
502 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
503 auto t_stress =
504 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
505
506 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
507 auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
508
509 auto t_diff_plastic_strain = diff_tensor(FTensor::Number<DIM>());
510 auto t_diff_deviator = diff_deviator(diff_tensor(FTensor::Number<DIM>()));
511
512 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstress;
513 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstrain;
514 t_flow_dir_dstress(i, j, k, l) =
515 1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l));
516 t_flow_dir_dstrain(i, j, k, l) =
517 t_flow_dir_dstress(i, j, m, n) * t_D_Op(m, n, k, l);
518
519
520 auto t_alpha_dir =
521 kinematic_hardening_dplastic_strain<DIM>(params[CommonData::C1_k]);
522
523 commonDataPtr->resC.resize(nb_gauss_pts, false);
524 commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
525 commonDataPtr->resCdStrain.resize(size_symm, nb_gauss_pts, false);
526 commonDataPtr->resCdPlasticStrain.resize(size_symm, nb_gauss_pts, false);
527 commonDataPtr->resFlow.resize(size_symm, nb_gauss_pts, false);
528 commonDataPtr->resFlowDtau.resize(size_symm, nb_gauss_pts, false);
529 commonDataPtr->resFlowDstrain.resize(size_symm * size_symm, nb_gauss_pts,
530 false);
531 commonDataPtr->resFlowDstrainDot.resize(size_symm * size_symm, nb_gauss_pts,
532 false);
533
534 commonDataPtr->resC.clear();
535 commonDataPtr->resCdTau.clear();
536 commonDataPtr->resCdStrain.clear();
537 commonDataPtr->resCdPlasticStrain.clear();
538 commonDataPtr->resFlow.clear();
539 commonDataPtr->resFlowDtau.clear();
540 commonDataPtr->resFlowDstrain.clear();
541 commonDataPtr->resFlowDstrainDot.clear();
542
543 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
544 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
545 auto t_res_c_dstrain =
546 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
547 auto t_res_c_plastic_strain =
548 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
549 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
550 auto t_res_flow_dtau =
551 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
552 auto t_res_flow_dstrain =
553 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
554 auto t_res_flow_dplastic_strain =
555 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
556
557 auto next = [&]() {
558 ++t_tau;
559 ++t_tau_dot;
560 ++t_f;
561 ++t_flow;
562 ++t_plastic_strain;
563 ++t_plastic_strain_dot;
564 ++t_stress;
565 ++t_res_c;
566 ++t_res_c_dtau;
567 ++t_res_c_dstrain;
568 ++t_res_c_plastic_strain;
569 ++t_res_flow;
570 ++t_res_flow_dtau;
571 ++t_res_flow_dstrain;
572 ++t_res_flow_dplastic_strain;
573 ++t_w;
574 };
575
576 auto get_avtive_pts = [&]() {
577 int nb_points_avtive_on_elem = 0;
578 int nb_points_on_elem = 0;
579
580 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
581 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
582 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
583 auto t_plastic_strain_dot =
584 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
585
586 for (auto &f : commonDataPtr->plasticSurface) {
587 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
588 const auto ww = w(
589 eqiv, t_tau_dot, t_f,
590 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
591 params[CommonData::BISO], params[CommonData::SIGMA_Y]),
592 params[CommonData::SIGMA_Y]);
593 const auto sign_ww = constrian_sign(ww);
594
595 ++nb_points_on_elem;
596 if (sign_ww > 0) {
597 ++nb_points_avtive_on_elem;
598 }
599
600 ++t_tau;
601 ++t_tau_dot;
602 ++t_f;
603 ++t_plastic_strain_dot;
604 }
605
606 int &active_points = PlasticOps::CommonData::activityData[0];
607 int &avtive_full_elems = PlasticOps::CommonData::activityData[1];
608 int &avtive_elems = PlasticOps::CommonData::activityData[2];
609 int &nb_points = PlasticOps::CommonData::activityData[3];
610 int &nb_elements = PlasticOps::CommonData::activityData[4];
611
612 ++nb_elements;
613 nb_points += nb_points_on_elem;
614 if (nb_points_avtive_on_elem > 0) {
615 ++avtive_elems;
616 active_points += nb_points_avtive_on_elem;
617 if (nb_points_avtive_on_elem == nb_points_on_elem) {
618 ++avtive_full_elems;
619 }
620 }
621
622 if (nb_points_avtive_on_elem != nb_points_on_elem)
623 return 1;
624 else
625 return 0;
626 };
627
628 if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
629 get_avtive_pts();
630 }
631
632 for (auto &f : commonDataPtr->plasticSurface) {
633
634 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
635 auto t_diff_eqiv = diff_equivalent_strain_dot(eqiv, t_plastic_strain_dot,
636 t_diff_plastic_strain,
638
639 const auto sigma_y =
640 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
641 params[CommonData::BISO], params[CommonData::SIGMA_Y]);
642 const auto d_sigma_y =
643 iso_hardening_dtau(t_tau, params[CommonData::H],
644 params[CommonData::QINF], params[CommonData::BISO]);
645
646 auto ww = w(eqiv, t_tau_dot, t_f, sigma_y, params[CommonData::SIGMA_Y]);
647 auto abs_ww = constrain_abs(ww);
648 auto sign_ww = constrian_sign(ww);
649
650 auto c = constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
651 params[CommonData::VIS_H], params[CommonData::SIGMA_Y]);
652 auto c_dot_tau = diff_constrain_ddot_tau(sign_ww, eqiv, t_tau_dot,
653 params[CommonData::VIS_H],
654 params[CommonData::SIGMA_Y]);
655 auto c_equiv = diff_constrain_deqiv(sign_ww, eqiv, t_tau_dot,
656 params[CommonData::SIGMA_Y]);
657 auto c_sigma_y = diff_constrain_dsigma_y(sign_ww);
658 auto c_f = diff_constrain_df(sign_ww);
659
660 auto t_dev_stress = deviator(
661
662 t_stress, trace(t_stress),
663
664 kinematic_hardening(t_plastic_strain, params[CommonData::C1_k])
665
666 );
667
669 t_flow_dir(k, l) = 1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
671 t_flow_dstrain(i, j) = t_flow(k, l) * t_D_Op(k, l, i, j);
672
673 auto get_res_c = [&]() { return c; };
674
675 auto get_res_c_dstrain = [&](auto &t_diff_res) {
676 t_diff_res(i, j) = c_f * t_flow_dstrain(i, j);
677 };
678
679 auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
680 t_diff_res(i, j) = (this->getTSa() * c_equiv) * t_diff_eqiv(i, j);
681 t_diff_res(k, l) -= c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
682 };
683
684 auto get_res_c_dtau = [&]() {
685 return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
686 };
687
688 auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
689 t_diff_res(k, l) = -c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
690 };
691
692 auto get_res_flow = [&](auto &t_res_flow) {
693 const auto a = sigma_y;
694 const auto b = t_tau_dot;
695 t_res_flow(k, l) = a * t_plastic_strain_dot(k, l) - b * t_flow_dir(k, l);
696 };
697
698 auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
699 const auto da = d_sigma_y;
700 const auto db = this->getTSa();
701 t_res_flow_dtau(k, l) =
702 da * t_plastic_strain_dot(k, l) - db * t_flow_dir(k, l);
703 };
704
705 auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
706 const auto b = t_tau_dot;
707 t_res_flow_dstrain(m, n, k, l) = -t_flow_dir_dstrain(m, n, k, l) * b;
708 };
709
710 auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
711 const auto a = sigma_y;
712 t_res_flow_dplastic_strain(m, n, k, l) =
713 (a * this->getTSa()) * t_diff_plastic_strain(m, n, k, l);
714 const auto b = t_tau_dot;
715 t_res_flow_dplastic_strain(m, n, i, j) +=
716 (t_flow_dir_dstrain(m, n, k, l) * t_alpha_dir(k, l, i, j)) * b;
717 };
718
719 t_res_c = get_res_c();
720 get_res_flow(t_res_flow);
721
722 if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
723 t_res_c_dtau = get_res_c_dtau();
724 get_res_c_dstrain(t_res_c_dstrain);
725 get_res_c_dplastic_strain(t_res_c_plastic_strain);
726 get_res_flow_dtau(t_res_flow_dtau);
727 get_res_flow_dstrain(t_res_flow_dstrain);
728 get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
729 }
730
731 next();
732 }
733
735}
736
737template <int DIM, typename DomainEleOp>
738struct OpPlasticStressImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
739 OpPlasticStressImpl(const std::string field_name,
740 boost::shared_ptr<CommonData> common_data_ptr,
741 boost::shared_ptr<MatrixDouble> mDPtr);
742 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
743
744private:
745 boost::shared_ptr<MatrixDouble> mDPtr;
746 boost::shared_ptr<CommonData> commonDataPtr;
747};
748
749template <int DIM, typename DomainEleOp>
751 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
752 boost::shared_ptr<MatrixDouble> m_D_ptr)
754 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
755 // Operator is only executed for vertices
756 std::fill(&DomainEleOp::doEntities[MBEDGE],
757 &DomainEleOp::doEntities[MBMAXTYPE], false);
758}
759
760//! [Calculate stress]
761template <int DIM, typename DomainEleOp>
762MoFEMErrorCode
764 EntData &data) {
766
771
772 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
773 commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
774 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
775 auto t_strain =
776 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
777 auto t_plastic_strain =
778 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
779 auto t_stress =
780 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
781
782 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
783 t_stress(i, j) =
784 t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
785 ++t_strain;
786 ++t_plastic_strain;
787 ++t_stress;
788 }
789
791}
792//! [Calculate stress]
793
794template <int DIM, typename AssemblyDomainEleOp>
796 : public AssemblyDomainEleOp {
798 boost::shared_ptr<CommonData> common_data_ptr,
799 boost::shared_ptr<MatrixDouble> m_D_ptr);
800 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
801
802private:
803 boost::shared_ptr<CommonData> commonDataPtr;
804 boost::shared_ptr<MatrixDouble> mDPtr;
805};
806
807template <int DIM, typename AssemblyDomainEleOp>
810 boost::shared_ptr<CommonData> common_data_ptr,
811 boost::shared_ptr<MatrixDouble> m_D_ptr)
813 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
814
815template <int DIM, typename AssemblyDomainEleOp>
816MoFEMErrorCode
818 EntitiesFieldData::EntData &data) {
820
825 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
827
828 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
829 const auto nb_base_functions = data.getN().size2();
830
831 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
832
834
835 auto next = [&]() { ++t_res_flow; };
836
837 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
838 auto t_base = data.getFTensor0N();
839 auto &nf = AssemblyDomainEleOp::locF;
840 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
841 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
842 ++t_w;
843
845 t_rhs(L) = alpha * (t_res_flow(i, j) * t_L(i, j, L));
846 next();
847
848 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
849 size_t bb = 0;
850 for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
851 t_nf(L) += t_base * t_rhs(L);
852 ++t_base;
853 ++t_nf;
854 }
855 for (; bb < nb_base_functions; ++bb)
856 ++t_base;
857 }
858
860}
861
862template <typename AssemblyDomainEleOp>
864 : public AssemblyDomainEleOp {
866 boost::shared_ptr<CommonData> common_data_ptr,
867 boost::shared_ptr<MatrixDouble> m_D_ptr);
868 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
869
870private:
871 boost::shared_ptr<CommonData> commonDataPtr;
872 boost::shared_ptr<MatrixDouble> mDPtr;
873};
874
875template <typename AssemblyDomainEleOp>
878 boost::shared_ptr<CommonData> common_data_ptr,
879 boost::shared_ptr<MatrixDouble> m_D_ptr)
881 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
882
883template <typename AssemblyDomainEleOp>
884MoFEMErrorCode
886 EntitiesFieldData::EntData &data) {
888
889 const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
890 const size_t nb_base_functions = data.getN().size2();
891
892 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
893
894 auto next = [&]() { ++t_res_c; };
895
896 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
897 auto &nf = AssemblyDomainEleOp::locF;
898 auto t_base = data.getFTensor0N();
899 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
900 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
901 ++t_w;
902 const auto res = alpha * t_res_c;
903 next();
904
905 size_t bb = 0;
906 for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
907 nf[bb] += t_base * res;
908 ++t_base;
909 }
910 for (; bb < nb_base_functions; ++bb)
911 ++t_base;
912 }
913
915}
916
917template <int DIM, typename AssemblyDomainEleOp>
919 : public AssemblyDomainEleOp {
921 const std::string row_field_name, const std::string col_field_name,
922 boost::shared_ptr<CommonData> common_data_ptr,
923 boost::shared_ptr<MatrixDouble> m_D_ptr);
924 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
925 EntitiesFieldData::EntData &col_data);
926
927private:
928 boost::shared_ptr<CommonData> commonDataPtr;
929 boost::shared_ptr<MatrixDouble> mDPtr;
930};
931
932template <int DIM, typename AssemblyDomainEleOp>
935 const std::string row_field_name, const std::string col_field_name,
936 boost::shared_ptr<CommonData> common_data_ptr,
937 boost::shared_ptr<MatrixDouble> m_D_ptr)
938 : AssemblyDomainEleOp(row_field_name, col_field_name,
939 AssemblyDomainEleOp::OPROWCOL),
940 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
941 AssemblyDomainEleOp::sYmm = false;
942}
943
944static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
947 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
948 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
949 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
950}
951
952static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
955 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
956 &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
957 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
958 &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
959 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
960 &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
961 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
962 &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
963 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
964 &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
965 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
966 &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
967}
968
969template <int DIM, typename AssemblyDomainEleOp>
970MoFEMErrorCode
972 EntitiesFieldData::EntData &row_data,
973 EntitiesFieldData::EntData &col_data) {
975
980 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
983
984 auto &locMat = AssemblyDomainEleOp::locMat;
985
986 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
987 const auto nb_row_base_functions = row_data.getN().size2();
988
989 auto t_res_flow_dstrain =
990 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
991 auto t_res_flow_dplastic_strain =
992 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
994
995 auto next = [&]() {
996 ++t_res_flow_dstrain;
997 ++t_res_flow_dplastic_strain;
998 };
999
1000 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1001 auto t_row_base = row_data.getFTensor0N();
1002 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1003 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1004 ++t_w;
1005
1007 t_res_mat(O, L) =
1008 alpha * (t_L(i, j, O) * ((t_res_flow_dplastic_strain(i, j, k, l) -
1009 t_res_flow_dstrain(i, j, k, l)) *
1010 t_L(k, l, L)));
1011 next();
1012
1013 size_t rr = 0;
1014 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1015 auto t_mat = get_mat_tensor_sym_dtensor_sym(rr, locMat,
1017 auto t_col_base = col_data.getFTensor0N(gg, 0);
1018 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
1019 t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
1020 ++t_mat;
1021 ++t_col_base;
1022 }
1023
1024 ++t_row_base;
1025 }
1026
1027 for (; rr < nb_row_base_functions; ++rr)
1028 ++t_row_base;
1029 }
1030
1032}
1033
1034template <int DIM, typename AssemblyDomainEleOp>
1036 : public AssemblyDomainEleOp {
1038 const std::string row_field_name, const std::string col_field_name,
1039 boost::shared_ptr<CommonData> common_data_ptr,
1040 boost::shared_ptr<MatrixDouble> m_D_ptr);
1041 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1042 EntitiesFieldData::EntData &col_data);
1043
1044private:
1045 boost::shared_ptr<CommonData> commonDataPtr;
1046 boost::shared_ptr<MatrixDouble> mDPtr;
1047};
1048
1049template <int DIM, typename AssemblyDomainEleOp>
1051 : public AssemblyDomainEleOp {
1053 const std::string row_field_name, const std::string col_field_name,
1054 boost::shared_ptr<CommonData> common_data_ptr,
1055 boost::shared_ptr<MatrixDouble> m_D_ptr);
1056 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1057 EntitiesFieldData::EntData &col_data);
1058
1059private:
1060 boost::shared_ptr<CommonData> commonDataPtr;
1061 boost::shared_ptr<MatrixDouble> mDPtr;
1062};
1063
1064template <int DIM, typename AssemblyDomainEleOp>
1067 const std::string row_field_name, const std::string col_field_name,
1068 boost::shared_ptr<CommonData> common_data_ptr,
1069 boost::shared_ptr<MatrixDouble> m_D_ptr)
1070 : AssemblyDomainEleOp(row_field_name, col_field_name,
1071 DomainEleOp::OPROWCOL),
1072 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1073 AssemblyDomainEleOp::sYmm = false;
1074}
1075
1076static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1079 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
1080}
1081
1082static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1085 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
1086 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
1087}
1088
1089template <int DIM, typename AssemblyDomainEleOp>
1090MoFEMErrorCode
1092 EntitiesFieldData::EntData &row_data,
1093 EntitiesFieldData::EntData &col_data) {
1095
1098 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1100
1101 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1102 const size_t nb_row_base_functions = row_data.getN().size2();
1103 auto &locMat = AssemblyDomainEleOp::locMat;
1104
1105 const auto type = AssemblyDomainEleOp::getFEType();
1106 const auto nb_nodes = moab::CN::VerticesPerEntity(type);
1107
1108 auto t_res_flow_dtau =
1109 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
1110
1111 auto t_L = symm_L_tensor(FTensor::Number<DIM>());
1112
1113 auto next = [&]() { ++t_res_flow_dtau; };
1114
1115 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1116 auto t_row_base = row_data.getFTensor0N();
1117 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1118 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1119 ++t_w;
1121 t_res_vec(L) = alpha * (t_res_flow_dtau(i, j) * t_L(i, j, L));
1122 next();
1123
1124 size_t rr = 0;
1125 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1126 auto t_mat =
1128 auto t_col_base = col_data.getFTensor0N(gg, 0);
1129 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
1130 t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
1131 ++t_mat;
1132 ++t_col_base;
1133 }
1134 ++t_row_base;
1135 }
1136 for (; rr != nb_row_base_functions; ++rr)
1137 ++t_row_base;
1138 }
1139
1141}
1142
1143template <int DIM, typename AssemblyDomainEleOp>
1145 : public AssemblyDomainEleOp {
1147 const std::string row_field_name, const std::string col_field_name,
1148 boost::shared_ptr<CommonData> common_data_ptr,
1149 boost::shared_ptr<MatrixDouble> mat_D_ptr);
1150 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1151 EntitiesFieldData::EntData &col_data);
1152
1153private:
1154 boost::shared_ptr<CommonData> commonDataPtr;
1155 boost::shared_ptr<MatrixDouble> mDPtr;
1156};
1157
1158template <int DIM, typename AssemblyDomainEleOp>
1161 const std::string row_field_name, const std::string col_field_name,
1162 boost::shared_ptr<CommonData> common_data_ptr,
1163 boost::shared_ptr<MatrixDouble> m_D_ptr)
1164 : AssemblyDomainEleOp(row_field_name, col_field_name,
1165 DomainEleOp::OPROWCOL),
1166 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1167 AssemblyDomainEleOp::sYmm = false;
1168}
1169
1172 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
1173}
1174
1177 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
1178}
1179
1180template <int DIM, typename AssemblyDomainEleOp>
1181MoFEMErrorCode
1183 EntitiesFieldData::EntData &row_data,
1184 EntitiesFieldData::EntData &col_data) {
1186
1189 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1191
1192 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1193 const auto nb_row_base_functions = row_data.getN().size2();
1194
1195 auto t_c_dstrain =
1196 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
1197 auto t_c_dplastic_strain =
1198 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
1199
1200 auto next = [&]() {
1201 ++t_c_dstrain;
1202 ++t_c_dplastic_strain;
1203 };
1204
1206
1207 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1208 auto t_row_base = row_data.getFTensor0N();
1209 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1210 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1211 ++t_w;
1212
1214 t_res_vec(L) =
1215 t_L(i, j, L) * (t_c_dplastic_strain(i, j) - t_c_dstrain(i, j));
1216 next();
1217
1220 size_t rr = 0;
1221 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1222 const auto row_base = alpha * t_row_base;
1223 auto t_col_base = col_data.getFTensor0N(gg, 0);
1224 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
1225 t_mat(L) += (row_base * t_col_base) * t_res_vec(L);
1226 ++t_mat;
1227 ++t_col_base;
1228 }
1229 ++t_row_base;
1230 }
1231 for (; rr != nb_row_base_functions; ++rr)
1232 ++t_row_base;
1233 }
1234
1236}
1237
1238template <typename AssemblyDomainEleOp>
1240 : public AssemblyDomainEleOp {
1242 const std::string row_field_name, const std::string col_field_name,
1243 boost::shared_ptr<CommonData> common_data_ptr);
1244 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1245 EntitiesFieldData::EntData &col_data);
1246
1247private:
1248 boost::shared_ptr<CommonData> commonDataPtr;
1249};
1250
1251template <typename AssemblyDomainEleOp>
1254 const std::string row_field_name, const std::string col_field_name,
1255 boost::shared_ptr<CommonData> common_data_ptr)
1256 : AssemblyDomainEleOp(row_field_name, col_field_name,
1257 DomainEleOp::OPROWCOL),
1258 commonDataPtr(common_data_ptr) {
1259 AssemblyDomainEleOp::sYmm = false;
1260}
1261
1262template <typename AssemblyDomainEleOp>
1263MoFEMErrorCode
1265 EntitiesFieldData::EntData &row_data,
1266 EntitiesFieldData::EntData &col_data) {
1268
1269 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1270 const auto nb_row_base_functions = row_data.getN().size2();
1271
1272 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
1273 auto next = [&]() { ++t_res_c_dtau; };
1274
1275 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1276 auto t_row_base = row_data.getFTensor0N();
1277 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1278 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1279 ++t_w;
1280
1281 const auto res = alpha * (t_res_c_dtau);
1282 next();
1283
1284 auto mat_ptr = AssemblyDomainEleOp::locMat.data().begin();
1285 size_t rr = 0;
1286 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1287 auto t_col_base = col_data.getFTensor0N(gg, 0);
1288 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
1289 *mat_ptr += t_row_base * t_col_base * res;
1290 ++t_col_base;
1291 ++mat_ptr;
1292 }
1293 ++t_row_base;
1294 }
1295 for (; rr < nb_row_base_functions; ++rr)
1296 ++t_row_base;
1297 }
1298
1300}
1301}; // namespace PlasticOps
static Index< 'L', 3 > L
constexpr double third
constexpr double a
Kronecker Delta class symmetric.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, DIM > &t_flow, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
double constrian_sign(double x)
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
auto symm_L_tensor(FTensor::Number< DIM >)
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
auto diff_constrain_dstress(double diff_constrain_df, FTensor::Tensor2_symmetric< T, DIM > &t_plastic_flow)
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain, FTensor::Number< DIM >)
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w, double vis_H, double sigma_Y)
auto diff_constrain_df(double sign)
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:118
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
[Common data]
Definition: PlasticOps.hpp:115
double constrain_abs(double x)
auto diff_symmetrize(FTensor::Number< DIM >)
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:116
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
auto diff_deviator(FTensor::Ddg< double, DIM, DIM > &&t_diff_stress, FTensor::Number< DIM >)
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, DIM, DIM > &t_D, FTensor::Ddg< double, DIM, DIM > &&t_diff_plastic_flow_dstress)
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
auto deviator(FTensor::Tensor2_symmetric< T, DIM > &t_stress, double trace, FTensor::Tensor2_symmetric< double, DIM > &t_alpha, FTensor::Number< DIM >)
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
constexpr AssemblyType A
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:142
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:128
constexpr auto size_symm
Definition: plastic.cpp:42
double zeta
Viscous hardening.
Definition: plastic.cpp:173
double cn0
Definition: plastic.cpp:178
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:123
double cn1
Definition: plastic.cpp:179
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107