v0.14.0
Loading...
Searching...
No Matches
HenckyOps.hpp
Go to the documentation of this file.
1/**
2 * \file HenkyOps.hpp
3 * \example HenkyOps.hpp
4 *
5 * @copyright Copyright (c) 2023
6 */
7
8#ifndef __HENKY_OPS_HPP__
9#define __HENKY_OPS_HPP__
10
11namespace HenckyOps {
12
13constexpr double eps = std::numeric_limits<float>::epsilon();
14
15auto f = [](double v) { return 0.5 * std::log(v); };
16auto d_f = [](double v) { return 0.5 / v; };
17auto dd_f = [](double v) { return -0.5 / (v * v); };
18
19struct isEq {
20 static inline auto check(const double &a, const double &b) {
21 return std::abs(a - b) / absMax < eps;
22 }
23 static double absMax;
24};
25
26double isEq::absMax = 1;
27
28inline auto is_eq(const double &a, const double &b) {
29 return isEq::check(a, b);
30};
31
32template <int DIM> inline auto get_uniq_nb(double *ptr) {
33 std::array<double, DIM> tmp;
34 std::copy(ptr, &ptr[DIM], tmp.begin());
35 std::sort(tmp.begin(), tmp.end());
36 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
37 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
38};
39
40template <int DIM>
43
45 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
46
47 int i = 0, j = 1, k = 2;
48
49 if (is_eq(eig(0), eig(1))) {
50 i = 0;
51 j = 2;
52 k = 1;
53 } else if (is_eq(eig(0), eig(2))) {
54 i = 0;
55 j = 1;
56 k = 2;
57 } else if (is_eq(eig(1), eig(2))) {
58 i = 1;
59 j = 0;
60 k = 2;
61 }
62
64 eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
65
66 eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
67
68 eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
69
70 FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
71
72 {
75 eig(i) = eig_c(i);
76 eigen_vec(i, j) = eigen_vec_c(i, j);
77 }
78};
79
80struct CommonData : public boost::enable_shared_from_this<CommonData> {
81 boost::shared_ptr<MatrixDouble> matGradPtr;
82 boost::shared_ptr<MatrixDouble> matDPtr;
83 boost::shared_ptr<MatrixDouble> matLogCPlastic;
84
85 MatrixDouble matEigVal;
86 MatrixDouble matEigVec;
87 MatrixDouble matLogC;
88 MatrixDouble matLogCdC;
89 MatrixDouble matFirstPiolaStress;
91 MatrixDouble matHenckyStress;
92 MatrixDouble matTangent;
93
94 inline auto getMatFirstPiolaStress() {
95 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
97 }
98
99 inline auto getMatHenckyStress() {
100 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
102 }
103
104 inline auto getMatLogC() {
105 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matLogC);
106 }
107
108 inline auto getMatTangent() {
109 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matTangent);
110 }
111};
112
113template <int DIM, IntegrationType I, typename DomainEleOp>
115
116template <int DIM, IntegrationType I, typename DomainEleOp>
118
119template <int DIM, IntegrationType I, typename DomainEleOp>
121
122template <int DIM, IntegrationType I, typename DomainEleOp>
124
125template <int DIM, IntegrationType I, typename DomainEleOp>
127
128template <int DIM, IntegrationType I, typename DomainEleOp>
130
131template <int DIM, IntegrationType I, typename DomainEleOp>
133
134template <int DIM, typename DomainEleOp>
136
138 boost::shared_ptr<CommonData> common_data)
140 commonDataPtr(common_data) {
141 std::fill(&DomainEleOp::doEntities[MBEDGE],
142 &DomainEleOp::doEntities[MBMAXTYPE], false);
143 }
144
145 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
147
151
152 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
153 // const size_t nb_gauss_pts = matGradPtr->size2();
154 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
155
156 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
157
158 commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts, false);
159 commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts, false);
160 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
161 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
162
163 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
164
169
170 F(i, j) = t_grad(i, j) + t_kd(i, j);
171 C(i, j) = F(k, i) ^ F(k, j);
172
173 for (int ii = 0; ii != DIM; ii++)
174 for (int jj = 0; jj != DIM; jj++)
175 eigen_vec(ii, jj) = C(ii, jj);
176
177 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
178 for (auto ii = 0; ii != DIM; ++ii)
179 eig(ii) = std::max(std::numeric_limits<double>::epsilon(), eig(ii));
180
181 // rare case when two eigen values are equal
182 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
183 if constexpr (DIM == 3) {
184 if (nb_uniq == 2) {
185 sort_eigen_vals<DIM>(eig, eigen_vec);
186 }
187 }
188
189 t_eig_val(i) = eig(i);
190 t_eig_vec(i, j) = eigen_vec(i, j);
191
192 ++t_grad;
193 ++t_eig_val;
194 ++t_eig_vec;
195 }
196
198 }
199
200private:
201 boost::shared_ptr<CommonData> commonDataPtr;
202};
203
204template <int DIM, typename DomainEleOp>
205struct OpCalculateLogCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
206
208 boost::shared_ptr<CommonData> common_data)
210 commonDataPtr(common_data) {
211 std::fill(&DomainEleOp::doEntities[MBEDGE],
212 &DomainEleOp::doEntities[MBMAXTYPE], false);
213 }
214
215 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
217
220
221 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
222 // const size_t nb_gauss_pts = matGradPtr->size2();
223 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
224 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
225 commonDataPtr->matLogC.resize(size_symm, nb_gauss_pts, false);
226
227 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
228 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
229
230 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
231
232 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
233
234 // rare case when two eigen values are equal
235 auto nb_uniq = get_uniq_nb<DIM>(&(t_eig_val(0)));
236
239 eig(i) = t_eig_val(i);
240 eigen_vec(i, j) = t_eig_vec(i, j);
241 auto logC = EigenMatrix::getMat(eig, eigen_vec, f);
242 t_logC(i, j) = logC(i, j);
243
244 ++t_eig_val;
245 ++t_eig_vec;
246 ++t_logC;
247 }
248
250 }
251
252private:
253 boost::shared_ptr<CommonData> commonDataPtr;
254};
255
256template <int DIM, typename DomainEleOp>
257struct OpCalculateLogC_dCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
258
260 boost::shared_ptr<CommonData> common_data)
262 commonDataPtr(common_data) {
263 std::fill(&DomainEleOp::doEntities[MBEDGE],
264 &DomainEleOp::doEntities[MBMAXTYPE], false);
265 }
266
267 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
269
274
275 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
276 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
277 commonDataPtr->matLogCdC.resize(size_symm * size_symm, nb_gauss_pts, false);
278 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
279 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
280 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
281
282 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
283
286 eig(i) = t_eig_val(i);
287 eigen_vec(i, j) = t_eig_vec(i, j);
288
289 // rare case when two eigen values are equal
290 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
291 auto dlogC_dC = EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
292 dlogC_dC(i, j, k, l) *= 2;
293
294 t_logC_dC(i, j, k, l) = dlogC_dC(i, j, k, l);
295
296 ++t_logC_dC;
297 ++t_eig_val;
298 ++t_eig_vec;
299 }
300
302 }
303
304private:
305 boost::shared_ptr<CommonData> commonDataPtr;
306};
307
308template <int DIM, typename DomainEleOp>
310 : public DomainEleOp {
311
313 boost::shared_ptr<CommonData> common_data)
315 commonDataPtr(common_data) {
316 std::fill(&DomainEleOp::doEntities[MBEDGE],
317 &DomainEleOp::doEntities[MBMAXTYPE], false);
318 }
319
320 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
322
327
328 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
329
330 // const size_t nb_gauss_pts = matGradPtr->size2();
331 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
332 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
333 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
334 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
335 commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
336 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
337
338 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
339 t_T(i, j) = t_D(i, j, k, l) * t_logC(k, l);
340 ++t_logC;
341 ++t_T;
342 ++t_D;
343 }
344
346 }
347
348private:
349 boost::shared_ptr<CommonData> commonDataPtr;
350};
351
352template <int DIM, typename DomainEleOp>
354 : public DomainEleOp {
355
357 boost::shared_ptr<CommonData> common_data,
358 boost::shared_ptr<MatrixDouble> mat_D_ptr,
359 const double scale = 1)
360 : DomainEleOp(field_name, DomainEleOp::OPROW), commonDataPtr(common_data),
361 scaleStress(scale), matDPtr(mat_D_ptr) {
362 std::fill(&DomainEleOp::doEntities[MBEDGE],
363 &DomainEleOp::doEntities[MBMAXTYPE], false);
364
365 matLogCPlastic = commonDataPtr->matLogCPlastic;
366 }
367
368 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
370
375
376 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
377
378 // const size_t nb_gauss_pts = matGradPtr->size2();
379 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
380 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
381 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
382 auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
383 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
384 commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
385 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
386
387 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
388 t_T(i, j) = t_D(i, j, k, l) * (t_logC(k, l) - t_logCPlastic(k, l));
389 t_T(i, j) /= scaleStress;
390 ++t_logC;
391 ++t_T;
392 ++t_D;
393 ++t_logCPlastic;
394 }
395
397 }
398
399private:
400 boost::shared_ptr<CommonData> commonDataPtr;
401 boost::shared_ptr<MatrixDouble> matDPtr;
402 boost::shared_ptr<MatrixDouble> matLogCPlastic;
403 const double scaleStress;
404};
405
406template <int DIM, typename DomainEleOp>
408 : public DomainEleOp {
409
411 boost::shared_ptr<CommonData> common_data)
413 commonDataPtr(common_data) {
414 std::fill(&DomainEleOp::doEntities[MBEDGE],
415 &DomainEleOp::doEntities[MBMAXTYPE], false);
416 }
417
418 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
420
425
426 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
427
428 // const size_t nb_gauss_pts = matGradPtr->size2();
429 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
430 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
431 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
432 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
433 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
434 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts, false);
435 commonDataPtr->matSecondPiolaStress.resize(size_symm, nb_gauss_pts, false);
436 auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
437 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
438 auto t_S =
439 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
440 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
441
442 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
443
444#ifdef HENCKY_SMALL_STRAIN
445 t_P(i, j) = t_D(i, j, k, l) * t_grad(k, l);
446#else
448 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
449 t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
450 t_P(i, l) = t_F(i, k) * t_S(k, l);
451#endif
452
453 ++t_grad;
454 ++t_logC;
455 ++t_logC_dC;
456 ++t_P;
457 ++t_T;
458 ++t_S;
459 ++t_D;
460 }
461
463 }
464
465private:
466 boost::shared_ptr<CommonData> commonDataPtr;
467};
468
469template <int DIM, typename DomainEleOp>
470struct OpHenckyTangentImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
472 boost::shared_ptr<CommonData> common_data,
473 boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
475 commonDataPtr(common_data) {
476 std::fill(&DomainEleOp::doEntities[MBEDGE],
477 &DomainEleOp::doEntities[MBMAXTYPE], false);
478 if (mat_D_ptr)
479 matDPtr = mat_D_ptr;
480 else
481 matDPtr = commonDataPtr->matDPtr;
482 }
483
484 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
486
495
496 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
497 // const size_t nb_gauss_pts = matGradPtr->size2();
498 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
499 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
500 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
501 auto dP_dF =
502 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
503
504 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
505 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
506 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
507 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
508 auto t_S =
509 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
510 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
511 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
512
513 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
514
515#ifdef HENCKY_SMALL_STRAIN
516 dP_dF(i, j, k, l) = t_D(i, j, k, l);
517#else
518
520 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
521
525 eig(i) = t_eig_val(i);
526 eigen_vec(i, j) = t_eig_vec(i, j);
527 T(i, j) = t_T(i, j);
528
529 // rare case when two eigen values are equal
530 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
531
533 dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
534
535 auto TL =
536 EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, T, nb_uniq);
537
538 TL(i, j, k, l) *= 4;
539 FTensor::Ddg<double, DIM, DIM> P_D_P_plus_TL;
540 P_D_P_plus_TL(i, j, k, l) =
541 TL(i, j, k, l) +
542 (t_logC_dC(i, j, o, p) * t_D(o, p, m, n)) * t_logC_dC(m, n, k, l);
543 P_D_P_plus_TL(i, j, k, l) *= 0.5;
544 dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_S(k, j));
545 dP_dF(i, j, m, n) +=
546 t_F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
547
548#endif
549
550 ++dP_dF;
551
552 ++t_grad;
553 ++t_eig_val;
554 ++t_eig_vec;
555 ++t_logC_dC;
556 ++t_S;
557 ++t_T;
558 ++t_D;
559 }
560
562 }
563
564private:
565 boost::shared_ptr<CommonData> commonDataPtr;
566 boost::shared_ptr<MatrixDouble> matDPtr;
567};
568
569template <typename DomainEleOp> struct HenkyIntegrators {
570 template <int DIM, IntegrationType I>
572
573 template <int DIM, IntegrationType I>
575
576 template <int DIM, IntegrationType I>
578
579 template <int DIM, IntegrationType I>
582
583 template <int DIM, IntegrationType I>
586
587 template <int DIM, IntegrationType I>
590
591 template <int DIM, IntegrationType I>
593};
594
595template <int DIM>
596MoFEMErrorCode
598 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
599 std::string block_name,
600 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
601 double scale = 1) {
603
604 struct OpMatBlocks : public DomainEleOp {
605 OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
606 double shear_modulus_G, MoFEM::Interface &m_field, Sev sev,
607 std::vector<const CubitMeshSets *> meshset_vec_ptr,
608 double scale)
610 bulkModulusKDefault(bulk_modulus_K),
611 shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
612 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
613 "Can not get data from block");
614 }
615
616 MoFEMErrorCode doWork(int side, EntityType type,
617 EntitiesFieldData::EntData &data) {
619
620 for (auto &b : blockData) {
621
622 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
623 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
624 b.shearModulusG * scaleYoungModulus);
626 }
627 }
628
629 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
630 shearModulusGDefault * scaleYoungModulus);
632 }
633
634 private:
635 boost::shared_ptr<MatrixDouble> matDPtr;
636 const double scaleYoungModulus;
637
638 struct BlockData {
639 double bulkModulusK;
640 double shearModulusG;
641 Range blockEnts;
642 };
643
644 double bulkModulusKDefault;
645 double shearModulusGDefault;
646 std::vector<BlockData> blockData;
647
648 MoFEMErrorCode
649 extractBlockData(MoFEM::Interface &m_field,
650 std::vector<const CubitMeshSets *> meshset_vec_ptr,
651 Sev sev) {
653
654 for (auto m : meshset_vec_ptr) {
655 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
656 std::vector<double> block_data;
657 CHKERR m->getAttributes(block_data);
658 if (block_data.size() != 2) {
659 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
660 "Expected that block has two attribute");
661 }
662 auto get_block_ents = [&]() {
663 Range ents;
664 CHKERR
665 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
666 return ents;
667 };
668
669 double young_modulus = block_data[0];
670 double poisson_ratio = block_data[1];
671 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
672 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
673
674 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
675 << "E = " << young_modulus << " nu = " << poisson_ratio;
676
677 blockData.push_back(
678 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
679 }
680 MOFEM_LOG_CHANNEL("WORLD");
682 }
683
684 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
685 double bulk_modulus_K, double shear_modulus_G) {
687 //! [Calculate elasticity tensor]
688 auto set_material_stiffness = [&]() {
694 double A = (DIM == 2)
695 ? 2 * shear_modulus_G /
696 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
697 : 1;
698 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
699 t_D(i, j, k, l) =
700 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
701 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
702 t_kd(k, l);
703 };
704 //! [Calculate elasticity tensor]
705 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
706 mat_D_ptr->resize(size_symm * size_symm, 1);
707 set_material_stiffness();
709 }
710 };
711
712 double E = young_modulus;
713 double nu = poisson_ratio;
714
715 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
716 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
717 PETSC_NULL);
718 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
719 PETSC_NULL);
720 ierr = PetscOptionsEnd();
721
722 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
723 double shear_modulus_G = E / (2 * (1 + nu));
724 pip.push_back(new OpMatBlocks(
725 mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
726
727 // Get blockset using regular expression
728 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
729
730 (boost::format("%s(.*)") % block_name).str()
731
732 )),
733 scale
734
735 ));
736
738}
739
740template <int DIM, IntegrationType I, typename DomainEleOp>
742 MoFEM::Interface &m_field,
743 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
744 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
745
746 auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
747 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
748 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
749
750 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
751 common_ptr->matDPtr, sev, scale),
752 "addMatBlockOps");
753
755
757 field_name, common_ptr->matGradPtr));
758 pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
759 field_name, common_ptr));
760 pip.push_back(
761 new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
762 pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
763 field_name, common_ptr));
764 pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I>(
765 field_name, common_ptr));
766 pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
767 field_name, common_ptr));
768
769 return common_ptr;
770}
771
772template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
773MoFEMErrorCode opFactoryDomainRhs(
774 MoFEM::Interface &m_field,
775 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
776 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
777 Sev sev) {
779
780 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
781 A>::template LinearForm<I>;
783 typename B::template OpGradTimesTensor<1, DIM, DIM>;
784 pip.push_back(
785 new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
786
788}
789
790template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
791MoFEMErrorCode opFactoryDomainRhs(
792 MoFEM::Interface &m_field,
793 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
794 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
796
797 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
798 m_field, pip, field_name, block_name, sev, scale);
799 CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
800 common_ptr, sev);
801
803}
804
805template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
806MoFEMErrorCode opFactoryDomainLhs(
807 MoFEM::Interface &m_field,
808 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
809 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
810 Sev sev) {
812
813 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
814 A>::template BiLinearForm<I>;
815 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
816
818 pip.push_back(
819 new typename H::template OpHenckyTangent<DIM, I>(field_name, common_ptr));
820 pip.push_back(
821 new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
822
824}
825
826template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
827MoFEMErrorCode opFactoryDomainLhs(
828 MoFEM::Interface &m_field,
829 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
830 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
832
833 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
834 m_field, pip, field_name, block_name, sev, scale);
835 CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
836 common_ptr, sev);
837
839}
840} // namespace HenckyOps
841
842#endif // __HENKY_OPS_HPP__
static Index< 'o', 3 > o
static Index< 'p', 3 > p
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
constexpr double a
static PetscErrorCode ierr
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
Kronecker Delta class.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
double bulk_modulus_K
double shear_modulus_G
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
constexpr auto t_kd
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto get_uniq_nb(double *ptr)
Definition: HenckyOps.hpp:32
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: HenckyOps.hpp:41
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
Definition: HenckyOps.hpp:806
constexpr double eps
Definition: HenckyOps.hpp:13
auto dd_f
Definition: HenckyOps.hpp:17
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
Definition: HenckyOps.hpp:773
auto d_f
Definition: HenckyOps.hpp:16
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
Definition: HenckyOps.hpp:597
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
Definition: HenckyOps.hpp:741
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:28
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr AssemblyType A
double young_modulus
Young modulus.
Definition: plastic.cpp:172
double scale
Definition: plastic.cpp:170
constexpr auto size_symm
Definition: plastic.cpp:42
double H
Hardening.
Definition: plastic.cpp:175
constexpr auto field_name
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
auto getMatFirstPiolaStress()
Definition: HenckyOps.hpp:94
MatrixDouble matEigVal
Definition: HenckyOps.hpp:85
MatrixDouble matLogCdC
Definition: HenckyOps.hpp:88
MatrixDouble matEigVec
Definition: HenckyOps.hpp:86
MatrixDouble matTangent
Definition: HenckyOps.hpp:92
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:83
MatrixDouble matFirstPiolaStress
Definition: HenckyOps.hpp:89
MatrixDouble matLogC
Definition: HenckyOps.hpp:87
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:82
MatrixDouble matSecondPiolaStress
Definition: HenckyOps.hpp:90
MatrixDouble matHenckyStress
Definition: HenckyOps.hpp:91
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: HenckyOps.hpp:81
OpCalculateEigenValsImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:137
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:145
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:368
OpCalculateHenckyPlasticStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr, const double scale=1)
Definition: HenckyOps.hpp:356
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:320
OpCalculateHenckyStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:312
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:267
OpCalculateLogC_dCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:259
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:215
OpCalculateLogCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:207
OpCalculatePiolaStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:410
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:418
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:484
OpHenckyTangentImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
Definition: HenckyOps.hpp:471
static double absMax
Definition: HenckyOps.hpp:23
static auto check(const double &a, const double &b)
Definition: HenckyOps.hpp:20
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.