v0.14.0
Loading...
Searching...
No Matches
PlasticOps.hpp
Go to the documentation of this file.
1
2
3/** \file PlasticOps.hpp
4 * \example PlasticOps.hpp
5
6\f[
7\left\{
8\begin{array}{ll}
9\frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
10\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
11\frac{\partial u_j}{\partial x_i} \right)\\
12\sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
13\dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
14\sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
15f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
16u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
17\sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
18\Omega_u \cup \Omega_\sigma = \Omega \\
19\Omega_u \cap \Omega_\sigma = \emptyset
20\end{array}
21\right.
22\f]
23
24\f[
25\left\{
26\begin{array}{ll}
27\left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
28u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
29\delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
30\dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
31& \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
32\left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
33(f(\pmb\sigma,\tau) - \sigma_y) +
34\| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
35\forall \delta\tau \in L^2(\Omega) \end{array} \right.
36\f]
37
38*/
39
40namespace PlasticOps {
41
42//! [Common data]
43struct CommonData : public boost::enable_shared_from_this<CommonData> {
44
56
57 using BlockParams = std::array<double, LAST_PARAM>;
59
60 inline auto getParamsPtr() {
61 return boost::shared_ptr<BlockParams>(shared_from_this(), &blockParams);
62 };
63
64 //! [Common data set externally]
65 boost::shared_ptr<MatrixDouble> mDPtr;
66 boost::shared_ptr<MatrixDouble> mGradPtr;
67 boost::shared_ptr<MatrixDouble> mStrainPtr;
68 boost::shared_ptr<MatrixDouble> mStressPtr;
69 //! [Common data set externally]
70
71 VectorDouble plasticSurface;
72 MatrixDouble plasticFlow;
73 VectorDouble plasticTau;
74 VectorDouble plasticTauDot;
75 MatrixDouble plasticStrain;
76 MatrixDouble plasticStrainDot;
77
78 VectorDouble resC;
79 VectorDouble resCdTau;
80 MatrixDouble resCdStrain;
81 MatrixDouble resCdPlasticStrain;
82 MatrixDouble resFlow;
83 MatrixDouble resFlowDtau;
84 MatrixDouble resFlowDstrain;
85 MatrixDouble resFlowDstrainDot;
86
87 inline auto getPlasticSurfacePtr() {
88 return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticSurface);
89 }
90 inline auto getPlasticTauPtr() {
91 return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTau);
92 }
93 inline auto getPlasticTauDotPtr() {
94 return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTauDot);
95 }
96 inline auto getPlasticStrainPtr() {
97 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticStrain);
98 }
99 inline auto getPlasticStrainDotPtr() {
100 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
102 }
103 inline auto getPlasticFlowPtr() {
104 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticFlow);
105 }
106
107 static std::array<int, 5> activityData;
108};
109
110std::array<int, 5> CommonData::activityData = {0, 0, 0, 0, 0};
111
112//! [Common data]
113
114
119
120template <int DIM, IntegrationType I, typename DomainEleOp>
122
123template <int DIM, IntegrationType I, typename DomainEleOp>
125
126template <int DIM, IntegrationType I, typename DomainEleOp>
128
129template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
131
132template <IntegrationType I, typename AssemblyDomainEleOp>
134
135template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
137
138template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
140
141template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
143
144template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
146
147template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
149
150template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
152
153template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
155
156template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
158
159template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
161
162template <IntegrationType I, typename AssemblyDomainEleOp>
164
165template <typename DomainEleOp> struct PlasticityIntegrators {
166 template <int DIM, IntegrationType I>
169
170 template <int DIM, IntegrationType I>
172
173 template <int DIM, IntegrationType I>
175
176 template <AssemblyType A> struct Assembly {
177
179 typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
180
181 template <int DIM, IntegrationType I>
184
185 template <IntegrationType I>
188
189 template <int DIM, IntegrationType I>
192
193 template <int DIM, IntegrationType I>
196 DIM, I, AssemblyDomainEleOp>;
197
198 template <int DIM, IntegrationType I>
201
202 template <int DIM, IntegrationType I>
205
206 template <int DIM, IntegrationType I>
209
210 template <int DIM, IntegrationType I>
213
214 template <int DIM, IntegrationType I>
217
218 template <int DIM, IntegrationType I>
221
222 template <int DIM, IntegrationType I>
225
226 template <IntegrationType I>
229
230 };
231};
232
233}; // namespace PlasticOps
234
235#include <PlasticOpsGeneric.hpp>
238#include <PlasticOpsMonitor.hpp>
239
240namespace PlasticOps {
241
242using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
243using CommonPlasticPtr = boost::shared_ptr<PlasticOps::CommonData>;
244using CommonHenkyPtr = boost::shared_ptr<HenckyOps::CommonData>;
245
246template <int DIM>
247MoFEMErrorCode
248addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
249 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
250 boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
251 double scale_value, Sev sev) {
253
254 struct OpMatBlocks : public DomainEleOp {
255 OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
256 boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
257 double scale_value, MoFEM::Interface &m_field, Sev sev,
258 std::vector<const CubitMeshSets *> meshset_vec_ptr)
259 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m_D_ptr),
260 matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
261 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
262 "Can not get data from block");
263 }
264
265 MoFEMErrorCode doWork(int side, EntityType type,
266 EntitiesFieldData::EntData &data) {
268
269 auto getK = [](auto &p) {
270 auto young_modulus = p[CommonData::YOUNG_MODULUS];
271 auto poisson_ratio = p[CommonData::POISSON_RATIO];
272 return young_modulus / (3 * (1 - 2 * poisson_ratio));
273 };
274
275 auto getG = [](auto &p) {
276 auto young_modulus = p[CommonData::YOUNG_MODULUS];
277 auto poisson_ratio = p[CommonData::POISSON_RATIO];
278 return young_modulus / (2 * (1 + poisson_ratio));
279 };
280
281 auto scale_fun = [this](auto &p) {
282 p[CommonData::YOUNG_MODULUS] *= scaleVal;
283 p[CommonData::SIGMA_Y] *= scaleVal;
284 p[CommonData::H] *= scaleVal;
285 p[CommonData::VIS_H] *= scaleVal;
286 p[CommonData::QINF] *= scaleVal;
287 p[CommonData::C1_k] *= scaleVal;
288 };
289
290 for (auto &b : blockData) {
291 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
292 *matParamsPtr = b.bParams;
293 scale_fun(*matParamsPtr);
294 CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
296 }
297 }
298
299 (*matParamsPtr) = {young_modulus, poisson_ratio, sigmaY, H,
300 visH, Qinf, b_iso, C1_k};
301 scale_fun(*matParamsPtr);
302 CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
303
305 }
306
307 private:
308 boost::shared_ptr<MatrixDouble> matDPtr;
309 boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
310 const double scaleVal;
311
312 struct BlockData {
313 std::array<double, CommonData::LAST_PARAM> bParams;
314 Range blockEnts;
315 };
316 std::vector<BlockData> blockData;
317
318 /**
319 * @brief Extract block data from meshsets
320 *
321 * @param m_field
322 * @param meshset_vec_ptr
323 * @param sev
324 * @return MoFEMErrorCode
325 */
326 MoFEMErrorCode
327 extractBlockData(MoFEM::Interface &m_field,
328 std::vector<const CubitMeshSets *> meshset_vec_ptr,
329 Sev sev) {
331
332 for (auto m : meshset_vec_ptr) {
333 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
334 std::vector<double> block_data;
335 CHKERR m->getAttributes(block_data);
336 if (block_data.size() != CommonData::LAST_PARAM) {
337 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338 "Wrong number of block attribute");
339 }
340 auto get_block_ents = [&]() {
341 Range ents;
342 CHKERR m_field.get_moab().get_entities_by_handle(m->meshset, ents,
343 true);
344 return ents;
345 };
346
347 CommonData::BlockParams block_params;
348 for (auto i = 0; i != CommonData::LAST_PARAM; ++i) {
349 block_params[i] = block_data[i];
350 }
351
352 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
353 << "E = " << block_params[CommonData::YOUNG_MODULUS]
354 << " nu = " << block_params[CommonData::POISSON_RATIO];
355 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
356 << std::endl
357 << "sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
358 << "h = " << block_params[CommonData::H] << std::endl
359 << "vis_h = " << block_params[CommonData::VIS_H] << std::endl
360 << "qinf = " << block_params[CommonData::QINF] << std::endl
361 << "biso = " << block_params[CommonData::BISO] << std::endl
362 << "C1_k = " << block_params[CommonData::C1_k] << std::endl;
363
364 blockData.push_back({block_params, get_block_ents()});
365 }
366 MOFEM_LOG_CHANNEL("WORLD");
368 }
369
370 /**
371 * @brief Get elasticity tensor
372 *
373 * Calculate elasticity tensor for given material parameters
374 *
375 * @param mat_D_ptr
376 * @param bulk_modulus_K
377 * @param shear_modulus_G
378 * @return MoFEMErrorCode
379 *
380 */
381 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
382 double bulk_modulus_K, double shear_modulus_G) {
384 //! [Calculate elasticity tensor]
385 auto set_material_stiffness = [&]() {
386 FTensor::Index<'i', DIM> i;
387 FTensor::Index<'j', DIM> j;
388 FTensor::Index<'k', DIM> k;
389 FTensor::Index<'l', DIM> l;
391 double A = (DIM == 2)
392 ? 2 * shear_modulus_G /
393 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
394 : 1;
395 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
396 t_D(i, j, k, l) =
397 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
398 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
399 t_kd(k, l);
400 };
401 //! [Calculate elasticity tensor]
402 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
403 mat_D_ptr->resize(size_symm * size_symm, 1);
404 set_material_stiffness();
406 }
407 };
408
409 // push operator to calculate material stiffness matrix for each block
410 pip.push_back(new OpMatBlocks(
411 mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
412
413 // Get blockset using regular expression
414 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
415
416 (boost::format("%s(.*)") % block_name).str()
417
418 ))
419
420 ));
421
423}
424
425template <int DIM, IntegrationType I, typename DomainEleOp>
427 MoFEM::Interface &m_field, std::string block_name,
428 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
429 std::string u, std::string ep, std::string tau, double scale, Sev sev) {
430
432
433 auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
434 common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
435
436 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
437 auto make_d_mat = []() {
438 return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
439 };
440
441 common_plastic_ptr->mDPtr = make_d_mat();
442 common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
443 common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
444 common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
445
446 auto m_D_ptr = common_plastic_ptr->mDPtr;
447
448 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, block_name, pip, m_D_ptr,
449 common_plastic_ptr->getParamsPtr(),
450 scale, sev),
451 "add mat block ops");
452
453 pip.push_back(new OpCalculateScalarFieldValues(
454 tau, common_plastic_ptr->getPlasticTauPtr()));
455 pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
456 ep, common_plastic_ptr->getPlasticStrainPtr()));
458 u, common_plastic_ptr->mGradPtr));
459
460 CommonHenkyPtr common_henky_ptr;
461
462 if (is_large_strains) {
463 common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
464 common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
465 common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
466 common_henky_ptr->matLogCPlastic =
467 common_plastic_ptr->getPlasticStrainPtr();
468 common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
469 common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
470
472
473 pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
474 u, common_henky_ptr));
475 pip.push_back(
476 new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
477 pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
478 u, common_henky_ptr));
479 pip.push_back(new
480 typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
481 u, common_henky_ptr, m_D_ptr));
482 pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
483 u, common_henky_ptr));
484 } else {
485
486 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
487 common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
488 pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
489 common_plastic_ptr, m_D_ptr));
490 }
491
492 pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
493 u, common_plastic_ptr));
494
495 return std::make_tuple(common_plastic_ptr, common_henky_ptr);
496}
497
498template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
499MoFEMErrorCode
500opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
501 std::string u, std::string ep, std::string tau) {
503
504 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
505 A>::template LinearForm<I>;
507 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
509 typename B::template OpGradTimesTensor<1, DIM, DIM>;
510
512
513 auto [common_plastic_ptr, common_henky_ptr] =
514 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
515 ep, tau, scale, Sev::inform);
516
517 auto m_D_ptr = common_plastic_ptr->mDPtr;
518
519 pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
520 ep, common_plastic_ptr->getPlasticStrainDotPtr()));
521 pip.push_back(new OpCalculateScalarFieldValuesDot(
522 tau, common_plastic_ptr->getPlasticTauDotPtr()));
523 pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
524 u, common_plastic_ptr, m_D_ptr));
525
526 // Calculate internal forces
527 if (common_henky_ptr) {
528 pip.push_back(new OpInternalForcePiola(
529 u, common_henky_ptr->getMatFirstPiolaStress()));
530 } else {
531 pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
532 }
533
534 pip.push_back(
535 new
536 typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
537 tau, common_plastic_ptr, m_D_ptr));
538 pip.push_back(
539 new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
540 DIM, I>(ep, common_plastic_ptr, m_D_ptr));
541
543}
544
545template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
546MoFEMErrorCode
547opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
548 std::string u, std::string ep, std::string tau) {
550
551 using namespace HenckyOps;
552
553 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
554 A>::template BiLinearForm<I>;
555 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
556 using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
557
559
560 auto [common_plastic_ptr, common_henky_ptr] =
561 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
562 ep, tau, scale, Sev::verbose);
563
564 auto m_D_ptr = common_plastic_ptr->mDPtr;
565
566 pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
567 ep, common_plastic_ptr->getPlasticStrainDotPtr()));
568 pip.push_back(new OpCalculateScalarFieldValuesDot(
569 tau, common_plastic_ptr->getPlasticTauDotPtr()));
570 pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
571 u, common_plastic_ptr, m_D_ptr));
572
573 if (common_henky_ptr) {
575 pip.push_back(new typename H::template OpHenckyTangent<DIM, I>(
576 u, common_henky_ptr, m_D_ptr));
577 pip.push_back(new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
578 pip.push_back(
579 new typename P::template Assembly<A>::
580 template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
581 u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
582 } else {
583 pip.push_back(new OpKCauchy(u, u, m_D_ptr));
584 pip.push_back(new typename P::template Assembly<A>::
585 template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
586 u, ep, common_plastic_ptr, m_D_ptr));
587 }
588
589 if (common_henky_ptr) {
590 pip.push_back(
591 new typename P::template Assembly<A>::
592 template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
593 tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
594 pip.push_back(
595 new typename P::template Assembly<A>::
596 template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
597 ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
598 } else {
599 pip.push_back(
600 new
601 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
602 DIM, I>(tau, u, common_plastic_ptr, m_D_ptr));
603 pip.push_back(
604 new
605 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
606 DIM, I>(ep, u, common_plastic_ptr, m_D_ptr));
607 }
608
609 pip.push_back(
610 new
611 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
612 DIM, I>(ep, ep, common_plastic_ptr, m_D_ptr));
613 pip.push_back(
614 new
615 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
616 DIM, I>(ep, tau, common_plastic_ptr, m_D_ptr));
617 pip.push_back(
618 new
619 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
620 DIM, I>(tau, ep, common_plastic_ptr, m_D_ptr));
621 pip.push_back(
622 new
623 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
624 I>(tau, tau, common_plastic_ptr));
625
627}
628
629template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
631 std::string block_name, Pip &pip,
632 std::string u, std::string ep,
633 std::string tau) {
635
636 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
637 A>::template LinearForm<I>;
639 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
641 typename B::template OpGradTimesTensor<1, DIM, DIM>;
642
643 auto [common_plastic_ptr, common_henky_ptr] =
644 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
645 ep, tau, 1., Sev::inform);
646
647 // Calculate internal forces
648 if (common_henky_ptr) {
649 pip.push_back(new OpInternalForcePiola(
650 u, common_henky_ptr->getMatFirstPiolaStress()));
651 } else {
652 pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
653 }
654
656}
657
658} // namespace PlasticOps
static Index< 'p', 3 > p
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double bulk_modulus_K
double shear_modulus_G
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'M', 3 > M
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
FTensor::Index< 'N', 3 > N
FTensor::Index< 'I', 3 > I
[Common data]
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'J', 3 > J
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
double young_modulus
Young modulus.
Definition plastic.cpp:172
double C1_k
Kinematic hardening.
Definition plastic.cpp:180
double Qinf
Saturation yield stress.
Definition plastic.cpp:178
double visH
Viscous hardening.
Definition plastic.cpp:176
double scale
Definition plastic.cpp:170
constexpr auto size_symm
Definition plastic.cpp:42
double H
Hardening.
Definition plastic.cpp:175
double b_iso
Saturation exponent.
Definition plastic.cpp:179
PetscBool is_large_strains
Large strains.
Definition plastic.cpp:167
double sigmaY
Yield stress.
Definition plastic.cpp:174
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition seepage.cpp:63
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition seepage.cpp:65
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
@ OPSPACE
operator do Work is execute on space data
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< MatrixDouble > mStrainPtr
MatrixDouble resFlowDstrainDot
MatrixDouble resCdStrain
boost::shared_ptr< MatrixDouble > mGradPtr
static std::array< int, 5 > activityData
MatrixDouble resCdPlasticStrain
VectorDouble plasticTauDot
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
VectorDouble plasticSurface
[Common data set externally]
MatrixDouble resFlowDtau
MatrixDouble plasticStrain
boost::shared_ptr< MatrixDouble > mStressPtr
MatrixDouble plasticFlow
MatrixDouble plasticStrainDot
VectorDouble plasticTau
MatrixDouble resFlowDstrain
std::array< double, LAST_PARAM > BlockParams
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy