v0.14.0
Loading...
Searching...
No Matches
EshelbianOperators.cpp
Go to the documentation of this file.
1/**
2 * \file EshelbianOperators.cpp
3 * \example EshelbianOperators.cpp
4 *
5 * \brief Implementation of operators
6 */
7
8#include <MoFEM.hpp>
9using namespace MoFEM;
10
12#include <boost/math/constants/constants.hpp>
13
14#include <lapack_wrap.h>
15
16#include <MatrixFunction.hpp>
17
18namespace EshelbianPlasticity {
19
20constexpr static auto size_symm = (3 * (3 + 1)) / 2;
21
22inline auto diff_tensor() {
29 t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
30 return t_diff;
31};
32
33inline auto symm_L_tensor() {
38 t_L(i, j, L) = 0;
39 t_L(0, 0, 0) = 1;
40 t_L(1, 1, 3) = 1;
41 t_L(2, 2, 5) = 1;
42 t_L(1, 0, 1) = 1;
43 t_L(2, 0, 2) = 1;
44 t_L(2, 1, 4) = 1;
45 return t_L;
46}
47
48template <class T> struct TensorTypeExtractor {
49 typedef typename std::remove_pointer<T>::type Type;
50};
51template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
52 typedef typename std::remove_pointer<T>::type Type;
53};
54
55template <typename T>
57 RotSelector rotSelector = LARGE_ROT) {
58
59 using D = typename TensorTypeExtractor<T>::Type;
60
65
66 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
67 t_R(i, j) = t_kd(i, j);
68
70 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
71
72 const auto angle =
73 sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
74 if (rotSelector == SMALL_ROT) {
75 t_R(i, j) += t_Omega(i, j);
76 return t_R;
77 }
78
79 const auto a = sin(angle) / angle;
80 t_R(i, j) += a * t_Omega(i, j);
81 if (rotSelector == MODERATE_ROT)
82 return t_R;
83
84 const auto ss_2 = sin(angle / 2.);
85 const auto b = 2. * ss_2 * ss_2 / (angle * angle);
86 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
87
88 return t_R;
89}
90
91template <typename T>
93 RotSelector rotSelector = LARGE_ROT) {
94
95 using D = typename TensorTypeExtractor<T>::Type;
96
101
103
104 const auto angle =
105 sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
106 if (rotSelector == SMALL_ROT) {
107 t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
108 return t_diff_R;
109 }
110
111 const auto ss = sin(angle);
112 const auto a = ss / angle;
113 t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
114
116 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
117
118 const auto angle2 = angle * angle;
119 const auto cc = cos(angle);
120 const auto diff_a = (angle * cc - ss) / angle2;
121 t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
122 if (rotSelector == MODERATE_ROT)
123 return t_diff_R;
124
125 const auto ss_2 = sin(angle / 2.);
126 const auto cc_2 = cos(angle / 2.);
127 const auto b = 2. * ss_2 * ss_2 / angle2;
128 t_diff_R(i, j, k) +=
129 b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
130 FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
131 const auto diff_b =
132 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
133 t_diff_R(i, j, k) +=
134 diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
135
136 return t_diff_R;
137}
138
139template <typename T>
141 RotSelector rotSelector = LARGE_ROT) {
142
143 using D = typename TensorTypeExtractor<T>::Type;
144
147
148 constexpr double eps = 1e-10;
149 for (int l = 0; l != 3; ++l) {
151 t_omega_c(i) = t_omega(i);
152 t_omega_c(l) += std::complex<double>(0, eps);
153 auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
154 for (int i = 0; i != 3; ++i) {
155 for (int j = 0; j != 3; ++j) {
156 for (int k = 0; k != 3; ++k) {
157 t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
158 }
159 }
160 }
161 }
162
163 return t_diff2_R;
164}
165
166struct isEq {
167 static inline auto check(const double &a, const double &b) {
168 constexpr double eps = std::numeric_limits<float>::epsilon();
169 return std::abs(a - b) / absMax < eps;
170 }
171 static double absMax;
172};
173
174double isEq::absMax = 1;
175
176inline auto is_eq(const double &a, const double &b) {
177 return isEq::check(a, b);
178};
179
180template <int DIM> inline auto get_uniq_nb(double *ptr) {
181 std::array<double, DIM> tmp;
182 std::copy(ptr, &ptr[DIM], tmp.begin());
183 std::sort(tmp.begin(), tmp.end());
184 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
185 constexpr double eps = std::numeric_limits<float>::epsilon();
186 isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
187 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
188}
189
190template <int DIM>
193
195 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
196
197 int i = 0, j = 1, k = 2;
198
199 if (is_eq(eig(0), eig(1))) {
200 i = 0;
201 j = 2;
202 k = 1;
203 } else if (is_eq(eig(0), eig(2))) {
204 i = 0;
205 j = 1;
206 k = 2;
207 } else if (is_eq(eig(1), eig(2))) {
208 i = 1;
209 j = 0;
210 k = 2;
211 }
212
214 eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
215
216 eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
217
218 eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
219
220 FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
221
222 {
225 eig(i) = eig_c(i);
226 eigen_vec(i, j) = eigen_vec_c(i, j);
227 }
228}
229
231 EntityType type,
232 EntData &data) {
234
236 auto t_L = symm_L_tensor();
237
238 int nb_integration_pts = getGaussPts().size2();
245
246 dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
247 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
248 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
249 dataAtPts->leviKirchoffAtPts.resize(3, nb_integration_pts, false);
250 dataAtPts->leviKirchoffdPAtPts.resize(9 * 3, nb_integration_pts, false);
251 dataAtPts->leviKirchoffdOmegaAtPts.resize(9, nb_integration_pts, false);
252
253 dataAtPts->adjontPdstretchAtPts.resize(9, nb_integration_pts, false);
254 dataAtPts->adjontPdUAtPts.resize(size_symm, nb_integration_pts, false);
255 dataAtPts->adjontPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
256 dataAtPts->adjontPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
257 false);
258 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
259 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
260 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
261 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
262 dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
263 dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
264 dataAtPts->nbUniq.resize(nb_integration_pts, false);
265
266 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
267
268 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
269 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
270 auto t_h_dlog_u = getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
271 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchoffAtPts);
272 auto t_levi_kirchoff_dP =
273 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchoffdPAtPts);
274 auto t_levi_kirchoff_domega =
275 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchoffdOmegaAtPts);
276 auto t_approx_P_adjont_dstretch =
277 getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
278 auto t_approx_P_adjont_log_du =
279 getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
280 auto t_approx_P_adjont_log_du_dP =
281 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
282 auto t_approx_P_adjont_log_du_domega =
283 getFTensor2FromMat<3, size_symm>(dataAtPts->adjontPdUdOmegaAtPts);
284 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
285 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
286 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
287 auto t_log_u =
288 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
289 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
290 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
291 auto t_diff_u =
292 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
293 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
294 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
295
296 auto &nbUniq = dataAtPts->nbUniq;
297 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
298 auto t_log_stretch_total =
299 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
301
302 for (int gg = 0; gg != nb_integration_pts; ++gg) {
303
309 FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiata;
310
311 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
312
313 auto calulate_rotation = [&]() {
314 auto t0_diff =
317 t_diff_R(i, j, k) = t0_diff(i, j, k);
318 t_R(i, j) = t0(i, j);
319 };
320
321 auto calulate_streach = [&]() {
322 eigen_vec(i, j) = t_log_u(i, j);
323 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
324 // rare case when two eigen values are equal
325 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
326 if (nbUniq[gg] < 3) {
327 sort_eigen_vals<3>(eig, eigen_vec);
328 }
329 t_eigen_vals(i) = eig(i);
330 t_eigen_vecs(i, j) = eigen_vec(i, j);
331 auto t_u_tmp =
332 EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
333 t_u(i, j) = t_u_tmp(i, j);
334 auto t_diff_u_tmp =
335 EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f,
336 EshelbianCore::d_f, nbUniq[gg]);
337 t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
338 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
339 };
340
341 calulate_rotation();
342 calulate_streach();
343
345 case LARGE_ROT:
346
347 t_Ru(i, m) = t_R(i, l) * t_u(l, m);
348 t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
349 t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
350 t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
351
352 t_approx_P_intermidiata(i, m) = t_approx_P(i, j) * t_h1(m, j);
353 t_approx_P_adjont_dstretch(l, m) =
354 t_approx_P_intermidiata(i, m) * t_R(i, l);
355
356 t_levi_kirchoff(k) =
357 levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
358 t_levi_kirchoff_dP(i, j, k) =
359 (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
360 t_levi_kirchoff_domega(k, n) =
361 levi_civita(l, m, k) *
362 (t_approx_P_intermidiata(i, m) * t_diff_R(i, l, n));
363
364 t_approx_P_adjont_log_du(L) =
365 t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
366 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
367 t_approx_P_adjont_log_du_domega(n, L) =
368 t_Ldiff_u(l, m, L) *
369 (t_approx_P_intermidiata(i, m) * t_diff_R(i, l, n));
370
371 break;
372 case MODERATE_ROT:
373
374 t_Ru(i, m) =
375 t_kd(i, m) + (t_R(i, m) - t_kd(i, m)) + (t_u(i, m) - t_kd(i, m));
376 t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
377 t_h_domega(i, j, k) = t_diff_R(i, m, k) * t_h1(m, j);
378 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
379
380 t_approx_P_intermidiata(i, m) = t_approx_P(i, j) * t_h1(m, j);
381 t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiata(i, m);
382
383 t_levi_kirchoff(k) =
384 levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
385 t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
386 t_levi_kirchoff_domega(k, n) = 0;
387
388 t_approx_P_adjont_log_du(L) =
389 t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
390 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
391 t_approx_P_adjont_log_du_domega(n, L) = 0;
392
393 break;
394
395 case SMALL_ROT:
396
397 t_h(i, j) =
398 t_kd(i, j) + (t_R(i, j) - t_kd(i, j)) + (t_u(i, j) - t_kd(i, j));
399 t_h_domega(i, j, k) = t_diff_R(i, j, k);
400 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
401
402 t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
403 t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
404 t_levi_kirchoff_domega(k, l) = 0;
405
406 t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
407 t_approx_P_adjont_log_du(L) =
408 t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
409 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
410 t_approx_P_adjont_log_du_domega(m, L) = 0;
411
412 break;
413 }
414
416 t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
417 eigen_vec(i, j) = t_C_h1(i, j);
418 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
421 break;
423 for (int ii = 0; ii != 3; ++ii) {
424 eig(ii) = std::max(eig(ii),
426 std::numeric_limits<double>::min_exponent)));
427 }
428 break;
429 }
430
431 auto t_log_u2_h1_tmp =
433
435 case LARGE_ROT:
436 case MODERATE_ROT:
437 t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
438 break;
439 case SMALL_ROT:
440 t_log_stretch_total(i, j) = t_log_u(i, j);
441 break;
442 };
443
444 ++t_h;
445 ++t_h_domega;
446 ++t_h_dlog_u;
447 ++t_levi_kirchoff;
448 ++t_levi_kirchoff_dP;
449 ++t_levi_kirchoff_domega;
450 ++t_approx_P_adjont_dstretch;
451 ++t_approx_P_adjont_log_du;
452 ++t_approx_P_adjont_log_du_dP;
453 ++t_approx_P_adjont_log_du_domega;
454 ++t_approx_P;
455 ++t_R;
456 ++t_diff_R;
457 ++t_log_u;
458 ++t_u;
459 ++t_diff_u;
460 ++t_eigen_vals;
461 ++t_eigen_vecs;
462 ++t_omega;
463 ++t_grad_h1;
464 ++t_log_stretch_total;
465 }
466
468}
469
472 int nb_dofs = data.getIndices().size();
473 int nb_integration_pts = data.getN().size1();
474 auto v = getVolume();
475 auto t_w = getFTensor0IntegrationWeight();
476 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
477 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
478 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
479 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
480 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
481 dataAtPts->wL2DotDotAtPts.clear();
482 }
483 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
484
485 int nb_base_functions = data.getN().size2();
486 auto t_row_base_fun = data.getFTensor0N();
487
489 auto get_ftensor1 = [](auto &v) {
491 &v[2]);
492 };
493
494 for (int gg = 0; gg != nb_integration_pts; ++gg) {
495 double a = v * t_w;
496 auto t_nf = get_ftensor1(nF);
497 int bb = 0;
498 for (; bb != nb_dofs / 3; ++bb) {
499 t_nf(i) += a * t_row_base_fun * t_div_P(i);
500 t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
501 t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
502 ++t_nf;
503 ++t_row_base_fun;
504 }
505 for (; bb != nb_base_functions; ++bb)
506 ++t_row_base_fun;
507 ++t_w;
508 ++t_div_P;
509 ++t_s_dot_w;
510 ++t_s_dot_dot_w;
511 }
512
514}
515
518 int nb_dofs = data.getIndices().size();
519 int nb_integration_pts = getGaussPts().size2();
520 auto v = getVolume();
521 auto t_w = getFTensor0IntegrationWeight();
522 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchoffAtPts);
523 int nb_base_functions = data.getN().size2();
524 auto t_row_base_fun = data.getFTensor0N();
528 auto get_ftensor1 = [](auto &v) {
530 &v[2]);
531 };
532 for (int gg = 0; gg != nb_integration_pts; ++gg) {
533 double a = v * t_w;
534 auto t_nf = get_ftensor1(nF);
535 int bb = 0;
536 for (; bb != nb_dofs / 3; ++bb) {
537 t_nf(k) += (a * t_row_base_fun) * t_levi_kirchoff(k);
538 ++t_nf;
539 ++t_row_base_fun;
540 }
541 for (; bb != nb_base_functions; ++bb) {
542 ++t_row_base_fun;
543 }
544 ++t_w;
545 ++t_levi_kirchoff;
546 }
548}
549
552
553 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
555 } else {
557 }
558
560}
561
564
566 auto t_L = symm_L_tensor();
567
568 int nb_dofs = data.getIndices().size();
569 int nb_integration_pts = data.getN().size1();
570 auto v = getVolume();
571 auto t_w = getFTensor0IntegrationWeight();
572 auto t_approx_P_adjont_log_du =
573 getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
574 auto t_log_streach_h1 =
575 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
576 auto t_dot_log_u =
577 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
578
579 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
580
585 auto get_ftensor2 = [](auto &v) {
587 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
588 };
589
590 int nb_base_functions = data.getN().size2();
591 auto t_row_base_fun = data.getFTensor0N();
592 for (int gg = 0; gg != nb_integration_pts; ++gg) {
593 double a = v * t_w;
594 auto t_nf = get_ftensor2(nF);
595
597 t_T(i, j) =
598 t_D(i, j, k, l) * (t_log_streach_h1(k, l) + alphaU * t_dot_log_u(k, l));
600 t_residual(L) = a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
601
602 int bb = 0;
603 for (; bb != nb_dofs / 6; ++bb) {
604 t_nf(L) += t_row_base_fun * t_residual(L);
605 ++t_nf;
606 ++t_row_base_fun;
607 }
608 for (; bb != nb_base_functions; ++bb)
609 ++t_row_base_fun;
610
611 ++t_w;
612 ++t_approx_P_adjont_log_du;
613 ++t_dot_log_u;
614 ++t_log_streach_h1;
615 }
617}
618
621
623 auto t_L = symm_L_tensor();
624
625 int nb_dofs = data.getIndices().size();
626 int nb_integration_pts = data.getN().size1();
627 auto v = getVolume();
628 auto t_w = getFTensor0IntegrationWeight();
629 auto t_approx_P_adjont_log_du =
630 getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
631 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
632 auto t_dot_log_u =
633 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
634 auto t_diff_u =
635 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
636
641 auto get_ftensor2 = [](auto &v) {
643 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
644 };
645
646 int nb_base_functions = data.getN().size2();
647 auto t_row_base_fun = data.getFTensor0N();
648 for (int gg = 0; gg != nb_integration_pts; ++gg) {
649 double a = v * t_w;
650 auto t_nf = get_ftensor2(nF);
651
653 t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
654
656 t_viscous_P(i, j) =
657 alphaU *
658 t_dot_log_u(i,
659 j); // That is chit, should be split on axiator and deviator
660
662 t_residual(L) =
663 a * (t_approx_P_adjont_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
664 t_L(i, j, L) * t_viscous_P(i, j));
665
666 int bb = 0;
667 for (; bb != nb_dofs / 6; ++bb) {
668 t_nf(L) += t_row_base_fun * t_residual(L);
669 ++t_nf;
670 ++t_row_base_fun;
671 }
672 for (; bb != nb_base_functions; ++bb)
673 ++t_row_base_fun;
674
675 ++t_w;
676 ++t_approx_P_adjont_log_du;
677 ++t_P;
678 ++t_dot_log_u;
679 ++t_diff_u;
680 }
682}
683
686 int nb_dofs = data.getIndices().size();
687 int nb_integration_pts = data.getN().size1();
688 auto v = getVolume();
689 auto t_w = getFTensor0IntegrationWeight();
690 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
691 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
692
693 int nb_base_functions = data.getN().size2() / 3;
694 auto t_row_base_fun = data.getFTensor1N<3>();
699
700 auto get_ftensor1 = [](auto &v) {
702 &v[2]);
703 };
704
705 for (int gg = 0; gg != nb_integration_pts; ++gg) {
706 double a = v * t_w;
707 auto t_nf = get_ftensor1(nF);
708
709 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
711 t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
712
713 int bb = 0;
714 for (; bb != nb_dofs / 3; ++bb) {
715 t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
716 ++t_nf;
717 ++t_row_base_fun;
718 }
719
720 for (; bb != nb_base_functions; ++bb)
721 ++t_row_base_fun;
722
723 ++t_w;
724 ++t_h;
725 ++t_omega_dot;
726 }
727
729}
730
733 int nb_dofs = data.getIndices().size();
734 int nb_integration_pts = data.getN().size1();
735 auto v = getVolume();
736 auto t_w = getFTensor0IntegrationWeight();
737 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
738 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
739
740 int nb_base_functions = data.getN().size2() / 9;
741 auto t_row_base_fun = data.getFTensor2N<3, 3>();
746
747 auto get_ftensor0 = [](auto &v) {
749 };
750
751 for (int gg = 0; gg != nb_integration_pts; ++gg) {
752 double a = v * t_w;
753 auto t_nf = get_ftensor0(nF);
754
755 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
757 t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
758
759 int bb = 0;
760 for (; bb != nb_dofs; ++bb) {
761 t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
762 ++t_nf;
763 ++t_row_base_fun;
764 }
765 for (; bb != nb_base_functions; ++bb) {
766 ++t_row_base_fun;
767 }
768 ++t_w;
769 ++t_h;
770 ++t_omega_dot;
771 }
772
774}
775
778 int nb_dofs = data.getIndices().size();
779 int nb_integration_pts = data.getN().size1();
780 auto v = getVolume();
781 auto t_w = getFTensor0IntegrationWeight();
782 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
783 int nb_base_functions = data.getN().size2() / 3;
784 auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
786 auto get_ftensor1 = [](auto &v) {
788 &v[2]);
789 };
790
791 for (int gg = 0; gg != nb_integration_pts; ++gg) {
792 double a = v * t_w;
793 auto t_nf = get_ftensor1(nF);
794 int bb = 0;
795 for (; bb != nb_dofs / 3; ++bb) {
796 double div_row_base = t_row_diff_base_fun(i, i);
797 t_nf(i) += a * div_row_base * t_w_l2(i);
798 ++t_nf;
799 ++t_row_diff_base_fun;
800 }
801 for (; bb != nb_base_functions; ++bb) {
802 ++t_row_diff_base_fun;
803 }
804 ++t_w;
805 ++t_w_l2;
806 }
807
809}
810
813 // get entity of face
814 EntityHandle fe_ent = getFEEntityHandle();
815 // interate over all boundary data
816 for (auto &bc : (*bcDispPtr)) {
817 // check if finite element entity is part of boundary condition
818 if (bc.faces.find(fe_ent) != bc.faces.end()) {
819 int nb_dofs = data.getIndices().size();
820 int nb_integration_pts = data.getN().size1();
821 auto t_normal = getFTensor1Normal();
822 auto t_w = getFTensor0IntegrationWeight();
823 int nb_base_functions = data.getN().size2() / 3;
824 auto t_row_base_fun = data.getFTensor1N<3>();
825
828 auto get_ftensor1 = [](auto &v) {
830 &v[2]);
831 };
832
833 // get bc data
834 FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
835 t_bc_disp(i) *= getFEMethod()->ts_t;
836
837 for (int gg = 0; gg != nb_integration_pts; ++gg) {
839 t_bc_res(i) = t_bc_disp(i);
840 auto t_nf = get_ftensor1(nF);
841 int bb = 0;
842 for (; bb != nb_dofs / 3; ++bb) {
843 t_nf(i) -=
844 t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
845 ++t_nf;
846 ++t_row_base_fun;
847 }
848 for (; bb != nb_base_functions; ++bb)
849 ++t_row_base_fun;
850
851 ++t_w;
852 }
853 }
854 }
856}
857
860 // get entity of face
861 EntityHandle fe_ent = getFEEntityHandle();
862 // interate over all boundary data
863 for (auto &bc : (*bcRotPtr)) {
864 // check if finite element entity is part of boundary condition
865 if (bc.faces.find(fe_ent) != bc.faces.end()) {
866 int nb_dofs = data.getIndices().size();
867 int nb_integration_pts = data.getN().size1();
868 auto t_normal = getFTensor1Normal();
869 auto t_w = getFTensor0IntegrationWeight();
870
871 int nb_base_functions = data.getN().size2() / 3;
872 auto t_row_base_fun = data.getFTensor1N<3>();
876 auto get_ftensor1 = [](auto &v) {
878 &v[2]);
879 };
880
881 // get bc data
882 FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
883
884 double theta = bc.theta;
885 theta *= getFEMethod()->ts_t;
886
888 const double a = sqrt(t_normal(i) * t_normal(i));
889 t_omega(i) = t_normal(i) * (theta / a);
890 auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
891 auto t_coords = getFTensor1CoordsAtGaussPts();
892
893 for (int gg = 0; gg != nb_integration_pts; ++gg) {
895 t_delta(i) = t_center(i) - t_coords(i);
897 t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
898
899 auto t_nf = get_ftensor1(nF);
900 int bb = 0;
901 for (; bb != nb_dofs / 3; ++bb) {
902 t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
903 ++t_nf;
904 ++t_row_base_fun;
905 }
906 for (; bb != nb_base_functions; ++bb)
907 ++t_row_base_fun;
908
909 ++t_w;
910 ++t_coords;
911 }
912 }
913 }
915}
916
919
920 struct FaceRule {
921 int operator()(int p_row, int p_col, int p_data) const {
922 return 2 * (p_data + 1);
923 }
924 };
925
926 if (ts_ctx == CTX_TSSETIFUNCTION) {
927
928 // Loop boundary elements with traction boundary conditions
931 {HDIV});
932 fe.getOpPtrVector().push_back(
933 new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
934 fe.getRuleHook = FaceRule();
935 fe.ts_t = ts_t;
936 CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
937 fe, nullptr, MF_ZERO,
938 this->getCacheWeakPtr());
939 }
940
942}
943
946
949
950 auto t_normal = getFTensor1Normal();
951 const double nrm2 = sqrt(t_normal(i) * t_normal(i));
952 FTensor::Tensor1<double, 3> t_unit_normal;
953 t_unit_normal(i) = t_normal(i) / nrm2;
954 int nb_dofs = data.getFieldData().size();
955 int nb_integration_pts = data.getN().size1();
956 int nb_base_functions = data.getN().size2() / 3;
957 double ts_t = getFEMethod()->ts_t;
958
959 auto integrate_matrix = [&]() {
961
962 auto t_w = getFTensor0IntegrationWeight();
963 matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
964 matPP.clear();
965
966 auto t_row_base_fun = data.getFTensor1N<3>();
967 for (int gg = 0; gg != nb_integration_pts; ++gg) {
968
969 int rr = 0;
970 for (; rr != nb_dofs / 3; ++rr) {
971 const double a = t_row_base_fun(i) * t_unit_normal(i);
972 auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
973 for (int cc = 0; cc != nb_dofs / 3; ++cc) {
974 const double b = t_col_base_fun(i) * t_unit_normal(i);
975 matPP(rr, cc) += t_w * a * b;
976 ++t_col_base_fun;
977 }
978 ++t_row_base_fun;
979 }
980
981 for (; rr != nb_base_functions; ++rr)
982 ++t_row_base_fun;
983
984 ++t_w;
985 }
986
988 };
989
990 auto integrate_rhs = [&](auto &bc) {
992
993 auto t_w = getFTensor0IntegrationWeight();
994 vecPv.resize(3, nb_dofs / 3, false);
995 vecPv.clear();
996
997 auto t_row_base_fun = data.getFTensor1N<3>();
998 double ts_t = getFEMethod()->ts_t;
999
1000 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1001 int rr = 0;
1002 for (; rr != nb_dofs / 3; ++rr) {
1003 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1004 for (int dd = 0; dd != 3; ++dd)
1005 if (bc.flags[dd])
1006 vecPv(dd, rr) += t * bc.vals[dd];
1007 ++t_row_base_fun;
1008 }
1009 for (; rr != nb_base_functions; ++rr)
1010 ++t_row_base_fun;
1011 ++t_w;
1012 }
1014 };
1015
1016 auto integrate_rhs_cook = [&](auto &bc) {
1018
1019 vecPv.resize(3, nb_dofs / 3, false);
1020 vecPv.clear();
1021
1022 auto t_w = getFTensor0IntegrationWeight();
1023 auto t_row_base_fun = data.getFTensor1N<3>();
1024 auto t_coords = getFTensor1CoordsAtGaussPts();
1025
1026 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1027
1028 auto calc_tau = [](double y) {
1029 y -= 44;
1030 y /= (60 - 44);
1031 return -y * (y - 1) / 0.25;
1032 };
1033
1034 const double tau = calc_tau(t_coords(1));
1035
1036 int rr = 0;
1037 for (; rr != nb_dofs / 3; ++rr) {
1038 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1039 for (int dd = 0; dd != 3; ++dd)
1040 if (bc.flags[dd])
1041 vecPv(dd, rr) += t * tau * bc.vals[dd];
1042 ++t_row_base_fun;
1043 }
1044
1045 for (; rr != nb_base_functions; ++rr)
1046 ++t_row_base_fun;
1047 ++t_w;
1048 ++t_coords;
1049 }
1051 };
1052
1053 // get entity of face
1054 EntityHandle fe_ent = getFEEntityHandle();
1055 for (auto &bc : *(bcFe.bcData)) {
1056 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1057
1058 int nb_dofs = data.getFieldData().size();
1059 if (nb_dofs) {
1060
1061 CHKERR integrate_matrix();
1062 if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
1063 CHKERR integrate_rhs_cook(bc);
1064 else
1065 CHKERR integrate_rhs(bc);
1066
1067 const auto info =
1068 lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1069 nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1070 if (info > 0)
1071 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1072 "The leading minor of order %i is "
1073 "not positive; definite;\nthe "
1074 "solution could not be computed",
1075 info);
1076
1077 for (int dd = 0; dd != 3; ++dd)
1078 if (bc.flags[dd])
1079 for (int rr = 0; rr != nb_dofs / 3; ++rr)
1080 data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1081 }
1082 }
1083 }
1084
1086}
1087
1089 EntData &col_data) {
1091 int nb_integration_pts = row_data.getN().size1();
1092 int row_nb_dofs = row_data.getIndices().size();
1093 int col_nb_dofs = col_data.getIndices().size();
1094 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1096 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
1097 };
1099 auto v = getVolume();
1100 auto t_w = getFTensor0IntegrationWeight();
1101 int row_nb_base_functions = row_data.getN().size2();
1102 auto t_row_base_fun = row_data.getFTensor0N();
1103 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1104 double a = v * t_w;
1105 int rr = 0;
1106 for (; rr != row_nb_dofs / 3; ++rr) {
1107 auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
1108 auto t_m = get_ftensor1(K, 3 * rr, 0);
1109 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1110 double div_col_base = t_col_diff_base_fun(i, i);
1111 t_m(i) += a * t_row_base_fun * div_col_base;
1112 ++t_m;
1113 ++t_col_diff_base_fun;
1114 }
1115 ++t_row_base_fun;
1116 }
1117 for (; rr != row_nb_base_functions; ++rr)
1118 ++t_row_base_fun;
1119 ++t_w;
1120 }
1122}
1123
1125 EntData &col_data) {
1127
1128 if (alphaW < std::numeric_limits<double>::epsilon() &&
1129 alphaRho < std::numeric_limits<double>::epsilon())
1131
1132 const int nb_integration_pts = row_data.getN().size1();
1133 const int row_nb_dofs = row_data.getIndices().size();
1134 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1136 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1137
1138 );
1139 };
1141
1142 auto v = getVolume();
1143 auto t_w = getFTensor0IntegrationWeight();
1144
1145 int row_nb_base_functions = row_data.getN().size2();
1146 auto t_row_base_fun = row_data.getFTensor0N();
1147
1148 double ts_scale = alphaW * getTSa();
1149 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1150 ts_scale += alphaRho * getTSaa();
1151
1152 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1153 double a = v * t_w * ts_scale;
1154
1155 int rr = 0;
1156 for (; rr != row_nb_dofs / 3; ++rr) {
1157
1158 auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1159 auto t_m = get_ftensor1(K, 3 * rr, 0);
1160 for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1161 const double b = a * t_row_base_fun * t_col_base_fun;
1162 t_m(i) -= b;
1163 ++t_m;
1164 ++t_col_base_fun;
1165 }
1166
1167 ++t_row_base_fun;
1168 }
1169
1170 for (; rr != row_nb_base_functions; ++rr)
1171 ++t_row_base_fun;
1172
1173 ++t_w;
1174 }
1175
1177}
1178
1180 EntData &col_data) {
1182
1184
1185 int nb_integration_pts = row_data.getN().size1();
1186 int row_nb_dofs = row_data.getIndices().size();
1187 int col_nb_dofs = col_data.getIndices().size();
1188 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1190
1191 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1192
1193 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1194
1195 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1196
1197 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1198
1199 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1200
1201 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1202 };
1203
1206
1207 auto v = getVolume();
1208 auto t_w = getFTensor0IntegrationWeight();
1209 auto t_approx_P_adjont_log_du_dP =
1210 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
1211 int row_nb_base_functions = row_data.getN().size2();
1212 auto t_row_base_fun = row_data.getFTensor0N();
1213 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1214 double a = v * t_w;
1215 int rr = 0;
1216 for (; rr != row_nb_dofs / 6; ++rr) {
1217
1218 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1219 auto t_m = get_ftensor3(K, 6 * rr, 0);
1220
1221 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1222 t_m(L, i) += a *
1223 (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(j)) *
1224 t_row_base_fun;
1225 ++t_col_base_fun;
1226 ++t_m;
1227 }
1228
1229 ++t_row_base_fun;
1230 }
1231 for (; rr != row_nb_base_functions; ++rr)
1232 ++t_row_base_fun;
1233 ++t_w;
1234 ++t_approx_P_adjont_log_du_dP;
1235 }
1237}
1238
1240 EntData &col_data) {
1242
1244
1245 int nb_integration_pts = row_data.getN().size1();
1246 int row_nb_dofs = row_data.getIndices().size();
1247 int col_nb_dofs = col_data.getIndices().size();
1248 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1250 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1251 &m(r + 5, c));
1252 };
1253
1256
1257 auto v = getVolume();
1258 auto t_w = getFTensor0IntegrationWeight();
1259 auto t_approx_P_adjont_log_du_dP =
1260 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
1261
1262 int row_nb_base_functions = row_data.getN().size2();
1263 auto t_row_base_fun = row_data.getFTensor0N();
1264 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1265 double a = v * t_w;
1266 int rr = 0;
1267 for (; rr != row_nb_dofs / 6; ++rr) {
1268 auto t_m = get_ftensor2(K, 6 * rr, 0);
1269 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1270 for (int cc = 0; cc != col_nb_dofs; ++cc) {
1271 t_m(L) += a *
1272 (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
1273 t_row_base_fun;
1274 ++t_m;
1275 ++t_col_base_fun;
1276 }
1277 ++t_row_base_fun;
1278 }
1279 for (; rr != row_nb_base_functions; ++rr)
1280 ++t_row_base_fun;
1281 ++t_w;
1282 ++t_approx_P_adjont_log_du_dP;
1283 }
1285}
1286
1288 EntData &col_data) {
1290 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1291 CHKERR integratePiola(row_data, col_data);
1292 } else {
1293 CHKERR integrateHencky(row_data, col_data);
1294 }
1296}
1297
1299 EntData &col_data) {
1301
1304 auto t_L = symm_L_tensor();
1305 auto t_diff = diff_tensor();
1306
1307 int nb_integration_pts = row_data.getN().size1();
1308 int row_nb_dofs = row_data.getIndices().size();
1309 int col_nb_dofs = col_data.getIndices().size();
1310
1311 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1313 size_symm>(
1314
1315 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1316 &m(r + 0, c + 4), &m(r + 0, c + 5),
1317
1318 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1319 &m(r + 1, c + 4), &m(r + 1, c + 5),
1320
1321 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1322 &m(r + 2, c + 4), &m(r + 2, c + 5),
1323
1324 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1325 &m(r + 3, c + 4), &m(r + 3, c + 5),
1326
1327 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1328 &m(r + 4, c + 4), &m(r + 4, c + 5),
1329
1330 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1331 &m(r + 5, c + 4), &m(r + 5, c + 5)
1332
1333 );
1334 };
1341
1345
1346 auto v = getVolume();
1347 auto t_w = getFTensor0IntegrationWeight();
1348
1349 auto t_approx_P_adjont_dstretch =
1350 getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
1351 auto t_dot_log_u =
1352 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
1353 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1354 auto t_diff_u =
1355 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
1356 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1357 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1358 auto &nbUniq = dataAtPts->nbUniq;
1359
1360 int row_nb_base_functions = row_data.getN().size2();
1361 auto t_row_base_fun = row_data.getFTensor0N();
1362
1363 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1364
1365 const double ts_a = getTSa();
1366 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1367 double a = v * t_w;
1368
1370
1372 // Work of symmetric tensor on undefined tensor is equal to the work of
1373 // the symmetric part of it
1375 t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
1376 t_approx_P_adjont_dstretch(j, i));
1377 t_sym(i, j) /= 2.0;
1378 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1379 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1380 EshelbianCore::dd_f, t_sym, nbUniq[gg]);
1381 t_dP(L, J) =
1382 t_L(i, j, L) *
1383 ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) * t_L(k, l, J)) /
1384 2.;
1385 } else {
1386 t_dP(L, J) = 0;
1387 }
1388
1389 t_dP(L, J) -= (1 + alphaU * ts_a) *
1390 (t_L(i, j, L) *
1391 ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
1392
1393 int rr = 0;
1394 for (; rr != row_nb_dofs / 6; ++rr) {
1395 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1396 auto t_m = get_ftensor2(K, 6 * rr, 0);
1397 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1398 const double b = a * t_row_base_fun * t_col_base_fun;
1399 t_m(L, J) += b * t_dP(L, J);
1400 ++t_m;
1401 ++t_col_base_fun;
1402 }
1403 ++t_row_base_fun;
1404 }
1405
1406 for (; rr != row_nb_base_functions; ++rr) {
1407 ++t_row_base_fun;
1408 }
1409
1410 ++t_w;
1411 ++t_approx_P_adjont_dstretch;
1412 ++t_dot_log_u;
1413 ++t_u;
1414 ++t_diff_u;
1415 ++t_eigen_vals;
1416 ++t_eigen_vecs;
1417 }
1419}
1420
1422 EntData &col_data) {
1424
1427 auto t_L = symm_L_tensor();
1428 auto t_diff = diff_tensor();
1429
1430 int nb_integration_pts = row_data.getN().size1();
1431 int row_nb_dofs = row_data.getIndices().size();
1432 int col_nb_dofs = col_data.getIndices().size();
1433
1434 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1436 size_symm>(
1437
1438 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1439 &m(r + 0, c + 4), &m(r + 0, c + 5),
1440
1441 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1442 &m(r + 1, c + 4), &m(r + 1, c + 5),
1443
1444 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1445 &m(r + 2, c + 4), &m(r + 2, c + 5),
1446
1447 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1448 &m(r + 3, c + 4), &m(r + 3, c + 5),
1449
1450 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1451 &m(r + 4, c + 4), &m(r + 4, c + 5),
1452
1453 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1454 &m(r + 5, c + 4), &m(r + 5, c + 5)
1455
1456 );
1457 };
1464
1465 auto v = getVolume();
1466 auto t_w = getFTensor0IntegrationWeight();
1467 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
1468
1469 auto t_approx_P_adjont_dstretch =
1470 getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
1471 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1472 auto t_dot_log_u =
1473 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
1474 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1475 auto t_diff_u =
1476 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
1477 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1478 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1479 auto &nbUniq = dataAtPts->nbUniq;
1480
1481 int row_nb_base_functions = row_data.getN().size2();
1482 auto t_row_base_fun = row_data.getFTensor0N();
1483
1484 const double ts_a = getTSa();
1485 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1486 double a = v * t_w;
1487
1489 t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
1490
1491 // Work of symmetric tensor on undefined tensor is equal to the work of the
1492 // symmetric part of it
1494
1497 t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
1498 t_deltaP_sym(i, j) /= 2.0;
1499 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1500 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1501 EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
1502 t_dP(L, J) = t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
1503 } else {
1504 t_dP(L, J) = 0;
1505 }
1506
1508 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
1509 t_dP(L, J) -=
1510 t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
1511 // viscous stress derivative
1512 t_dP(L, J) -=
1513 (alphaU * ts_a) * (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
1514
1515 int rr = 0;
1516 for (; rr != row_nb_dofs / 6; ++rr) {
1517 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1518 auto t_m = get_ftensor2(K, 6 * rr, 0);
1519 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1520 const double b = a * t_row_base_fun * t_col_base_fun;
1521 t_m(L, J) += b * t_dP(L, J);
1522 ++t_m;
1523 ++t_col_base_fun;
1524 }
1525 ++t_row_base_fun;
1526 }
1527
1528 for (; rr != row_nb_base_functions; ++rr) {
1529 ++t_row_base_fun;
1530 }
1531
1532 ++t_w;
1533 ++r_P_du;
1534 ++t_P;
1535 ++t_approx_P_adjont_dstretch;
1536 ++t_dot_log_u;
1537 ++t_u;
1538 ++t_diff_u;
1539 ++t_eigen_vals;
1540 ++t_eigen_vecs;
1541 }
1543}
1544
1546 EntData &col_data) {
1548
1550 auto t_L = symm_L_tensor();
1551
1552 int row_nb_dofs = row_data.getIndices().size();
1553 int col_nb_dofs = col_data.getIndices().size();
1554 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1556
1557 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1558
1559 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1560
1561 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1562
1563 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1564
1565 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1566
1567 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1568
1569 );
1570 };
1576
1577 auto v = getVolume();
1578 auto t_w = getFTensor0IntegrationWeight();
1579 auto t_approx_P_adjont_log_du_domega =
1580 getFTensor2FromMat<3, size_symm>(dataAtPts->adjontPdUdOmegaAtPts);
1581
1582 int row_nb_base_functions = row_data.getN().size2();
1583 auto t_row_base_fun = row_data.getFTensor0N();
1584
1585 int nb_integration_pts = row_data.getN().size1();
1586
1587 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1588 double a = v * t_w;
1589
1590 int rr = 0;
1591 for (; rr != row_nb_dofs / 6; ++rr) {
1592 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1593 auto t_m = get_ftensor3(K, 6 * rr, 0);
1594 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1595 double v = a * t_row_base_fun * t_col_base_fun;
1596 t_m(L, k) += v * t_approx_P_adjont_log_du_domega(k, L);
1597 ++t_m;
1598 ++t_col_base_fun;
1599 }
1600 ++t_row_base_fun;
1601 }
1602
1603 for (; rr != row_nb_base_functions; ++rr)
1604 ++t_row_base_fun;
1605
1606 ++t_w;
1607 ++t_approx_P_adjont_log_du_domega;
1608 }
1609
1611}
1612
1614 EntData &col_data) {
1616 int nb_integration_pts = getGaussPts().size2();
1617 int row_nb_dofs = row_data.getIndices().size();
1618 int col_nb_dofs = col_data.getIndices().size();
1619 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1621
1622 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1623
1624 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1625
1626 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1627
1628 );
1629 };
1636 auto v = getVolume();
1637 auto t_w = getFTensor0IntegrationWeight();
1638 auto t_levi_kirchoff_dP =
1639 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchoffdPAtPts);
1640
1641 int row_nb_base_functions = row_data.getN().size2();
1642 auto t_row_base_fun = row_data.getFTensor0N();
1643
1644 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1645 double a = v * t_w;
1646 int rr = 0;
1647 for (; rr != row_nb_dofs / 3; ++rr) {
1648 double b = a * t_row_base_fun;
1649 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1650 auto t_m = get_ftensor2(K, 3 * rr, 0);
1651 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1652 t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(l));
1653 ++t_m;
1654 ++t_col_base_fun;
1655 }
1656 ++t_row_base_fun;
1657 }
1658 for (; rr != row_nb_base_functions; ++rr) {
1659 ++t_row_base_fun;
1660 }
1661
1662 ++t_w;
1663 ++t_levi_kirchoff_dP;
1664 }
1666}
1667
1669 EntData &col_data) {
1671 int nb_integration_pts = getGaussPts().size2();
1672 int row_nb_dofs = row_data.getIndices().size();
1673 int col_nb_dofs = col_data.getIndices().size();
1674
1675 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1677 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1678 };
1679
1685
1686 auto v = getVolume();
1687 auto t_w = getFTensor0IntegrationWeight();
1688 auto t_levi_kirchoff_dP =
1689 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchoffdPAtPts);
1690
1691 int row_nb_base_functions = row_data.getN().size2();
1692 auto t_row_base_fun = row_data.getFTensor0N();
1693
1694 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1695 double a = v * t_w;
1696 int rr = 0;
1697 for (; rr != row_nb_dofs / 3; ++rr) {
1698 double b = a * t_row_base_fun;
1699 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1700 auto t_m = get_ftensor1(K, 3 * rr, 0);
1701 for (int cc = 0; cc != col_nb_dofs; ++cc) {
1702 t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(i, j));
1703 ++t_m;
1704 ++t_col_base_fun;
1705 }
1706 ++t_row_base_fun;
1707 }
1708
1709 for (; rr != row_nb_base_functions; ++rr) {
1710 ++t_row_base_fun;
1711 }
1712 ++t_w;
1713 ++t_levi_kirchoff_dP;
1714 }
1716}
1717
1719 EntData &col_data) {
1721 int nb_integration_pts = getGaussPts().size2();
1722 int row_nb_dofs = row_data.getIndices().size();
1723 int col_nb_dofs = col_data.getIndices().size();
1724 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1726
1727 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1728
1729 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1730
1731 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1732
1733 );
1734 };
1741
1742 auto v = getVolume();
1743 auto t_w = getFTensor0IntegrationWeight();
1744 auto t_levi_kirchoff_domega =
1745 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchoffdOmegaAtPts);
1746 int row_nb_base_functions = row_data.getN().size2();
1747 auto t_row_base_fun = row_data.getFTensor0N();
1748 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1749 double a = v * t_w;
1750 int rr = 0;
1751 for (; rr != row_nb_dofs / 3; ++rr) {
1752 auto t_m = get_ftensor2(K, 3 * rr, 0);
1753 const double b = a * t_row_base_fun;
1754 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1755 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1756 t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(k, l);
1757 ++t_m;
1758 ++t_col_base_fun;
1759 }
1760 ++t_row_base_fun;
1761 }
1762 for (; rr != row_nb_base_functions; ++rr) {
1763 ++t_row_base_fun;
1764 }
1765 ++t_w;
1766 ++t_levi_kirchoff_domega;
1767 }
1769}
1770
1772 EntData &col_data) {
1774
1775 int nb_integration_pts = row_data.getN().size1();
1776 int row_nb_dofs = row_data.getIndices().size();
1777 int col_nb_dofs = col_data.getIndices().size();
1778
1779 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1781
1782 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1783
1784 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1785
1786 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1787
1788 );
1789 };
1790
1797
1798 auto v = getVolume();
1799 auto t_w = getFTensor0IntegrationWeight();
1800 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1801 int row_nb_base_functions = row_data.getN().size2() / 3;
1802 auto t_row_base_fun = row_data.getFTensor1N<3>();
1803
1804 const double ts_a = getTSa();
1805
1806 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1807 double a = v * t_w;
1808
1809 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1810
1811 int rr = 0;
1812 for (; rr != row_nb_dofs / 3; ++rr) {
1813
1815 t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
1816
1817 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1818 auto t_m = get_ftensor2(K, 3 * rr, 0);
1819 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1820 t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1821 ++t_m;
1822 ++t_col_base_fun;
1823 }
1824
1825 ++t_row_base_fun;
1826 }
1827
1828 for (; rr != row_nb_base_functions; ++rr)
1829 ++t_row_base_fun;
1830 ++t_w;
1831 ++t_h_domega;
1832 }
1834}
1835
1838 EntData &col_data) {
1840
1841 int nb_integration_pts = row_data.getN().size1();
1842 int row_nb_dofs = row_data.getIndices().size();
1843 int col_nb_dofs = col_data.getIndices().size();
1844
1845 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1847 &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1848 };
1849
1856
1857 auto v = getVolume();
1858 auto t_w = getFTensor0IntegrationWeight();
1859 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1860 int row_nb_base_functions = row_data.getN().size2() / 9;
1861 auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1862
1863 const double ts_a = getTSa();
1864
1865 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1866 double a = v * t_w;
1867
1868 int rr = 0;
1869 for (; rr != row_nb_dofs; ++rr) {
1870
1872 t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
1873
1874 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1875 auto t_m = get_ftensor2(K, rr, 0);
1876 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1877 t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1878 ++t_m;
1879 ++t_col_base_fun;
1880 }
1881
1882 ++t_row_base_fun;
1883 }
1884
1885 for (; rr != row_nb_base_functions; ++rr)
1886 ++t_row_base_fun;
1887
1888 ++t_w;
1889 ++t_h_domega;
1890 }
1892}
1893
1895 EntData &data) {
1897
1898 auto create_tag = [this](const std::string tag_name, const int size) {
1899 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1900 Tag th;
1901 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
1902 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1903 def_VAL);
1904 return th;
1905 };
1906
1907 Tag th_w = create_tag("SpatialDisplacement", 3);
1908 Tag th_omega = create_tag("Omega", 3);
1909 Tag th_approxP = create_tag("Piola1Stress", 9);
1910 Tag th_sigma = create_tag("CauchyStress", 9);
1911 Tag th_log_u = create_tag("LogSpatialStretch", 9);
1912 Tag th_u = create_tag("SpatialStretch", 9);
1913 Tag th_h = create_tag("h", 9);
1914 Tag th_X = create_tag("X", 3);
1915 Tag th_detF = create_tag("detF", 1);
1916 Tag th_angular_momentum = create_tag("AngularMomentum", 3);
1917
1918 Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
1919 Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
1920 Tag th_traction = create_tag("traction", 3);
1921
1922 Tag th_disp = create_tag("U", 3);
1923 Tag th_disp_error = create_tag("U_ERROR", 1);
1924
1925 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1926 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
1927 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1928 auto t_log_u =
1929 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
1930 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1931 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1932 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1933 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchoffAtPts);
1934 auto t_coords = getFTensor1CoordsAtGaussPts();
1935 auto t_normal = getFTensor1NormalsAtGaussPts();
1936 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
1937
1942
1943 auto set_float_precision = [](const double x) {
1944 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1945 return 0.;
1946 else
1947 return x;
1948 };
1949
1950 // scalars
1951 auto save_scal_tag = [&](auto &th, auto v, const int gg) {
1953 v = set_float_precision(v);
1954 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
1956 };
1957
1958 // vectors
1959 VectorDouble3 v(3);
1960 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
1961 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
1963 t_v(i) = t_d(i);
1964 for (auto &a : v.data())
1965 a = set_float_precision(a);
1966 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1967 &*v.data().begin());
1969 };
1970
1971 // tensors
1972
1973 MatrixDouble3by3 m(3, 3);
1975 &m(0, 0), &m(0, 1), &m(0, 2),
1976
1977 &m(1, 0), &m(1, 1), &m(1, 2),
1978
1979 &m(2, 0), &m(2, 1), &m(2, 2));
1980
1981 auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
1983 t_m(i, j) = t_d(i, j);
1984 for (auto &v : m.data())
1985 v = set_float_precision(v);
1986 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1987 &*m.data().begin());
1989 };
1990
1991 const auto nb_gauss_pts = getGaussPts().size2();
1992 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1993
1994 FTensor::Tensor1<double, 3> t_traction;
1995 t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
1996
1997 // vectors
1998 CHKERR save_vec_tag(th_w, t_w, gg);
1999 CHKERR save_vec_tag(th_X, t_coords, gg);
2000 CHKERR save_vec_tag(th_omega, t_omega, gg);
2001 CHKERR save_vec_tag(th_traction, t_traction, gg);
2002
2003 // tensors
2004 CHKERR save_mat_tag(th_h, t_h, gg);
2005
2007 for (int ii = 0; ii != 3; ++ii)
2008 for (int jj = 0; jj != 3; ++jj)
2009 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2010
2011 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2012
2014 for (int ii = 0; ii != 3; ++ii)
2015 for (int jj = 0; jj != 3; ++jj)
2016 t_u_tmp(ii, jj) = t_u(ii, jj);
2017
2018 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2019 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2020 CHKERR save_vec_tag(th_disp, t_disp, gg);
2021
2022 double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
2023 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2024
2025 const double jac = determinantTensor3by3(t_h);
2027 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2028 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2029 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2030
2032 t_levi(k) = t_levi_kirchoff(k);
2033 CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2034 &t_levi(0));
2035
2036 auto get_eiegn_vector_symmetric = [&](auto &t_u) {
2038
2039 for (int ii = 0; ii != 3; ++ii)
2040 for (int jj = 0; jj != 3; ++jj)
2041 t_m(ii, jj) = t_u(ii, jj);
2042
2043 VectorDouble3 eigen_values(3);
2044 auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2045 CHKERR computeEigenValuesSymmetric(t_m, t_eigen_values);
2046
2047 CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2048 &*m.data().begin());
2049 CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2050 &*eigen_values.data().begin());
2051
2053 };
2054
2055 CHKERR get_eiegn_vector_symmetric(t_u);
2056
2057 ++t_w;
2058 ++t_h;
2059 ++t_log_u;
2060 ++t_u;
2061 ++t_omega;
2062 ++t_R;
2063 ++t_approx_P;
2064 ++t_levi_kirchoff;
2065 ++t_coords;
2066 ++t_normal;
2067 ++t_disp;
2068 }
2069
2071}
2072
2077 if (type == MBTET) {
2078 int nb_integration_pts = data.getN().size1();
2079 auto v = getVolume();
2080 auto t_w = getFTensor0IntegrationWeight();
2081 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2082 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2083
2086
2087 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2088 const double a = t_w * v;
2089 // FIXME: this is wrong, energy should be calculated in material
2090 (*energy) += a * t_P(i, J) * t_h(i, J);
2091 ++t_w;
2092 ++t_P;
2093 ++t_h;
2094 }
2095 }
2097}
2098
2099} // namespace EshelbianPlasticity
static Index< 'J', 3 > J
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
Eshelbian plasticity interface.
constexpr double a
static const double eps
Kronecker Delta class symmetric.
Kronecker Delta class.
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:211
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto get_uniq_nb(double *ptr)
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
static constexpr auto size_symm
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
auto is_eq(const double &a, const double &b)
auto get_diff2_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1469
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1557
constexpr double t
plate stiffness
Definition: plate.cpp:59
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> inv_f
static enum StretchSelector stretchSelector
static boost::function< double(const double)> d_f
static boost::function< double(const double)> f
boost::shared_ptr< TractionBcVec > bcData
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< BcDispVec > bcDispPtr
MoFEMErrorCode integrate(EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< BcRotVec > bcRotPtr
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrateHencky(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
static auto check(const double &a, const double &b)
Add operators pushing bases from local to physical configuration.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.
PetscReal ts_t
time