v0.14.0
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 
11 namespace HenckyOps {
12 
13 constexpr double eps = std::numeric_limits<float>::epsilon();
14 
15 auto f = [](double v) { return 0.5 * std::log(v); };
16 auto d_f = [](double v) { return 0.5 / v; };
17 auto dd_f = [](double v) { return -0.5 / (v * v); };
18 
19 struct 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 
26 double isEq::absMax = 1;
27 
28 inline auto is_eq(const double &a, const double &b) {
29  return isEq::check(a, b);
30 };
31 
32 template <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 
40 template <int DIM>
43 
44  isEq::absMax =
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 
80 struct 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 
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(),
101  &matHenckyStress);
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 
113 template <int DIM, IntegrationType I, typename DomainEleOp>
115 
116 template <int DIM, IntegrationType I, typename DomainEleOp>
118 
119 template <int DIM, IntegrationType I, typename DomainEleOp>
121 
122 template <int DIM, IntegrationType I, typename DomainEleOp>
124 
125 template <int DIM, IntegrationType I, typename DomainEleOp>
127 
128 template <int DIM, IntegrationType I, typename DomainEleOp>
130 
131 template <int DIM, IntegrationType I, typename DomainEleOp>
133 
134 template <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 
200 private:
201  boost::shared_ptr<CommonData> commonDataPtr;
202 };
203 
204 template <int DIM, typename DomainEleOp>
206 
207  OpCalculateLogCImpl(const std::string field_name,
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 
252 private:
253  boost::shared_ptr<CommonData> commonDataPtr;
254 };
255 
256 template <int DIM, typename 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 
304 private:
305  boost::shared_ptr<CommonData> commonDataPtr;
306 };
307 
308 template <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 
348 private:
349  boost::shared_ptr<CommonData> commonDataPtr;
350 };
351 
352 template <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 
399 private:
400  boost::shared_ptr<CommonData> commonDataPtr;
401  boost::shared_ptr<MatrixDouble> matDPtr;
402  boost::shared_ptr<MatrixDouble> matLogCPlastic;
403  const double scaleStress;
404 };
405 
406 template <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 
465 private:
466  boost::shared_ptr<CommonData> commonDataPtr;
467 };
468 
469 template <int DIM, typename DomainEleOp>
471  OpHenckyTangentImpl(const std::string field_name,
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 
564 private:
565  boost::shared_ptr<CommonData> commonDataPtr;
566  boost::shared_ptr<MatrixDouble> matDPtr;
567 };
568 
569 template <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>
588  using OpCalculatePiolaStress =
590 
591  template <int DIM, IntegrationType I>
593 };
594 
595 template <int DIM>
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)
609  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
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,
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 
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 
740 template <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 
772 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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>;
782  using OpInternalForcePiola =
783  typename B::template OpGradTimesTensor<1, DIM, DIM>;
784  pip.push_back(
785  new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
786 
788 }
789 
790 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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 
805 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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 
826 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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__
NOSPACE
@ NOSPACE
Definition: definitions.h:83
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
HenckyOps::OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:565
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:400
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
HenckyOps::CommonData::matLogCPlastic
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:83
HenckyOps::OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >::OpCalculatePiolaStressImpl
OpCalculatePiolaStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:410
EshelbianPlasticity::size_symm
constexpr static auto size_symm
Definition: EshelbianOperators.cpp:47
HenckyOps::CommonData::getMatFirstPiolaStress
auto getMatFirstPiolaStress()
Definition: HenckyOps.hpp:94
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
HenckyOps::OpCalculateHenckyStressImpl
Definition: HenckyOps.hpp:123
HenckyOps::OpCalculateHenckyStressImpl< DIM, GAUSS, DomainEleOp >::OpCalculateHenckyStressImpl
OpCalculateHenckyStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:312
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
HenckyOps::isEq::absMax
static double absMax
Definition: HenckyOps.hpp:23
HenckyOps::CommonData::matEigVec
MatrixDouble matEigVec
Definition: HenckyOps.hpp:86
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HenckyOps::OpCalculateLogC_dCImpl< DIM, GAUSS, DomainEleOp >::OpCalculateLogC_dCImpl
OpCalculateLogC_dCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:259
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
HenckyOps::CommonData::getMatLogC
auto getMatLogC()
Definition: HenckyOps.hpp:104
HenckyOps::opFactoryDomainLhs
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
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
HenckyOps::OpCalculateLogCImpl< DIM, GAUSS, DomainEleOp >::OpCalculateLogCImpl
OpCalculateLogCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:207
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::OpCalculateEigenValsImpl
OpCalculateEigenValsImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:137
HenckyOps::CommonData
Definition: HenckyOps.hpp:80
HenckyOps::CommonData::matEigVal
MatrixDouble matEigVal
Definition: HenckyOps.hpp:85
HenckyOps
Definition: HenckyOps.hpp:11
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
HenckyOps::CommonData::matLogC
MatrixDouble matLogC
Definition: HenckyOps.hpp:87
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
HenckyOps::OpCalculateHenckyStressImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:349
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
HenckyOps::OpCalculateHenckyPlasticStressImpl
Definition: HenckyOps.hpp:126
HenckyOps::OpCalculatePiolaStressImpl
Definition: HenckyOps.hpp:129
HenckyOps::isEq::check
static auto check(const double &a, const double &b)
Definition: HenckyOps.hpp:20
HenckyOps::is_eq
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:28
HenckyOps::OpCalculateLogC_dCImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:305
HenckyOps::OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:418
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
HenckyOps::OpCalculateLogC_dCImpl
Definition: HenckyOps.hpp:120
HenckyOps::OpHenckyTangentImpl
Definition: HenckyOps.hpp:132
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
HenckyOps::OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:466
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
HenckyOps::OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >
Definition: HenckyOps.hpp:470
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition: child_and_parent.cpp:36
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
HenckyOps::OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >::OpHenckyTangentImpl
OpHenckyTangentImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
Definition: HenckyOps.hpp:471
BlockData
Definition: gel_analysis.cpp:74
a
constexpr double a
Definition: approx_sphere.cpp:30
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:401
convert.type
type
Definition: convert.py:64
HenckyOps::OpCalculateLogCImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:215
HenckyOps::OpCalculateLogCImpl
Definition: HenckyOps.hpp:117
HenckyOps::eps
constexpr double eps
Definition: HenckyOps.hpp:13
HenckyOps::OpCalculateHenckyStressImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:320
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
HenckyOps::CommonData::getMatHenckyStress
auto getMatHenckyStress()
Definition: HenckyOps.hpp:99
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::matLogCPlastic
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:402
HenckyOps::CommonData::matFirstPiolaStress
MatrixDouble matFirstPiolaStress
Definition: HenckyOps.hpp:89
HenckyOps::HenkyIntegrators
Definition: HenckyOps.hpp:569
HenckyOps::CommonData::matSecondPiolaStress
MatrixDouble matSecondPiolaStress
Definition: HenckyOps.hpp:90
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::OpCalculateHenckyPlasticStressImpl
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
HenckyOps::OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:484
HenckyOps::OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >
Definition: HenckyOps.hpp:407
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:145
HenckyOps::OpCalculateEigenValsImpl
Definition: HenckyOps.hpp:114
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
HenckyOps::isEq
Definition: HenckyOps.hpp:19
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', DIM >
convert.n
n
Definition: convert.py:82
HenckyOps::CommonData::matLogCdC
MatrixDouble matLogCdC
Definition: HenckyOps.hpp:88
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
DomainEleOp
HenckyOps::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: HenckyOps.hpp:41
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
HenckyOps::CommonData::matGradPtr
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: HenckyOps.hpp:81
HenckyOps::CommonData::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:82
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
HenckyOps::CommonData::getMatTangent
auto getMatTangent()
Definition: HenckyOps.hpp:108
HenckyOps::OpCalculateLogCImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:253
HenckyOps::opFactoryDomainRhs
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
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:368
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EigenMatrix::getDiffDiffMat
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)
Definition: MatrixFunction.cpp:78
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:201
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
HenckyOps::OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >::scaleStress
const double scaleStress
Definition: HenckyOps.hpp:403
HenckyOps::OpCalculateLogC_dCImpl< DIM, GAUSS, DomainEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: HenckyOps.hpp:267
HenckyOps::commonDataFactory
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
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
HenckyOps::addMatBlockOps
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
HenckyOps::OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:566
HenckyOps::CommonData::matTangent
MatrixDouble matTangent
Definition: HenckyOps.hpp:92
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EigenMatrix::getDiffMat
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.
Definition: MatrixFunction.cpp:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HenckyOps::CommonData::matHenckyStress
MatrixDouble matHenckyStress
Definition: HenckyOps.hpp:91
HenckyOps::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: HenckyOps.hpp:32
H
double H
Hardening.
Definition: plastic.cpp:124
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1430
F
@ F
Definition: free_surface.cpp:394