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) {
444 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
445 t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
446 t_P(i, l) = t_F(i, k) * t_S(k, l);
447
448 ++t_grad;
449 ++t_logC;
450 ++t_logC_dC;
451 ++t_P;
452 ++t_T;
453 ++t_S;
454 ++t_D;
455 }
456
458 }
459
460private:
461 boost::shared_ptr<CommonData> commonDataPtr;
462};
463
464template <int DIM, typename DomainEleOp>
465struct OpHenckyTangentImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
467 boost::shared_ptr<CommonData> common_data,
468 boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
470 commonDataPtr(common_data) {
471 std::fill(&DomainEleOp::doEntities[MBEDGE],
472 &DomainEleOp::doEntities[MBMAXTYPE], false);
473 if (mat_D_ptr)
474 matDPtr = mat_D_ptr;
475 else
476 matDPtr = commonDataPtr->matDPtr;
477 }
478
479 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
481
490
491 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
492 // const size_t nb_gauss_pts = matGradPtr->size2();
493 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
494 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
495 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
496 auto dP_dF =
497 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
498
499 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
500 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
501 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
502 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
503 auto t_S =
504 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
505 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
506 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
507
508 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
509
511 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
512
516 eig(i) = t_eig_val(i);
517 eigen_vec(i, j) = t_eig_vec(i, j);
518 T(i, j) = t_T(i, j);
519
520 // rare case when two eigen values are equal
521 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
522
524 dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
525
526 auto TL =
527 EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, T, nb_uniq);
528
529 TL(i, j, k, l) *= 4;
530 FTensor::Ddg<double, DIM, DIM> P_D_P_plus_TL;
531 P_D_P_plus_TL(i, j, k, l) =
532 TL(i, j, k, l) +
533 (t_logC_dC(i, j, o, p) * t_D(o, p, m, n)) * t_logC_dC(m, n, k, l);
534 P_D_P_plus_TL(i, j, k, l) *= 0.5;
535 dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_S(k, j));
536 dP_dF(i, j, m, n) +=
537 t_F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
538
539 ++dP_dF;
540
541 ++t_grad;
542 ++t_eig_val;
543 ++t_eig_vec;
544 ++t_logC_dC;
545 ++t_S;
546 ++t_T;
547 ++t_D;
548 }
549
551 }
552
553private:
554 boost::shared_ptr<CommonData> commonDataPtr;
555 boost::shared_ptr<MatrixDouble> matDPtr;
556};
557
558template <typename DomainEleOp> struct HenkyIntegrators {
559 template <int DIM, IntegrationType I>
561
562 template <int DIM, IntegrationType I>
564
565 template <int DIM, IntegrationType I>
567
568 template <int DIM, IntegrationType I>
571
572 template <int DIM, IntegrationType I>
575
576 template <int DIM, IntegrationType I>
579
580 template <int DIM, IntegrationType I>
582};
583
584template <int DIM>
585MoFEMErrorCode
587 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
588 std::string field_name, std::string block_name,
589 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
590 double scale = 1) {
592
593 struct OpMatBlocks : public DomainEleOp {
594 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
595 double bulk_modulus_K, double shear_modulus_G,
596 MoFEM::Interface &m_field, Sev sev,
597 std::vector<const CubitMeshSets *> meshset_vec_ptr,
598 double scale)
600 bulkModulusKDefault(bulk_modulus_K),
601 shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
602 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
603 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
604 "Can not get data from block");
605 }
606
607 MoFEMErrorCode doWork(int side, EntityType type,
608 EntitiesFieldData::EntData &data) {
610
611 for (auto &b : blockData) {
612
613 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
614 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
615 b.shearModulusG * scaleYoungModulus);
617 }
618 }
619
620 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
621 shearModulusGDefault * scaleYoungModulus);
623 }
624
625 private:
626 boost::shared_ptr<MatrixDouble> matDPtr;
627 const double scaleYoungModulus;
628
629 struct BlockData {
630 double bulkModulusK;
631 double shearModulusG;
632 Range blockEnts;
633 };
634
635 double bulkModulusKDefault;
636 double shearModulusGDefault;
637 std::vector<BlockData> blockData;
638
639 MoFEMErrorCode
640 extractBlockData(MoFEM::Interface &m_field,
641 std::vector<const CubitMeshSets *> meshset_vec_ptr,
642 Sev sev) {
644
645 for (auto m : meshset_vec_ptr) {
646 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
647 std::vector<double> block_data;
648 CHKERR m->getAttributes(block_data);
649 if (block_data.size() != 2) {
650 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
651 "Expected that block has two attribute");
652 }
653 auto get_block_ents = [&]() {
654 Range ents;
655 CHKERR
656 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
657 return ents;
658 };
659
660 double young_modulus = block_data[0];
661 double poisson_ratio = block_data[1];
662 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
663 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
664
665 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
666 << "E = " << young_modulus << " nu = " << poisson_ratio;
667
668 blockData.push_back(
669 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
670 }
671 MOFEM_LOG_CHANNEL("WORLD");
673 }
674
675 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
676 double bulk_modulus_K, double shear_modulus_G) {
678 //! [Calculate elasticity tensor]
679 auto set_material_stiffness = [&]() {
685 double A = (DIM == 2)
686 ? 2 * shear_modulus_G /
687 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
688 : 1;
689 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
690 t_D(i, j, k, l) =
691 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
692 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
693 t_kd(k, l);
694 };
695 //! [Calculate elasticity tensor]
696 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
697 mat_D_ptr->resize(size_symm * size_symm, 1);
698 set_material_stiffness();
700 }
701 };
702
703 double E = young_modulus;
704 double nu = poisson_ratio;
705
706 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
707 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
708 PETSC_NULL);
709 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
710 PETSC_NULL);
711 ierr = PetscOptionsEnd();
712
713 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
714 double shear_modulus_G = E / (2 * (1 + nu));
715 pip.push_back(new OpMatBlocks(
716 field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
717
718 // Get blockset using regular expression
719 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
720
721 (boost::format("%s(.*)") % block_name).str()
722
723 )),
724 scale
725
726 ));
727
729}
730
731template <int DIM, IntegrationType I, typename DomainEleOp>
733 MoFEM::Interface &m_field,
734 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
735 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
736
737 auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
738 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
739 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
740
741 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, field_name, block_name,
742 common_ptr->matDPtr, sev, scale),
743 "addMatBlockOps");
744
746
748 field_name, common_ptr->matGradPtr));
749 pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
750 field_name, common_ptr));
751 pip.push_back(
752 new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
753 pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
754 field_name, common_ptr));
755 pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I>(
756 field_name, common_ptr));
757 pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
758 field_name, common_ptr));
759
760 return common_ptr;
761}
762
763template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
764MoFEMErrorCode opFactoryDomainRhs(
765 MoFEM::Interface &m_field,
766 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
767 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
768 Sev sev) {
770
771 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
772 A>::template LinearForm<I>;
773 using OpInternalForcePiola =
774 typename B::template OpGradTimesTensor<1, DIM, DIM>;
775 pip.push_back(
776 new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
777
779}
780
781template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
782MoFEMErrorCode opFactoryDomainRhs(
783 MoFEM::Interface &m_field,
784 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
785 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
787
788 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
789 m_field, pip, field_name, block_name, sev, scale);
790 CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
791 common_ptr, sev);
792
794}
795
796template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
797MoFEMErrorCode opFactoryDomainLhs(
798 MoFEM::Interface &m_field,
799 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
800 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
801 Sev sev) {
803
804 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
805 A>::template BiLinearForm<I>;
806 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
807
809 pip.push_back(
810 new typename H::template OpHenckyTangent<DIM, I>(field_name, common_ptr));
811 pip.push_back(
812 new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
813
815}
816
817template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
818MoFEMErrorCode opFactoryDomainLhs(
819 MoFEM::Interface &m_field,
820 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
821 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
823
824 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
825 m_field, pip, field_name, block_name, sev, scale);
826 CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
827 common_ptr, sev);
828
830}
831} // namespace HenckyOps
832
833#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
#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:797
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:764
auto d_f
Definition: HenckyOps.hpp:16
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:732
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:28
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
Definition: HenckyOps.hpp:586
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr AssemblyType A
double young_modulus
Young modulus.
Definition: plastic.cpp:168
double scale
Definition: plastic.cpp:166
constexpr auto size_symm
Definition: plastic.cpp:42
double H
Hardening.
Definition: plastic.cpp:171
constexpr auto field_name
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:479
OpHenckyTangentImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
Definition: HenckyOps.hpp:466
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)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.