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
20inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
21 FTensor::Index<'i', 3> i;
22 FTensor::Index<'j', 3> j;
23 FTensor::Index<'k', 3> k;
24 FTensor::Index<'l', 3> l;
25
26 FTensor::Ddg<double, 3, 3> t_diff_deviator;
27 t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
28
29 constexpr double third = boost::math::constants::third<double>();
30
31 t_diff_deviator(0, 0, 0, 0) -= third;
32 t_diff_deviator(0, 0, 1, 1) -= third;
33
34 t_diff_deviator(1, 1, 0, 0) -= third;
35 t_diff_deviator(1, 1, 1, 1) -= third;
36
37 t_diff_deviator(2, 2, 0, 0) -= third;
38 t_diff_deviator(2, 2, 1, 1) -= third;
39
40 t_diff_deviator(0, 0, 2, 2) -= third;
41 t_diff_deviator(1, 1, 2, 2) -= third;
42 t_diff_deviator(2, 2, 2, 2) -= third;
43
44 return t_diff_deviator;
45}
46
47constexpr static auto size_symm = (3 * (3 + 1)) / 2;
48
49inline auto diff_tensor() {
50 FTensor::Index<'i', 3> i;
51 FTensor::Index<'j', 3> j;
52 FTensor::Index<'k', 3> k;
53 FTensor::Index<'l', 3> l;
56 t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
57 return t_diff;
58};
59
60inline auto symm_L_tensor() {
61 FTensor::Index<'i', 3> i;
62 FTensor::Index<'j', 3> j;
65 t_L(i, j, L) = 0;
66 t_L(0, 0, 0) = 1;
67 t_L(1, 1, 3) = 1;
68 t_L(2, 2, 5) = 1;
69 t_L(1, 0, 1) = 1;
70 t_L(2, 0, 2) = 1;
71 t_L(2, 1, 4) = 1;
72 return t_L;
73}
74
75template <class T> struct TensorTypeExtractor {
76 typedef typename std::remove_pointer<T>::type Type;
77};
78template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
79 typedef typename std::remove_pointer<T>::type Type;
80};
81
82template <typename T>
84 RotSelector rotSelector = LARGE_ROT) {
85
86 using D = typename TensorTypeExtractor<T>::Type;
87
88 FTensor::Index<'i', 3> i;
89 FTensor::Index<'j', 3> j;
90 FTensor::Index<'k', 3> k;
92
93 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
94 t_R(i, j) = t_kd(i, j);
95
97 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
98
99 if (rotSelector == SMALL_ROT) {
100 t_R(i, j) += t_Omega(i, j);
101 return t_R;
102 }
103
104 const auto angle =
105 sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
106 const auto a = sin(angle) / angle;
107 t_R(i, j) += a * t_Omega(i, j);
108 if (rotSelector == MODERATE_ROT)
109 return t_R;
110
111 const auto ss_2 = sin(angle / 2.) / angle;
112 const auto b = 2. * ss_2 * ss_2;
113 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
114
115 return t_R;
116}
117
118template <typename T>
120 RotSelector rotSelector = LARGE_ROT) {
121
122 using D = typename TensorTypeExtractor<T>::Type;
123
124 FTensor::Index<'i', 3> i;
125 FTensor::Index<'j', 3> j;
126 FTensor::Index<'k', 3> k;
127 FTensor::Index<'l', 3> l;
128
130
131 const auto angle =
132 sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
133 if (rotSelector == SMALL_ROT) {
134 t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
135 return t_diff_R;
136 }
137
138 const auto ss = sin(angle);
139 const auto a = ss / angle;
140 t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
141
143 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
144
145 const auto angle2 = angle * angle;
146 const auto cc = cos(angle);
147 const auto diff_a = (angle * cc - ss) / angle2;
148 t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
149 if (rotSelector == MODERATE_ROT)
150 return t_diff_R;
151
152 const auto ss_2 = sin(angle / 2.);
153 const auto cc_2 = cos(angle / 2.);
154 const auto b = 2. * ss_2 * ss_2 / angle2;
155 t_diff_R(i, j, k) +=
156 b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
157 FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
158 const auto diff_b =
159 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
160 t_diff_R(i, j, k) +=
161 diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
162
163 return t_diff_R;
164}
165
166template <typename T>
168 RotSelector rotSelector = LARGE_ROT) {
169
170 using D = typename TensorTypeExtractor<T>::Type;
171
173 FTensor::Index<'i', 3> i;
174
175 constexpr double eps = 1e-10;
176 for (int l = 0; l != 3; ++l) {
178 t_omega_c(i) = t_omega(i);
179 t_omega_c(l) += std::complex<double>(0, eps);
180 auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
181 for (int i = 0; i != 3; ++i) {
182 for (int j = 0; j != 3; ++j) {
183 for (int k = 0; k != 3; ++k) {
184 t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
185 }
186 }
187 }
188 }
189
190 return t_diff2_R;
191}
192
193struct isEq {
194 static inline auto check(const double &a, const double &b) {
195 constexpr double eps = std::numeric_limits<float>::epsilon();
196 return std::abs(a - b) / absMax < eps;
197 }
198 static double absMax;
199};
200
201double isEq::absMax = 1;
202
203inline auto is_eq(const double &a, const double &b) {
204 return isEq::check(a, b);
205};
206
207template <int DIM> inline auto get_uniq_nb(double *ptr) {
208 std::array<double, DIM> tmp;
209 std::copy(ptr, &ptr[DIM], tmp.begin());
210 std::sort(tmp.begin(), tmp.end());
211 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
212 constexpr double eps = std::numeric_limits<float>::epsilon();
213 isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
214 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
215}
216
217template <int DIM>
220
222 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
223
224 int i = 0, j = 1, k = 2;
225
226 if (is_eq(eig(0), eig(1))) {
227 i = 0;
228 j = 2;
229 k = 1;
230 } else if (is_eq(eig(0), eig(2))) {
231 i = 0;
232 j = 1;
233 k = 2;
234 } else if (is_eq(eig(1), eig(2))) {
235 i = 1;
236 j = 0;
237 k = 2;
238 }
239
241 eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
242
243 eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
244
245 eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
246
247 FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
248
249 {
250 FTensor::Index<'i', DIM> i;
251 FTensor::Index<'j', DIM> j;
252 eig(i) = eig_c(i);
253 eigen_vec(i, j) = eigen_vec_c(i, j);
254 }
255}
256
258 EntityType type,
259 EntData &data) {
261
263 auto t_L = symm_L_tensor();
264
265 int nb_integration_pts = getGaussPts().size2();
266 FTensor::Index<'i', 3> i;
267 FTensor::Index<'j', 3> j;
268 FTensor::Index<'k', 3> k;
269 FTensor::Index<'l', 3> l;
270 FTensor::Index<'m', 3> m;
271 FTensor::Index<'l', 3> n;
272
273 dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
274 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
275 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
276 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
277 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
278 dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts, false);
279
280 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
281 dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
282 dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
283 dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
284 false);
285 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
286 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
287 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
288 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
289 dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
290 dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
291 dataAtPts->nbUniq.resize(nb_integration_pts, false);
292
293 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
294
295 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
296 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
297 auto t_h_dlog_u =
298 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
299 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
300 auto t_levi_kirchoff_dP =
301 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
302 auto t_levi_kirchoff_domega =
303 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
304 auto t_approx_P_adjont_dstretch =
305 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
306 auto t_approx_P_adjont_log_du =
307 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
308 auto t_approx_P_adjont_log_du_dP =
309 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
310 auto t_approx_P_adjont_log_du_domega =
311 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
312 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
313 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
314 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
315 auto t_log_u =
316 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
317 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
318 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
319 auto t_diff_u =
320 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
321 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
322 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
323
324 auto &nbUniq = dataAtPts->nbUniq;
325 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
326 auto t_log_stretch_total =
327 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
329
330 for (int gg = 0; gg != nb_integration_pts; ++gg) {
331
337 FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
338
339 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
340
341 auto calulate_rotation = [&]() {
342 auto t0_diff =
345 t_diff_R(i, j, k) = t0_diff(i, j, k);
346 t_R(i, j) = t0(i, j);
347 };
348
349 auto calulate_streach = [&]() {
350 eigen_vec(i, j) = t_log_u(i, j);
351 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
352 // rare case when two eigen values are equal
353 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
354 if (nbUniq[gg] < 3) {
355 sort_eigen_vals<3>(eig, eigen_vec);
356 }
357 t_eigen_vals(i) = eig(i);
358 t_eigen_vecs(i, j) = eigen_vec(i, j);
359 auto t_u_tmp =
360 EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
361 t_u(i, j) = t_u_tmp(i, j);
362 auto t_diff_u_tmp =
363 EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f,
364 EshelbianCore::d_f, nbUniq[gg]);
365 t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
366 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
367 };
368
369 calulate_rotation();
370 calulate_streach();
371
373 case LARGE_ROT:
374
375 // rotation can be large, moderate or small
376 t_Ru(i, m) = t_R(i, l) * t_u(l, m);
377 t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
378 t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
379 t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
380
381 // conservation of angular momentum at current configuration
382 t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
383 t_approx_P_adjont_dstretch(l, m) =
384 t_approx_P_intermidiate(i, m) * t_R(i, l);
385
386 t_levi_kirchoff(k) =
387 levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
388 t_levi_kirchoff_dP(i, j, k) =
389 (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
390 t_levi_kirchoff_domega(k, n) =
391 levi_civita(l, m, k) *
392 (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
393
394 t_approx_P_adjont_log_du(L) =
395 t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
396 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
397 t_approx_P_adjont_log_du_domega(n, L) =
398 t_Ldiff_u(l, m, L) *
399 (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
400
401 break;
402 case MODERATE_ROT:
403
404 {
405
406 // assume small current rotation
408 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
409 t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
410 t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
411 t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
412 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
413
414 // conservation of angular momentum at intermediate configuration
415 t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
416 t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
417
418 t_levi_kirchoff(k) =
419 levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
420 t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
421 t_levi_kirchoff_domega(k, n) = 0;
422
423 t_approx_P_adjont_log_du(L) =
424 t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
425 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
426 t_approx_P_adjont_log_du_domega(n, L) = 0;
427
428 }
429
430 break;
431
432 case SMALL_ROT:
433
434 {
435
436 // assume small current rotation
438 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
439 t_h(i, j) = t_Omega(i, j) + t_u(i, j);
440 t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
441 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
442
443 // conservation of angular momentum at reference condition
444 t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
445 t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
446 t_levi_kirchoff_domega(k, l) = 0;
447
448 t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
449 t_approx_P_adjont_log_du(L) =
450 t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
451 t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
452 t_approx_P_adjont_log_du_domega(m, L) = 0;
453
454 }
455
456 break;
457 }
458
460 t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
461 eigen_vec(i, j) = t_C_h1(i, j);
462 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
465 break;
467 for (int ii = 0; ii != 3; ++ii) {
468 eig(ii) = std::max(eig(ii),
470 std::numeric_limits<double>::min_exponent)));
471 }
472 break;
473 }
474
475 auto t_log_u2_h1_tmp =
477
479 case LARGE_ROT:
480 case MODERATE_ROT:
481 t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
482 break;
483 case SMALL_ROT:
484 t_log_stretch_total(i, j) = t_log_u(i, j);
485 break;
486 };
487
488 ++t_h;
489 ++t_h_domega;
490 ++t_h_dlog_u;
491 ++t_levi_kirchoff;
492 ++t_levi_kirchoff_dP;
493 ++t_levi_kirchoff_domega;
494 ++t_approx_P_adjont_dstretch;
495 ++t_approx_P_adjont_log_du;
496 ++t_approx_P_adjont_log_du_dP;
497 ++t_approx_P_adjont_log_du_domega;
498 ++t_approx_P;
499 ++t_R;
500 ++t_diff_R;
501 ++t_log_u;
502 ++t_u;
503 ++t_diff_u;
504 ++t_eigen_vals;
505 ++t_eigen_vecs;
506 ++t_omega;
507 ++t_grad_h1;
508 ++t_log_stretch_total;
509 }
510
512}
513
516 int nb_dofs = data.getIndices().size();
517 int nb_integration_pts = data.getN().size1();
518 auto v = getVolume();
519 auto t_w = getFTensor0IntegrationWeight();
520 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
521 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
522 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
523 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
524 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
525 dataAtPts->wL2DotDotAtPts.clear();
526 }
527 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
528
529 int nb_base_functions = data.getN().size2();
530 auto t_row_base_fun = data.getFTensor0N();
531
532 FTensor::Index<'i', 3> i;
533 auto get_ftensor1 = [](auto &v) {
535 &v[2]);
536 };
537
538 for (int gg = 0; gg != nb_integration_pts; ++gg) {
539 double a = v * t_w;
540 auto t_nf = get_ftensor1(nF);
541 int bb = 0;
542 for (; bb != nb_dofs / 3; ++bb) {
543 t_nf(i) += a * t_row_base_fun * t_div_P(i);
544 t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
545 t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
546 ++t_nf;
547 ++t_row_base_fun;
548 }
549 for (; bb != nb_base_functions; ++bb)
550 ++t_row_base_fun;
551 ++t_w;
552 ++t_div_P;
553 ++t_s_dot_w;
554 ++t_s_dot_dot_w;
555 }
556
558}
559
562 int nb_dofs = data.getIndices().size();
563 int nb_integration_pts = getGaussPts().size2();
564 auto v = getVolume();
565 auto t_w = getFTensor0IntegrationWeight();
566 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
567 int nb_base_functions = data.getN().size2();
568 auto t_row_base_fun = data.getFTensor0N();
569 FTensor::Index<'i', 3> i;
570 FTensor::Index<'j', 3> j;
571 FTensor::Index<'k', 3> k;
572 auto get_ftensor1 = [](auto &v) {
574 &v[2]);
575 };
576 for (int gg = 0; gg != nb_integration_pts; ++gg) {
577 double a = v * t_w;
578 auto t_nf = get_ftensor1(nF);
579 int bb = 0;
580 for (; bb != nb_dofs / 3; ++bb) {
581 t_nf(k) += (a * t_row_base_fun) * t_levi_kirchoff(k);
582 ++t_nf;
583 ++t_row_base_fun;
584 }
585 for (; bb != nb_base_functions; ++bb) {
586 ++t_row_base_fun;
587 }
588 ++t_w;
589 ++t_levi_kirchoff;
590 }
592}
593
596
597 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
599 } else {
600 if (polyConvex) {
602 } else {
604 }
605 }
606
608}
609
612
614 auto t_L = symm_L_tensor();
615
616 int nb_dofs = data.getIndices().size();
617 int nb_integration_pts = data.getN().size1();
618 auto v = getVolume();
619 auto t_w = getFTensor0IntegrationWeight();
620 auto t_approx_P_adjont_log_du =
621 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
622 auto t_log_streach_h1 =
623 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
624 auto t_dot_log_u =
625 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
626
627 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
628
629 FTensor::Index<'i', 3> i;
630 FTensor::Index<'j', 3> j;
631 FTensor::Index<'k', 3> k;
632 FTensor::Index<'l', 3> l;
633 auto get_ftensor2 = [](auto &v) {
635 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
636 };
637
638 int nb_base_functions = data.getN().size2();
639 auto t_row_base_fun = data.getFTensor0N();
640 for (int gg = 0; gg != nb_integration_pts; ++gg) {
641 double a = v * t_w;
642 auto t_nf = get_ftensor2(nF);
643
645 t_T(i, j) =
646 t_D(i, j, k, l) * (t_log_streach_h1(k, l) + alphaU * t_dot_log_u(k, l));
648 t_residual(L) =
649 a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
650
651 int bb = 0;
652 for (; bb != nb_dofs / 6; ++bb) {
653 t_nf(L) += t_row_base_fun * t_residual(L);
654 ++t_nf;
655 ++t_row_base_fun;
656 }
657 for (; bb != nb_base_functions; ++bb)
658 ++t_row_base_fun;
659
660 ++t_w;
661 ++t_approx_P_adjont_log_du;
662 ++t_dot_log_u;
663 ++t_log_streach_h1;
664 }
666}
667
670
672 auto t_L = symm_L_tensor();
673
674 int nb_dofs = data.getIndices().size();
675 int nb_integration_pts = data.getN().size1();
676 auto v = getVolume();
677 auto t_w = getFTensor0IntegrationWeight();
678 auto t_approx_P_adjont_log_du =
679 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
680 auto t_log_streach_h1 =
681 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
682 auto t_dot_log_u =
683 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
684
685 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
686
687 FTensor::Index<'i', 3> i;
688 FTensor::Index<'j', 3> j;
689 FTensor::Index<'k', 3> k;
690 FTensor::Index<'l', 3> l;
691 auto get_ftensor2 = [](auto &v) {
693 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
694 };
695
696 constexpr double nohat_k = 1. / 4;
697 constexpr double hat_k = 1. / 8;
698 double mu = dataAtPts->mu;
699 double lambda = dataAtPts->lambda;
700
701 constexpr double third = boost::math::constants::third<double>();
703 auto t_diff_deviator = diff_deviator(diff_tensor());
704
705 int nb_base_functions = data.getN().size2();
706 auto t_row_base_fun = data.getFTensor0N();
707 for (int gg = 0; gg != nb_integration_pts; ++gg) {
708 double a = v * t_w;
709 auto t_nf = get_ftensor2(nF);
710
711 double log_det = t_log_streach_h1(i, i);
712 double log_det2 = log_det * log_det;
714 t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
715 double dev_norm2 = t_dev(i, j) * t_dev(i, j);
716
718 auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
719 auto B = lambda * std::exp(hat_k * log_det2) * log_det;
720 t_T(i, j) =
721
722 A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
723
724 +
725
726 B * t_kd(i, j)
727
728 +
729
730 alphaU * t_D(i, j, k, l) * t_log_streach_h1(k, l);
731
733 t_residual(L) =
734 a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
735
736 int bb = 0;
737 for (; bb != nb_dofs / 6; ++bb) {
738 t_nf(L) += t_row_base_fun * t_residual(L);
739 ++t_nf;
740 ++t_row_base_fun;
741 }
742 for (; bb != nb_base_functions; ++bb)
743 ++t_row_base_fun;
744
745 ++t_w;
746 ++t_approx_P_adjont_log_du;
747 ++t_dot_log_u;
748 ++t_log_streach_h1;
749 }
751}
752
755
757 auto t_L = symm_L_tensor();
758
759 int nb_dofs = data.getIndices().size();
760 int nb_integration_pts = data.getN().size1();
761 auto v = getVolume();
762 auto t_w = getFTensor0IntegrationWeight();
763 auto t_approx_P_adjont_log_du =
764 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
765 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
766 auto t_dot_log_u =
767 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
768 auto t_diff_u =
769 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
770
771 FTensor::Index<'i', 3> i;
772 FTensor::Index<'j', 3> j;
773 FTensor::Index<'k', 3> k;
774 FTensor::Index<'l', 3> l;
775 auto get_ftensor2 = [](auto &v) {
777 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
778 };
779
780 int nb_base_functions = data.getN().size2();
781 auto t_row_base_fun = data.getFTensor0N();
782 for (int gg = 0; gg != nb_integration_pts; ++gg) {
783 double a = v * t_w;
784 auto t_nf = get_ftensor2(nF);
785
787 t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
788
790 t_viscous_P(i, j) =
791 alphaU *
792 t_dot_log_u(i,
793 j); // That is chit, should be split on axiator and deviator
794
796 t_residual(L) =
797 a * (t_approx_P_adjont_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
798 t_L(i, j, L) * t_viscous_P(i, j));
799
800 int bb = 0;
801 for (; bb != nb_dofs / 6; ++bb) {
802 t_nf(L) += t_row_base_fun * t_residual(L);
803 ++t_nf;
804 ++t_row_base_fun;
805 }
806 for (; bb != nb_base_functions; ++bb)
807 ++t_row_base_fun;
808
809 ++t_w;
810 ++t_approx_P_adjont_log_du;
811 ++t_P;
812 ++t_dot_log_u;
813 ++t_diff_u;
814 }
816}
817
820 int nb_dofs = data.getIndices().size();
821 int nb_integration_pts = data.getN().size1();
822 auto v = getVolume();
823 auto t_w = getFTensor0IntegrationWeight();
824 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
825 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
826
827 int nb_base_functions = data.getN().size2() / 3;
828 auto t_row_base_fun = data.getFTensor1N<3>();
829 FTensor::Index<'i', 3> i;
830 FTensor::Index<'j', 3> j;
831 FTensor::Index<'k', 3> k;
832 FTensor::Index<'m', 3> m;
833
834 auto get_ftensor1 = [](auto &v) {
836 &v[2]);
837 };
838
839 for (int gg = 0; gg != nb_integration_pts; ++gg) {
840 double a = v * t_w;
841 auto t_nf = get_ftensor1(nF);
842
843 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
845 t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
846
847 int bb = 0;
848 for (; bb != nb_dofs / 3; ++bb) {
849 t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
850 ++t_nf;
851 ++t_row_base_fun;
852 }
853
854 for (; bb != nb_base_functions; ++bb)
855 ++t_row_base_fun;
856
857 ++t_w;
858 ++t_h;
859 ++t_omega_dot;
860 }
861
863}
864
867 int nb_dofs = data.getIndices().size();
868 int nb_integration_pts = data.getN().size1();
869 auto v = getVolume();
870 auto t_w = getFTensor0IntegrationWeight();
871 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
872 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
873
874 int nb_base_functions = data.getN().size2() / 9;
875 auto t_row_base_fun = data.getFTensor2N<3, 3>();
876 FTensor::Index<'i', 3> i;
877 FTensor::Index<'j', 3> j;
878 FTensor::Index<'k', 3> k;
879 FTensor::Index<'m', 3> m;
880
881 auto get_ftensor0 = [](auto &v) {
883 };
884
885 for (int gg = 0; gg != nb_integration_pts; ++gg) {
886 double a = v * t_w;
887 auto t_nf = get_ftensor0(nF);
888
889 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
891 t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
892
893 int bb = 0;
894 for (; bb != nb_dofs; ++bb) {
895 t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
896 ++t_nf;
897 ++t_row_base_fun;
898 }
899 for (; bb != nb_base_functions; ++bb) {
900 ++t_row_base_fun;
901 }
902 ++t_w;
903 ++t_h;
904 ++t_omega_dot;
905 }
906
908}
909
912 int nb_dofs = data.getIndices().size();
913 int nb_integration_pts = data.getN().size1();
914 auto v = getVolume();
915 auto t_w = getFTensor0IntegrationWeight();
916 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
917 int nb_base_functions = data.getN().size2() / 3;
918 auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
919 FTensor::Index<'i', 3> i;
920 auto get_ftensor1 = [](auto &v) {
922 &v[2]);
923 };
924
925 for (int gg = 0; gg != nb_integration_pts; ++gg) {
926 double a = v * t_w;
927 auto t_nf = get_ftensor1(nF);
928 int bb = 0;
929 for (; bb != nb_dofs / 3; ++bb) {
930 double div_row_base = t_row_diff_base_fun(i, i);
931 t_nf(i) += a * div_row_base * t_w_l2(i);
932 ++t_nf;
933 ++t_row_diff_base_fun;
934 }
935 for (; bb != nb_base_functions; ++bb) {
936 ++t_row_diff_base_fun;
937 }
938 ++t_w;
939 ++t_w_l2;
940 }
941
943}
944
947 // get entity of face
948 EntityHandle fe_ent = getFEEntityHandle();
949 // interate over all boundary data
950 for (auto &bc : (*bcDispPtr)) {
951 // check if finite element entity is part of boundary condition
952 if (bc.faces.find(fe_ent) != bc.faces.end()) {
953 int nb_dofs = data.getIndices().size();
954 int nb_integration_pts = data.getN().size1();
955 auto t_normal = getFTensor1Normal();
956 auto t_w = getFTensor0IntegrationWeight();
957 int nb_base_functions = data.getN().size2() / 3;
958 auto t_row_base_fun = data.getFTensor1N<3>();
959
960 FTensor::Index<'i', 3> i;
961 FTensor::Index<'j', 3> j;
962 auto get_ftensor1 = [](auto &v) {
964 &v[2]);
965 };
966
967 double scale = 1;
968 for (auto &sm : scalingMethodsVec) {
969 scale *= sm->getScale(getFEMethod()->ts_t);
970 }
971
972 // get bc data
973 FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
974 t_bc_disp(i) *= scale;
975
976 for (int gg = 0; gg != nb_integration_pts; ++gg) {
978 t_bc_res(i) = t_bc_disp(i);
979 auto t_nf = get_ftensor1(nF);
980 int bb = 0;
981 for (; bb != nb_dofs / 3; ++bb) {
982 t_nf(i) -=
983 t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
984 ++t_nf;
985 ++t_row_base_fun;
986 }
987 for (; bb != nb_base_functions; ++bb)
988 ++t_row_base_fun;
989
990 ++t_w;
991 }
992 }
993 }
995}
996
999 // get entity of face
1000 EntityHandle fe_ent = getFEEntityHandle();
1001 // interate over all boundary data
1002 for (auto &bc : (*bcRotPtr)) {
1003 // check if finite element entity is part of boundary condition
1004 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1005 int nb_dofs = data.getIndices().size();
1006 int nb_integration_pts = data.getN().size1();
1007 auto t_normal = getFTensor1Normal();
1008 auto t_w = getFTensor0IntegrationWeight();
1009
1010 int nb_base_functions = data.getN().size2() / 3;
1011 auto t_row_base_fun = data.getFTensor1N<3>();
1012 FTensor::Index<'i', 3> i;
1013 FTensor::Index<'j', 3> j;
1014 FTensor::Index<'k', 3> k;
1015 auto get_ftensor1 = [](auto &v) {
1017 &v[2]);
1018 };
1019
1020 // get bc data
1021 FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
1022
1023 double theta = bc.theta;
1024 theta *= getFEMethod()->ts_t;
1025
1027 const double a = sqrt(t_normal(i) * t_normal(i));
1028 t_omega(i) = t_normal(i) * (theta / a);
1029 auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
1030 auto t_coords = getFTensor1CoordsAtGaussPts();
1031
1032 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1034 t_delta(i) = t_center(i) - t_coords(i);
1036 t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
1037
1038 auto t_nf = get_ftensor1(nF);
1039 int bb = 0;
1040 for (; bb != nb_dofs / 3; ++bb) {
1041 t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
1042 ++t_nf;
1043 ++t_row_base_fun;
1044 }
1045 for (; bb != nb_base_functions; ++bb)
1046 ++t_row_base_fun;
1047
1048 ++t_w;
1049 ++t_coords;
1050 }
1051 }
1052 }
1054}
1055
1058
1059 struct FaceRule {
1060 int operator()(int p_row, int p_col, int p_data) const {
1061 return 2 * (p_data + 1);
1062 }
1063 };
1064
1065 if (ts_ctx == CTX_TSSETIFUNCTION) {
1066
1067 // Loop boundary elements with traction boundary conditions
1070 {HDIV});
1071 fe.getOpPtrVector().push_back(
1072 new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
1073 fe.getRuleHook = FaceRule();
1074 fe.ts_t = ts_t;
1075 CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
1076 fe, nullptr, MF_ZERO,
1077 this->getCacheWeakPtr());
1078 }
1079
1081}
1082
1085
1086 FTensor::Index<'i', 3> i;
1087 FTensor::Index<'j', 3> j;
1088
1089 auto t_normal = getFTensor1Normal();
1090 const double nrm2 = sqrt(t_normal(i) * t_normal(i));
1091 FTensor::Tensor1<double, 3> t_unit_normal;
1092 t_unit_normal(i) = t_normal(i) / nrm2;
1093 int nb_dofs = data.getFieldData().size();
1094 int nb_integration_pts = data.getN().size1();
1095 int nb_base_functions = data.getN().size2() / 3;
1096 double ts_t = getFEMethod()->ts_t;
1097
1098 auto integrate_matrix = [&]() {
1100
1101 auto t_w = getFTensor0IntegrationWeight();
1102 matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
1103 matPP.clear();
1104
1105 auto t_row_base_fun = data.getFTensor1N<3>();
1106 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1107
1108 int rr = 0;
1109 for (; rr != nb_dofs / 3; ++rr) {
1110 const double a = t_row_base_fun(i) * t_unit_normal(i);
1111 auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
1112 for (int cc = 0; cc != nb_dofs / 3; ++cc) {
1113 const double b = t_col_base_fun(i) * t_unit_normal(i);
1114 matPP(rr, cc) += t_w * a * b;
1115 ++t_col_base_fun;
1116 }
1117 ++t_row_base_fun;
1118 }
1119
1120 for (; rr != nb_base_functions; ++rr)
1121 ++t_row_base_fun;
1122
1123 ++t_w;
1124 }
1125
1127 };
1128
1129 auto integrate_rhs = [&](auto &bc) {
1131
1132 auto t_w = getFTensor0IntegrationWeight();
1133 vecPv.resize(3, nb_dofs / 3, false);
1134 vecPv.clear();
1135
1136 auto t_row_base_fun = data.getFTensor1N<3>();
1137 double ts_t = getFEMethod()->ts_t;
1138
1139 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1140 int rr = 0;
1141 for (; rr != nb_dofs / 3; ++rr) {
1142 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1143 for (int dd = 0; dd != 3; ++dd)
1144 if (bc.flags[dd])
1145 vecPv(dd, rr) += t * bc.vals[dd];
1146 ++t_row_base_fun;
1147 }
1148 for (; rr != nb_base_functions; ++rr)
1149 ++t_row_base_fun;
1150 ++t_w;
1151 }
1153 };
1154
1155 auto integrate_rhs_cook = [&](auto &bc) {
1157
1158 vecPv.resize(3, nb_dofs / 3, false);
1159 vecPv.clear();
1160
1161 auto t_w = getFTensor0IntegrationWeight();
1162 auto t_row_base_fun = data.getFTensor1N<3>();
1163 auto t_coords = getFTensor1CoordsAtGaussPts();
1164
1165 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1166
1167 auto calc_tau = [](double y) {
1168 y -= 44;
1169 y /= (60 - 44);
1170 return -y * (y - 1) / 0.25;
1171 };
1172
1173 const double tau = calc_tau(t_coords(1));
1174
1175 int rr = 0;
1176 for (; rr != nb_dofs / 3; ++rr) {
1177 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1178 for (int dd = 0; dd != 3; ++dd)
1179 if (bc.flags[dd])
1180 vecPv(dd, rr) += t * tau * bc.vals[dd];
1181 ++t_row_base_fun;
1182 }
1183
1184 for (; rr != nb_base_functions; ++rr)
1185 ++t_row_base_fun;
1186 ++t_w;
1187 ++t_coords;
1188 }
1190 };
1191
1192 // get entity of face
1193 EntityHandle fe_ent = getFEEntityHandle();
1194 for (auto &bc : *(bcFe.bcData)) {
1195 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1196
1197 int nb_dofs = data.getFieldData().size();
1198 if (nb_dofs) {
1199
1200 CHKERR integrate_matrix();
1201 if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
1202 CHKERR integrate_rhs_cook(bc);
1203 else
1204 CHKERR integrate_rhs(bc);
1205
1206 const auto info =
1207 lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1208 nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1209 if (info > 0)
1210 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1211 "The leading minor of order %i is "
1212 "not positive; definite;\nthe "
1213 "solution could not be computed",
1214 info);
1215
1216 for (int dd = 0; dd != 3; ++dd)
1217 if (bc.flags[dd])
1218 for (int rr = 0; rr != nb_dofs / 3; ++rr)
1219 data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1220 }
1221 }
1222 }
1223
1225}
1226
1228 EntData &col_data) {
1230 int nb_integration_pts = row_data.getN().size1();
1231 int row_nb_dofs = row_data.getIndices().size();
1232 int col_nb_dofs = col_data.getIndices().size();
1233 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1235 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
1236 };
1237 FTensor::Index<'i', 3> i;
1238 auto v = getVolume();
1239 auto t_w = getFTensor0IntegrationWeight();
1240 int row_nb_base_functions = row_data.getN().size2();
1241 auto t_row_base_fun = row_data.getFTensor0N();
1242 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1243 double a = v * t_w;
1244 int rr = 0;
1245 for (; rr != row_nb_dofs / 3; ++rr) {
1246 auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
1247 auto t_m = get_ftensor1(K, 3 * rr, 0);
1248 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1249 double div_col_base = t_col_diff_base_fun(i, i);
1250 t_m(i) += a * t_row_base_fun * div_col_base;
1251 ++t_m;
1252 ++t_col_diff_base_fun;
1253 }
1254 ++t_row_base_fun;
1255 }
1256 for (; rr != row_nb_base_functions; ++rr)
1257 ++t_row_base_fun;
1258 ++t_w;
1259 }
1261}
1262
1264 EntData &col_data) {
1266
1267 if (alphaW < std::numeric_limits<double>::epsilon() &&
1268 alphaRho < std::numeric_limits<double>::epsilon())
1270
1271 const int nb_integration_pts = row_data.getN().size1();
1272 const int row_nb_dofs = row_data.getIndices().size();
1273 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1275 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1276
1277 );
1278 };
1279 FTensor::Index<'i', 3> i;
1280
1281 auto v = getVolume();
1282 auto t_w = getFTensor0IntegrationWeight();
1283
1284 int row_nb_base_functions = row_data.getN().size2();
1285 auto t_row_base_fun = row_data.getFTensor0N();
1286
1287 double ts_scale = alphaW * getTSa();
1288 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1289 ts_scale += alphaRho * getTSaa();
1290
1291 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1292 double a = v * t_w * ts_scale;
1293
1294 int rr = 0;
1295 for (; rr != row_nb_dofs / 3; ++rr) {
1296
1297 auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1298 auto t_m = get_ftensor1(K, 3 * rr, 0);
1299 for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1300 const double b = a * t_row_base_fun * t_col_base_fun;
1301 t_m(i) -= b;
1302 ++t_m;
1303 ++t_col_base_fun;
1304 }
1305
1306 ++t_row_base_fun;
1307 }
1308
1309 for (; rr != row_nb_base_functions; ++rr)
1310 ++t_row_base_fun;
1311
1312 ++t_w;
1313 }
1314
1316}
1317
1319 EntData &col_data) {
1321
1323
1324 int nb_integration_pts = row_data.getN().size1();
1325 int row_nb_dofs = row_data.getIndices().size();
1326 int col_nb_dofs = col_data.getIndices().size();
1327 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1329
1330 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1331
1332 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1333
1334 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1335
1336 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1337
1338 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1339
1340 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1341 };
1342
1343 FTensor::Index<'i', 3> i;
1344 FTensor::Index<'j', 3> j;
1345
1346 auto v = getVolume();
1347 auto t_w = getFTensor0IntegrationWeight();
1348 auto t_approx_P_adjont_log_du_dP =
1349 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1350 int row_nb_base_functions = row_data.getN().size2();
1351 auto t_row_base_fun = row_data.getFTensor0N();
1352 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1353 double a = v * t_w;
1354 int rr = 0;
1355 for (; rr != row_nb_dofs / 6; ++rr) {
1356
1357 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1358 auto t_m = get_ftensor3(K, 6 * rr, 0);
1359
1360 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1361 t_m(L, i) +=
1362 a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(j)) *
1363 t_row_base_fun;
1364 ++t_col_base_fun;
1365 ++t_m;
1366 }
1367
1368 ++t_row_base_fun;
1369 }
1370 for (; rr != row_nb_base_functions; ++rr)
1371 ++t_row_base_fun;
1372 ++t_w;
1373 ++t_approx_P_adjont_log_du_dP;
1374 }
1376}
1377
1379 EntData &col_data) {
1381
1383
1384 int nb_integration_pts = row_data.getN().size1();
1385 int row_nb_dofs = row_data.getIndices().size();
1386 int col_nb_dofs = col_data.getIndices().size();
1387 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1389 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1390 &m(r + 5, c));
1391 };
1392
1393 FTensor::Index<'i', 3> i;
1394 FTensor::Index<'j', 3> j;
1395
1396 auto v = getVolume();
1397 auto t_w = getFTensor0IntegrationWeight();
1398 auto t_approx_P_adjont_log_du_dP =
1399 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1400
1401 int row_nb_base_functions = row_data.getN().size2();
1402 auto t_row_base_fun = row_data.getFTensor0N();
1403 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1404 double a = v * t_w;
1405 int rr = 0;
1406 for (; rr != row_nb_dofs / 6; ++rr) {
1407 auto t_m = get_ftensor2(K, 6 * rr, 0);
1408 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1409 for (int cc = 0; cc != col_nb_dofs; ++cc) {
1410 t_m(L) +=
1411 a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
1412 t_row_base_fun;
1413 ++t_m;
1414 ++t_col_base_fun;
1415 }
1416 ++t_row_base_fun;
1417 }
1418 for (; rr != row_nb_base_functions; ++rr)
1419 ++t_row_base_fun;
1420 ++t_w;
1421 ++t_approx_P_adjont_log_du_dP;
1422 }
1424}
1425
1427 EntData &col_data) {
1429 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1430 CHKERR integratePiola(row_data, col_data);
1431 } else {
1432 if (polyConvex) {
1433 CHKERR integratePolyconvexHencky(row_data, col_data);
1434 } else {
1435 CHKERR integrateHencky(row_data, col_data);
1436 }
1437 }
1439}
1440
1442 EntData &col_data) {
1444
1447 auto t_L = symm_L_tensor();
1448 auto t_diff = diff_tensor();
1449
1450 int nb_integration_pts = row_data.getN().size1();
1451 int row_nb_dofs = row_data.getIndices().size();
1452 int col_nb_dofs = col_data.getIndices().size();
1453
1454 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1456 size_symm>(
1457
1458 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1459 &m(r + 0, c + 4), &m(r + 0, c + 5),
1460
1461 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1462 &m(r + 1, c + 4), &m(r + 1, c + 5),
1463
1464 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1465 &m(r + 2, c + 4), &m(r + 2, c + 5),
1466
1467 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1468 &m(r + 3, c + 4), &m(r + 3, c + 5),
1469
1470 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1471 &m(r + 4, c + 4), &m(r + 4, c + 5),
1472
1473 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1474 &m(r + 5, c + 4), &m(r + 5, c + 5)
1475
1476 );
1477 };
1478 FTensor::Index<'i', 3> i;
1479 FTensor::Index<'j', 3> j;
1480 FTensor::Index<'k', 3> k;
1481 FTensor::Index<'l', 3> l;
1482 FTensor::Index<'m', 3> m;
1483 FTensor::Index<'n', 3> n;
1484
1485 auto v = getVolume();
1486 auto t_w = getFTensor0IntegrationWeight();
1487
1488 auto t_approx_P_adjont_dstretch =
1489 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1490 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1491 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1492 auto &nbUniq = dataAtPts->nbUniq;
1493
1494 int row_nb_base_functions = row_data.getN().size2();
1495 auto t_row_base_fun = row_data.getFTensor0N();
1496
1497 auto get_dP = [&]() {
1498 dP.resize(size_symm * size_symm, nb_integration_pts, false);
1499 auto ts_a = getTSa();
1500
1501 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1503 t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
1504 (t_L(i, j, L) *
1505 ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
1506
1507
1509 auto t_approx_P_adjont_dstretch =
1510 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1511 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1512 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1513 auto &nbUniq = dataAtPts->nbUniq;
1514
1515 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1516 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1517
1518 // Work of symmetric tensor on undefined tensor is equal to the work
1519 // of the symmetric part of it
1521 t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
1522 t_approx_P_adjont_dstretch(j, i));
1523 t_sym(i, j) /= 2.0;
1524 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1525 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1526 EshelbianCore::dd_f, t_sym, nbUniq[gg]);
1527 t_dP(L, J) = t_L(i, j, L) *
1528 ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
1529 t_L(k, l, J)) /
1530 2. +
1531 t_dP_tmp(L, J);
1532
1533 ++t_dP;
1534 ++t_approx_P_adjont_dstretch;
1535 ++t_eigen_vals;
1536 ++t_eigen_vecs;
1537 }
1538 } else {
1539 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1540 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1541 t_dP(L, J) = t_dP_tmp(L, J);
1542 ++t_dP;
1543 }
1544 }
1545
1546 return getFTensor2FromMat<size_symm, size_symm>(dP);
1547 };
1548
1549 auto t_dP = get_dP();
1550 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1551 double a = v * t_w;
1552
1553 int rr = 0;
1554 for (; rr != row_nb_dofs / 6; ++rr) {
1555 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1556 auto t_m = get_ftensor2(K, 6 * rr, 0);
1557 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1558 const double b = a * t_row_base_fun * t_col_base_fun;
1559 t_m(L, J) += b * t_dP(L, J);
1560 ++t_m;
1561 ++t_col_base_fun;
1562 }
1563 ++t_row_base_fun;
1564 }
1565
1566 for (; rr != row_nb_base_functions; ++rr) {
1567 ++t_row_base_fun;
1568 }
1569
1570 ++t_w;
1571 ++t_dP;
1572 }
1574}
1575
1578 EntData &col_data) {
1580
1583 auto t_L = symm_L_tensor();
1584 auto t_diff = diff_tensor();
1585
1586 int nb_integration_pts = row_data.getN().size1();
1587 int row_nb_dofs = row_data.getIndices().size();
1588 int col_nb_dofs = col_data.getIndices().size();
1589
1590 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1592 size_symm>(
1593
1594 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1595 &m(r + 0, c + 4), &m(r + 0, c + 5),
1596
1597 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1598 &m(r + 1, c + 4), &m(r + 1, c + 5),
1599
1600 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1601 &m(r + 2, c + 4), &m(r + 2, c + 5),
1602
1603 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1604 &m(r + 3, c + 4), &m(r + 3, c + 5),
1605
1606 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1607 &m(r + 4, c + 4), &m(r + 4, c + 5),
1608
1609 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1610 &m(r + 5, c + 4), &m(r + 5, c + 5)
1611
1612 );
1613 };
1614 FTensor::Index<'i', 3> i;
1615 FTensor::Index<'j', 3> j;
1616 FTensor::Index<'k', 3> k;
1617 FTensor::Index<'l', 3> l;
1618 FTensor::Index<'m', 3> m;
1619 FTensor::Index<'n', 3> n;
1620
1621 auto v = getVolume();
1622 auto t_w = getFTensor0IntegrationWeight();
1623
1624 int row_nb_base_functions = row_data.getN().size2();
1625 auto t_row_base_fun = row_data.getFTensor0N();
1626
1627 auto get_dP = [&]() {
1628 dP.resize(size_symm * size_symm, nb_integration_pts, false);
1629 auto ts_a = getTSa();
1630
1631 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1632
1633 constexpr double nohat_k = 1. / 4;
1634 constexpr double hat_k = 1. / 8;
1635 double mu = dataAtPts->mu;
1636 double lambda = dataAtPts->lambda;
1637
1638 constexpr double third = boost::math::constants::third<double>();
1640 auto t_diff_deviator = diff_deviator(diff_tensor());
1641
1642 auto t_approx_P_adjont_dstretch =
1643 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1644 auto t_log_streach_h1 =
1645 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
1646 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1647 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1648 auto &nbUniq = dataAtPts->nbUniq;
1649
1650 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1651 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1652
1653 double log_det = t_log_streach_h1(i, i);
1654 double log_det2 = log_det * log_det;
1656 t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
1657 double dev_norm2 = t_dev(i, j) * t_dev(i, j);
1658
1659 auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
1660 auto B = lambda * std::exp(hat_k * log_det2) * log_det;
1661
1662 FTensor::Tensor2_symmetric<double, 3> t_A_diff, t_B_diff;
1663 t_A_diff(i, j) =
1664 (A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
1665 t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
1666 lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
1668 t_dT(i, j, k, l) =
1669 t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
1670
1671 +
1672
1673 A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
1674
1675 +
1676
1677 t_B_diff(i, j) * t_kd(k, l);
1678
1679 t_dP(L, J) = -t_L(i, j, L) *
1680 ((
1681
1682 t_dT(i, j, k, l)
1683
1684 +
1685
1686 (alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
1687
1688 )) *
1689 t_L(k, l, J));
1690
1691 // Work of symmetric tensor on undefined tensor is equal to the work
1692 // of the symmetric part of it
1695 t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
1696 t_approx_P_adjont_dstretch(j, i));
1697 t_sym(i, j) /= 2.0;
1698 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1699 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1700 EshelbianCore::dd_f, t_sym, nbUniq[gg]);
1701 t_dP(L, J) += t_L(i, j, L) *
1702 ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
1703 t_L(k, l, J)) /
1704 2.;
1705 }
1706
1707 ++t_dP;
1708 ++t_approx_P_adjont_dstretch;
1709 ++t_log_streach_h1;
1710 ++t_eigen_vals;
1711 ++t_eigen_vecs;
1712 }
1713
1714 return getFTensor2FromMat<size_symm, size_symm>(dP);
1715 };
1716
1717 auto t_dP = get_dP();
1718 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1719 double a = v * t_w;
1720
1721 int rr = 0;
1722 for (; rr != row_nb_dofs / 6; ++rr) {
1723 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1724 auto t_m = get_ftensor2(K, 6 * rr, 0);
1725 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1726 const double b = a * t_row_base_fun * t_col_base_fun;
1727 t_m(L, J) += b * t_dP(L, J);
1728 ++t_m;
1729 ++t_col_base_fun;
1730 }
1731 ++t_row_base_fun;
1732 }
1733
1734 for (; rr != row_nb_base_functions; ++rr) {
1735 ++t_row_base_fun;
1736 }
1737
1738 ++t_w;
1739 ++t_dP;
1740 }
1742}
1743
1745 EntData &col_data) {
1747
1750 auto t_L = symm_L_tensor();
1751 auto t_diff = diff_tensor();
1752
1753 int nb_integration_pts = row_data.getN().size1();
1754 int row_nb_dofs = row_data.getIndices().size();
1755 int col_nb_dofs = col_data.getIndices().size();
1756
1757 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1759 size_symm>(
1760
1761 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1762 &m(r + 0, c + 4), &m(r + 0, c + 5),
1763
1764 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1765 &m(r + 1, c + 4), &m(r + 1, c + 5),
1766
1767 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1768 &m(r + 2, c + 4), &m(r + 2, c + 5),
1769
1770 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1771 &m(r + 3, c + 4), &m(r + 3, c + 5),
1772
1773 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1774 &m(r + 4, c + 4), &m(r + 4, c + 5),
1775
1776 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1777 &m(r + 5, c + 4), &m(r + 5, c + 5)
1778
1779 );
1780 };
1781 FTensor::Index<'i', 3> i;
1782 FTensor::Index<'j', 3> j;
1783 FTensor::Index<'k', 3> k;
1784 FTensor::Index<'l', 3> l;
1785 FTensor::Index<'m', 3> m;
1786 FTensor::Index<'n', 3> n;
1787
1788 auto v = getVolume();
1789 auto t_w = getFTensor0IntegrationWeight();
1790
1791 auto get_dP = [&]() {
1792 dP.resize(size_symm * size_symm, nb_integration_pts, false);
1793 auto ts_a = getTSa();
1794
1795 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1796 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
1797 auto t_approx_P_adjont_dstretch =
1798 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1799 auto t_dot_log_u =
1800 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
1801 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1802 auto t_diff_u =
1803 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
1804 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1805 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1806 auto &nbUniq = dataAtPts->nbUniq;
1807
1808 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1809 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1810
1812 t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
1813
1815 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
1816 t_dP(L, J) =
1817 -t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
1818 // viscous stress derivative
1819 t_dP(L, J) -= (alphaU * ts_a) *
1820 (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
1821
1824 t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
1825 t_deltaP_sym(i, j) /= 2.0;
1826 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1827 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1828 EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
1829 t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
1830 }
1831
1832 ++t_P;
1833 ++t_dP;
1834 ++r_P_du;
1835 ++t_approx_P_adjont_dstretch;
1836 ++t_dot_log_u;
1837 ++t_u;
1838 ++t_diff_u;
1839 ++t_eigen_vals;
1840 ++t_eigen_vecs;
1841 }
1842
1843 return getFTensor2FromMat<size_symm, size_symm>(dP);
1844 };
1845
1846 int row_nb_base_functions = row_data.getN().size2();
1847 auto t_row_base_fun = row_data.getFTensor0N();
1848
1849 auto t_dP = get_dP();
1850 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1851 double a = v * t_w;
1852
1853 int rr = 0;
1854 for (; rr != row_nb_dofs / 6; ++rr) {
1855 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1856 auto t_m = get_ftensor2(K, 6 * rr, 0);
1857 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1858 const double b = a * t_row_base_fun * t_col_base_fun;
1859 t_m(L, J) += b * t_dP(L, J);
1860 ++t_m;
1861 ++t_col_base_fun;
1862 }
1863 ++t_row_base_fun;
1864 }
1865
1866 for (; rr != row_nb_base_functions; ++rr) {
1867 ++t_row_base_fun;
1868 }
1869
1870 ++t_w;
1871 ++t_dP;
1872 }
1874}
1875
1877 EntData &col_data) {
1879
1881 auto t_L = symm_L_tensor();
1882
1883 int row_nb_dofs = row_data.getIndices().size();
1884 int col_nb_dofs = col_data.getIndices().size();
1885 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1887
1888 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1889
1890 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1891
1892 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1893
1894 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1895
1896 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1897
1898 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1899
1900 );
1901 };
1902 FTensor::Index<'i', 3> i;
1903 FTensor::Index<'j', 3> j;
1904 FTensor::Index<'k', 3> k;
1905 FTensor::Index<'m', 3> m;
1906 FTensor::Index<'n', 3> n;
1907
1908 auto v = getVolume();
1909 auto t_w = getFTensor0IntegrationWeight();
1910 auto t_approx_P_adjont_log_du_domega =
1911 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
1912
1913 int row_nb_base_functions = row_data.getN().size2();
1914 auto t_row_base_fun = row_data.getFTensor0N();
1915
1916 int nb_integration_pts = row_data.getN().size1();
1917
1918 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1919 double a = v * t_w;
1920
1921 int rr = 0;
1922 for (; rr != row_nb_dofs / 6; ++rr) {
1923 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1924 auto t_m = get_ftensor3(K, 6 * rr, 0);
1925 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1926 double v = a * t_row_base_fun * t_col_base_fun;
1927 t_m(L, k) += v * t_approx_P_adjont_log_du_domega(k, L);
1928 ++t_m;
1929 ++t_col_base_fun;
1930 }
1931 ++t_row_base_fun;
1932 }
1933
1934 for (; rr != row_nb_base_functions; ++rr)
1935 ++t_row_base_fun;
1936
1937 ++t_w;
1938 ++t_approx_P_adjont_log_du_domega;
1939 }
1940
1942}
1943
1945 EntData &col_data) {
1947 int nb_integration_pts = getGaussPts().size2();
1948 int row_nb_dofs = row_data.getIndices().size();
1949 int col_nb_dofs = col_data.getIndices().size();
1950 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1952
1953 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1954
1955 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1956
1957 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1958
1959 );
1960 };
1961 FTensor::Index<'i', 3> i;
1962 FTensor::Index<'j', 3> j;
1963 FTensor::Index<'k', 3> k;
1964 FTensor::Index<'l', 3> l;
1965 FTensor::Index<'m', 3> m;
1966 FTensor::Index<'n', 3> n;
1967 auto v = getVolume();
1968 auto t_w = getFTensor0IntegrationWeight();
1969 auto t_levi_kirchoff_dP =
1970 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1971
1972 int row_nb_base_functions = row_data.getN().size2();
1973 auto t_row_base_fun = row_data.getFTensor0N();
1974
1975 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1976 double a = v * t_w;
1977 int rr = 0;
1978 for (; rr != row_nb_dofs / 3; ++rr) {
1979 double b = a * t_row_base_fun;
1980 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1981 auto t_m = get_ftensor2(K, 3 * rr, 0);
1982 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1983 t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(l));
1984 ++t_m;
1985 ++t_col_base_fun;
1986 }
1987 ++t_row_base_fun;
1988 }
1989 for (; rr != row_nb_base_functions; ++rr) {
1990 ++t_row_base_fun;
1991 }
1992
1993 ++t_w;
1994 ++t_levi_kirchoff_dP;
1995 }
1997}
1998
2000 EntData &col_data) {
2002 int nb_integration_pts = getGaussPts().size2();
2003 int row_nb_dofs = row_data.getIndices().size();
2004 int col_nb_dofs = col_data.getIndices().size();
2005
2006 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2008 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
2009 };
2010
2011 FTensor::Index<'i', 3> i;
2012 FTensor::Index<'j', 3> j;
2013 FTensor::Index<'k', 3> k;
2014 FTensor::Index<'l', 3> l;
2015 FTensor::Index<'m', 3> m;
2016
2017 auto v = getVolume();
2018 auto t_w = getFTensor0IntegrationWeight();
2019 auto t_levi_kirchoff_dP =
2020 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2021
2022 int row_nb_base_functions = row_data.getN().size2();
2023 auto t_row_base_fun = row_data.getFTensor0N();
2024
2025 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2026 double a = v * t_w;
2027 int rr = 0;
2028 for (; rr != row_nb_dofs / 3; ++rr) {
2029 double b = a * t_row_base_fun;
2030 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
2031 auto t_m = get_ftensor1(K, 3 * rr, 0);
2032 for (int cc = 0; cc != col_nb_dofs; ++cc) {
2033 t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(i, j));
2034 ++t_m;
2035 ++t_col_base_fun;
2036 }
2037 ++t_row_base_fun;
2038 }
2039
2040 for (; rr != row_nb_base_functions; ++rr) {
2041 ++t_row_base_fun;
2042 }
2043 ++t_w;
2044 ++t_levi_kirchoff_dP;
2045 }
2047}
2048
2050 EntData &col_data) {
2052 int nb_integration_pts = getGaussPts().size2();
2053 int row_nb_dofs = row_data.getIndices().size();
2054 int col_nb_dofs = col_data.getIndices().size();
2055 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2057
2058 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2059
2060 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2061
2062 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2063
2064 );
2065 };
2066 FTensor::Index<'i', 3> i;
2067 FTensor::Index<'j', 3> j;
2068 FTensor::Index<'k', 3> k;
2069 FTensor::Index<'l', 3> l;
2070 FTensor::Index<'m', 3> m;
2071 FTensor::Index<'n', 3> n;
2072
2073 auto v = getVolume();
2074 auto t_w = getFTensor0IntegrationWeight();
2075 auto t_levi_kirchoff_domega =
2076 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
2077 int row_nb_base_functions = row_data.getN().size2();
2078 auto t_row_base_fun = row_data.getFTensor0N();
2079 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2080 double a = v * t_w;
2081 int rr = 0;
2082 for (; rr != row_nb_dofs / 3; ++rr) {
2083 auto t_m = get_ftensor2(K, 3 * rr, 0);
2084 const double b = a * t_row_base_fun;
2085 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2086 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2087 t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(k, l);
2088 ++t_m;
2089 ++t_col_base_fun;
2090 }
2091 ++t_row_base_fun;
2092 }
2093 for (; rr != row_nb_base_functions; ++rr) {
2094 ++t_row_base_fun;
2095 }
2096 ++t_w;
2097 ++t_levi_kirchoff_domega;
2098 }
2100}
2101
2103 EntData &col_data) {
2105
2106 int nb_integration_pts = row_data.getN().size1();
2107 int row_nb_dofs = row_data.getIndices().size();
2108 int col_nb_dofs = col_data.getIndices().size();
2109
2110 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2112
2113 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2114
2115 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2116
2117 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2118
2119 );
2120 };
2121
2122 FTensor::Index<'i', 3> i;
2123 FTensor::Index<'j', 3> j;
2124 FTensor::Index<'k', 3> k;
2125 FTensor::Index<'l', 3> l;
2126 FTensor::Index<'m', 3> m;
2127 FTensor::Index<'n', 3> n;
2128
2129 auto v = getVolume();
2130 auto t_w = getFTensor0IntegrationWeight();
2131 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2132 int row_nb_base_functions = row_data.getN().size2() / 3;
2133 auto t_row_base_fun = row_data.getFTensor1N<3>();
2134
2135 const double ts_a = getTSa();
2136
2137 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2138 double a = v * t_w;
2139
2140 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
2141
2142 int rr = 0;
2143 for (; rr != row_nb_dofs / 3; ++rr) {
2144
2146 t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
2147
2148 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2149 auto t_m = get_ftensor2(K, 3 * rr, 0);
2150 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2151 t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
2152 ++t_m;
2153 ++t_col_base_fun;
2154 }
2155
2156 ++t_row_base_fun;
2157 }
2158
2159 for (; rr != row_nb_base_functions; ++rr)
2160 ++t_row_base_fun;
2161 ++t_w;
2162 ++t_h_domega;
2163 }
2165}
2166
2169 EntData &col_data) {
2171
2172 int nb_integration_pts = row_data.getN().size1();
2173 int row_nb_dofs = row_data.getIndices().size();
2174 int col_nb_dofs = col_data.getIndices().size();
2175
2176 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2178 &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
2179 };
2180
2181 FTensor::Index<'i', 3> i;
2182 FTensor::Index<'j', 3> j;
2183 FTensor::Index<'k', 3> k;
2184 FTensor::Index<'l', 3> l;
2185 FTensor::Index<'m', 3> m;
2186 FTensor::Index<'n', 3> n;
2187
2188 auto v = getVolume();
2189 auto t_w = getFTensor0IntegrationWeight();
2190 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2191 int row_nb_base_functions = row_data.getN().size2() / 9;
2192 auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
2193
2194 const double ts_a = getTSa();
2195
2196 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2197 double a = v * t_w;
2198
2199 int rr = 0;
2200 for (; rr != row_nb_dofs; ++rr) {
2201
2203 t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
2204
2205 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2206 auto t_m = get_ftensor2(K, rr, 0);
2207 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2208 t_m(j) += (a * t_col_base_fun) * t_PRT(j);
2209 ++t_m;
2210 ++t_col_base_fun;
2211 }
2212
2213 ++t_row_base_fun;
2214 }
2215
2216 for (; rr != row_nb_base_functions; ++rr)
2217 ++t_row_base_fun;
2218
2219 ++t_w;
2220 ++t_h_domega;
2221 }
2223}
2224
2226 EntData &data) {
2228
2229 auto create_tag = [this](const std::string tag_name, const int size) {
2230 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2231 Tag th;
2232 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2233 th, MB_TAG_CREAT | MB_TAG_SPARSE,
2234 def_VAL);
2235 return th;
2236 };
2237
2238 Tag th_w = create_tag("SpatialDisplacement", 3);
2239 Tag th_omega = create_tag("Omega", 3);
2240 Tag th_approxP = create_tag("Piola1Stress", 9);
2241 Tag th_sigma = create_tag("CauchyStress", 9);
2242 Tag th_log_u = create_tag("LogSpatialStretch", 9);
2243 Tag th_u = create_tag("SpatialStretch", 9);
2244 Tag th_h = create_tag("h", 9);
2245 Tag th_X = create_tag("X", 3);
2246 Tag th_detF = create_tag("detF", 1);
2247 Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2248
2249 Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
2250 Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
2251 Tag th_traction = create_tag("traction", 3);
2252
2253 Tag th_disp = create_tag("H1Displacement", 3);
2254 Tag th_disp_error = create_tag("DisplacementError", 1);
2255 Tag th_lambda_disp = create_tag("ContactDisplacement", 3);
2256
2257 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
2258 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
2259 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2260 auto t_log_u =
2261 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
2262 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
2263 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
2264 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2265 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
2266 auto t_coords = getFTensor1CoordsAtGaussPts();
2267 auto t_normal = getFTensor1NormalsAtGaussPts();
2268 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
2269 auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
2270
2271 FTensor::Index<'i', 3> i;
2272 FTensor::Index<'j', 3> j;
2273 FTensor::Index<'k', 3> k;
2274 FTensor::Index<'l', 3> l;
2275
2276 auto set_float_precision = [](const double x) {
2277 if (std::abs(x) < std::numeric_limits<float>::epsilon())
2278 return 0.;
2279 else
2280 return x;
2281 };
2282
2283 // scalars
2284 auto save_scal_tag = [&](auto &th, auto v, const int gg) {
2286 v = set_float_precision(v);
2287 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
2289 };
2290
2291 // vectors
2292 VectorDouble3 v(3);
2293 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2294 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2296 t_v(i) = t_d(i);
2297 for (auto &a : v.data())
2298 a = set_float_precision(a);
2299 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2300 &*v.data().begin());
2302 };
2303
2304 // tensors
2305
2306 MatrixDouble3by3 m(3, 3);
2308 &m(0, 0), &m(0, 1), &m(0, 2),
2309
2310 &m(1, 0), &m(1, 1), &m(1, 2),
2311
2312 &m(2, 0), &m(2, 1), &m(2, 2));
2313
2314 auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2316 t_m(i, j) = t_d(i, j);
2317 for (auto &v : m.data())
2318 v = set_float_precision(v);
2319 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2320 &*m.data().begin());
2322 };
2323
2324 const auto nb_gauss_pts = getGaussPts().size2();
2325 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
2326
2327 FTensor::Tensor1<double, 3> t_traction;
2328 t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
2329
2330 // vectors
2331 CHKERR save_vec_tag(th_w, t_w, gg);
2332 CHKERR save_vec_tag(th_X, t_coords, gg);
2333 CHKERR save_vec_tag(th_omega, t_omega, gg);
2334 CHKERR save_vec_tag(th_traction, t_traction, gg);
2335
2336 // tensors
2337 CHKERR save_mat_tag(th_h, t_h, gg);
2338
2340 for (int ii = 0; ii != 3; ++ii)
2341 for (int jj = 0; jj != 3; ++jj)
2342 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2343
2344 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2345
2347 for (int ii = 0; ii != 3; ++ii)
2348 for (int jj = 0; jj != 3; ++jj)
2349 t_u_tmp(ii, jj) = t_u(ii, jj);
2350
2351 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2352 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2353 CHKERR save_vec_tag(th_disp, t_disp, gg);
2354 CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
2355
2356 double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
2357 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2358
2359 const double jac = determinantTensor3by3(t_h);
2361 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2362 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2363 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2364
2366 t_levi(k) = t_levi_kirchoff(k);
2367 CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2368 &t_levi(0));
2369
2370 auto get_eiegn_vector_symmetric = [&](auto &t_u) {
2372
2373 for (int ii = 0; ii != 3; ++ii)
2374 for (int jj = 0; jj != 3; ++jj)
2375 t_m(ii, jj) = t_u(ii, jj);
2376
2377 VectorDouble3 eigen_values(3);
2378 auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2379 CHKERR computeEigenValuesSymmetric(t_m, t_eigen_values);
2380
2381 CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2382 &*m.data().begin());
2383 CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2384 &*eigen_values.data().begin());
2385
2387 };
2388
2389 CHKERR get_eiegn_vector_symmetric(t_u);
2390
2391 ++t_w;
2392 ++t_h;
2393 ++t_log_u;
2394 ++t_u;
2395 ++t_omega;
2396 ++t_R;
2397 ++t_approx_P;
2398 ++t_levi_kirchoff;
2399 ++t_coords;
2400 ++t_normal;
2401 ++t_disp;
2402 ++t_lambda_disp;
2403 }
2404
2406}
2407
2412 if (type == MBTET) {
2413 int nb_integration_pts = data.getN().size1();
2414 auto v = getVolume();
2415 auto t_w = getFTensor0IntegrationWeight();
2416 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2417 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2418
2419 FTensor::Index<'i', 3> i;
2420 FTensor::Index<'J', 3> J;
2421
2422 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2423 const double a = t_w * v;
2424 // FIXME: this is wrong, energy should be calculated in material
2425 (*energy) += a * t_P(i, J) * t_h(i, J);
2426 ++t_w;
2427 ++t_P;
2428 ++t_h;
2429 }
2430 }
2432}
2433
2434} // namespace EshelbianPlasticity
static Index< 'J', 3 > J
constexpr double third
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()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
static double lambda
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)
MoFEM::TsCtx * ts_ctx
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 diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
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.
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
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
double scale
Definition plastic.cpp:170
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
MoFEMErrorCode integrate(EntData &data)
std::vector< boost::shared_ptr< ScalingMethod > > scalingMethodsVec
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(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 &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 integratePolyconvexHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(EntData &data)
MoFEMErrorCode integratePolyconvexHencky(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.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
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