v0.13.0
HenckyOps.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 namespace HenckyOps {
16 
17 constexpr double eps = std::numeric_limits<float>::epsilon();
18 
19 auto f = [](double v) { return 0.5 * std::log(static_cast<long double>(v)); };
20 auto d_f = [](double v) { return 0.5 / static_cast<long double>(v); };
21 auto dd_f = [](double v) {
22  return -0.5 / (static_cast<long double>(v) * static_cast<long double>(v));
23 };
24 
25 inline auto is_eq(const double &a, const double &b) {
26  return std::abs(a - b) < eps;
27 };
28 
29 template <int DIM> inline auto get_uniq_nb(double *ptr) {
30  std::array<double, DIM> tmp;
31  std::copy(ptr, &ptr[DIM], tmp.begin());
32  std::sort(tmp.begin(), tmp.end());
33  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
34 };
35 
36 template <int DIM>
39  std::array<std::pair<double, size_t>, DIM> tmp;
40  auto is_eq_pair = [](const auto &a, const auto &b) {
41  return a.first < b.first;
42  };
43 
44  for (size_t i = 0; i != DIM; ++i)
45  tmp[i] = std::make_pair(eig(i), i);
46  std::sort(tmp.begin(), tmp.end(), is_eq_pair);
47 
48  int i, j, k;
49  if (is_eq_pair(tmp[0], tmp[1])) {
50  i = tmp[0].second;
51  j = tmp[2].second;
52  k = tmp[1].second;
53  } else {
54  i = tmp[0].second;
55  j = tmp[1].second;
56  k = tmp[2].second;
57  }
58 
60  eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
61 
62  eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
63 
64  eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
65 
66  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
67 
68  {
71  eig(i) = eig_c(i);
72  eigen_vec(i, j) = eigen_vec_c(i, j);
73  }
74 };
75 
76 struct CommonData : public boost::enable_shared_from_this<CommonData> {
77  boost::shared_ptr<MatrixDouble> matGradPtr;
78  boost::shared_ptr<MatrixDouble> matDPtr;
79  boost::shared_ptr<MatrixDouble> matLogCPlastic;
80 
89 
90  inline auto getMatFirstPiolaStress() {
91  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
93  }
94 
95  inline auto getMatHenckyStress() {
96  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
98  }
99 
100  inline auto getMatLogC() {
101  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matLogC);
102  }
103 
104  inline auto getMatTangent() {
105  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matTangent);
106  }
107 };
108 
109 template <int DIM> struct OpCalculateEigenVals : public DomainEleOp {
110 
111  OpCalculateEigenVals(const std::string field_name,
112  boost::shared_ptr<CommonData> common_data)
113  : DomainEleOp(field_name, DomainEleOp::OPROW),
114  commonDataPtr(common_data) {
115  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
116  }
117 
118  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
120 
124 
125  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
126  // const size_t nb_gauss_pts = matGradPtr->size2();
127  const size_t nb_gauss_pts = getGaussPts().size2();
128 
129  auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
130 
131  commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts, false);
132  commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts, false);
133  auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
134  auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
135 
136  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
137 
142 
143  F(i, j) = t_grad(i, j) + t_kd(i, j);
144  C(i, j) = F(k, i) ^ F(k, j);
145 
146  for (int ii = 0; ii != DIM; ii++)
147  for (int jj = 0; jj != DIM; jj++)
148  eigen_vec(ii, jj) = C(ii, jj);
149 
150  CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
151 
152  // rare case when two eigen values are equal
153  auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
154  if (DIM == 3 && nb_uniq == 2)
155  sort_eigen_vals<DIM>(eig, eigen_vec);
156 
157  t_eig_val(i) = eig(i);
158  t_eig_vec(i, j) = eigen_vec(i, j);
159 
160  ++t_grad;
161  ++t_eig_val;
162  ++t_eig_vec;
163  }
164 
166  }
167 
168 private:
169  boost::shared_ptr<CommonData> commonDataPtr;
170 };
171 
172 template <int DIM> struct OpCalculateLogC : public DomainEleOp {
173 
174  OpCalculateLogC(const std::string field_name,
175  boost::shared_ptr<CommonData> common_data)
176  : DomainEleOp(field_name, DomainEleOp::OPROW),
177  commonDataPtr(common_data) {
178  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
179  }
180 
181  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
183 
186 
187  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
188  // const size_t nb_gauss_pts = matGradPtr->size2();
189  const size_t nb_gauss_pts = getGaussPts().size2();
190  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
191  commonDataPtr->matLogC.resize(size_symm, nb_gauss_pts, false);
192 
193  auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
194  auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
195 
196  auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
197 
198  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
199 
200  // rare case when two eigen values are equal
201  auto nb_uniq = get_uniq_nb<DIM>(&(t_eig_val(0)));
202 
205  eig(i) = t_eig_val(i);
206  eigen_vec(i, j) = t_eig_vec(i, j);
207  auto logC = EigenMatrix::getMat(eig, eigen_vec, f);
208  t_logC(i, j) = logC(i, j);
209 
210  ++t_eig_val;
211  ++t_eig_vec;
212  ++t_logC;
213  }
214 
216  }
217 
218 private:
219  boost::shared_ptr<CommonData> commonDataPtr;
220 };
221 
222 template <int DIM> struct OpCalculateLogC_dC : public DomainEleOp {
223 
224  OpCalculateLogC_dC(const std::string field_name,
225  boost::shared_ptr<CommonData> common_data)
226  : DomainEleOp(field_name, DomainEleOp::OPROW),
227  commonDataPtr(common_data) {
228  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
229  }
230 
231  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
233 
238 
239  // const size_t nb_gauss_pts = matGradPtr->size2();
240  const size_t nb_gauss_pts = getGaussPts().size2();
241  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
242  commonDataPtr->matLogCdC.resize(size_symm * size_symm, nb_gauss_pts, false);
243  auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
244  auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
245  auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
246 
247  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
248 
251  eig(i) = t_eig_val(i);
252  eigen_vec(i, j) = t_eig_vec(i, j);
253 
254  // rare case when two eigen values are equal
255  auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
256  auto dlogC_dC = EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
257  dlogC_dC(i, j, k, l) *= 2;
258 
259  t_logC_dC(i, j, k, l) = dlogC_dC(i, j, k, l);
260 
261  ++t_logC_dC;
262  ++t_eig_val;
263  ++t_eig_vec;
264  }
265 
267  }
268 
269 private:
270  boost::shared_ptr<CommonData> commonDataPtr;
271 };
272 
273 template <int DIM> struct OpCalculateHenckyStress : public DomainEleOp {
274 
275  OpCalculateHenckyStress(const std::string field_name,
276  boost::shared_ptr<CommonData> common_data)
277  : DomainEleOp(field_name, DomainEleOp::OPROW),
278  commonDataPtr(common_data) {
279  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
280  }
281 
282  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
284 
289 
290  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
291 
292  // const size_t nb_gauss_pts = matGradPtr->size2();
293  const size_t nb_gauss_pts = getGaussPts().size2();
294  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
295  auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
296  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
297  commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
298  auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
299 
300  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
301  t_T(i, j) = t_D(i, j, k, l) * t_logC(k, l);
302  ++t_logC;
303  ++t_T;
304  ++t_D;
305  }
306 
308  }
309 
310 private:
311  boost::shared_ptr<CommonData> commonDataPtr;
312 };
313 
314 template <int DIM> struct OpCalculateHenckyPlasticStress : public DomainEleOp {
315 
316  OpCalculateHenckyPlasticStress(const std::string field_name,
317  boost::shared_ptr<CommonData> common_data,
318  boost::shared_ptr<MatrixDouble> mat_D_ptr,
319  const double scale = 1)
320  : DomainEleOp(field_name, DomainEleOp::OPROW), commonDataPtr(common_data),
321  scaleStress(scale), matDPtr(mat_D_ptr) {
322  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
323 
324  matLogCPlastic = commonDataPtr->matLogCPlastic;
325  }
326 
327  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
329 
334 
335  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
336 
337  // const size_t nb_gauss_pts = matGradPtr->size2();
338  const size_t nb_gauss_pts = getGaussPts().size2();
339  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
340  auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
341  auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
342  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
343  commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
344  auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
345 
346  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
347  t_T(i, j) = t_D(i, j, k, l) * (t_logC(k, l) - t_logCPlastic(k, l));
348  t_T(i, j) /= scaleStress;
349  ++t_logC;
350  ++t_T;
351  ++t_D;
352  ++t_logCPlastic;
353  }
354 
356  }
357 
358 private:
359  boost::shared_ptr<CommonData> commonDataPtr;
360  boost::shared_ptr<MatrixDouble> matDPtr;
361  boost::shared_ptr<MatrixDouble> matLogCPlastic;
362  const double scaleStress;
363 };
364 
365 template <int DIM> struct OpCalculatePiolaStress : public DomainEleOp {
366 
367  OpCalculatePiolaStress(const std::string field_name,
368  boost::shared_ptr<CommonData> common_data)
369  : DomainEleOp(field_name, DomainEleOp::OPROW),
370  commonDataPtr(common_data) {
371  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
372  }
373 
374  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
376 
381 
382  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
383 
384  // const size_t nb_gauss_pts = matGradPtr->size2();
385  const size_t nb_gauss_pts = getGaussPts().size2();
386  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
387  auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
388  auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
389  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
390  commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts, false);
391  commonDataPtr->matSecondPiolaStress.resize(size_symm, nb_gauss_pts, false);
392  auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
393  auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
394  auto t_S =
395  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
396  auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
397 
398  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
400  t_F(i, j) = t_grad(i, j) + t_kd(i, j);
401  t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
402  t_P(i, l) = t_F(i, k) * t_S(k, l);
403 
404  ++t_grad;
405  ++t_logC;
406  ++t_logC_dC;
407  ++t_P;
408  ++t_T;
409  ++t_S;
410  ++t_D;
411  }
412 
414  }
415 
416 private:
417  boost::shared_ptr<CommonData> commonDataPtr;
418 };
419 
420 template <int DIM> struct OpHenckyTangent : public DomainEleOp {
421  OpHenckyTangent(const std::string field_name,
422  boost::shared_ptr<CommonData> common_data,
423  boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
424  : DomainEleOp(field_name, DomainEleOp::OPROW),
425  commonDataPtr(common_data) {
426  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
427  if (mat_D_ptr)
428  matDPtr = mat_D_ptr;
429  else
430  matDPtr = commonDataPtr->matDPtr;
431  }
432 
433  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
435 
444 
445  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
446  // const size_t nb_gauss_pts = matGradPtr->size2();
447  const size_t nb_gauss_pts = getGaussPts().size2();
448  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
449  commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
450  auto dP_dF =
451  getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
452 
453  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
454  auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
455  auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
456  auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
457  auto t_S =
458  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
459  auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
460  auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
461 
462  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
463 
465  t_F(i, j) = t_grad(i, j) + t_kd(i, j);
466 
470  eig(i) = t_eig_val(i);
471  eigen_vec(i, j) = t_eig_vec(i, j);
472  T(i, j) = t_T(i, j);
473 
474  // rare case when two eigen values are equal
475  auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
476 
478  dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
479 
480  auto TL =
481  EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, T, nb_uniq);
482 
483  TL(i, j, k, l) *= 4;
484  FTensor::Ddg<double, DIM, DIM> P_D_P_plus_TL;
485  P_D_P_plus_TL(i, j, k, l) =
486  TL(i, j, k, l) +
487  (t_logC_dC(i, j, o, p) * t_D(o, p, m, n)) * t_logC_dC(m, n, k, l);
488  P_D_P_plus_TL(i, j, k, l) *= 0.5;
489  dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_S(k, j));
490  dP_dF(i, j, m, n) +=
491  t_F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
492 
493  ++dP_dF;
494 
495  ++t_grad;
496  ++t_eig_val;
497  ++t_eig_vec;
498  ++t_logC_dC;
499  ++t_S;
500  ++t_T;
501  ++t_D;
502  }
503 
505  }
506 
507 private:
508  boost::shared_ptr<CommonData> commonDataPtr;
509  boost::shared_ptr<MatrixDouble> matDPtr;
510 };
511 
512 template <int DIM> struct OpPostProcHencky : public DomainEleOp {
513  OpPostProcHencky(const std::string field_name,
514  moab::Interface &post_proc_mesh,
515  std::vector<EntityHandle> &map_gauss_pts,
516  boost::shared_ptr<CommonData> common_data_ptr);
517  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
518 
519 private:
521  std::vector<EntityHandle> &mapGaussPts;
522  boost::shared_ptr<CommonData> commonDataPtr;
523 };
524 
525 template <int DIM>
527  const std::string field_name, moab::Interface &post_proc_mesh,
528  std::vector<EntityHandle> &map_gauss_pts,
529  boost::shared_ptr<CommonData> common_data_ptr)
530  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
531  mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
532  // Opetor is only executed for vertices
533  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
534 }
535 
536 //! [Postprocessing]
537 template <int DIM>
539  EntData &data) {
541 
545 
546  std::array<double, 9> def;
547  std::fill(def.begin(), def.end(), 0);
548 
549  auto get_tag = [&](const std::string name, size_t size) {
550  Tag th;
551  CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
552  MB_TAG_CREAT | MB_TAG_SPARSE,
553  def.data());
554  return th;
555  };
556 
557  MatrixDouble3by3 mat(3, 3);
558 
559  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
560  mat.clear();
561  for (size_t r = 0; r != SPACE_DIM; ++r)
562  for (size_t c = 0; c != SPACE_DIM; ++c)
563  mat(r, c) = t(r, c);
564  return mat;
565  };
566 
567  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
568  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
569  &*mat.data().begin());
570  };
571 
572  auto th_grad = get_tag("GRAD", 9);
573  auto th_stress = get_tag("STRESS", 9);
574 
575  auto t_piola =
576  getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
577  auto t_grad = getFTensor2FromMat<DIM, DIM>(*commonDataPtr->matGradPtr);
578 
581 
582  for (int gg = 0; gg != commonDataPtr->matGradPtr->size2(); ++gg) {
583 
584  F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
585  double inv_t_detF;
586 
587  if (DIM == 2)
588  determinantTensor2by2(F, inv_t_detF);
589  else
590  determinantTensor3by3(F, inv_t_detF);
591 
592  inv_t_detF = 1. / inv_t_detF;
593 
594  cauchy_stress(i, j) = inv_t_detF * (t_piola(i, k) ^ F(j, k));
595 
596  CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
597  CHKERR set_tag(th_stress, gg, set_matrix(cauchy_stress));
598 
599  ++t_piola;
600  ++t_grad;
601  }
602 
604 }
605 
606 } // namespace HenckyOps
static Index< 'n', 3 > n
static Index< 'o', 3 > o
static Index< 'p', 3 > p
constexpr double a
constexpr int SPACE_DIM
Kronecker Delta class.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
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)
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
auto get_uniq_nb(double *ptr)
Definition: HenckyOps.hpp:29
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: HenckyOps.hpp:37
constexpr double eps
Definition: HenckyOps.hpp:17
auto dd_f
Definition: HenckyOps.hpp:21
auto d_f
Definition: HenckyOps.hpp:20
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:25
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:116
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1202
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1161
const double r
rate factor
double scale
Definition: plastic.cpp:103
constexpr double t
plate stiffness
Definition: plate.cpp:76
FTensor::Index< 'm', 3 > m
auto getMatFirstPiolaStress()
Definition: HenckyOps.hpp:90
MatrixDouble matEigVal
Definition: HenckyOps.hpp:81
MatrixDouble matLogCdC
Definition: HenckyOps.hpp:84
MatrixDouble matEigVec
Definition: HenckyOps.hpp:82
MatrixDouble matTangent
Definition: HenckyOps.hpp:88
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:79
MatrixDouble matFirstPiolaStress
Definition: HenckyOps.hpp:85
MatrixDouble matLogC
Definition: HenckyOps.hpp:83
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:78
MatrixDouble matSecondPiolaStress
Definition: HenckyOps.hpp:86
MatrixDouble matHenckyStress
Definition: HenckyOps.hpp:87
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: HenckyOps.hpp:77
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:169
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:118
OpCalculateEigenVals(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:111
OpCalculateHenckyPlasticStress(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:316
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:360
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:361
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:359
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:327
OpCalculateHenckyStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:275
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:311
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:282
OpCalculateLogC_dC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:224
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:231
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:270
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:219
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:181
OpCalculateLogC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:174
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:374
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:417
OpCalculatePiolaStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:367
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:508
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:433
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:509
OpHenckyTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
Definition: HenckyOps.hpp:421
Definition: HenckyOps.hpp:512
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
Definition: HenckyOps.hpp:538
OpPostProcHencky(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
Definition: HenckyOps.hpp:526
std::vector< EntityHandle > & mapGaussPts
Definition: HenckyOps.hpp:521
moab::Interface & postProcMesh
Definition: HenckyOps.hpp:520
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:522
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element