v0.15.0
Loading...
Searching...
No Matches
ADOLCPlasticity.cpp
Go to the documentation of this file.
1/**
2 * \file ADOLCPlasticity.cpp
3 * \ingroup adoc_plasticity
4 * \brief Operators and data structures for ADOL-C plasticity
5 *
6 */
7
8#include <MoFEM.hpp>
9using namespace MoFEM;
10
11#include <ADOLCPlasticity.hpp>
12
13namespace ADOLCPlasticity {
14
16
17template <typename T1, typename T2, int DIM1, int DIM2>
18double calcAveStrain(bool b_bar, const int nb_gauss_pts,
21 if constexpr (DIM1 != DIM2)
23 "Case of mixed dimension by gradient not implemented");
24
25
26 FTensor::Index<'i', DIM1> i;
27 double ave_tr_strain = 0;
28 if (b_bar) {
29 double v = 0;
30 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
31 v += t_w;
32 ave_tr_strain += t_w * t_grad(i, i);
33 ++t_grad;
34 ++t_w;
35 }
36 ave_tr_strain /= (DIM1 * v);
37 }
38 return ave_tr_strain;
39}
40
41template <typename T1, typename T2, int DIM1, int DIM2, int DIM3>
43calcStrain(bool b_bar, double ave_tr_strain,
45 FTensor::Tensor1<T2, DIM3> &&t_voight_strain) {
46
47 if constexpr (DIM1 != DIM2)
49 "Case of mixed dimension by gradient not implemented");
50
51 FTensor::Index<'i', DIM1> i;
52 FTensor::Index<'j', DIM2> j;
53 FTensor::Index<'Z', DIM3> Z;
54 FTensor::Tensor2_symmetric<double, 3> t_strain{0., 0., 0., 0., 0., 0.};
55 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
56 double tr_strain = t_grad(i, i) / DIM1;
57 if (b_bar) {
58 FTensor::Index<'i', DIM1> i;
59 FTensor::Index<'j', DIM2> j;
61 t_strain(i, j) += (ave_tr_strain - tr_strain) * t_kd(i, j);
62 }
63 {
64 FTensor::Index<'i', 3> i;
65 FTensor::Index<'j', 3> j;
66 t_voight_strain(Z) = tStrainToVoightOp(i, j, Z) * t_strain(i, j);
67 }
68 return t_voight_strain;
69};
70
71template <typename T1, typename T2, int DIM1, int DIM3>
73calcStrain(bool b_bar, double ave_tr_strain,
75 FTensor::Tensor1<T2, DIM3> &&t_voight_strain) {
76
77// TODO: implement this function to calculate b_bar strain
78
79 FTensor::Index<'Z', DIM3> Z;
80 {
81 FTensor::Index<'i', 3> i;
82 FTensor::Index<'j', 3> j;
83 t_voight_strain(Z) = tStrainToVoightOp(i, j, Z) * t_strain(i, j);
84 }
85 return t_voight_strain;
86};
87
88//! [BBar method]
89/**
90 * @brief Calculate base for B-bar method
91 *
92 * @tparam DIM Dimension of problem
93 * @param data access to base functions
94 * @param storage storage for base functions at integration points
95 * @param b_bar flag to sith on B-bar method
96 * @param nb_integration_pts number of integration points
97 * @param w_ptr integration weights
98 * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
99 */
100template <int DIM>
102getFTensor2SymmetricDiffBaseImpl(DataForcesAndSourcesCore::EntData &data,
103 MatrixDouble &storage, const bool b_bar,
104 const int nb_integration_pts, double *w_ptr,
106 const auto nb_dofs = data.getFieldData().size();
107 const auto nb_bases = data.getN().size2();
108
109 storage.resize(nb_integration_pts, 6 * nb_dofs, false);
110
111 FTensor::Index<'I', DIM> I;
112 FTensor::Index<'J', DIM> J;
113 FTensor::Index<'i', 3> i;
114 FTensor::Index<'j', 3> j;
118
119 // get tensorial base
120 auto get_ftensor2_symmetric = [&](const int gg, const int rr) {
121 return getFTensor2SymmetricFromPtr<3>(&storage(gg, 6 * rr));
122 };
123
124 // calculate base
125 auto calc_base = [&]() {
126 auto t_t2_diff = get_ftensor2_symmetric(0, 0);
127 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
128 auto t_diff = data.getFTensor1DiffN<DIM>(gg, 0);
129 for (auto b = 0; b != nb_dofs / DIM; ++b) {
131
132 t_grad(i, j) = 0;
133
134 t_grad(N0, J) = t_diff(J);
135 t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
136 ++t_t2_diff;
137
138 t_grad(i, j) = 0;
139 t_grad(N1, J) = t_diff(J);
140 t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
141 ++t_t2_diff;
142
143 t_grad(i, j) = 0;
144 if constexpr (DIM == 3) {
145 t_grad(N2, J) = t_diff(J);
146 t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
147 ++t_t2_diff;
148 }
149
150 ++t_diff;
151 }
152 }
153 return get_ftensor2_symmetric(0, 0);
154 };
155
156 // calculate volume average of trace
157 auto calc_vol = [&](auto &&t_t2_dff) {
158 std::vector<double> vol(nb_dofs, 0);
159 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
160 for (auto b = 0; b != nb_dofs; ++b) {
161 vol[b] += w_ptr[gg] * t_t2_dff(i, i) / DIM;
162 ++t_t2_dff;
163 }
164 }
165 double sum = 0;
166 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
167 sum += w_ptr[gg];
168 }
169 for (auto &v : vol) {
170 v /= sum;
171 }
172 return vol;
173 };
174
175 // modify base for B-bar method
176 auto make_b_bar = [&](auto &&vol) {
177 auto t_t2_diff = get_ftensor2_symmetric(0, 0);
179 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
180 for (auto b = 0; b != nb_dofs; ++b) {
181 const auto trace = t_t2_diff(J, J) / DIM;
182 t_t2_diff(I, J) += (vol[b] - trace) * t_kd(I, J);
183 ++t_t2_diff;
184 }
185 }
186 return get_ftensor2_symmetric(0, 0);
187 };
188
189 return b_bar ? make_b_bar(calc_vol(calc_base())) : calc_base();
190};
191//! [BBar method]
192
194MakeB::getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
195 MatrixDouble &storage, const bool b_bar,
196 const int nb_integration_pts, double *w_ptr,
199 data, storage, b_bar, nb_integration_pts, w_ptr, FTensor::Number<3>());
200}
201
203MakeB::getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
204 MatrixDouble &storage, const bool b_bar,
205 const int nb_integration_pts, double *w_ptr,
208 data, storage, b_bar, nb_integration_pts, w_ptr, FTensor::Number<2>());
209}
210
211struct OpUpdate : public ForcesAndSourcesCore::UserDataOperator {
212
213 OpUpdate(boost::shared_ptr<CommonData> common_data_ptr);
214
215 MoFEMErrorCode doWork(int side, EntityType type,
216 DataForcesAndSourcesCore::EntData &data);
217
218protected:
219 boost::shared_ptr<CommonData> commonDataPtr;
220};
221
222OpUpdate::OpUpdate(boost::shared_ptr<CommonData> common_data_ptr)
224 UserDataOperator::OPSPACE),
225 commonDataPtr(common_data_ptr) {}
226
227MoFEMErrorCode OpUpdate::doWork(int side, EntityType type,
228 DataForcesAndSourcesCore::EntData &data) {
230
231 int nb_gauss_pts = getGaussPts().size2();
232
233 for (int gg = 0; gg < nb_gauss_pts; gg++) {
234 if (commonDataPtr->deltaGamma[gg] > 0) {
235 auto cp_plastic_strain = commonDataPtr->getPlasticStrain(gg);
236 auto cp_internal_variables = commonDataPtr->getInternalVariables(gg);
237 const auto nb_internal_variables = cp_internal_variables.size();
238 auto plastic_strain =
239 getVectorAdaptor(&commonDataPtr->plasticStrainPtr[gg * 6], 6);
240 auto internal_variables = getVectorAdaptor(
241 &commonDataPtr->internalVariablesPtr[gg * nb_internal_variables],
242 nb_internal_variables);
243 plastic_strain = cp_plastic_strain;
244 internal_variables = cp_internal_variables;
245 }
246 }
247
249}
250
251template <int DIM, StrainType STRAIN>
252struct OpCalculateStressImpl : public ForcesAndSourcesCore::UserDataOperator {
253
255 boost::shared_ptr<CommonData> common_data_ptr,
256 boost::shared_ptr<ClosestPointProjection> cp_ptr,
257 bool calc_lhs);
258
259protected:
261 boost::shared_ptr<CommonData> commonDataPtr;
262 boost::shared_ptr<ClosestPointProjection> cpPtr;
264
267
269
270 MoFEMErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts,
271 const int nb_internal_variables);
272};
273
274template <int DIM, StrainType STRAIN>
275struct OpCalculateStress : public OpCalculateStressImpl<DIM, STRAIN> {
276
278
279 MoFEMErrorCode doWork(int side, EntityType type,
280 DataForcesAndSourcesCore::EntData &data) override;
281};
282template <int DIM>
283struct OpCalculateStress<DIM, LARGE_STRAIN> : public OpCalculateStressImpl<DIM, LARGE_STRAIN> {
284
286
287 MoFEMErrorCode doWork(int side, EntityType type,
288 DataForcesAndSourcesCore::EntData &data) override;
289};
290
291template <int DIM>
292struct OpCalculateStress<DIM, SMALL_STRAIN> : public OpCalculateStressImpl<DIM, SMALL_STRAIN> {
293
295
296 MoFEMErrorCode doWork(int side, EntityType type,
297 DataForcesAndSourcesCore::EntData &data) override;
298};
299
300template <int DIM, StrainType STRAIN>
302 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
303 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs)
305 NOSPACE, UserDataOperator::OPSPACE),
306 mField(m_field), commonDataPtr(common_data_ptr), cpPtr(cp_ptr),
307 calcLhs(calc_lhs) {
308 CHK_THROW_MESSAGE(getTags(), "get tags");
309}
310
311template <int DIM, StrainType STRAIN>
314 int def_length = 0;
315 CHKERR mField.get_moab().tag_get_handle(
316 "PLASTIC_STRAIN", def_length, MB_TYPE_DOUBLE, thPlasticStrain,
317 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, NULL);
318 CHKERR mField.get_moab().tag_get_handle(
319 "INTERNAL_VARIABLES", def_length, MB_TYPE_DOUBLE, thInternalVariables,
320 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, NULL);
322}
323
324template <int DIM, StrainType STRAIN>
326 const EntityHandle tet, const int nb_gauss_pts,
327 const int nb_internal_variables) {
329
330 VectorDouble v;
331 {
332 rval = mField.get_moab().tag_get_by_ptr(
333 thPlasticStrain, &tet, 1,
334 (const void **)&commonDataPtr->plasticStrainPtr,
335 &commonDataPtr->plasticStrainSize);
336 if (rval != MB_SUCCESS ||
337 commonDataPtr->plasticStrainSize != 6 * nb_gauss_pts) {
338 v.resize(6 * nb_gauss_pts, false);
339 v.clear();
340 int tag_size[1];
341 tag_size[0] = v.size();
342 void const *tag_data[] = {&v[0]};
343 CHKERR mField.get_moab().tag_set_by_ptr(thPlasticStrain, &tet, 1,
344 tag_data, tag_size);
345 CHKERR mField.get_moab().tag_get_by_ptr(
346 thPlasticStrain, &tet, 1,
347 (const void **)&commonDataPtr->plasticStrainPtr,
348 &commonDataPtr->plasticStrainSize);
349 }
350 }
351 if (nb_internal_variables > 0) {
352 rval = mField.get_moab().tag_get_by_ptr(
353 thInternalVariables, &tet, 1,
354 (const void **)&commonDataPtr->internalVariablesPtr,
355 &commonDataPtr->internalVariablesSize);
356 if (rval != MB_SUCCESS || commonDataPtr->internalVariablesSize !=
357 nb_internal_variables * nb_gauss_pts) {
358 v.resize(nb_internal_variables * nb_gauss_pts, false);
359 v.clear();
360 int tag_size[1];
361 tag_size[0] = v.size();
362 void const *tag_data[] = {&v[0]};
363 CHKERR mField.get_moab().tag_set_by_ptr(thInternalVariables, &tet, 1,
364 tag_data, tag_size);
365 CHKERR mField.get_moab().tag_get_by_ptr(
366 thInternalVariables, &tet, 1,
367 (const void **)&commonDataPtr->internalVariablesPtr,
368 &commonDataPtr->internalVariablesSize);
369 }
370 }
372}
373
374template <int DIM>
377 DataForcesAndSourcesCore::EntData &data) {
379
380 auto do_work_impl = [this](auto commonDataPtr, auto cpPtr) {
382
383 int nb_gauss_pts = this->getGaussPts().size2();
384
385 int nb_internal_variables = cpPtr->nbInternalVariables;
386 auto tet = this->getNumeredEntFiniteElementPtr()->getEnt();
387 CHKERR this->setTagsData(tet, nb_gauss_pts, nb_internal_variables);
388
389 commonDataPtr->activeVariablesW.resize(
390 nb_gauss_pts, 12 + cpPtr->nbInternalVariables, false);
391 commonDataPtr->gradientW.resize(nb_gauss_pts,
392 12 + cpPtr->nbInternalVariables, false);
393 commonDataPtr->materialTangentOperator.resize(36, nb_gauss_pts, false);
394 commonDataPtr->deltaGamma.resize(nb_gauss_pts);
395
396 FTensor::Index<'i', 3> i;
397 FTensor::Index<'j', 3> j;
398 FTensor::Index<'k', 3> k;
399 FTensor::Index<'l', 3> l;
400 FTensor::Index<'Z', 6> Z;
401 FTensor::Index<'Y', 6> Y;
402
403 // Code uses trial stress, until first solve of linear system of equations
404 int iter = 1;
405 if (this->getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
406 CHKERR SNESGetLinearSolveIterations(this->getFEMethod()->snes, &iter);
407 }
408
409 double ave_tr_strain = 0;
410
411 auto t_voight_stress_op = voight_to_stress_op();
412
413 auto t_grad =
414 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->gradAtGaussPts);
415
416 auto t_Cep =
417 getFTensor4DdgFromMat<3, 3, 1>(commonDataPtr->materialTangentOperator);
418 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
419
420 auto tmp_active =
421 getVectorAdaptor(&(commonDataPtr->activeVariablesW(gg, 0)),
422 12 + cpPtr->nbInternalVariables);
423 cpPtr->activeVariablesW.swap(tmp_active);
424 auto tmp_gradient = getVectorAdaptor(&(commonDataPtr->gradientW(gg, 0)),
425 12 + cpPtr->nbInternalVariables);
426 this->cpPtr->gradientW.swap(tmp_gradient);
427
428 auto t_voigt_strain =
429 calcStrain(commonDataPtr->bBar, ave_tr_strain, t_grad,
430 getFTensor1FromPtr<6>(&(cpPtr->activeVariablesW[6])));
431 ++t_grad;
432
433 // Copy plastic strains and internal variables from tags
434 auto copy_plastic_strain_and_internal_variables = [&]() {
436
437 // Get data stored on Tags
438 VectorAdaptor plastic_strain_tags =
439 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(
440 6, &commonDataPtr->plasticStrainPtr[gg * 6]));
441 VectorAdaptor internal_variables_tags = VectorAdaptor(
442 nb_internal_variables,
443 ublas::shallow_array_adaptor<double>(
444 nb_internal_variables,
445 &commonDataPtr
446 ->internalVariablesPtr[gg * nb_internal_variables]));
447
448 // Set values to closest point projection
449 cpPtr->plasticStrain0.resize(6, false);
450 noalias(cpPtr->plasticStrain0) = plastic_strain_tags;
451 cpPtr->internalVariables0.resize(nb_internal_variables, false);
452 noalias(cpPtr->internalVariables0) = internal_variables_tags;
453
454 auto cp_plastic_strain = cpPtr->getPlasticStrain();
455 auto cp_internal_variables = cpPtr->getInternalVariables();
456 noalias(cp_plastic_strain) = plastic_strain_tags;
457 noalias(cp_internal_variables) = internal_variables_tags;
459 };
460
461 CHKERR copy_plastic_strain_and_internal_variables();
462
463 auto closest_point_projection = [&](auto &recalculate_elastic_tangent) {
467 CHKERR cpPtr->setParams(t, recalculate_elastic_tangent);
468 }
469 CHKERR cpPtr->playW_NoHessian();
470 CHKERR cpPtr->playY_NoGradient();
471 double y = cpPtr->y;
472 cpPtr->deltaGamma = 0; // zero lagrange multiplier
473 if (iter > 0 && y > std::numeric_limits<double>::epsilon()) {
474 CHKERR cpPtr->solveClosestProjection();
475 }
476 commonDataPtr->deltaGamma[gg] = cpPtr->deltaGamma;
478 };
479
480 auto evaluate_consistent_tangent_matrix =
481 [&](auto &recalculate_elastic_tangent) {
483 if (recalculate_elastic_tangent)
484 CHKERR cpPtr->playW();
485
486 if (iter > 0 &&
487 cpPtr->deltaGamma > std::numeric_limits<double>::epsilon()) {
488
489 CHKERR cpPtr->consistentTangent();
490
491 auto &m = cpPtr->Cep;
492 // for generic case of non-associated plasticity consistent
493 // tangent matrix is non-symmetric
494 auto t_voight_cep =
496 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
497 &m(0, 5),
498
499 &m(1, 0), &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4),
500 &m(1, 5),
501
502 &m(2, 0), &m(2, 1), &m(2, 2), &m(2, 3), &m(2, 4),
503 &m(2, 5),
504
505 &m(3, 0), &m(3, 1), &m(3, 2), &m(3, 3), &m(3, 4),
506 &m(3, 5),
507
508 &m(4, 0), &m(4, 1), &m(4, 2), &m(4, 3), &m(4, 4),
509 &m(4, 5),
510
511 &m(5, 0), &m(5, 1), &m(5, 2), &m(5, 3), &m(5, 4), &m(5, 5)
512
513 );
514
515 t_Cep(i, j, k, l) =
516 t_voight_stress_op(i, j, Z) *
517 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
518
519 } else {
520
521 auto &m = cpPtr->C;
522 auto t_voight_cep =
524 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
525 &m(0, 5),
526
527 &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4), &m(1, 5),
528
529 &m(2, 2), &m(2, 3), &m(2, 4), &m(2, 5),
530
531 &m(3, 3), &m(3, 4), &m(3, 5),
532
533 &m(4, 4), &m(4, 5),
534
535 &m(5, 5));
536
537 t_Cep(i, j, k, l) =
538 t_voight_stress_op(i, j, Z) *
539 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
540 }
541
542 ++t_Cep;
544 };
545
546 // Always recalculate elastic tangent if first element
547 bool recalculate_elastic_tangent =
548 (this->getNinTheLoop() == 0) ? true : false;
549 CHKERR closest_point_projection(recalculate_elastic_tangent);
550 // always calculate consistent tangent matrix for large strain
551 CHKERR evaluate_consistent_tangent_matrix(recalculate_elastic_tangent);
552 }
553
554 auto calcs_stress_matrix = [&]() {
556 auto t_voight_stress_op = voight_to_stress_op();
557 commonDataPtr->stressMatrix.resize(
558 6, this->commonDataPtr->gradientW.size1(), false);
559 auto t_stess =
560 getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
561 auto t_voight_stress = commonDataPtr->getFTensor1StressAtGaussPts();
562 for (auto gg = 0; gg != commonDataPtr->gradientW.size1(); ++gg) {
563 t_stess(i, j) = t_voight_stress_op(i, j, Z) * t_voight_stress(Z);
564 ++t_voight_stress;
565 ++t_stess;
566 }
568 };
569
570 auto calcs_plastic_strain_matrix = [&]() {
572 auto t_voight_strain_op = voight_to_stress_op();
573 commonDataPtr->plasticStrainMatrix.resize(
574 6, commonDataPtr->activeVariablesW.size1(), false);
575 auto t_plastic_strain =
576 getFTensor2SymmetricFromMat<3>(commonDataPtr->plasticStrainMatrix);
577 auto t_voight_plastic_strain =
578 commonDataPtr->getFTensor1PlasticStrainAtGaussPts();
579 for (auto gg = 0; gg != commonDataPtr->activeVariablesW.size1(); ++gg) {
580 t_plastic_strain(i, j) =
581 t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
582 ++t_voight_plastic_strain;
583 ++t_plastic_strain;
584 }
586 };
587
588 CHKERR calcs_stress_matrix();
589 CHKERR calcs_plastic_strain_matrix();
591 };
592
593 CHKERR do_work_impl(this->commonDataPtr, this->cpPtr);
594
596}
597
598template <int DIM>
601 DataForcesAndSourcesCore::EntData &data) {
603
604 auto do_work_impl = [this](auto commonDataPtr, auto cpPtr) {
606
607 int nb_gauss_pts = this->getGaussPts().size2();
608
609 int nb_internal_variables = cpPtr->nbInternalVariables;
610 auto tet = this->getNumeredEntFiniteElementPtr()->getEnt();
611 CHKERR this->setTagsData(tet, nb_gauss_pts, nb_internal_variables);
612
613 commonDataPtr->activeVariablesW.resize(
614 nb_gauss_pts, 12 + cpPtr->nbInternalVariables, false);
615 commonDataPtr->gradientW.resize(nb_gauss_pts,
616 12 + cpPtr->nbInternalVariables, false);
617 commonDataPtr->materialTangentOperator.resize(36, nb_gauss_pts, false);
618 commonDataPtr->deltaGamma.resize(nb_gauss_pts);
619
620 FTensor::Index<'i', 3> i;
621 FTensor::Index<'j', 3> j;
622 FTensor::Index<'k', 3> k;
623 FTensor::Index<'l', 3> l;
624 FTensor::Index<'Z', 6> Z;
625 FTensor::Index<'Y', 6> Y;
626
627 // Code uses trial stress, until first solve of linear system of equations
628 int iter = 1;
629 if (this->getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
630 CHKERR SNESGetLinearSolveIterations(this->getFEMethod()->snes, &iter);
631 }
632
633 // b-bar
634 auto ave_tr_strain = calcAveStrain(
635 this->commonDataPtr->bBar, nb_gauss_pts,
636 getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts),
637 this->getFTensor0IntegrationWeight());
638
639 auto t_voight_stress_op = voight_to_stress_op();
640
641 auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts);
642
644 this->commonDataPtr->materialTangentOperator);
645 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
646
647 auto tmp_active =
648 getVectorAdaptor(&(commonDataPtr->activeVariablesW(gg, 0)),
649 12 + cpPtr->nbInternalVariables);
650 cpPtr->activeVariablesW.swap(tmp_active);
651 auto tmp_gradient = getVectorAdaptor(&(commonDataPtr->gradientW(gg, 0)),
652 12 + cpPtr->nbInternalVariables);
653 cpPtr->gradientW.swap(tmp_gradient);
654
655 auto t_voigt_strain =
656 calcStrain(commonDataPtr->bBar, ave_tr_strain, t_grad,
657 getFTensor1FromPtr<6>(&(cpPtr->activeVariablesW[6])));
658 ++t_grad;
659
660 // Copy plastic strains and internal variables from tags
661 auto copy_plastic_strain_and_internal_variables = [&]() {
663
664 // Get data stored on Tags
665 VectorAdaptor plastic_strain_tags =
666 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(
667 6, &commonDataPtr->plasticStrainPtr[gg * 6]));
668 VectorAdaptor internal_variables_tags = VectorAdaptor(
669 nb_internal_variables,
670 ublas::shallow_array_adaptor<double>(
671 nb_internal_variables,
672 &commonDataPtr
673 ->internalVariablesPtr[gg * nb_internal_variables]));
674
675 // Set values to closest point projection
676 cpPtr->plasticStrain0.resize(6, false);
677 noalias(cpPtr->plasticStrain0) = plastic_strain_tags;
678 cpPtr->internalVariables0.resize(nb_internal_variables, false);
679 noalias(cpPtr->internalVariables0) = internal_variables_tags;
680
681 auto cp_plastic_strain = cpPtr->getPlasticStrain();
682 auto cp_internal_variables = cpPtr->getInternalVariables();
683 noalias(cp_plastic_strain) = plastic_strain_tags;
684 noalias(cp_internal_variables) = internal_variables_tags;
686 };
687
688 CHKERR copy_plastic_strain_and_internal_variables();
689
690 auto closest_point_projection = [&](auto &recalculate_elastic_tangent) {
694 CHKERR cpPtr->setParams(t, recalculate_elastic_tangent);
695 }
696 CHKERR cpPtr->playW_NoHessian();
697 CHKERR cpPtr->playY_NoGradient();
698 double y = cpPtr->y;
699 cpPtr->deltaGamma = 0; // zero lagrange multiplier
700 if (iter > 0 && y > std::numeric_limits<double>::epsilon()) {
701 CHKERR cpPtr->solveClosestProjection();
702 }
703 commonDataPtr->deltaGamma[gg] = cpPtr->deltaGamma;
705 };
706
707 auto evaluate_consistent_tangent_matrix =
708 [&](auto &recalculate_elastic_tangent) {
710 if (recalculate_elastic_tangent)
711 CHKERR cpPtr->playW();
712
713 if (iter > 0 &&
714 cpPtr->deltaGamma > std::numeric_limits<double>::epsilon()) {
715
716 CHKERR cpPtr->consistentTangent();
717
718 auto &m = cpPtr->Cep;
719 // for generic case of non-associated plasticity consistent
720 // tangent matrix is non-symmetric
721 auto t_voight_cep =
723 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
724 &m(0, 5),
725
726 &m(1, 0), &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4),
727 &m(1, 5),
728
729 &m(2, 0), &m(2, 1), &m(2, 2), &m(2, 3), &m(2, 4),
730 &m(2, 5),
731
732 &m(3, 0), &m(3, 1), &m(3, 2), &m(3, 3), &m(3, 4),
733 &m(3, 5),
734
735 &m(4, 0), &m(4, 1), &m(4, 2), &m(4, 3), &m(4, 4),
736 &m(4, 5),
737
738 &m(5, 0), &m(5, 1), &m(5, 2), &m(5, 3), &m(5, 4), &m(5, 5)
739
740 );
741
742 t_Cep(i, j, k, l) =
743 t_voight_stress_op(i, j, Z) *
744 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
745
746 } else {
747
748 auto &m = cpPtr->C;
749 auto t_voight_cep =
751 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
752 &m(0, 5),
753
754 &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4), &m(1, 5),
755
756 &m(2, 2), &m(2, 3), &m(2, 4), &m(2, 5),
757
758 &m(3, 3), &m(3, 4), &m(3, 5),
759
760 &m(4, 4), &m(4, 5),
761
762 &m(5, 5));
763
764 t_Cep(i, j, k, l) =
765 t_voight_stress_op(i, j, Z) *
766 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
767 }
768
769 ++t_Cep;
771 };
772
773 // Always recalculate elastic tangent if first element
774 bool recalculate_elastic_tangent =
775 (this->getNinTheLoop() == 0) ? true : false;
776 CHKERR closest_point_projection(recalculate_elastic_tangent);
777 if (this->calcLhs)
778 CHKERR evaluate_consistent_tangent_matrix(recalculate_elastic_tangent);
779 }
780
781 auto calcs_stress_matrix = [&]() {
783 auto t_voight_stress_op = voight_to_stress_op();
784 commonDataPtr->stressMatrix.resize(6, commonDataPtr->gradientW.size1(),
785 false);
786 auto t_stess =
787 getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
788 auto t_voight_stress = commonDataPtr->getFTensor1StressAtGaussPts();
789 for (auto gg = 0; gg != commonDataPtr->gradientW.size1(); ++gg) {
790 t_stess(i, j) = t_voight_stress_op(i, j, Z) * t_voight_stress(Z);
791 ++t_voight_stress;
792 ++t_stess;
793 }
795 };
796
797 auto calcs_plastic_strain_matrix = [&]() {
799 auto t_voight_strain_op = voight_to_stress_op();
800 commonDataPtr->plasticStrainMatrix.resize(
801 6, commonDataPtr->activeVariablesW.size1(), false);
802 auto t_plastic_strain =
803 getFTensor2SymmetricFromMat<3>(commonDataPtr->plasticStrainMatrix);
804 auto t_voight_plastic_strain =
805 commonDataPtr->getFTensor1PlasticStrainAtGaussPts();
806 for (auto gg = 0; gg != commonDataPtr->activeVariablesW.size1(); ++gg) {
807 t_plastic_strain(i, j) =
808 t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
809 ++t_voight_plastic_strain;
810 ++t_plastic_strain;
811 }
813 };
814
815 CHKERR calcs_stress_matrix();
816 CHKERR calcs_plastic_strain_matrix();
818 };
819
820 CHKERR do_work_impl(this->commonDataPtr, this->cpPtr);
821
823}
824
825template <>
826ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, LARGE_STRAIN>(
827 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
828 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
829 return new OpCalculateStress<3, LARGE_STRAIN>(m_field, common_data_ptr, cp_ptr,
830 calc_lhs);
831}
832
833template <>
834ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, SMALL_STRAIN>(
835 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
836 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
837 return new OpCalculateStress<3, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr,
838 calc_lhs);
839}
840
841template <>
842ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, LARGE_STRAIN>(
843 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
844 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
845 return new OpCalculateStress<2, LARGE_STRAIN>(m_field, common_data_ptr, cp_ptr,
846 calc_lhs);
847}
848
849template <>
850ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, SMALL_STRAIN>(
851 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
852 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
853 return new OpCalculateStress<2, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr,
854 calc_lhs);
855}
856
857template <int DIM>
860 std::string block_name,
861 boost::shared_ptr<CommonData> common_data_ptr,
862 boost::shared_ptr<ClosestPointProjection> cp_ptr) {
864 pip.push_back(new OpUpdate(common_data_ptr));
866};
867
868template <>
871 std::string block_name,
872 boost::shared_ptr<CommonData> common_data_ptr,
873 boost::shared_ptr<ClosestPointProjection> cp_ptr) {
874 return opFactoryDomainUpdateImpl<3>(m_field, pip, block_name, common_data_ptr,
875 cp_ptr);
876};
877
878template <>
881 std::string block_name,
882 boost::shared_ptr<CommonData> common_data_ptr,
883 boost::shared_ptr<ClosestPointProjection> cp_ptr) {
884 return opFactoryDomainUpdateImpl<2>(m_field, pip, block_name, common_data_ptr,
885 cp_ptr);
886};
887
888
889//! [TSUpdateImpl]
890/**
891 * @brief TS update history variables
892 */
893struct TSUpdateImpl : public TSUpdate {
894
895 TSUpdateImpl(std::string fe_name, boost::shared_ptr<FEMethod> fe_ptr);
897
898private:
899 boost::shared_ptr<FEMethod> fePtr;
900 const std::string feName;
901};
902
903TSUpdateImpl::TSUpdateImpl(std::string fe_name,
904 boost::shared_ptr<FEMethod> fe_ptr)
905 : feName(fe_name), fePtr(fe_ptr) {}
906
909 DM dm;
910 CHKERR TSGetDM(ts, &dm);
911 // This is where update is preformed. Element list is iterated and on each
912 // element pipeline accessed throng fePtr is executed.
915}
916
917/**
918 * @brief Create test update object
919 *
920 * @param fe_name
921 * @param fe_ptr
922 * @return boost::shared_ptr<TSUpdate>
923 */
924boost::shared_ptr<TSUpdate> createTSUpdate(std::string fe_name,
925 boost::shared_ptr<FEMethod> fe_ptr) {
926 return boost::make_shared<TSUpdateImpl>(fe_name, fe_ptr);
927}
928//! [TSUpdateImpl]
929
930} // namespace ADOLCPlasticity
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
constexpr auto t_kd
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
constexpr int DIM2
Definition level_set.cpp:22
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
constexpr int DIM1
Definition level_set.cpp:21
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
FTensor::Dg< double, 3, 6 > strain_to_voight_op()
Op convert strain tensor to Vight strain vector.
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
MoFEMErrorCode opFactoryDomainUpdate< 2 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
static FTensor::Dg< double, 3, 6 > tStrainToVoightOp
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricDiffBaseImpl(DataForcesAndSourcesCore::EntData &data, MatrixDouble &storage, const bool b_bar, const int nb_integration_pts, double *w_ptr, FTensor::Number< DIM >)
[BBar method]
double calcAveStrain(bool b_bar, const int nb_gauss_pts, FTensor::Tensor2< T1, DIM1, DIM2 > &&t_grad, FTensor::Tensor0< T2 > &&t_w)
MoFEMErrorCode opFactoryDomainUpdateImpl(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
FTensor::Tensor1< T2, DIM3 > calcStrain(bool b_bar, double ave_tr_strain, FTensor::Tensor2< T1, DIM1, DIM2 > &t_grad, FTensor::Tensor1< T2, DIM3 > &&t_voight_strain)
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
MoFEMErrorCode opFactoryDomainUpdate< 3 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:31
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricFromPtr< 3 >(double *ptr)
constexpr IntegrationType I
constexpr double t
plate stiffness
Definition plate.cpp:58
FTensor::Index< 'm', 3 > m
static FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data, MatrixDouble &storage, const bool b_bar, const int nb_integration_pts, double *w_ptr, FTensor::Number< 2 >)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
boost::shared_ptr< ClosestPointProjection > cpPtr
OpCalculateStressImpl(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) override
boost::shared_ptr< CommonData > commonDataPtr
OpUpdate(boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode postProcess(TS ts)
boost::shared_ptr< FEMethod > fePtr
TSUpdateImpl(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
Update internal fluxes (update history variables)
Deprecated interface functions.
structure to get information form mofem into EntitiesFieldData