v0.15.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>) {
11 FTensor::Index<'i', DIM> i;
12 FTensor::Index<'j', DIM> j;
13 FTensor::Index<'k', DIM> k;
14 FTensor::Index<'l', DIM> l;
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>) {
22 FTensor::Index<'i', DIM> i;
23 FTensor::Index<'j', DIM> j;
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
45 FTensor::Index<'i', DIM> i;
46 FTensor::Index<'j', DIM> j;
47 FTensor::Index<'k', DIM> k;
48 FTensor::Index<'l', DIM> l;
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>
127 FTensor::Index<'k', DIM> k;
128 FTensor::Index<'l', DIM> l;
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) {
198 FTensor::Index<'k', DIM> k;
199 FTensor::Index<'l', DIM> l;
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) {
211 FTensor::Index<'i', DIM> i;
212 FTensor::Index<'j', DIM> j;
213 FTensor::Index<'k', DIM> k;
214 FTensor::Index<'l', DIM> l;
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) {
227 FTensor::Index<'i', DIM> i;
228 FTensor::Index<'j', DIM> j;
229 FTensor::Index<'k', DIM> k;
230 FTensor::Index<'l', DIM> l;
231 FTensor::Index<'m', DIM> m;
232 FTensor::Index<'n', DIM> n;
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, double dt) {
252 const auto y = x / (zeta / dt);
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, double dt) {
266 const auto y = -x / (zeta / dt);
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 / dt) * 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 (1. / cn1) * ((f - sigma_y) / sigma_Y) + dot_tau;
279};
280
281/**
282
283\f[
284v^H \cdot \dot{\tau} +
285\left(\frac{\sigma_Y}{2}\right) \cdot
286\left(
287 \left(
288 c_{\textrm{n0}} \cdot c_{\textrm{n1}} \cdot \left(\dot{\tau} - \dot{e_{\textrm{p}}}\right)
289 \right) +
290 \left(
291 c_{\textrm{n1}} \cdot \left(\dot{\tau} - \frac{1}{c_{\textrm{n1}}} \cdot \frac{f - \sigma_y}{\sigma_Y} - \text{abs\_w}\right)
292 \right)
293\right)
294\f]
295 */
296inline double constraint(double eqiv, double dot_tau, double f, double sigma_y,
297 double abs_w, double vis_H, double sigma_Y) {
298
299 return
300
301 vis_H * dot_tau +
302 (sigma_Y / 2) *
303 (
304
305 (cn0 * cn1 * ((dot_tau - eqiv)) +
306
307 cn1 * ((dot_tau) - (1. / cn1) * (f - sigma_y) / sigma_Y - abs_w))
308
309 );
310};
311
312inline double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau,
313 double vis_H, double sigma_Y) {
314 return vis_H + (sigma_Y / 2) * (cn0 * cn1 + cn1 * (1 - sign));
315};
316
317inline double diff_constrain_deqiv(double sign, double eqiv, double dot_tau,
318 double sigma_Y) {
319 return (sigma_Y / 2) * (-cn0 * cn1);
320};
321
322inline auto diff_constrain_df(double sign) { return (-1 - sign) / 2; };
323
324inline auto diff_constrain_dsigma_y(double sign) { return (1 + sign) / 2; }
325
326template <typename T, int DIM>
327inline auto
329 FTensor::Tensor2_symmetric<T, DIM> &t_plastic_flow) {
330 FTensor::Index<'i', DIM> i;
331 FTensor::Index<'j', DIM> j;
332 FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstress;
333 t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
334 return t_diff_constrain_dstress;
335};
336
337template <typename T1, typename T2, int DIM>
338inline auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress) {
339 FTensor::Index<'i', DIM> i;
340 FTensor::Index<'j', DIM> j;
341 FTensor::Index<'k', DIM> k;
342 FTensor::Index<'l', DIM> l;
343 FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstrain;
344 t_diff_constrain_dstrain(k, l) =
345 t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
346 return t_diff_constrain_dstrain;
347};
348
349template <typename T, int DIM>
351 FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain_dot) {
352 FTensor::Index<'i', DIM> i;
353 FTensor::Index<'j', DIM> j;
354 constexpr double A = 2. / 3;
355 return std::sqrt(A * t_plastic_strain_dot(i, j) *
356 t_plastic_strain_dot(i, j)) +
357 std::numeric_limits<double>::epsilon();
358};
359
360template <typename T1, typename T2, typename T3, int DIM>
361inline auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot,
362 T3 &t_diff_plastic_strain,
364 FTensor::Index<'i', DIM> i;
365 FTensor::Index<'j', DIM> j;
366 FTensor::Index<'k', DIM> k;
367 FTensor::Index<'l', DIM> l;
368 constexpr double A = 2. / 3;
370 t_diff_eqiv(i, j) = A * (t_plastic_strain_dot(k, l) / eqiv) *
371 t_diff_plastic_strain(k, l, i, j);
372 return t_diff_eqiv;
373};
374
375//! [Lambda functions]
376
377//! [Auxiliary functions functions
378static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
381 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
382 &mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
383}
384
385static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
388 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
389 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
390 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
391 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
392 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
393 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
394}
395
396//! [Auxiliary functions functions
397
398template <int DIM, typename DomainEleOp>
400 : public DomainEleOp {
402 boost::shared_ptr<CommonData> common_data_ptr);
403 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
404
405protected:
406 boost::shared_ptr<CommonData> commonDataPtr;
407};
408
409template <int DIM, typename DomainEleOp>
412 boost::shared_ptr<CommonData> common_data_ptr)
414 commonDataPtr(common_data_ptr) {
415 // Operator is only executed for vertices
416 std::fill(&DomainEleOp::doEntities[MBEDGE],
417 &DomainEleOp::doEntities[MBMAXTYPE], false);
418}
419
420template <int DIM, typename DomainEleOp>
422 int side, EntityType type, EntData &data) {
424
425 FTensor::Index<'i', DIM> i;
426 FTensor::Index<'j', DIM> j;
427
428 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
429 auto t_stress =
430 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
431 auto t_plastic_strain =
432 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
433
434 commonDataPtr->plasticSurface.resize(nb_gauss_pts, false);
435 commonDataPtr->plasticFlow.resize(size_symm, nb_gauss_pts, false);
436 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
437
438 auto &params = commonDataPtr->blockParams;
439
440 for (auto &f : commonDataPtr->plasticSurface) {
441
442 f = platsic_surface(
443
444 deviator(
445 t_stress, trace(t_stress),
446 kinematic_hardening(t_plastic_strain, params[CommonData::C1_k]))
447
448 );
449
450 auto t_flow_tmp =
451 plastic_flow(f,
452
453 deviator(t_stress, trace(t_stress),
454 kinematic_hardening(t_plastic_strain,
455 params[CommonData::C1_k])),
456
458 t_flow(i, j) = t_flow_tmp(i, j);
459
460 ++t_plastic_strain;
461 ++t_flow;
462 ++t_stress;
463 }
464
466}
467
468template <int DIM, typename DomainEleOp>
470 OpCalculatePlasticityImpl(const std::string field_name,
471 boost::shared_ptr<CommonData> common_data_ptr,
472 boost::shared_ptr<MatrixDouble> m_D_ptr);
473 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
474
475protected:
476 boost::shared_ptr<CommonData> commonDataPtr;
477 boost::shared_ptr<MatrixDouble> mDPtr;
478};
479
480template <int DIM, typename DomainEleOp>
482 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
483 boost::shared_ptr<MatrixDouble> m_D_ptr)
485 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
486 // Opetor is only executed for vertices
487 std::fill(&DomainEleOp::doEntities[MBEDGE],
488 &DomainEleOp::doEntities[MBMAXTYPE], false);
489}
490
491template <int DIM, typename DomainEleOp>
493 int side, EntityType type, EntData &data) {
495
496 FTensor::Index<'i', DIM> i;
497 FTensor::Index<'j', DIM> j;
498 FTensor::Index<'k', DIM> k;
499 FTensor::Index<'l', DIM> l;
500 FTensor::Index<'m', DIM> m;
501 FTensor::Index<'n', DIM> n;
502
503 auto &params = commonDataPtr->blockParams; ///< material parameters
504
505 auto nb_gauss_pts = DomainEleOp::getGaussPts().size2();
506 auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
507 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
508 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
509 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
510 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
511 auto t_plastic_strain =
512 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
513 auto t_plastic_strain_dot =
514 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
515 auto t_stress =
516 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
517
518 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
519 auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
520
521 auto t_diff_plastic_strain = diff_tensor(FTensor::Number<DIM>());
522 auto t_diff_deviator = diff_deviator(diff_tensor(FTensor::Number<DIM>()));
523
524 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstress;
525 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstrain;
526 t_flow_dir_dstress(i, j, k, l) =
527 1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l));
528 t_flow_dir_dstrain(i, j, k, l) =
529 t_flow_dir_dstress(i, j, m, n) * t_D_Op(m, n, k, l);
530
531
532 auto t_alpha_dir =
534
535 commonDataPtr->resC.resize(nb_gauss_pts, false);
536 commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
537 commonDataPtr->resCdStrain.resize(size_symm, nb_gauss_pts, false);
538 commonDataPtr->resCdPlasticStrain.resize(size_symm, nb_gauss_pts, false);
539 commonDataPtr->resFlow.resize(size_symm, nb_gauss_pts, false);
540 commonDataPtr->resFlowDtau.resize(size_symm, nb_gauss_pts, false);
541 commonDataPtr->resFlowDstrain.resize(size_symm * size_symm, nb_gauss_pts,
542 false);
543 commonDataPtr->resFlowDstrainDot.resize(size_symm * size_symm, nb_gauss_pts,
544 false);
545
546 commonDataPtr->resC.clear();
547 commonDataPtr->resCdTau.clear();
548 commonDataPtr->resCdStrain.clear();
549 commonDataPtr->resCdPlasticStrain.clear();
550 commonDataPtr->resFlow.clear();
551 commonDataPtr->resFlowDtau.clear();
552 commonDataPtr->resFlowDstrain.clear();
553 commonDataPtr->resFlowDstrainDot.clear();
554
555 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
556 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
557 auto t_res_c_dstrain =
558 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
559 auto t_res_c_plastic_strain =
560 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
561 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
562 auto t_res_flow_dtau =
563 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
564 auto t_res_flow_dstrain =
565 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
566 auto t_res_flow_dplastic_strain =
567 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
568
569 auto next = [&]() {
570 ++t_tau;
571 ++t_tau_dot;
572 ++t_f;
573 ++t_flow;
574 ++t_plastic_strain;
575 ++t_plastic_strain_dot;
576 ++t_stress;
577 ++t_res_c;
578 ++t_res_c_dtau;
579 ++t_res_c_dstrain;
580 ++t_res_c_plastic_strain;
581 ++t_res_flow;
582 ++t_res_flow_dtau;
583 ++t_res_flow_dstrain;
584 ++t_res_flow_dplastic_strain;
585 ++t_w;
586 };
587
588 auto get_avtive_pts = [&]() {
589 int nb_points_avtive_on_elem = 0;
590 int nb_points_on_elem = 0;
591
592 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
593 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
594 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
595 auto t_plastic_strain_dot =
596 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
597
598 auto dt = this->getTStimeStep();
599
600 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
601 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
602 const auto ww = w(
603 eqiv, t_tau_dot, t_f,
604 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
605 params[CommonData::BISO], params[CommonData::SIGMA_Y]),
606 params[CommonData::SIGMA_Y]);
607 const auto sign_ww = constrian_sign(ww, dt);
608
609 ++nb_points_on_elem;
610 if (sign_ww > 0) {
611 ++nb_points_avtive_on_elem;
612 }
613
614 ++t_tau;
615 ++t_tau_dot;
616 ++t_f;
617 ++t_plastic_strain_dot;
618 }
619
620 int &active_points = PlasticOps::CommonData::activityData[0];
621 int &avtive_full_elems = PlasticOps::CommonData::activityData[1];
622 int &avtive_elems = PlasticOps::CommonData::activityData[2];
623 int &nb_points = PlasticOps::CommonData::activityData[3];
624 int &nb_elements = PlasticOps::CommonData::activityData[4];
625
626 ++nb_elements;
627 nb_points += nb_points_on_elem;
628 if (nb_points_avtive_on_elem > 0) {
629 ++avtive_elems;
630 active_points += nb_points_avtive_on_elem;
631 if (nb_points_avtive_on_elem == nb_points_on_elem) {
632 ++avtive_full_elems;
633 }
634 }
635
636 if (nb_points_avtive_on_elem != nb_points_on_elem)
637 return 1;
638 else
639 return 0;
640 };
641
642 if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
643 get_avtive_pts();
644 }
645
646 auto dt = this->getTStimeStep();
647 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
648
649 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
650 auto t_diff_eqiv = diff_equivalent_strain_dot(eqiv, t_plastic_strain_dot,
651 t_diff_plastic_strain,
653
654 const auto sigma_y =
655 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
656 params[CommonData::BISO], params[CommonData::SIGMA_Y]);
657 const auto d_sigma_y =
658 iso_hardening_dtau(t_tau, params[CommonData::H],
659 params[CommonData::QINF], params[CommonData::BISO]);
660
661 auto ww = w(eqiv, t_tau_dot, t_f, sigma_y, params[CommonData::SIGMA_Y]);
662 auto abs_ww = constrain_abs(ww, dt);
663 auto sign_ww = constrian_sign(ww, dt);
664
665 auto c = constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
666 params[CommonData::VIS_H], params[CommonData::SIGMA_Y]);
667 auto c_dot_tau = diff_constrain_ddot_tau(sign_ww, eqiv, t_tau_dot,
668 params[CommonData::VIS_H],
669 params[CommonData::SIGMA_Y]);
670 auto c_equiv = diff_constrain_deqiv(sign_ww, eqiv, t_tau_dot,
671 params[CommonData::SIGMA_Y]);
672 auto c_sigma_y = diff_constrain_dsigma_y(sign_ww);
673 auto c_f = diff_constrain_df(sign_ww);
674
675 auto t_dev_stress = deviator(
676
677 t_stress, trace(t_stress),
678
679 kinematic_hardening(t_plastic_strain, params[CommonData::C1_k])
680
681 );
682
684 t_flow_dir(k, l) = 1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
686 t_flow_dstrain(i, j) = t_flow(k, l) * t_D_Op(k, l, i, j);
687
688 auto get_res_c = [&]() { return c; };
689
690 auto get_res_c_dstrain = [&](auto &t_diff_res) {
691 t_diff_res(i, j) = c_f * t_flow_dstrain(i, j);
692 };
693
694 auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
695 t_diff_res(i, j) = (this->getTSa() * c_equiv) * t_diff_eqiv(i, j);
696 t_diff_res(k, l) -= c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
697 };
698
699 auto get_res_c_dtau = [&]() {
700 return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
701 };
702
703 auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
704 t_diff_res(k, l) = -c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
705 };
706
707 auto get_res_flow = [&](auto &t_res_flow) {
708 const auto a = sigma_y;
709 const auto b = t_tau_dot;
710 t_res_flow(k, l) = a * t_plastic_strain_dot(k, l) - b * t_flow_dir(k, l);
711 };
712
713 auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
714 const auto da = d_sigma_y;
715 const auto db = this->getTSa();
716 t_res_flow_dtau(k, l) =
717 da * t_plastic_strain_dot(k, l) - db * t_flow_dir(k, l);
718 };
719
720 auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
721 const auto b = t_tau_dot;
722 t_res_flow_dstrain(m, n, k, l) = -t_flow_dir_dstrain(m, n, k, l) * b;
723 };
724
725 auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
726 const auto a = sigma_y;
727 t_res_flow_dplastic_strain(m, n, k, l) =
728 (a * this->getTSa()) * t_diff_plastic_strain(m, n, k, l);
729 const auto b = t_tau_dot;
730 t_res_flow_dplastic_strain(m, n, i, j) +=
731 (t_flow_dir_dstrain(m, n, k, l) * t_alpha_dir(k, l, i, j)) * b;
732 };
733
734 t_res_c = get_res_c();
735 get_res_flow(t_res_flow);
736
737 if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
738 t_res_c_dtau = get_res_c_dtau();
739 get_res_c_dstrain(t_res_c_dstrain);
740 get_res_c_dplastic_strain(t_res_c_plastic_strain);
741 get_res_flow_dtau(t_res_flow_dtau);
742 get_res_flow_dstrain(t_res_flow_dstrain);
743 get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
744 }
745
746 next();
747 }
748
750}
751
752template <int DIM, typename DomainEleOp>
753struct OpPlasticStressImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
754
755 /**
756 * @deprecated do not use this constructor
757 */
758 DEPRECATED OpPlasticStressImpl(const std::string field_name,
759 boost::shared_ptr<CommonData> common_data_ptr,
760 boost::shared_ptr<MatrixDouble> mDPtr);
761 OpPlasticStressImpl(boost::shared_ptr<CommonData> common_data_ptr,
762 boost::shared_ptr<MatrixDouble> mDPtr);
763
764 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
765
766private:
767 boost::shared_ptr<MatrixDouble> mDPtr;
768 boost::shared_ptr<CommonData> commonDataPtr;
769};
770
771template <int DIM, typename DomainEleOp>
773 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
774 boost::shared_ptr<MatrixDouble> m_D_ptr)
776 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
777 // Operator is only executed for vertices
778 std::fill(&DomainEleOp::doEntities[MBEDGE],
779 &DomainEleOp::doEntities[MBMAXTYPE], false);
780}
781
782template <int DIM, typename DomainEleOp>
784 boost::shared_ptr<CommonData> common_data_ptr,
785 boost::shared_ptr<MatrixDouble> m_D_ptr)
786 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
787 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
788
789//! [Calculate stress]
790template <int DIM, typename DomainEleOp>
791MoFEMErrorCode
793 EntData &data) {
795
796 FTensor::Index<'i', DIM> i;
797 FTensor::Index<'j', DIM> j;
798 FTensor::Index<'k', DIM> k;
799 FTensor::Index<'l', DIM> l;
800
801 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
802 commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
803 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
804 auto t_strain =
805 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
806 auto t_plastic_strain =
807 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
808 auto t_stress =
809 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
810
811 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
812 t_stress(i, j) =
813 t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
814 ++t_strain;
815 ++t_plastic_strain;
816 ++t_stress;
817 }
818
820}
821//! [Calculate stress]
822
823template <int DIM, typename AssemblyDomainEleOp>
825 : public AssemblyDomainEleOp {
827 boost::shared_ptr<CommonData> common_data_ptr,
828 boost::shared_ptr<MatrixDouble> m_D_ptr);
829 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
830
831private:
832 boost::shared_ptr<CommonData> commonDataPtr;
833 boost::shared_ptr<MatrixDouble> mDPtr;
834};
835
836template <int DIM, typename AssemblyDomainEleOp>
839 boost::shared_ptr<CommonData> common_data_ptr,
840 boost::shared_ptr<MatrixDouble> m_D_ptr)
842 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
843
844template <int DIM, typename AssemblyDomainEleOp>
845MoFEMErrorCode
847 EntitiesFieldData::EntData &data) {
849
850 FTensor::Index<'i', DIM> i;
851 FTensor::Index<'j', DIM> j;
852 FTensor::Index<'k', DIM> k;
853 FTensor::Index<'l', DIM> l;
854 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
856
857 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
858 const auto nb_base_functions = data.getN().size2();
859
860 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
861
863
864 auto next = [&]() { ++t_res_flow; };
865
866 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
867 auto t_base = data.getFTensor0N();
868 auto &nf = AssemblyDomainEleOp::locF;
869 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
870 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
871 ++t_w;
872
874 t_rhs(L) = alpha * (t_res_flow(i, j) * t_L(i, j, L));
875 next();
876
877 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
878 size_t bb = 0;
879 for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
880 t_nf(L) += t_base * t_rhs(L);
881 ++t_base;
882 ++t_nf;
883 }
884 for (; bb < nb_base_functions; ++bb)
885 ++t_base;
886 }
887
889}
890
891template <typename AssemblyDomainEleOp>
893 : public AssemblyDomainEleOp {
895 boost::shared_ptr<CommonData> common_data_ptr,
896 boost::shared_ptr<MatrixDouble> m_D_ptr);
897 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
898
899private:
900 boost::shared_ptr<CommonData> commonDataPtr;
901 boost::shared_ptr<MatrixDouble> mDPtr;
902};
903
904template <typename AssemblyDomainEleOp>
907 boost::shared_ptr<CommonData> common_data_ptr,
908 boost::shared_ptr<MatrixDouble> m_D_ptr)
910 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
911
912template <typename AssemblyDomainEleOp>
913MoFEMErrorCode
915 EntitiesFieldData::EntData &data) {
917
918 const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
919 const size_t nb_base_functions = data.getN().size2();
920
921 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
922
923 auto next = [&]() { ++t_res_c; };
924
925 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
926 auto &nf = AssemblyDomainEleOp::locF;
927 auto t_base = data.getFTensor0N();
928 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
929 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
930 ++t_w;
931 const auto res = alpha * t_res_c;
932 next();
933
934 size_t bb = 0;
935 for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
936 nf[bb] += t_base * res;
937 ++t_base;
938 }
939 for (; bb < nb_base_functions; ++bb)
940 ++t_base;
941 }
942
944}
945
946template <int DIM, typename AssemblyDomainEleOp>
948 : public AssemblyDomainEleOp {
950 const std::string row_field_name, const std::string col_field_name,
951 boost::shared_ptr<CommonData> common_data_ptr,
952 boost::shared_ptr<MatrixDouble> m_D_ptr);
953 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
954 EntitiesFieldData::EntData &col_data);
955
956private:
957 boost::shared_ptr<CommonData> commonDataPtr;
958 boost::shared_ptr<MatrixDouble> mDPtr;
959};
960
961template <int DIM, typename AssemblyDomainEleOp>
964 const std::string row_field_name, const std::string col_field_name,
965 boost::shared_ptr<CommonData> common_data_ptr,
966 boost::shared_ptr<MatrixDouble> m_D_ptr)
967 : AssemblyDomainEleOp(row_field_name, col_field_name,
968 AssemblyDomainEleOp::OPROWCOL),
969 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
970 AssemblyDomainEleOp::sYmm = false;
971}
972
973static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
976 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
977 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
978 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
979}
980
981static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
984 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
985 &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
986 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
987 &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
988 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
989 &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
990 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
991 &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
992 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
993 &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
994 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
995 &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
996}
997
998template <int DIM, typename AssemblyDomainEleOp>
999MoFEMErrorCode
1001 EntitiesFieldData::EntData &row_data,
1002 EntitiesFieldData::EntData &col_data) {
1004
1005 FTensor::Index<'i', DIM> i;
1006 FTensor::Index<'j', DIM> j;
1007 FTensor::Index<'k', DIM> k;
1008 FTensor::Index<'l', DIM> l;
1009 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1010 FTensor::Index<'L', size_symm> L;
1011 FTensor::Index<'O', size_symm> O;
1012
1013 auto &locMat = AssemblyDomainEleOp::locMat;
1014
1015 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1016 const auto nb_row_base_functions = row_data.getN().size2();
1017
1018 auto t_res_flow_dstrain =
1019 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
1020 auto t_res_flow_dplastic_strain =
1021 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
1022 auto t_L = symm_L_tensor(FTensor::Number<DIM>());
1023
1024 auto next = [&]() {
1025 ++t_res_flow_dstrain;
1026 ++t_res_flow_dplastic_strain;
1027 };
1028
1029 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1030 auto t_row_base = row_data.getFTensor0N();
1031 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1032 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1033 ++t_w;
1034
1036 t_res_mat(O, L) =
1037 alpha * (t_L(i, j, O) * ((t_res_flow_dplastic_strain(i, j, k, l) -
1038 t_res_flow_dstrain(i, j, k, l)) *
1039 t_L(k, l, L)));
1040 next();
1041
1042 size_t rr = 0;
1043 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1044 auto t_mat = get_mat_tensor_sym_dtensor_sym(rr, locMat,
1046 auto t_col_base = col_data.getFTensor0N(gg, 0);
1047 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
1048 t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
1049 ++t_mat;
1050 ++t_col_base;
1051 }
1052
1053 ++t_row_base;
1054 }
1055
1056 for (; rr < nb_row_base_functions; ++rr)
1057 ++t_row_base;
1058 }
1059
1061}
1062
1063template <int DIM, typename AssemblyDomainEleOp>
1065 : public 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 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1071 EntitiesFieldData::EntData &col_data);
1072
1073private:
1074 boost::shared_ptr<CommonData> commonDataPtr;
1075 boost::shared_ptr<MatrixDouble> mDPtr;
1076};
1077
1078template <int DIM, typename AssemblyDomainEleOp>
1080 : public AssemblyDomainEleOp {
1082 const std::string row_field_name, const std::string col_field_name,
1083 boost::shared_ptr<CommonData> common_data_ptr,
1084 boost::shared_ptr<MatrixDouble> m_D_ptr);
1085 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1086 EntitiesFieldData::EntData &col_data);
1087
1088private:
1089 boost::shared_ptr<CommonData> commonDataPtr;
1090 boost::shared_ptr<MatrixDouble> mDPtr;
1091};
1092
1093template <int DIM, typename AssemblyDomainEleOp>
1096 const std::string row_field_name, const std::string col_field_name,
1097 boost::shared_ptr<CommonData> common_data_ptr,
1098 boost::shared_ptr<MatrixDouble> m_D_ptr)
1099 : AssemblyDomainEleOp(row_field_name, col_field_name,
1100 DomainEleOp::OPROWCOL),
1101 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1102 AssemblyDomainEleOp::sYmm = false;
1103}
1104
1105static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1108 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
1109}
1110
1111static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1114 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
1115 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
1116}
1117
1118template <int DIM, typename AssemblyDomainEleOp>
1119MoFEMErrorCode
1121 EntitiesFieldData::EntData &row_data,
1122 EntitiesFieldData::EntData &col_data) {
1124
1125 FTensor::Index<'i', DIM> i;
1126 FTensor::Index<'j', DIM> j;
1127 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1128 FTensor::Index<'L', size_symm> L;
1129
1130 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1131 const size_t nb_row_base_functions = row_data.getN().size2();
1132 auto &locMat = AssemblyDomainEleOp::locMat;
1133
1134 auto t_res_flow_dtau =
1135 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
1136
1137 auto t_L = symm_L_tensor(FTensor::Number<DIM>());
1138
1139 auto next = [&]() { ++t_res_flow_dtau; };
1140
1141 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1142 auto t_row_base = row_data.getFTensor0N();
1143 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1144 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1145 ++t_w;
1147 t_res_vec(L) = alpha * (t_res_flow_dtau(i, j) * t_L(i, j, L));
1148 next();
1149
1150 size_t rr = 0;
1151 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1152 auto t_mat =
1154 auto t_col_base = col_data.getFTensor0N(gg, 0);
1155 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
1156 t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
1157 ++t_mat;
1158 ++t_col_base;
1159 }
1160 ++t_row_base;
1161 }
1162 for (; rr != nb_row_base_functions; ++rr)
1163 ++t_row_base;
1164 }
1165
1167}
1168
1169template <int DIM, typename AssemblyDomainEleOp>
1171 : public AssemblyDomainEleOp {
1173 const std::string row_field_name, const std::string col_field_name,
1174 boost::shared_ptr<CommonData> common_data_ptr,
1175 boost::shared_ptr<MatrixDouble> mat_D_ptr);
1176 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1177 EntitiesFieldData::EntData &col_data);
1178
1179private:
1180 boost::shared_ptr<CommonData> commonDataPtr;
1181 boost::shared_ptr<MatrixDouble> mDPtr;
1182};
1183
1184template <int DIM, typename AssemblyDomainEleOp>
1187 const std::string row_field_name, const std::string col_field_name,
1188 boost::shared_ptr<CommonData> common_data_ptr,
1189 boost::shared_ptr<MatrixDouble> m_D_ptr)
1190 : AssemblyDomainEleOp(row_field_name, col_field_name,
1191 DomainEleOp::OPROWCOL),
1192 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1193 AssemblyDomainEleOp::sYmm = false;
1194}
1195
1198 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
1199}
1200
1203 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
1204}
1205
1206template <int DIM, typename AssemblyDomainEleOp>
1207MoFEMErrorCode
1209 EntitiesFieldData::EntData &row_data,
1210 EntitiesFieldData::EntData &col_data) {
1212
1213 FTensor::Index<'i', DIM> i;
1214 FTensor::Index<'j', DIM> j;
1215 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1216 FTensor::Index<'L', size_symm> L;
1217
1218 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1219 const auto nb_row_base_functions = row_data.getN().size2();
1220
1221 auto t_c_dstrain =
1222 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
1223 auto t_c_dplastic_strain =
1224 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
1225
1226 auto next = [&]() {
1227 ++t_c_dstrain;
1228 ++t_c_dplastic_strain;
1229 };
1230
1232
1233 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1234 auto t_row_base = row_data.getFTensor0N();
1235 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1236 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1237 ++t_w;
1238
1240 t_res_vec(L) =
1241 t_L(i, j, L) * (t_c_dplastic_strain(i, j) - t_c_dstrain(i, j));
1242 next();
1243
1244 auto t_mat = get_mat_scalar_dtensor_sym(AssemblyDomainEleOp::locMat,
1246 size_t rr = 0;
1247 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1248 const auto row_base = alpha * t_row_base;
1249 auto t_col_base = col_data.getFTensor0N(gg, 0);
1250 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
1251 t_mat(L) += (row_base * t_col_base) * t_res_vec(L);
1252 ++t_mat;
1253 ++t_col_base;
1254 }
1255 ++t_row_base;
1256 }
1257 for (; rr != nb_row_base_functions; ++rr)
1258 ++t_row_base;
1259 }
1260
1262}
1263
1264template <typename AssemblyDomainEleOp>
1266 : public AssemblyDomainEleOp {
1268 const std::string row_field_name, const std::string col_field_name,
1269 boost::shared_ptr<CommonData> common_data_ptr);
1270 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1271 EntitiesFieldData::EntData &col_data);
1272
1273private:
1274 boost::shared_ptr<CommonData> commonDataPtr;
1275};
1276
1277template <typename AssemblyDomainEleOp>
1280 const std::string row_field_name, const std::string col_field_name,
1281 boost::shared_ptr<CommonData> common_data_ptr)
1282 : AssemblyDomainEleOp(row_field_name, col_field_name,
1283 DomainEleOp::OPROWCOL),
1284 commonDataPtr(common_data_ptr) {
1285 AssemblyDomainEleOp::sYmm = false;
1286}
1287
1288template <typename AssemblyDomainEleOp>
1289MoFEMErrorCode
1291 EntitiesFieldData::EntData &row_data,
1292 EntitiesFieldData::EntData &col_data) {
1294
1295 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1296 const auto nb_row_base_functions = row_data.getN().size2();
1297
1298 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
1299 auto next = [&]() { ++t_res_c_dtau; };
1300
1301 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1302 auto t_row_base = row_data.getFTensor0N();
1303 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1304 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1305 ++t_w;
1306
1307 const auto res = alpha * (t_res_c_dtau);
1308 next();
1309
1310 auto mat_ptr = AssemblyDomainEleOp::locMat.data().begin();
1311 size_t rr = 0;
1312 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1313 auto t_col_base = col_data.getFTensor0N(gg, 0);
1314 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
1315 *mat_ptr += t_row_base * t_col_base * res;
1316 ++t_col_base;
1317 ++mat_ptr;
1318 }
1319 ++t_row_base;
1320 }
1321 for (; rr < nb_row_base_functions; ++rr)
1322 ++t_row_base;
1323 }
1324
1326}
1327}; // namespace PlasticOps
constexpr double third
constexpr double a
EntitiesFieldData::EntData EntData
Data on entities.
Kronecker Delta class symmetric.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define DEPRECATED
Definition definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
double dt
const double c
speed of light (cm/ns)
const double n
refractive index of diffusive medium
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 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)
double constrian_sign(double x, double dt)
FTensor::Index< 'M', 3 > M
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
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
[Common data]
auto diff_symmetrize(FTensor::Number< DIM >)
FTensor::Index< 'J', 3 > J
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)
double constrain_abs(double x, double dt)
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition plastic.cpp:92
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition plastic.cpp:78
constexpr auto size_symm
Definition plastic.cpp:42
double zeta
Viscous hardening.
Definition plastic.cpp:130
double cn0
Definition plastic.cpp:135
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition plastic.cpp:73
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition plastic.cpp:106
double cn1
Definition plastic.cpp:136
constexpr auto field_name
FTensor::Index< 'm', 3 > m
Data on single entity (This is passed as argument to DataOperator::doWork)
static std::array< int, 5 > activityData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)