v0.15.0
Loading...
Searching...
No Matches
HenckyOps.hpp
Go to the documentation of this file.
1/**
2 * \file HenckyOps.hpp
3 * \example HenckyOps.hpp
4 *
5 * @copyright Copyright (c) 2023
6 */
7
8#ifndef __HENCKY_OPS_HPP__
9#define __HENCKY_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 {
73 FTensor::Index<'i', DIM> i;
74 FTensor::Index<'j', DIM> j;
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, int S>
124
125template <int DIM, IntegrationType I, typename DomainEleOp, int S>
127
128template <int DIM, IntegrationType I, typename DomainEleOp, int S>
130
131template <int DIM, IntegrationType I, typename DomainEleOp, int S>
133
134template <int DIM, IntegrationType I, typename DomainEleOp, int S>
136
137template <int DIM, IntegrationType I, typename DomainEleOp, int S>
139
140template <int DIM, IntegrationType I, typename DomainEleOp, int S>
142
143template <int DIM, typename DomainEleOp>
145
147 boost::shared_ptr<CommonData> common_data)
149 commonDataPtr(common_data) {
150 std::fill(&DomainEleOp::doEntities[MBEDGE],
151 &DomainEleOp::doEntities[MBMAXTYPE], false);
152 }
153
154 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
156
157 FTensor::Index<'i', DIM> i;
158 FTensor::Index<'j', DIM> j;
159 FTensor::Index<'k', DIM> k;
160
161 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
162 // const size_t nb_gauss_pts = matGradPtr->size2();
163 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
164
165 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
166
167 commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts, false);
168 commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts, false);
169 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
170 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
171
172 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
173
178
179 F(i, j) = t_grad(i, j) + t_kd(i, j);
180 C(i, j) = F(k, i) ^ F(k, j);
181
182 for (int ii = 0; ii != DIM; ii++)
183 for (int jj = 0; jj != DIM; jj++)
184 eigen_vec(ii, jj) = C(ii, jj);
185
186 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
187 for (auto ii = 0; ii != DIM; ++ii)
188 eig(ii) = std::max(std::numeric_limits<double>::epsilon(), eig(ii));
189
190 // rare case when two eigen values are equal
191 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
192 if constexpr (DIM == 3) {
193 if (nb_uniq == 2) {
194 sort_eigen_vals<DIM>(eig, eigen_vec);
195 }
196 }
197
198 t_eig_val(i) = eig(i);
199 t_eig_vec(i, j) = eigen_vec(i, j);
200
201 ++t_grad;
202 ++t_eig_val;
203 ++t_eig_vec;
204 }
205
207 }
208
209private:
210 boost::shared_ptr<CommonData> commonDataPtr;
211};
212
213template <int DIM, typename DomainEleOp>
214struct OpCalculateLogCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
215
217 boost::shared_ptr<CommonData> common_data)
219 commonDataPtr(common_data) {
220 std::fill(&DomainEleOp::doEntities[MBEDGE],
221 &DomainEleOp::doEntities[MBMAXTYPE], false);
222 }
223
224 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
226
227 FTensor::Index<'i', DIM> i;
228 FTensor::Index<'j', DIM> j;
229
230 // const size_t nb_gauss_pts = matGradPtr->size2();
231 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
232 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
233 commonDataPtr->matLogC.resize(size_symm, nb_gauss_pts, false);
234
235 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
236 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
237
238 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
239
240 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
241 t_logC(i, j) = EigenMatrix::getMat(t_eig_val, t_eig_vec, f)(i, j);
242 ++t_eig_val;
243 ++t_eig_vec;
244 ++t_logC;
245 }
246
248 }
249
250private:
251 boost::shared_ptr<CommonData> commonDataPtr;
252};
253
254template <int DIM, typename DomainEleOp>
255struct OpCalculateLogC_dCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
256
258 boost::shared_ptr<CommonData> common_data)
260 commonDataPtr(common_data) {
261 std::fill(&DomainEleOp::doEntities[MBEDGE],
262 &DomainEleOp::doEntities[MBMAXTYPE], false);
263 }
264
265 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
267
268 FTensor::Index<'i', DIM> i;
269 FTensor::Index<'j', DIM> j;
270 FTensor::Index<'k', DIM> k;
271 FTensor::Index<'l', DIM> l;
272
273 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
274 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
275 commonDataPtr->matLogCdC.resize(size_symm * size_symm, nb_gauss_pts, false);
276 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
277 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
278 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
279
280 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
281 // rare case when two eigen values are equal
282 auto nb_uniq = get_uniq_nb<DIM>(&t_eig_val(0));
283 t_logC_dC(i, j, k, l) =
284 2 * EigenMatrix::getDiffMat(t_eig_val, t_eig_vec, f, d_f,
285 nb_uniq)(i, j, k, l);
286
287 ++t_logC_dC;
288 ++t_eig_val;
289 ++t_eig_vec;
290 }
291
293 }
294
295private:
296 boost::shared_ptr<CommonData> commonDataPtr;
297};
298
299template <int DIM, typename DomainEleOp, int S>
301 : public DomainEleOp {
302
304 boost::shared_ptr<CommonData> common_data)
306 commonDataPtr(common_data) {
307 std::fill(&DomainEleOp::doEntities[MBEDGE],
308 &DomainEleOp::doEntities[MBMAXTYPE], false);
309 }
310
311 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
313
314 FTensor::Index<'i', DIM> i;
315 FTensor::Index<'j', DIM> j;
316 FTensor::Index<'k', DIM> k;
317 FTensor::Index<'l', DIM> l;
318
319 // const size_t nb_gauss_pts = matGradPtr->size2();
320 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
321 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
322 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
323 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
324 commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
325 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
326
327 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
328 t_T(i, j) = t_D(i, j, k, l) * t_logC(k, l);
329 ++t_logC;
330 ++t_T;
331 ++t_D;
332 }
333
335 }
336
337private:
338 boost::shared_ptr<CommonData> commonDataPtr;
339};
340
341template <int DIM, typename DomainEleOp, int S>
343 : public DomainEleOp {
344
346 const std::string field_name, boost::shared_ptr<VectorDouble> temperature,
347 boost::shared_ptr<CommonData> common_data,
348 boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
349 boost::shared_ptr<double> ref_temp_ptr)
350 : DomainEleOp(field_name, DomainEleOp::OPROW), tempPtr(temperature),
351 commonDataPtr(common_data), coeffExpansionPtr(coeff_expansion_ptr),
352 refTempPtr(ref_temp_ptr) {
353 std::fill(&DomainEleOp::doEntities[MBEDGE],
354 &DomainEleOp::doEntities[MBMAXTYPE], false);
355 }
356
357 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
359
360 FTensor::Index<'i', DIM> i;
361 FTensor::Index<'j', DIM> j;
362 FTensor::Index<'k', DIM> k;
363 FTensor::Index<'l', DIM> l;
364
365 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
366
367 // const size_t nb_gauss_pts = matGradPtr->size2();
368 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
369 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
370 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
371 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
372 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
373 commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
374 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts, false);
375 commonDataPtr->matSecondPiolaStress.resize(size_symm, nb_gauss_pts, false);
376 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
377 auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
378 auto t_S =
379 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
380 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
381 auto t_temp = getFTensor0FromVec(*tempPtr);
382
384 t_coeff_exp(i, j) = 0;
385 for (auto d = 0; d != SPACE_DIM; ++d) {
386 t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
387 }
388
389 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
390#ifdef HENCKY_SMALL_STRAIN
391 t_P(i, j) = t_D(i, j, k, l) *
392 (t_grad(k, l) - t_coeff_exp(k, l) * (t_temp - (*refTempPtr)));
393#else
394 t_T(i, j) = t_D(i, j, k, l) *
395 (t_logC(k, l) - t_coeff_exp(k, l) * (t_temp - (*refTempPtr)));
397 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
398 t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
399 t_P(i, l) = t_F(i, k) * t_S(k, l);
400#endif
401 ++t_grad;
402 ++t_logC;
403 ++t_logC_dC;
404 ++t_P;
405 ++t_T;
406 ++t_S;
407 ++t_D;
408 ++t_temp;
409 }
410
412 }
413
414private:
415 boost::shared_ptr<CommonData> commonDataPtr;
416 boost::shared_ptr<VectorDouble> tempPtr;
417 boost::shared_ptr<VectorDouble> coeffExpansionPtr;
418 boost::shared_ptr<double> refTempPtr;
419};
420
421template <int DIM, typename DomainEleOp, int S>
423 : public DomainEleOp {
424
426 boost::shared_ptr<CommonData> common_data,
427 boost::shared_ptr<MatrixDouble> mat_D_ptr,
428 const double scale = 1)
429 : DomainEleOp(field_name, DomainEleOp::OPROW), commonDataPtr(common_data),
430 scaleStress(scale), matDPtr(mat_D_ptr) {
431 std::fill(&DomainEleOp::doEntities[MBEDGE],
432 &DomainEleOp::doEntities[MBMAXTYPE], false);
433
434 matLogCPlastic = commonDataPtr->matLogCPlastic;
435 }
436
437 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
439
440 FTensor::Index<'i', DIM> i;
441 FTensor::Index<'j', DIM> j;
442 FTensor::Index<'k', DIM> k;
443 FTensor::Index<'l', DIM> l;
444
445 // const size_t nb_gauss_pts = matGradPtr->size2();
446 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
447 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
448 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
449 auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
450 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
451 commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
452 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
453
454 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
455 t_T(i, j) = t_D(i, j, k, l) * (t_logC(k, l) - t_logCPlastic(k, l));
456 t_T(i, j) /= scaleStress;
457 ++t_logC;
458 ++t_T;
459 ++t_D;
460 ++t_logCPlastic;
461 }
462
464 }
465
466private:
467 boost::shared_ptr<CommonData> commonDataPtr;
468 boost::shared_ptr<MatrixDouble> matDPtr;
469 boost::shared_ptr<MatrixDouble> matLogCPlastic;
470 const double scaleStress;
471};
472
473template <int DIM, typename DomainEleOp, int S>
475 : public DomainEleOp {
476
478 boost::shared_ptr<CommonData> common_data)
480 commonDataPtr(common_data) {
481 std::fill(&DomainEleOp::doEntities[MBEDGE],
482 &DomainEleOp::doEntities[MBMAXTYPE], false);
483 }
484
485 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
487
488 FTensor::Index<'i', DIM> i;
489 FTensor::Index<'j', DIM> j;
490 FTensor::Index<'k', DIM> k;
491 FTensor::Index<'l', DIM> l;
492
493 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
494
495 // const size_t nb_gauss_pts = matGradPtr->size2();
496 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
497#ifdef HENCKY_SMALL_STRAIN
498 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
499#endif
500 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
501 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
502 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
503 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts, false);
504 commonDataPtr->matSecondPiolaStress.resize(size_symm, nb_gauss_pts, false);
505 auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
506 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
507 auto t_S =
508 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
509 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
510
511 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
512
513#ifdef HENCKY_SMALL_STRAIN
514 t_P(i, j) = t_D(i, j, k, l) * t_grad(k, l);
515#else
517 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
518 t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
519 t_P(i, l) = t_F(i, k) * t_S(k, l);
520#endif
521
522 ++t_grad;
523 ++t_logC;
524 ++t_logC_dC;
525 ++t_P;
526 ++t_T;
527 ++t_S;
528#ifdef HENCKY_SMALL_STRAIN
529 ++t_D;
530#endif
531 }
532
534 }
535
536private:
537 boost::shared_ptr<CommonData> commonDataPtr;
538};
539
540template <int DIM, typename DomainEleOp, int S>
541struct OpHenckyTangentImpl<DIM, GAUSS, DomainEleOp, S> : public DomainEleOp {
543 boost::shared_ptr<CommonData> common_data,
544 boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
546 commonDataPtr(common_data) {
547 std::fill(&DomainEleOp::doEntities[MBEDGE],
548 &DomainEleOp::doEntities[MBMAXTYPE], false);
549 if (mat_D_ptr)
550 matDPtr = mat_D_ptr;
551 else
552 matDPtr = commonDataPtr->matDPtr;
553 }
554
555 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
557
558 FTensor::Index<'i', DIM> i;
559 FTensor::Index<'j', DIM> j;
560 FTensor::Index<'k', DIM> k;
561 FTensor::Index<'l', DIM> l;
562 FTensor::Index<'m', DIM> m;
563 FTensor::Index<'n', DIM> n;
564 FTensor::Index<'o', DIM> o;
565 FTensor::Index<'p', DIM> p;
566
567 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
568 // const size_t nb_gauss_pts = matGradPtr->size2();
569 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
570 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
571 auto dP_dF =
572 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
573
574 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
575 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
576 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
577 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
578 auto t_S =
579 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
580 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
581 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
582
583 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
584
585#ifdef HENCKY_SMALL_STRAIN
586 dP_dF(i, j, k, l) = t_D(i, j, k, l);
587#else
588
590 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
591
592 // rare case when two eigen values are equal
593 auto nb_uniq = get_uniq_nb<DIM>(&t_eig_val(0));
595 dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
596
597 auto TL = EigenMatrix::getDiffDiffMat(t_eig_val, t_eig_vec, f, d_f, dd_f,
598 t_T, nb_uniq);
599 TL(i, j, k, l) *= 4;
600 FTensor::Ddg<double, DIM, DIM> P_D_P_plus_TL;
601 P_D_P_plus_TL(i, j, k, l) =
602 TL(i, j, k, l) +
603 (t_logC_dC(i, j, o, p) * t_D(o, p, m, n)) * t_logC_dC(m, n, k, l);
604 P_D_P_plus_TL(i, j, k, l) *= 0.5;
605 dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_S(k, j));
606 dP_dF(i, j, m, n) +=
607 t_F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
608
609#endif
610
611 ++dP_dF;
612
613 ++t_grad;
614 ++t_eig_val;
615 ++t_eig_vec;
616 ++t_logC_dC;
617 ++t_S;
618 ++t_T;
619 ++t_D;
620 }
621
623 }
624
625private:
626 boost::shared_ptr<CommonData> commonDataPtr;
627 boost::shared_ptr<MatrixDouble> matDPtr;
628};
629
630template <int DIM, typename AssemblyDomainEleOp, int S>
632 : public AssemblyDomainEleOp {
634 const std::string row_field_name, const std::string col_field_name,
635 boost::shared_ptr<CommonData> elastic_common_data_ptr,
636 boost::shared_ptr<VectorDouble> coeff_expansion_ptr);
637
638 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
639 EntitiesFieldData::EntData &col_data);
640
641private:
642 boost::shared_ptr<CommonData> elasticCommonDataPtr;
643 boost::shared_ptr<VectorDouble> coeffExpansionPtr;
644};
645
646template <int DIM, typename AssemblyDomainEleOp, int S>
649 const std::string row_field_name, const std::string col_field_name,
650 boost::shared_ptr<CommonData> elastic_common_data_ptr,
651 boost::shared_ptr<VectorDouble> coeff_expansion_ptr)
652 : AssemblyDomainEleOp(row_field_name, col_field_name,
653 AssemblyDomainEleOp::OPROWCOL),
654 elasticCommonDataPtr(elastic_common_data_ptr),
655 coeffExpansionPtr(coeff_expansion_ptr) {
656 this->sYmm = false;
657}
658
659template <int DIM, typename AssemblyDomainEleOp, int S>
660MoFEMErrorCode
662 iNtegrate(EntitiesFieldData::EntData &row_data,
663 EntitiesFieldData::EntData &col_data) {
665
666 auto &locMat = AssemblyDomainEleOp::locMat;
667
668 const auto nb_integration_pts = row_data.getN().size1();
669 const auto nb_row_base_functions = row_data.getN().size2();
670 auto t_w = this->getFTensor0IntegrationWeight();
671
672 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
673 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
674 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*elasticCommonDataPtr->matDPtr);
675 auto t_grad =
676 getFTensor2FromMat<DIM, DIM>(*(elasticCommonDataPtr->matGradPtr));
677 auto t_logC_dC =
678 getFTensor4DdgFromMat<DIM, DIM>(elasticCommonDataPtr->matLogCdC);
681
682 FTensor::Index<'i', DIM> i;
683 FTensor::Index<'j', DIM> j;
684 FTensor::Index<'k', DIM> k;
685 FTensor::Index<'l', DIM> l;
686 FTensor::Index<'m', DIM> m;
687 FTensor::Index<'n', DIM> n;
688 FTensor::Index<'o', DIM> o;
689
691 t_coeff_exp(i, j) = 0;
692 for (auto d = 0; d != SPACE_DIM; ++d) {
693 t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
694 }
695
696 t_eigen_strain(i, j) = (t_D(i, j, k, l) * t_coeff_exp(k, l));
697
698 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
699
700 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
701
702 double alpha = this->getMeasure() * t_w;
703 auto rr = 0;
704 for (; rr != AssemblyDomainEleOp::nbRows / DIM; ++rr) {
705 auto t_mat = getFTensor1FromMat<DIM, 1>(locMat, rr * DIM);
706 auto t_col_base = col_data.getFTensor0N(gg, 0);
707 for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
708#ifdef HENCKY_SMALL_STRAIN
709 t_mat(i) -=
710 (t_row_diff_base(j) * t_eigen_strain(i, j)) * (t_col_base * alpha);
711#else
712 t_mat(i) -= (t_row_diff_base(j) *
713 (t_F(i, o) * ((t_D(m, n, k, l) * t_coeff_exp(k, l)) *
714 t_logC_dC(m, n, o, j)))) *
715 (t_col_base * alpha);
716#endif
717
718 ++t_mat;
719 ++t_col_base;
720 }
721
722 ++t_row_diff_base;
723 }
724 for (; rr != nb_row_base_functions; ++rr)
725 ++t_row_diff_base;
726
727 ++t_w;
728 ++t_grad;
729 ++t_logC_dC;
730 ++t_D;
731 }
732
734}
735
736template <typename DomainEleOp> struct HenckyIntegrators {
737 template <int DIM, IntegrationType I>
739
740 template <int DIM, IntegrationType I>
742
743 template <int DIM, IntegrationType I>
745
746 template <int DIM, IntegrationType I, int S>
749
750 template <int DIM, IntegrationType I, int S>
753
754 template <int DIM, IntegrationType I, int S>
757
758 template <int DIM, IntegrationType I, int S>
761
762 template <int DIM, IntegrationType I, int S>
764
765 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp, int S>
768};
769
770template <int DIM>
771MoFEMErrorCode
773 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
774 std::string block_name,
775 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
776 double scale = 1) {
778
779 PetscBool plane_strain_flag = PETSC_FALSE;
780 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
781 &plane_strain_flag, PETSC_NULLPTR);
782
783 struct OpMatBlocks : public DomainEleOp {
784 OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
785 double shear_modulus_G, MoFEM::Interface &m_field, Sev sev,
786 std::vector<const CubitMeshSets *> meshset_vec_ptr,
787 double scale, PetscBool plane_strain_flag)
788 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
789 bulkModulusKDefault(bulk_modulus_K),
790 shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale),
791 planeStrainFlag(plane_strain_flag) {
792 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
793 "Can not get data from block");
794 }
795
796 MoFEMErrorCode doWork(int side, EntityType type,
797 EntitiesFieldData::EntData &data) {
799
800 for (auto &b : blockData) {
801
802 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
803 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
804 b.shearModulusG * scaleYoungModulus,
805 planeStrainFlag);
807 }
808 }
809
810 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
811 shearModulusGDefault * scaleYoungModulus,
812 planeStrainFlag);
814 }
815
816 private:
817 boost::shared_ptr<MatrixDouble> matDPtr;
818 const double scaleYoungModulus;
819 const PetscBool planeStrainFlag;
820
821 struct BlockData {
822 double bulkModulusK;
823 double shearModulusG;
824 Range blockEnts;
825 };
826
827 double bulkModulusKDefault;
828 double shearModulusGDefault;
829 std::vector<BlockData> blockData;
830
831 MoFEMErrorCode
832 extractBlockData(MoFEM::Interface &m_field,
833 std::vector<const CubitMeshSets *> meshset_vec_ptr,
834 Sev sev) {
836
837 for (auto m : meshset_vec_ptr) {
838 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
839 std::vector<double> block_data;
840 CHKERR m->getAttributes(block_data);
841 if (block_data.size() != 2) {
842 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
843 "Expected that block has two attribute");
844 }
845 auto get_block_ents = [&]() {
846 Range ents;
847 CHKERR
848 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
849 return ents;
850 };
851
852 double young_modulus = block_data[0];
853 double poisson_ratio = block_data[1];
854 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
855 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
856
857 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
858 << "E = " << young_modulus << " nu = " << poisson_ratio;
859
860 blockData.push_back(
861 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
862 }
863 MOFEM_LOG_CHANNEL("WORLD");
865 }
866
867 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
868 double bulk_modulus_K, double shear_modulus_G,
869 PetscBool is_plane_strain) {
871 //! [Calculate elasticity tensor]
872 auto set_material_stiffness = [&]() {
873 FTensor::Index<'i', DIM> i;
874 FTensor::Index<'j', DIM> j;
875 FTensor::Index<'k', DIM> k;
876 FTensor::Index<'l', DIM> l;
878 double A = (SPACE_DIM == 2 && !is_plane_strain)
879 ? 2 * shear_modulus_G /
880 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
881 : 1;
882 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
883 t_D(i, j, k, l) =
884 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
885 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
886 t_kd(k, l);
887 };
888 //! [Calculate elasticity tensor]
889 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
890 mat_D_ptr->resize(size_symm * size_symm, 1);
891 set_material_stiffness();
893 }
894 };
895
896 double E = 1.0;
897 double nu = 0.3;
898
899 PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
900 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
901 PETSC_NULLPTR);
902 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
903 PETSC_NULLPTR);
904 PetscOptionsEnd();
905
906 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
907 double shear_modulus_G = E / (2 * (1 + nu));
908 pip.push_back(new OpMatBlocks(
909 mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
910
911 // Get blockset using regular expression
912 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
913
914 (boost::format("%s(.*)") % block_name).str()
915
916 )),
917 scale, plane_strain_flag
918
919 ));
920
922}
923
924template <int DIM, IntegrationType I, typename DomainEleOp>
926 MoFEM::Interface &m_field,
927 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
928 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
929
930 auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
931 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
932 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
933
934 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
935 common_ptr->matDPtr, sev, scale),
936 "addMatBlockOps");
937
939
940 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
941 field_name, common_ptr->matGradPtr));
942 pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
943 field_name, common_ptr));
944 pip.push_back(
945 new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
946 pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
947 field_name, common_ptr));
948 // Assumes constant D matrix per entity
949 pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I, 0>(
950 field_name, common_ptr));
951 pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
952 field_name, common_ptr));
953
954 return common_ptr;
955}
956
957template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
958MoFEMErrorCode opFactoryDomainRhs(
959 MoFEM::Interface &m_field,
960 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
961 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
962 Sev sev) {
964
965 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
966 A>::template LinearForm<I>;
968 typename B::template OpGradTimesTensor<1, DIM, DIM>;
969 pip.push_back(
970 new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
971
973}
974
975template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
976MoFEMErrorCode opFactoryDomainRhs(
977 MoFEM::Interface &m_field,
978 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
979 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
981
983 m_field, pip, field_name, block_name, sev, scale);
985 common_ptr, sev);
986
988}
989
990template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
991MoFEMErrorCode opFactoryDomainLhs(
992 MoFEM::Interface &m_field,
993 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
994 std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
995 Sev sev) {
997
998 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
999 A>::template BiLinearForm<I>;
1000 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
1001
1003 // Assumes constant D matrix per entity
1004 pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
1005 field_name, common_ptr));
1006 pip.push_back(
1007 new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
1008
1010}
1011
1012template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1013MoFEMErrorCode opFactoryDomainLhs(
1014 MoFEM::Interface &m_field,
1015 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1016 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
1018
1020 m_field, pip, field_name, block_name, sev, scale);
1022 common_ptr, sev);
1023
1025}
1026} // namespace HenckyOps
1027
1028#endif // __HENCKY_OPS_HPP__
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
constexpr double a
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
Kronecker Delta class.
PetscBool is_plane_strain
Definition contact.cpp:95
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double bulk_modulus_K
double shear_modulus_G
@ F
constexpr auto t_kd
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
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)
constexpr double eps
Definition HenckyOps.hpp:13
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)
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)
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)
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
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
double scale
Definition plastic.cpp:123
constexpr auto size_symm
Definition plastic.cpp:42
double H
Hardening.
Definition plastic.cpp:128
constexpr auto field_name
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition seepage.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition seepage.cpp:64
FTensor::Index< 'm', 3 > m
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)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateHenckyPlasticStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr, const double scale=1)
OpCalculateHenckyStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateHenckyThermalStressImpl(const std::string field_name, boost::shared_ptr< VectorDouble > temperature, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< VectorDouble > coeff_expansion_ptr, boost::shared_ptr< double > ref_temp_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateLogCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateLogC_dCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
OpCalculatePiolaStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpHenckyTangentImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.