v0.14.0
Loading...
Searching...
No Matches
PlasticOperators.cpp
Go to the documentation of this file.
1/** \file PlasticOperators.cpp
2 */
3
4/* This file is part of MoFEM.
5 * MoFEM is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU Lesser General Public License as published by the
7 * Free Software Foundation, either version 3 of the License, or (at your
8 * option) any later version.
9 *
10 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13 * License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17
18#include <MoFEM.hpp>
19using namespace MoFEM;
20using namespace FTensor;
21
24#include <RigidBodies.hpp>
25#include <BlockMatData.hpp>
27#include <BasicFeTools.hpp>
28#include <MatrixFunction.hpp>
29
30#include <ElasticOperators.hpp>
31#include <PlasticOperators.hpp>
32
33namespace OpPlasticTools {
34
37 auto diff_dev = diff_deviator(diff_tensor());
38 this->diffDeviator(I, J, k, l) = diff_dev(I, J, k, l);
40}
41
43 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
45 commonDataPtr(common_data_ptr) {
46 // Operator is only executed for vertices
47 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
48}
49
51 EntityType type,
52 EntData &data) {
54
55 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
56 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
57 auto &diff_dev = commonDataPtr->diffDeviator;
58 commonDataPtr->plasticSurfacePtr->resize(nb_gauss_pts, false);
59 commonDataPtr->plasticFlowPtr->resize(6, nb_gauss_pts, false);
60 auto t_flow = getFTensor2SymmetricFromMat<3>(*commonDataPtr->plasticFlowPtr);
61 auto t_plastic_strain =
62 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
63
64 for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
65
66 auto stress_eff = get_effective_stress(t_stress, t_plastic_strain);
67 f = plastic_surface(deviator(stress_eff));
68
69 auto t_flow_tmp = plastic_flow(f, deviator(stress_eff), diff_dev);
70 t_flow(i, j) = t_flow_tmp(i, j);
71
72 ++t_plastic_strain;
73 ++t_flow;
74 ++t_stress;
75 }
76
78}
79
81 boost::shared_ptr<CommonData> common_data_ptr,
82 boost::shared_ptr<MatrixDouble> mat_d_ptr)
84 commonDataPtr(common_data_ptr), matDPtr(mat_d_ptr) {
85 // Operator is only executed for vertices
86 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
87}
88
89//! [Calculate stress]
91 EntData &data) {
93 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
94 commonDataPtr->mStressPtr->resize(6, nb_gauss_pts);
95
96 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
97 auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
98 auto t_plastic_strain =
99 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
100
101 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
102
103 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
104 t_stress(i, j) =
105 t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
106
107 ++t_strain;
108 ++t_plastic_strain;
109 ++t_stress;
110 }
111
113}
114//! [Calculate stress]
115
117 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
119 commonDataPtr(common_data_ptr) {}
120
123 const size_t nb_dofs = data.getIndices().size();
124 if (nb_dofs) {
125
126 auto get_dt = [&]() {
127 double dt;
128 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
129 dt = 1;
130 return dt;
131 };
132 const auto dt = get_dt();
133
134 auto t_flow =
135 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
136
137 auto t_plastic_strain_dot =
138 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainDotPtr));
139 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
140
141 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
142
143 const size_t nb_integration_pts = data.getN().size1();
144 const size_t nb_base_functions = data.getN().size2();
145 auto t_w = getFTensor0IntegrationWeight();
146 auto t_base = data.getFTensor0N();
147
148 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
149 double alpha = dt * getMeasure() * t_w;
150
151 Tensor2<PackPtr<double *, 6>, 3, 3> t_nf{&locF(0), &locF(1), &locF(2),
152 &locF(1), &locF(3), &locF(4),
153 &locF(2), &locF(4), &locF(5)};
154
155 Tensor2_symmetric<double, 3> t_flow_stress_diff;
156 t_flow_stress_diff(i, j) = t_D(i, j, k, l) * (t_plastic_strain_dot(k, l) -
157 t_tau_dot * t_flow(k, l));
158
159 size_t bb = 0;
160 for (; bb != nb_dofs / 6; ++bb) {
161 t_nf(i, j) += alpha * t_base * t_flow_stress_diff(i, j);
162 ++t_base;
163 ++t_nf;
164 }
165 for (; bb < nb_base_functions; ++bb)
166 ++t_base;
167
168 ++t_flow;
169 ++t_plastic_strain_dot;
170 ++t_tau_dot;
171 ++t_w;
172 }
173 }
174
176}
177
179 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
181 commonDataPtr(common_data_ptr) {}
182
185
186 const size_t nb_dofs = data.getIndices().size();
187 if (nb_dofs) {
188
189 auto get_dt = [&]() {
190 double dt;
191 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
192 dt = 1;
193 return dt;
194 };
195 const auto dt = get_dt();
196
197 auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
198 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
199 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
200 auto t_w = getFTensor0IntegrationWeight();
201
202 auto t_base = data.getFTensor0N();
203 const size_t nb_integration_pts = data.getN().size1();
204 const size_t nb_base_functions = data.getN().size2();
205 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
206 const double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
207 const double beta = alpha * constraint(t_tau_dot, t_f, hardening(t_tau));
208
209 size_t bb = 0;
210 for (; bb != nb_dofs; ++bb) {
211 locF(bb) += beta * t_base;
212 ++t_base;
213 }
214 for (; bb < nb_base_functions; ++bb)
215 ++t_base;
216
217 ++t_tau;
218 ++t_tau_dot;
219 ++t_f;
220 ++t_w;
221 }
222 }
223
225}
226
227template <bool LOGSTRAIN>
230 const std::string row_field_name, const std::string col_field_name,
231 boost::shared_ptr<CommonData> common_data_ptr,
232 boost::shared_ptr<MatrixDouble> mat_D_ptr)
233 : DomainEleOpAssembly(row_field_name, col_field_name,
234 DomainEleOpAssembly::OPROWCOL),
235 commonDataPtr(common_data_ptr), matDPtr(mat_D_ptr) {
236 sYmm = false;
237}
238template <bool LOGSTRAIN>
240 EntData &row_data, EntData &col_data) {
242
243 const size_t nb_row_dofs = row_data.getIndices().size();
244 const size_t nb_col_dofs = col_data.getIndices().size();
245 if (nb_row_dofs && nb_col_dofs) {
246
247 const size_t nb_integration_pts = row_data.getN().size1();
248 const size_t nb_row_base_functions = row_data.getN().size2();
249 auto t_w = getFTensor0IntegrationWeight();
250 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
251 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
252
253 MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
254 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
256 auto t_L = symm_L_tensor();
257 auto t_diff_symmetrize = diff_symmetrize();
258
260 Dg<double, 3, 6> t_DL;
261 if (!LOGSTRAIN) {
262 t_DL(i, j, L) = 0;
263 for (int ii = 0; ii != 3; ++ii)
264 for (int jj = ii; jj != 3; ++jj)
265 for (int kk = 0; kk != 3; ++kk)
266 for (int ll = 0; ll != 3; ++ll)
267 for (int LL = 0; LL != 6; ++LL)
268 t_DL(ii, jj, LL) += t_D(ii, jj, kk, ll) * t_L(kk, ll, LL);
269 }
270
271 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
272 double alpha = getMeasure() * t_w;
273
274 if (LOGSTRAIN) {
276 // D2_symm(i, j, k, l) = 0.;
277
278 // FIXME: that transpose is not necessary
279 // for (int i = 0; i != 3; ++i)
280 // for (int jj = 0; jj != 3; ++jj)
281 // for (int k = 0; k != 3; ++k)
282 // for (int l = 0; l != 3; ++l)
283 // for (int m = 0; m != 3; ++m)
284 // for (int n = 0; n != 3; ++n)
285 // D2_symm(k, l, m, n) +=
286 // t_D(i, jj, m, n) * t_dE_dF(i, jj, k, l);
287 D2_symm(k, l, m, n) = t_D(m, n, i, j) * t_dE_dF(i, j, k, l);
288
289 t_DLs(i, j, L) = 0;
290 for (int ii = 0; ii != 3; ++ii)
291 for (int jj = 0; jj != 3; ++jj)
292 for (int kk = 0; kk != 3; ++kk)
293 for (int ll = 0; ll != 3; ++ll)
294 for (int LL = 0; LL != 6; ++LL)
295 t_DLs(ii, jj, LL) +=
296 D2_symm(ii, jj, kk, ll) * t_L(kk, ll, LL);
297 }
298
299 size_t rr = 0;
300 for (; rr != nb_row_dofs / 3; ++rr) {
301
302 Tensor2<PackPtr<double *, 6>, 3, 6> t_mat{
303 &locMat(3 * rr + 0, 0), &locMat(3 * rr + 0, 1),
304 &locMat(3 * rr + 0, 2), &locMat(3 * rr + 0, 3),
305 &locMat(3 * rr + 0, 4), &locMat(3 * rr + 0, 5),
306 &locMat(3 * rr + 1, 0), &locMat(3 * rr + 1, 1),
307 &locMat(3 * rr + 1, 2), &locMat(3 * rr + 1, 3),
308 &locMat(3 * rr + 1, 4), &locMat(3 * rr + 1, 5),
309 &locMat(3 * rr + 2, 0), &locMat(3 * rr + 2, 1),
310 &locMat(3 * rr + 2, 2), &locMat(3 * rr + 2, 3),
311 &locMat(3 * rr + 2, 4), &locMat(3 * rr + 2, 5)};
312
313 auto t_col_base = col_data.getFTensor0N(gg, 0);
315
316 if (LOGSTRAIN) {
317 t_tmp(i, L) = (t_DLs(i, j, L) * (alpha * t_row_diff_base(j)));
318
319 } else
320 t_tmp(i, L) = (t_DL(i, j, L) * (alpha * t_row_diff_base(j)));
321
322 for (size_t cc = 0; cc != nb_col_dofs / 6; ++cc) {
323
324 t_mat(i, L) -= (t_col_base * t_tmp(i, L));
325
326 ++t_mat;
327 ++t_col_base;
328 }
329
330 ++t_row_diff_base;
331 }
332
333 if (LOGSTRAIN)
334 ++t_dE_dF;
335
336 for (; rr < nb_row_base_functions; ++rr)
337 ++t_row_diff_base;
338 ++t_w;
339 }
340 }
341
343}
344
345template <bool LOGSTRAIN>
347 const std::string row_field_name, const std::string col_field_name,
348 boost::shared_ptr<CommonData> common_data_ptr)
349 : DomainEleOpAssembly(row_field_name, col_field_name,
350 DomainEleOpAssembly::OPROWCOL),
351 commonDataPtr(common_data_ptr) {
352 sYmm = false;
353}
354template <bool LOGSTRAIN>
356
357 EntData &row_data, EntData &col_data) {
359
360 const size_t nb_row_dofs = row_data.getIndices().size();
361 const size_t nb_col_dofs = col_data.getIndices().size();
362 if (nb_row_dofs && nb_col_dofs) {
363
364 auto get_dt = [&]() {
365 double dt;
366 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
367 dt = 1;
368 return dt;
369 };
370 const auto dt = get_dt();
371
372 const size_t nb_integration_pts = row_data.getN().size1();
373 const size_t nb_row_base_functions = row_data.getN().size2();
374 auto t_w = getFTensor0IntegrationWeight();
375
376 auto get_row_base = [&]() {
377 if (commonDataPtr->isDualBase) {
378 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
379 return Tensor0<PackPtr<double *, 1>>(base_ptr);
380 } else {
381 return row_data.getFTensor0N();
382 }
383 };
384 auto t_row_base = get_row_base();
385
386 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
387 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
388 auto &diff_dev = commonDataPtr->diffDeviator;
389 auto t_flow =
390 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
391 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
392 auto t_D_Deviator =
393 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
394
395 MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
396 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
397
398 auto t_diff_symmetrize = diff_symmetrize();
399
400 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
401
402 double alpha = dt * getMeasure() * t_w;
403 auto t_diff_plastic_flow_dstrain = diff_plastic_flow_dstrain(
404 t_D_Deviator, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
405 Ddg<double, 3, 3> t_flow_stress_dstrain;
406 t_flow_stress_dstrain(i, j, k, l) =
407 t_D(i, j, m, n) * t_diff_plastic_flow_dstrain(m, n, k, l);
408 Tensor4<double, 3, 3, 3, 3> t_diff_plastic_flow_stress_dgrad;
409
410 if (!LOGSTRAIN)
411 t_diff_plastic_flow_stress_dgrad(i, j, k, l) =
412 t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
413 else {
414
415 Tensor4<double, 3, 3, 3, 3> t_diff_plastic_flow_stress_dgrads;
416 t_diff_plastic_flow_stress_dgrads(i, j, k, l) =
417 t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
418 t_diff_plastic_flow_stress_dgrad(i, j, k, l) = 0.;
419 for (int i = 0; i != 3; ++i)
420 for (int j = 0; j != 3; ++j)
421 for (int k = 0; k != 3; ++k)
422 for (int l = 0; l != 3; ++l)
423 for (int m = 0; m != 3; ++m)
424 for (int n = 0; n != 3; ++n)
425 t_diff_plastic_flow_stress_dgrad(i, j, k, l) +=
426 t_diff_plastic_flow_stress_dgrads(i, j, m, n) *
427 t_dE_dF(m, n, k, l);
428 }
429
430 size_t rr = 0;
431 for (; rr != nb_row_dofs / 6; ++rr) {
432
433 // t3(T d000, T d001, T d002, T d010, T d011, T d012, T d020, T d021,
434 // T d022, T d100, T d101, T d102, T d110, T d111, T d112, T d120,
435 // T d121, T d122, T d200, T d201, T d202, T d210, T d211, T d212,
436 // T d220, T d221, T d222);
437
438 Tensor3<PackPtr<double *, 3>, 3, 3, 3> t_mat{
439 &locMat(6 * rr + 0, 0), // 000
440 &locMat(6 * rr + 0, 1), // 001
441 &locMat(6 * rr + 0, 2), // 002
442
443 &locMat(6 * rr + 1, 0), // 010
444 &locMat(6 * rr + 1, 1), // 011
445 &locMat(6 * rr + 1, 2), // 012
446
447 &locMat(6 * rr + 2, 0), // 020
448 &locMat(6 * rr + 2, 1), // 021
449 &locMat(6 * rr + 2, 2), // 022
450
451 &locMat(6 * rr + 1, 0), // 100
452 &locMat(6 * rr + 1, 1), // 101
453 &locMat(6 * rr + 1, 2), // 102
454
455 &locMat(6 * rr + 3, 0), // 110
456 &locMat(6 * rr + 3, 1), // 111
457 &locMat(6 * rr + 3, 2), // 112
458
459 &locMat(6 * rr + 4, 0), // 120
460 &locMat(6 * rr + 4, 1), // 121
461 &locMat(6 * rr + 4, 2), // 122
462
463 &locMat(6 * rr + 2, 0), // 200
464 &locMat(6 * rr + 2, 1), // 201
465 &locMat(6 * rr + 2, 2), // 202
466
467 &locMat(6 * rr + 4, 0), // 210
468 &locMat(6 * rr + 4, 1), // 211
469 &locMat(6 * rr + 4, 2), // 212
470
471 &locMat(6 * rr + 5, 0), // 220
472 &locMat(6 * rr + 5, 1), // 221
473 &locMat(6 * rr + 5, 2) // 222
474 };
475
476 const double c0 = alpha * t_row_base * t_tau_dot;
477
478 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
479 for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
480
481 t_mat(i, j, l) -= c0 * (t_diff_plastic_flow_stress_dgrad(i, j, l, k) *
482 t_col_diff_base(k));
483
484 ++t_mat;
485 ++t_col_diff_base;
486 }
487
488 ++t_row_base;
489 }
490 for (; rr < nb_row_base_functions; ++rr)
491 ++t_row_base;
492
493 if (LOGSTRAIN)
494 ++t_dE_dF;
495
496 ++t_w;
497 ++t_f;
498 ++t_flow;
499 ++t_tau_dot;
500 }
501 }
502
504}
505
507 const std::string row_field_name, const std::string col_field_name,
508 boost::shared_ptr<CommonData> common_data_ptr)
509 : DomainEleOpAssembly(row_field_name, col_field_name,
510 DomainEleOpAssembly::OPROWCOL),
511 commonDataPtr(common_data_ptr) {
512 sYmm = false;
513}
514
516 EntData &col_data) {
518
519 const size_t nb_row_dofs = row_data.getIndices().size();
520 const size_t nb_col_dofs = col_data.getIndices().size();
521 if (nb_row_dofs && nb_col_dofs) {
522
523 auto get_dt = [&]() {
524 double dt;
525 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
526 dt = 1;
527 return dt;
528 };
529 const auto dt = get_dt();
530
531 const size_t nb_integration_pts = row_data.getN().size1();
532 const size_t nb_row_base_functions = row_data.getN().size2();
533 auto t_w = getFTensor0IntegrationWeight();
534
535 auto get_row_base = [&]() {
536 if (commonDataPtr->isDualBase) {
537 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
538 return Tensor0<PackPtr<double *, 1>>(base_ptr);
539 } else {
540 return row_data.getFTensor0N();
541 }
542 };
543 auto t_row_base = get_row_base();
544
545 auto &diff_dev = commonDataPtr->diffDeviator;
546 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
547 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
548 auto t_flow =
549 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
550
551 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
552 auto t_D_Deviator =
553 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
554
555 auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
556 bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
557
558 auto t_diff_plastic_strain = diff_tensor();
559
560 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
561 double alpha = dt * getMeasure() * t_w;
562 double c1 = alpha * getTSa();
563 const double c0 = alpha * t_tau_dot;
564
565
566 auto t_diff_plastic_flow_dstrain = diff_plastic_flow_kin_hard_dstrain(
567 t_D_Deviator, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
568 Ddg<double, 3, 3> my_tensor;
569 my_tensor(i, j, k, l) =
570 (t_D(i, j, m, n) * (c1 * t_diff_plastic_strain(m, n, k, l) +
571 c0 * t_diff_plastic_flow_dstrain(m, n, k, l)));
572
573 size_t rr = 0;
574 for (; rr != nb_row_dofs / 6; ++rr) {
575
576 // t4(T d0000, T d0001, T d0002, T d0010, T d0011, T d0012, T d0020,
577 // T d0021, T d0022, T d0100, T d0101, T d0102, T d0110, T d0111,
578 // T d0112, T d0120, T d0121, T d0122, T d0200, T d0201, T d0202,
579 // T d0210, T d0211, T d0212, T d0220, T d0221, T d0222, T d1000,
580 // T d1001, T d1002, T d1010, T d1011, T d1012, T d1020, T d1021,
581 // T d1022, T d1100, T d1101, T d1102, T d1110, T d1111, T d1112,
582 // T d1120, T d1121, T d1122, T d1200, T d1201, T d1202, T d1210,
583 // T d1211, T d1212, T d1220, T d1221, T d1222, T d2000, T d2001,
584 // T d2002, T d2010, T d2011, T d2012, T d2020, T d2021, T d2022,
585 // T d2100, T d2101, T d2102, T d2110, T d2111, T d2112, T d2120,
586 // T d2121, T d2122, T d2200, T d2201, T d2202, T d2210, T d2211,
587 // T d2212, T d2220, T d2221, T d2222);
588
589 Tensor4<PackPtr<double *, 6>, 3, 3, 3, 3> t_mat{
590
591 &locMat(6 * rr + 0, 0), // 0000
592 &locMat(6 * rr + 0, 1), // 0001
593 &locMat(6 * rr + 0, 2), // 0002
594 &locMat(6 * rr + 0, 1), // 0010
595 &locMat(6 * rr + 0, 3), // 0011
596 &locMat(6 * rr + 0, 4), // 0012
597 &locMat(6 * rr + 0, 2), // 0020
598 &locMat(6 * rr + 0, 4), // 0021
599 &locMat(6 * rr + 0, 5), // 0022
600
601 &locMat(6 * rr + 1, 0), // 0100
602 &locMat(6 * rr + 1, 1), // 0101
603 &locMat(6 * rr + 1, 2), // 0102
604 &locMat(6 * rr + 1, 1), // 0110
605 &locMat(6 * rr + 1, 3), // 0111
606 &locMat(6 * rr + 1, 4), // 0112
607 &locMat(6 * rr + 1, 2), // 0120
608 &locMat(6 * rr + 1, 4), // 0121
609 &locMat(6 * rr + 1, 5), // 0122
610
611 &locMat(6 * rr + 2, 0), // 0200
612 &locMat(6 * rr + 2, 1), // 0201
613 &locMat(6 * rr + 2, 2), // 0202
614 &locMat(6 * rr + 2, 1), // 0210
615 &locMat(6 * rr + 2, 3), // 0211
616 &locMat(6 * rr + 2, 4), // 0212
617 &locMat(6 * rr + 2, 2), // 0220
618 &locMat(6 * rr + 2, 4), // 0221
619 &locMat(6 * rr + 2, 5), // 0222
620
621 &locMat(6 * rr + 1, 0), // 1000
622 &locMat(6 * rr + 1, 1), // 1001
623 &locMat(6 * rr + 1, 2), // 1002
624 &locMat(6 * rr + 1, 1), // 1010
625 &locMat(6 * rr + 1, 3), // 1011
626 &locMat(6 * rr + 1, 4), // 1012
627 &locMat(6 * rr + 1, 2), // 1020
628 &locMat(6 * rr + 1, 4), // 1021
629 &locMat(6 * rr + 1, 5), // 1022
630
631 &locMat(6 * rr + 3, 0), // 1100
632 &locMat(6 * rr + 3, 1), // 1101
633 &locMat(6 * rr + 3, 2), // 1102
634 &locMat(6 * rr + 3, 1), // 1110
635 &locMat(6 * rr + 3, 3), // 1111
636 &locMat(6 * rr + 3, 4), // 1112
637 &locMat(6 * rr + 3, 2), // 1120
638 &locMat(6 * rr + 3, 4), // 1121
639 &locMat(6 * rr + 3, 5), // 1122
640
641 &locMat(6 * rr + 4, 0), // 1200
642 &locMat(6 * rr + 4, 1), // 1201
643 &locMat(6 * rr + 4, 2), // 1202
644 &locMat(6 * rr + 4, 1), // 1210
645 &locMat(6 * rr + 4, 3), // 1211
646 &locMat(6 * rr + 4, 4), // 1212
647 &locMat(6 * rr + 4, 2), // 1220
648 &locMat(6 * rr + 4, 4), // 1221
649 &locMat(6 * rr + 4, 5), // 1222
650
651 &locMat(6 * rr + 2, 0), // 2000
652 &locMat(6 * rr + 2, 1), // 2001
653 &locMat(6 * rr + 2, 2), // 2002
654 &locMat(6 * rr + 2, 1), // 2010
655 &locMat(6 * rr + 2, 3), // 2011
656 &locMat(6 * rr + 2, 4), // 2012
657 &locMat(6 * rr + 2, 2), // 2020
658 &locMat(6 * rr + 2, 4), // 2021
659 &locMat(6 * rr + 2, 5), // 2022
660
661 &locMat(6 * rr + 4, 0), // 2100
662 &locMat(6 * rr + 4, 1), // 2101
663 &locMat(6 * rr + 4, 2), // 2102
664 &locMat(6 * rr + 4, 1), // 2110
665 &locMat(6 * rr + 4, 3), // 2111
666 &locMat(6 * rr + 4, 4), // 2112
667 &locMat(6 * rr + 4, 2), // 2120
668 &locMat(6 * rr + 4, 4), // 2121
669 &locMat(6 * rr + 4, 5), // 2122
670
671 &locMat(6 * rr + 5, 0), // 2200
672 &locMat(6 * rr + 5, 1), // 2201
673 &locMat(6 * rr + 5, 2), // 2202
674 &locMat(6 * rr + 5, 1), // 2210
675 &locMat(6 * rr + 5, 3), // 2211
676 &locMat(6 * rr + 5, 4), // 2212
677 &locMat(6 * rr + 5, 2), // 2220
678 &locMat(6 * rr + 5, 4), // 2221
679 &locMat(6 * rr + 5, 5) // 2222
680
681 };
682
683 auto t_col_base = col_data.getFTensor0N(gg, 0);
684 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
685 for (size_t cc = 0; cc != nb_col_dofs / 6; ++cc) {
686 t_mat(i, j, k, l) += my_tensor(i, j, k, l) * t_col_base * t_row_base;
687
688 if (is_rotating)
689 t_mat(i, j, k, l) +=
690 t_row_base *
691 (t_D(i, j, m, n) * (alpha * t_diff_plastic_strain(m, n, k, l) *
692 (t_col_diff_base(i) * t_omega(i))));
693
694 ++t_mat;
695 ++t_col_diff_base;
696 ++t_col_base;
697 }
698
699 ++t_row_base;
700 }
701 for (; rr < nb_row_base_functions; ++rr)
702 ++t_row_base;
703
704 if (is_rotating)
705 ++t_omega;
706
707 ++t_w;
708 ++t_f;
709 ++t_flow;
710 ++t_tau_dot;
711 }
712 }
713
715}
716
718 const std::string row_field_name, const std::string col_field_name,
719 boost::shared_ptr<CommonData> common_data_ptr)
720 : DomainEleOpAssembly(row_field_name, col_field_name,
721 DomainEleOp::OPROWCOL),
722 commonDataPtr(common_data_ptr) {
723 sYmm = false;
724}
725
727 EntData &col_data) {
729
730 const size_t nb_row_dofs = row_data.getIndices().size();
731 const size_t nb_col_dofs = col_data.getIndices().size();
732 if (nb_row_dofs && nb_col_dofs) {
733
734 auto get_dt = [&]() {
735 double dt;
736 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
737 dt = 1;
738 return dt;
739 };
740 const auto dt = get_dt();
741
742 const size_t nb_integration_pts = row_data.getN().size1();
743 auto t_w = getFTensor0IntegrationWeight();
744 auto t_flow =
745 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
746
747 auto get_row_base = [&]() {
748 if (commonDataPtr->isDualBase) {
749 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
750 return Tensor0<PackPtr<double *, 1>>(base_ptr);
751 } else {
752 return row_data.getFTensor0N();
753 }
754 };
755 auto t_row_base = get_row_base();
756
757 auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
758 bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
759
760 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
761
762 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
763 double alpha = dt * t_w * getMeasure();
764 double beta = alpha * getTSa();
765
766 Tensor2_symmetric<double, 3> t_flow_stress;
767 t_flow_stress(i, j) = t_D(i, j, m, n) * t_flow(m, n);
768
769 for (size_t rr = 0; rr != nb_row_dofs / 6; ++rr) {
770
771 Tensor2<PackPtr<double *, 1>, 3, 3> t_mat{
772 &locMat(6 * rr + 0, 0), &locMat(6 * rr + 1, 0),
773 &locMat(6 * rr + 2, 0), &locMat(6 * rr + 1, 0),
774 &locMat(6 * rr + 3, 0), &locMat(6 * rr + 4, 0),
775 &locMat(6 * rr + 2, 0), &locMat(6 * rr + 4, 0),
776 &locMat(6 * rr + 5, 0)};
777
778 auto t_col_base = col_data.getFTensor0N(gg, 0);
779 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
780 for (size_t cc = 0; cc != nb_col_dofs; cc++) {
781 t_mat(i, j) -= beta * t_row_base * t_col_base * t_flow_stress(i, j);
782 if (is_rotating)
783 t_mat(i, j) -= alpha * t_row_base *
784 t_flow_stress(i, j) *
785 (t_col_diff_base(k) * t_omega(k));
786 ++t_mat;
787 ++t_col_base;
788 ++t_col_diff_base;
789 }
790
791 ++t_row_base;
792 }
793 if (is_rotating)
794 ++t_omega;
795 ++t_w;
796 ++t_flow;
797 }
798 }
799
801}
802
805 const std::string row_field_name, const std::string col_field_name,
806 boost::shared_ptr<CommonData> common_data_ptr)
807 : DomainEleOpAssembly(row_field_name, col_field_name,
808 DomainEleOp::OPROWCOL),
809 commonDataPtr(common_data_ptr) {
810 sYmm = false;
811}
812
815 EntData &col_data) {
816 return 0;
817}
818
819template <bool LOGSTRAIN>
821 const std::string row_field_name, const std::string col_field_name,
822 boost::shared_ptr<CommonData> common_data_ptr)
823 : DomainEleOpAssembly(row_field_name, col_field_name,
824 DomainEleOp::OPROWCOL),
825 commonDataPtr(common_data_ptr) {
826 sYmm = false;
827}
828template <bool LOGSTRAIN>
830
831 EntData &row_data, EntData &col_data) {
833
834 const size_t nb_row_dofs = row_data.getIndices().size();
835 const size_t nb_col_dofs = col_data.getIndices().size();
836 if (nb_row_dofs && nb_col_dofs) {
837
838 auto get_dt = [&]() {
839 double dt;
840 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
841 dt = 1;
842 return dt;
843 };
844 const auto dt = get_dt();
845
846 const size_t nb_integration_pts = row_data.getN().size1();
847 const size_t nb_row_base_functions = row_data.getN().size2();
848 auto t_w = getFTensor0IntegrationWeight();
849
850 auto get_row_base = [&]() {
851 if (commonDataPtr->isDualBase) {
852 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
853 return Tensor0<PackPtr<double *, 1>>(base_ptr);
854 } else {
855 return row_data.getFTensor0N();
856 }
857 };
858 auto t_row_base = get_row_base();
859
860 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
861 auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
862 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
863 auto t_flow =
864 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
865 auto t_stress =
866 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
867 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
868 auto t_D_Deviator =
869 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
870 auto t_diff_symmetrize = diff_symmetrize();
871
872 MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
873 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
874
875 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
876 double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
877
878 auto t_diff_constrain_dstrain = diff_constrain_dstrain(
879 t_D_Deviator,
881 diff_constrain_df(t_tau_dot, t_f, hardening(t_tau)), t_flow));
882
883 Tensor2<double, 3, 3> t_diff_constrain_dgrad;
884 if (!LOGSTRAIN) {
885
886 t_diff_constrain_dgrad(k, l) =
887 t_diff_constrain_dstrain(i, j) * t_diff_symmetrize(i, j, k, l);
888 } else {
889
890 Tensor2<double, 3, 3> t_diff_constrain_dgrads;
891 t_diff_constrain_dgrad(k, l) =
892 t_diff_constrain_dstrain(i, j) * t_dE_dF(i, j, k, l);
893 }
894
895 Tensor1<PackPtr<double *, 3>, 3> t_mat{&locMat(0, 0), &locMat(0, 1),
896 &locMat(0, 2)};
897
898 size_t rr = 0;
899 for (; rr != nb_row_dofs; ++rr) {
900
901 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
902 for (size_t cc = 0; cc != nb_col_dofs / 3; cc++) {
903
904 t_mat(i) += alpha * t_row_base * t_diff_constrain_dgrad(i, j) *
905 t_col_diff_base(j);
906
907 ++t_mat;
908 ++t_col_diff_base;
909 }
910
911 ++t_row_base;
912 }
913 for (; rr != nb_row_base_functions; ++rr)
914 ++t_row_base;
915
916 if (LOGSTRAIN)
917 ++t_dE_dF;
918
919 ++t_f;
920 ++t_tau;
921 ++t_tau_dot;
922 ++t_flow;
923 ++t_stress;
924 ++t_w;
925 }
926 }
927
929}
930
932 const std::string row_field_name, const std::string col_field_name,
933 boost::shared_ptr<CommonData> common_data_ptr)
934 : DomainEleOpAssembly(row_field_name, col_field_name,
935 DomainEleOp::OPROWCOL),
936 commonDataPtr(common_data_ptr) {
937 sYmm = false;
938}
939
941 EntData &col_data) {
943
944 const size_t nb_row_dofs = row_data.getIndices().size();
945 const size_t nb_col_dofs = col_data.getIndices().size();
946 if (nb_row_dofs && nb_col_dofs) {
947
948 auto get_dt = [&]() {
949 double dt;
950 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
951 dt = 1;
952 return dt;
953 };
954 const auto dt = get_dt();
955
956 const size_t nb_integration_pts = row_data.getN().size1();
957 const size_t nb_row_base_functions = row_data.getN().size2();
958 auto t_w = getFTensor0IntegrationWeight();
959
960 auto get_row_base = [&]() {
961 if (commonDataPtr->isDualBase) {
962 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
963 return Tensor0<PackPtr<double *, 1>>(base_ptr);
964 } else {
965 return row_data.getFTensor0N();
966 }
967 };
968 auto t_row_base = get_row_base();
969
970 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
971 auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
972 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
973 auto t_flow =
974 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
975
976 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
977 auto t_D_Deviator =
978 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
979
980 auto t_diff_symmetrize = diff_symmetrize();
982 auto t_L = symm_L_tensor();
983
984 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
985 double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
986
987 auto mat_ptr = locMat.data().begin();
988 auto t_diff_constrain_dstrain = diff_constrain_kin_hard_dstrain(
989 t_D_Deviator,
991 diff_constrain_df(t_tau_dot, t_f, hardening(t_tau)), t_flow));
992
993 Tensor1<PackPtr<double *, 6>, 6> t_mat(&locMat(0, 0), &locMat(0, 1),
994 &locMat(0, 2), &locMat(0, 3),
995 &locMat(0, 4), &locMat(0, 5));
996
997 size_t rr = 0;
998 for (; rr != nb_row_dofs; ++rr) {
999
1000 auto t_col_base = col_data.getFTensor0N(gg, 0);
1001 for (size_t cc = 0; cc != nb_col_dofs / 6; cc++) {
1002
1003 t_mat(L) -= (alpha * t_row_base * t_col_base) *
1004 (t_diff_constrain_dstrain(i, j) * t_L(i, j, L));
1005
1006 ++t_mat;
1007 ++t_col_base;
1008 }
1009
1010 ++t_row_base;
1011 }
1012 for (; rr != nb_row_base_functions; ++rr)
1013 ++t_row_base;
1014
1015 ++t_f;
1016 ++t_tau;
1017 ++t_tau_dot;
1018 ++t_flow;
1019 ++t_w;
1020 }
1021 }
1022
1024}
1025
1027 const std::string row_field_name, const std::string col_field_name,
1028 boost::shared_ptr<CommonData> common_data_ptr)
1029 : DomainEleOpAssembly(row_field_name, col_field_name,
1030 DomainEleOp::OPROWCOL),
1031 commonDataPtr(common_data_ptr) {
1032 sYmm = false;
1033}
1034
1036 EntData &col_data) {
1038
1039 const size_t nb_row_dofs = row_data.getIndices().size();
1040 const size_t nb_col_dofs = col_data.getIndices().size();
1041 if (nb_row_dofs && nb_col_dofs) {
1042
1043 auto get_dt = [&]() {
1044 double dt;
1045 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1046 dt = 1;
1047 return dt;
1048 };
1049 const auto dt = get_dt();
1050
1051 const size_t nb_integration_pts = row_data.getN().size1();
1052 const size_t nb_row_base_functions = row_data.getN().size2();
1053 auto t_w = getFTensor0IntegrationWeight();
1054 auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1055 auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1056 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1057
1058 auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
1059 bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
1060
1061 auto &cn = (*cache).cn_pl;
1062
1063 auto get_row_base = [&]() {
1064 if (commonDataPtr->isDualBase) {
1065 double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
1066 return Tensor0<PackPtr<double *, 1>>(base_ptr);
1067 } else {
1068 return row_data.getFTensor0N();
1069 }
1070 };
1071 auto t_row_base = get_row_base();
1072
1073 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1074 const double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
1075 const double t_a = getTSa();
1076 const double c0 =
1077 alpha * /* t_a * */
1078 diff_constrain_dtau_dot(t_tau_dot, t_f, hardening(t_tau));
1079 const double c1 =
1080 alpha * diff_constrain_dsigma_y(t_tau_dot, t_f, hardening(t_tau)) *
1081 hardening_dtau(t_tau);
1082
1083 auto mat_ptr = locMat.data().begin();
1084
1085 size_t rr = 0;
1086 for (; rr != nb_row_dofs; ++rr) {
1087
1088 auto t_col_base = col_data.getFTensor0N(gg, 0);
1089 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1090 for (size_t cc = 0; cc != nb_col_dofs; ++cc) {
1091 if (!is_rotating) {
1092 *mat_ptr += (c0 * t_a + c1) * t_row_base * t_col_base;
1093
1094 } else {
1095 const double c0_p =
1096 c0 * (t_col_base * t_a + (t_col_diff_base(i) * t_omega(i)));
1097 *mat_ptr += (c0_p + c1 * t_col_base) * t_row_base;
1098 }
1099
1100 ++mat_ptr;
1101 ++t_col_base;
1102 ++t_col_diff_base;
1103 }
1104 ++t_row_base;
1105 }
1106 for (; rr < nb_row_base_functions; ++rr)
1107 ++t_row_base;
1108
1109 if (is_rotating)
1110 ++t_omega;
1111
1112 ++t_w;
1113 ++t_f;
1114 ++t_tau;
1115 ++t_tau_dot;
1116 }
1117 }
1118
1120}
1121
1123 const std::string field_name, moab::Interface &post_proc_mesh,
1124 std::vector<EntityHandle> &map_gauss_pts,
1125 boost::shared_ptr<CommonData> common_data_ptr)
1127 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1128 commonDataPtr(common_data_ptr) {
1129 // Operator is only executed for vertices
1130 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1131}
1132
1133//! [Postprocessing]
1135 EntData &data) {
1137
1138 std::array<double, 9> def;
1139 std::fill(def.begin(), def.end(), 0);
1140
1141 auto get_tag = [&](const std::string name, size_t size) {
1142 Tag th;
1143 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
1144 MB_TAG_CREAT | MB_TAG_SPARSE,
1145 def.data());
1146 return th;
1147 };
1148
1149 MatrixDouble3by3 mat(3, 3);
1150
1151 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
1152 mat.clear();
1153 for (size_t r = 0; r != 3; ++r)
1154 for (size_t c = 0; c != 3; ++c)
1155 mat(r, c) = t(r, c);
1156 return mat;
1157 };
1158
1159 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
1160 mat.clear();
1161 mat(0, 0) = t;
1162 return mat;
1163 };
1164
1165 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
1166 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1167 &*mat.data().begin());
1168 };
1169
1170 auto th_plastic_surface = get_tag("PLASTIC_SURFACE", 1);
1171 auto th_tau = get_tag("PLASTIC_MULTIPLIER", 1);
1172
1173 auto th_eqv_ep = get_tag("EQUIV_EP", 1);
1174 auto th_hardening = get_tag("HARDENING", 1);
1175
1176 auto th_plastic_flow = get_tag("PLASTIC_FLOW", 9);
1177 string strain_name = "PLASTIC_STRAIN";
1178
1179 auto th_plastic_strain = get_tag(strain_name, 9);
1180
1181 auto t_flow =
1182 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
1183 auto t_plastic_strain =
1184 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
1185 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
1186
1187 size_t gg = 0;
1188 for (; gg != commonDataPtr->plasticSurfacePtr->size(); ++gg) {
1189 const double f = (*(commonDataPtr->plasticSurfacePtr))[gg];
1190 const double tau = (*(commonDataPtr->plasticTauPtr))[gg];
1191
1192 auto plast_surface = f / (*cache).young_modulus_inv;
1193 CHKERR set_tag(th_plastic_surface, gg, set_scalar(plast_surface));
1194 CHKERR set_tag(th_tau, gg, set_scalar(tau));
1195 // FIXME: add option for more postprocessing
1196 // CHKERR set_tag(th_plastic_flow, gg, set_matrix_3d(t_flow));
1197
1198 // t_plastic_strain(i, j) *= 2.;
1199 // auto exp_plast_strain = exp_map(t_plastic_strain);
1200 // double det_g;
1201 // CHKERR determinantTensor3by3(exp_plast_strain, det_g);
1202 // CHKERR set_tag(th_det_g, gg, set_scalar(det_g));
1203
1204 CHKERR set_tag(th_plastic_strain, gg, set_matrix_3d(t_plastic_strain));
1205 const double eqv_ep =
1206 std::sqrt(t_plastic_strain(i, j) * t_plastic_strain(i, j));
1207 // CHKERR set_tag(th_eqv_ep, gg, set_scalar(eqv_ep));
1208 const double hard = hardening(tau) / (*cache).young_modulus_inv;
1209 CHKERR set_tag(th_hardening, gg, set_scalar(hard));
1210
1211 ++t_flow;
1212 ++t_stress;
1213 ++t_plastic_strain;
1214 }
1215
1217}
1218
1221
1224
1225template struct OpCalculateConstraintLhs_dU<true>;
1227
1228} // namespace OpPlasticTools
static Index< 'J', 3 > J
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:26
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
auto plastic_flow(long double f, Tensor2_symmetric< double, 3 > &&t_dev_stress, Ddg< double, 3, 3 > &t_diff_deviator)
double constraint(long double tau_dot, long double f, long double hardening)
auto deviator(Tensor2_symmetric< T1, 3 > &t_stress)
auto diff_plastic_flow_dstrain(Ddg< T, 3, 3 > &t_D, Ddg< double, 3, 3 > &&t_diff_plastic_flow_dstress)
auto diff_plastic_flow_kin_hard_dstrain(Ddg< T, 3, 3 > &t_D, Ddg< double, 3, 3 > &&t_diff_plastic_flow_dstress)
auto diff_constrain_dstrain(Ddg< T1, 3, 3 > &t_D, Tensor2_symmetric< T2, 3 > &&t_diff_constrain_dstress)
auto plastic_surface(Tensor2_symmetric< double, 3 > &&t_stress_deviator)
double diff_constrain_dtau_dot(long double tau_dot, long double f, long double hardening)
auto diff_constrain_dsigma_y(long double tau_dot, long double f, long double hardening)
auto diff_constrain_kin_hard_dstrain(Ddg< T1, 3, 3 > &t_D, Tensor2_symmetric< T2, 3 > &&t_diff_constrain_dstress)
auto get_effective_stress(Tensor2_symmetric< T, 3 > &t_stress, Tensor2_symmetric< T, 3 > &t_plastic_strain)
auto hardening_dtau(long double tau)
auto diff_deviator(Ddg< double, 3, 3 > &&t_diff_stress)
auto diff_plastic_flow_dstress(long double f, Tensor2_symmetric< T, 3 > &t_flow, Ddg< double, 3, 3 > &t_diff_deviator)
auto diff_constrain_df(long double tau_dot, long double f, long double hardening)
auto diff_tensor()
[Operators definitions]
auto hardening(long double tau)
auto diff_constrain_dstress(long double &&diff_constrain_df, Tensor2_symmetric< T, 3 > &t_plastic_flow)
constexpr IntegrationType I
constexpr double t
plate stiffness
Definition: plate.cpp:59
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)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode calculateDiffDeviator()
Ddg< double, 3, 3 > diffDeviator
OpCalculateConstraintLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculateConstraintLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpCalculateConstraintLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
OpCalculateConstraintRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePlasticFlowLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpCalculatePlasticInternalForceLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpCalculatePlasticInternalForceLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculatePlasticSurfaceAndFlow(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< MatrixDouble > matDPtr
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_d_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< CommonData > commonDataPtr
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
OpPostProcPlastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
moab::Interface & postProcMesh
boost::shared_ptr< CommonData > commonDataPtr