v0.13.2
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
14#include <boost/math/constants/constants.hpp>
15
16#include <lapack_wrap.h>
17
18#include <MatrixFunction.hpp>
19
20namespace EshelbianPlasticity {
21
22double f(const double v) { return exp(v); }
23double d_f(const double v) { return exp(v); }
24double dd_f(const double v) { return exp(v); }
25
26template <class T> struct TensorTypeExtractor {
27 typedef typename std::remove_pointer<T>::type Type;
28};
29template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
30 typedef typename std::remove_pointer<T>::type Type;
31};
32
33template <typename T>
36 &t_omega) {
41 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
42 t_R(i, j) = t_kd(i, j);
43
44 const typename TensorTypeExtractor<T>::Type angle =
45 sqrt(t_omega(i) * t_omega(i));
46
47 constexpr double tol = 1e-18;
48 if (std::abs(angle) < tol) {
49 return t_R;
50 }
51
53 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
54 const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
55 const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
56 const typename TensorTypeExtractor<T>::Type b =
57 2. * ss_2 * ss_2 / (angle * angle);
58 t_R(i, j) += a * t_Omega(i, j);
59 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
60
61 return t_R;
62}
63
64template <typename T>
67 &t_omega) {
68
73
74 const typename TensorTypeExtractor<T>::Type angle =
75 sqrt(t_omega(i) * t_omega(i));
76
77 constexpr double tol = 1e-18;
78 if (std::abs(angle) < tol) {
80 auto t_R = get_rotation_form_vector(t_omega);
81 t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
82 return t_diff_R;
83 }
84
87 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
88
89 const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
90 const typename TensorTypeExtractor<T>::Type ss = sin(angle);
91 const typename TensorTypeExtractor<T>::Type a = ss / angle;
92 const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
93 const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
94 t_diff_R(i, j, k) = 0;
95 t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
96 t_diff_R(i, j, k) +=
97 b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
98 FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
99
100 const typename TensorTypeExtractor<T>::Type cc = cos(angle);
101 const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
102 const typename TensorTypeExtractor<T>::Type diff_a =
103 (angle * cc - ss) / angle2;
104
105 const typename TensorTypeExtractor<T>::Type diff_b =
106 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
107 t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
108 t_diff_R(i, j, k) +=
109 diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
110
111 return t_diff_R;
112}
113
114template <typename T1, typename T2, typename T3>
116 T3 &eig_vec) {
118 for (int ii = 0; ii != 3; ii++)
119 for (int jj = 0; jj != 3; jj++)
120 eig_vec(ii, jj) = X(ii, jj);
121
122 int n = 3;
123 int lda = 3;
124 int lwork = (3 + 2) * 3;
125 std::array<double, (3 + 2) * 3> work;
126
127 if (lapack_dsyev('V', 'U', n, &(eig_vec(0, 0)), lda, &eig(0), work.data(),
128 lwork) > 0)
129 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
130 "The algorithm failed to compute eigenvalues.");
132}
133
134inline bool is_eq(const double &a, const double &b) {
135 constexpr double eps = 1e-8;
136 return abs(a - b) < eps;
137}
138
139template <typename T> inline int get_uniq_nb(T &ec) {
140 return distance(&ec(0), unique(&ec(0), &ec(0) + 3, is_eq));
141}
142
143template <typename T1, typename T2>
144void sort_eigen_vals(T1 &eig, T2 &eigen_vec) {
147 if (is_eq(eig(0), eig(1))) {
149 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
150 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
151 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
152 FTensor::Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
153 eig(i) = eig_c(i);
154 eigen_vec(i, j) = eigen_vec_c(i, j);
155 } else {
157 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
158 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
159 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
160 FTensor::Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
161 eig(i) = eig_c(i);
162 eigen_vec(i, j) = eigen_vec_c(i, j);
163 }
164}
165
167 EntityType type,
168 EntData &data) {
170 if (type != MBTET)
172 int nb_integration_pts = data.getN().size1();
177
178 dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
179 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
180 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
181 dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts, false);
182 dataAtPts->diffStreachTensorAtPts.resize(36, nb_integration_pts, false);
183 dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
184 dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
185 dataAtPts->nbUniq.resize(nb_integration_pts, false);
186
187 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
188 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
189 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
190 auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
191 auto t_log_u =
192 getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
193 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
194
195 auto t_diff_u =
196 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
197 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
198 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
199 auto &nbUniq = dataAtPts->nbUniq;
200
201 for (int gg = 0; gg != nb_integration_pts; ++gg) {
202
203 auto t0_diff = get_diff_rotation_form_vector(t_omega);
204 auto t0 = get_rotation_form_vector(t_omega);
205
206 t_diff_R(i, j, k) = t0_diff(i, j, k);
207 t_R(i, j) = t0(i, j);
208
211 CHKERR get_eigen_val_and_proj_lapack(t_log_u, eig, eigen_vec);
212
213 t_eigen_vals(i) = eig(i);
214 t_eigen_vecs(i, j) = eigen_vec(i, j);
215
216 // rare case when two eigen values are equal
217 nbUniq[gg] = get_uniq_nb(eig);
218 if (nbUniq[gg] == 2)
219 sort_eigen_vals(eig, eigen_vec);
220
221 auto t_u_tmp = EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, f);
222 auto t_diff_u_tmp =
223 EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, f, d_f, nbUniq[gg]);
224 for (int ii = 0; ii != 3; ++ii) {
225 for (int jj = ii; jj != 3; ++jj) {
226 for (int kk = 0; kk != 3; ++kk) {
227 for (int ll = 0; ll < kk; ++ll) {
228 t_diff_u_tmp(ii, jj, kk, ll) *= 2;
229 }
230 }
231 }
232 }
233
234 t_u(i, j) = t_u_tmp(i, j);
235 t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
236 t_h(i, j) = t_R(i, k) * t_u(k, j);
237
238 ++t_h;
239 ++t_R;
240 ++t_diff_R;
241 ++t_log_u;
242 ++t_u;
243 ++t_diff_u;
244 ++t_eigen_vals;
245 ++t_eigen_vecs;
246 ++t_omega;
247 }
248
250}
251
254 int nb_dofs = data.getIndices().size();
255 int nb_integration_pts = data.getN().size1();
256 auto v = getVolume();
257 auto t_w = getFTensor0IntegrationWeight();
258 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
259 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotAtPts);
260 if (dataAtPts->wDotDotAtPts.size1() == 0 &&
261 dataAtPts->wDotDotAtPts.size2() != nb_integration_pts) {
262 dataAtPts->wDotDotAtPts.resize(3, nb_integration_pts);
263 dataAtPts->wDotDotAtPts.clear();
264 }
265 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotDotAtPts);
266
267 int nb_base_functions = data.getN().size2();
268 auto t_row_base_fun = data.getFTensor0N();
269
271 auto get_ftensor1 = [](auto &v) {
273 &v[2]);
274 };
275
276 for (int gg = 0; gg != nb_integration_pts; ++gg) {
277 double a = v * t_w;
278 auto t_nf = get_ftensor1(nF);
279 int bb = 0;
280 for (; bb != nb_dofs / 3; ++bb) {
281 t_nf(i) += a * t_row_base_fun * t_div_P(i);
282 t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
283 t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
284 ++t_nf;
285 ++t_row_base_fun;
286 }
287 for (; bb != nb_base_functions; ++bb)
288 ++t_row_base_fun;
289 ++t_w;
290 ++t_div_P;
291 ++t_s_dot_w;
292 ++t_s_dot_dot_w;
293 }
294
296}
297
300 int nb_dofs = data.getIndices().size();
301 int nb_integration_pts = data.getN().size1();
302 auto v = getVolume();
303 auto t_w = getFTensor0IntegrationWeight();
304 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
305 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
306
307 int nb_base_functions = data.getN().size2();
308 auto t_row_base_fun = data.getFTensor0N();
314 auto get_ftensor1 = [](auto &v) {
316 &v[2]);
317 };
318
319 for (int gg = 0; gg != nb_integration_pts; ++gg) {
320 double a = v * t_w;
321 auto t_nf = get_ftensor1(nF);
323 t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
325 t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
326 int bb = 0;
327 for (; bb != nb_dofs / 3; ++bb) {
328 t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
329 ++t_nf;
330 ++t_row_base_fun;
331 }
332 for (; bb != nb_base_functions; ++bb)
333 ++t_row_base_fun;
334 ++t_w;
335 ++t_approx_P;
336 ++t_R;
337 }
339}
340
343 int nb_dofs = data.getIndices().size();
344 int nb_integration_pts = data.getN().size1();
345 auto v = getVolume();
346 auto t_w = getFTensor0IntegrationWeight();
347 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
348 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
349 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
350 auto t_dot_log_u =
351 getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
352 auto t_diff_u =
353 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
354
359 auto get_ftensor2 = [](auto &v) {
361 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
362 };
363
364 int nb_base_functions = data.getN().size2();
365 auto t_row_base_fun = data.getFTensor0N();
366 for (int gg = 0; gg != nb_integration_pts; ++gg) {
367 double a = v * t_w;
368 auto t_nf = get_ftensor2(nF);
369
371 t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
373 t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
374
376 t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
377
379 t1(i, j) =
380 a * (t_diff_u(k, l, i, j) * (t_deltaP(k, l) - t_viscous_stress(k, l)));
381
382 int bb = 0;
383 for (; bb != nb_dofs / 6; ++bb) {
384
385 t_nf(i, j) += t_row_base_fun * t1(i, j);
386
387 ++t_nf;
388 ++t_row_base_fun;
389 }
390 for (; bb != nb_base_functions; ++bb)
391 ++t_row_base_fun;
392
393 ++t_w;
394 ++t_approx_P;
395 ++t_P;
396 ++t_R;
397 ++t_dot_log_u;
398 ++t_diff_u;
399 }
401}
402
405 int nb_dofs = data.getIndices().size();
406 int nb_integration_pts = data.getN().size1();
407 auto v = getVolume();
408 auto t_w = getFTensor0IntegrationWeight();
409 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
410 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
411 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
412
413 int nb_base_functions = data.getN().size2() / 3;
414 auto t_row_base_fun = data.getFTensor1N<3>();
419
420 auto get_ftensor1 = [](auto &v) {
422 &v[2]);
423 };
424
425 for (int gg = 0; gg != nb_integration_pts; ++gg) {
426 double a = v * t_w;
427 auto t_nf = get_ftensor1(nF);
428
429 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
431 t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
432 FTensor::Tensor2<double, 3, 3> t_residuum_omega;
433 t_residuum_omega(i, j) =
434 (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, j);
435
436 int bb = 0;
437 for (; bb != nb_dofs / 3; ++bb) {
438 t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
439 t_nf(i) += a * t_row_base_fun(j) * t_residuum_omega(i, j);
440
441 ++t_nf;
442 ++t_row_base_fun;
443 }
444
445 for (; bb != nb_base_functions; ++bb)
446 ++t_row_base_fun;
447
448 ++t_w;
449 ++t_u;
450 ++t_R;
451 ++t_omega_dot;
452 }
453
455}
456
459 int nb_dofs = data.getIndices().size();
460 int nb_integration_pts = data.getN().size1();
461 auto v = getVolume();
462 auto t_w = getFTensor0IntegrationWeight();
463 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
464 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
465 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
466
467 int nb_base_functions = data.getN().size2() / 9;
468 auto t_row_base_fun = data.getFTensor2N<3, 3>();
473
474 auto get_ftensor0 = [](auto &v) {
476 };
477
478 for (int gg = 0; gg != nb_integration_pts; ++gg) {
479 double a = v * t_w;
480 auto t_nf = get_ftensor0(nF);
481
482 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
484 t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
485 FTensor::Tensor2<double, 3, 3> t_residuum_omega;
486 t_residuum_omega(i, j) =
487 (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, j);
488
489 int bb = 0;
490 for (; bb != nb_dofs; ++bb) {
491 t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
492 t_nf += a * t_row_base_fun(i, j) * t_residuum_omega(i, j);
493 ++t_nf;
494 ++t_row_base_fun;
495 }
496 for (; bb != nb_base_functions; ++bb) {
497 ++t_row_base_fun;
498 }
499 ++t_w;
500 ++t_R;
501 ++t_u;
502 ++t_omega_dot;
503 }
504
506}
507
510 int nb_dofs = data.getIndices().size();
511 int nb_integration_pts = data.getN().size1();
512 auto v = getVolume();
513 auto t_w = getFTensor0IntegrationWeight();
514 auto t_s_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
515 int nb_base_functions = data.getN().size2() / 3;
516 auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
518 auto get_ftensor1 = [](auto &v) {
520 &v[2]);
521 };
522
523 for (int gg = 0; gg != nb_integration_pts; ++gg) {
524 double a = v * t_w;
525 auto t_nf = get_ftensor1(nF);
526 int bb = 0;
527 for (; bb != nb_dofs / 3; ++bb) {
528 double div_row_base = t_row_diff_base_fun(i, i);
529 t_nf(i) += a * div_row_base * t_s_w(i);
530 ++t_nf;
531 ++t_row_diff_base_fun;
532 }
533 for (; bb != nb_base_functions; ++bb) {
534 ++t_row_diff_base_fun;
535 }
536 ++t_w;
537 ++t_s_w;
538 }
539
541}
542
545 // get entity of face
546 EntityHandle fe_ent = getFEEntityHandle();
547 // interate over all boundary data
548 for (auto &bc : (*bcDispPtr)) {
549 // check if finite element entity is part of boundary condition
550 if (bc.faces.find(fe_ent) != bc.faces.end()) {
551 int nb_dofs = data.getIndices().size();
552 int nb_integration_pts = data.getN().size1();
553 auto t_normal = getFTensor1Normal();
554 auto t_w = getFTensor0IntegrationWeight();
555 int nb_base_functions = data.getN().size2() / 3;
556 auto t_row_base_fun = data.getFTensor1N<3>();
557
560 auto get_ftensor1 = [](auto &v) {
562 &v[2]);
563 };
564
565 // get bc data
566 FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
567 t_bc_disp(i) *= getFEMethod()->ts_t;
568
569 for (int gg = 0; gg != nb_integration_pts; ++gg) {
571 t_bc_res(i) = t_bc_disp(i);
572 auto t_nf = get_ftensor1(nF);
573 int bb = 0;
574 for (; bb != nb_dofs / 3; ++bb) {
575 t_nf(i) -=
576 t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
577 ++t_nf;
578 ++t_row_base_fun;
579 }
580 for (; bb != nb_base_functions; ++bb)
581 ++t_row_base_fun;
582
583 ++t_w;
584 }
585 }
586 }
588}
589
592 // get entity of face
593 EntityHandle fe_ent = getFEEntityHandle();
594 // interate over all boundary data
595 for (auto &bc : (*bcRotPtr)) {
596 // check if finite element entity is part of boundary condition
597 if (bc.faces.find(fe_ent) != bc.faces.end()) {
598 int nb_dofs = data.getIndices().size();
599 int nb_integration_pts = data.getN().size1();
600 auto t_normal = getFTensor1Normal();
601 auto t_w = getFTensor0IntegrationWeight();
602
603 int nb_base_functions = data.getN().size2() / 3;
604 auto t_row_base_fun = data.getFTensor1N<3>();
608 auto get_ftensor1 = [](auto &v) {
610 &v[2]);
611 };
612
613 // get bc data
614 FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
615
616 double theta = bc.theta;
617 theta *= getFEMethod()->ts_t;
618
620 const double a = sqrt(t_normal(i) * t_normal(i));
621 t_omega(i) = t_normal(i) * (theta / a);
622 auto t_R = get_rotation_form_vector(t_omega);
623 auto t_coords = getFTensor1CoordsAtGaussPts();
624
625 for (int gg = 0; gg != nb_integration_pts; ++gg) {
627 t_delta(i) = t_center(i) - t_coords(i);
629 t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
630
631 auto t_nf = get_ftensor1(nF);
632 int bb = 0;
633 for (; bb != nb_dofs / 3; ++bb) {
634 t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
635 ++t_nf;
636 ++t_row_base_fun;
637 }
638 for (; bb != nb_base_functions; ++bb)
639 ++t_row_base_fun;
640
641 ++t_w;
642 ++t_coords;
643 }
644 }
645 }
647}
648
651
652 struct FaceRule {
653 int operator()(int p_row, int p_col, int p_data) const {
654 return 2 * (p_data + 1);
655 }
656 };
657
658 if (ts_ctx == CTX_TSSETIFUNCTION) {
659
660 // Loop boundary elements with traction boundary conditions
663 {HDIV});
664 fe.getOpPtrVector().push_back(
665 new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
666 fe.getRuleHook = FaceRule();
667 fe.ts_t = ts_t;
668 CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
669 fe, nullptr, MF_ZERO,
670 this->getCacheWeakPtr());
671 }
672
674}
675
678
681
682 auto t_normal = getFTensor1Normal();
683 const double nrm2 = sqrt(t_normal(i) * t_normal(i));
684 FTensor::Tensor1<double, 3> t_unit_normal;
685 t_unit_normal(i) = t_normal(i) / nrm2;
686 int nb_dofs = data.getFieldData().size();
687 int nb_integration_pts = data.getN().size1();
688 int nb_base_functions = data.getN().size2() / 3;
689 double ts_t = getFEMethod()->ts_t;
690
691 auto integrate_matrix = [&]() {
693
694 auto t_w = getFTensor0IntegrationWeight();
695 matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
696 matPP.clear();
697
698 auto t_row_base_fun = data.getFTensor1N<3>();
699 for (int gg = 0; gg != nb_integration_pts; ++gg) {
700
701 int rr = 0;
702 for (; rr != nb_dofs / 3; ++rr) {
703 const double a = t_row_base_fun(i) * t_unit_normal(i);
704 auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
705 for (int cc = 0; cc != nb_dofs / 3; ++cc) {
706 const double b = t_col_base_fun(i) * t_unit_normal(i);
707 matPP(rr, cc) += t_w * a * b;
708 ++t_col_base_fun;
709 }
710 ++t_row_base_fun;
711 }
712
713 for (; rr != nb_base_functions; ++rr)
714 ++t_row_base_fun;
715
716 ++t_w;
717 }
718
720 };
721
722 auto integrate_rhs = [&](auto &bc) {
724
725 auto t_w = getFTensor0IntegrationWeight();
726 vecPv.resize(3, nb_dofs / 3, false);
727 vecPv.clear();
728
729 auto t_row_base_fun = data.getFTensor1N<3>();
730 double ts_t = getFEMethod()->ts_t;
731
732 for (int gg = 0; gg != nb_integration_pts; ++gg) {
733 int rr = 0;
734 for (; rr != nb_dofs / 3; ++rr) {
735 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
736 for (int dd = 0; dd != 3; ++dd)
737 if (bc.flags[dd])
738 vecPv(dd, rr) += t * bc.vals[dd];
739 ++t_row_base_fun;
740 }
741 for (; rr != nb_base_functions; ++rr)
742 ++t_row_base_fun;
743 ++t_w;
744 }
746 };
747
748 auto integrate_rhs_cook = [&](auto &bc) {
750
751 vecPv.resize(3, nb_dofs / 3, false);
752 vecPv.clear();
753
754 auto t_w = getFTensor0IntegrationWeight();
755 auto t_row_base_fun = data.getFTensor1N<3>();
756 auto t_coords = getFTensor1CoordsAtGaussPts();
757
758 for (int gg = 0; gg != nb_integration_pts; ++gg) {
759
760 auto calc_tau = [](double y) {
761 y -= 44;
762 y /= (60 - 44);
763 return -y * (y - 1) / 0.25;
764 };
765
766 const double tau = calc_tau(t_coords(1));
767
768 int rr = 0;
769 for (; rr != nb_dofs / 3; ++rr) {
770 const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
771 for (int dd = 0; dd != 3; ++dd)
772 if (bc.flags[dd])
773 vecPv(dd, rr) += t * tau * bc.vals[dd];
774 ++t_row_base_fun;
775 }
776
777 for (; rr != nb_base_functions; ++rr)
778 ++t_row_base_fun;
779 ++t_w;
780 ++t_coords;
781 }
783 };
784
785 // get entity of face
786 EntityHandle fe_ent = getFEEntityHandle();
787 for (auto &bc : *(bcFe.bcData)) {
788 if (bc.faces.find(fe_ent) != bc.faces.end()) {
789
790 int nb_dofs = data.getFieldData().size();
791 if (nb_dofs) {
792
793 CHKERR integrate_matrix();
794 if (bc.blockName.compare(20, 4, "COOK") == 0)
795 CHKERR integrate_rhs_cook(bc);
796 else
797 CHKERR integrate_rhs(bc);
798
799 const auto info =
800 lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
801 nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
802 if (info > 0)
803 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
804 "The leading minor of order %i is "
805 "not positive; definite;\nthe "
806 "solution could not be computed",
807 info);
808
809 for (int dd = 0; dd != 3; ++dd)
810 if (bc.flags[dd])
811 for (int rr = 0; rr != nb_dofs / 3; ++rr)
812 data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
813 }
814 }
815 }
816
818}
819
821 EntData &col_data) {
823 int nb_integration_pts = row_data.getN().size1();
824 int row_nb_dofs = row_data.getIndices().size();
825 int col_nb_dofs = col_data.getIndices().size();
826 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
828 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
829 };
831 auto v = getVolume();
832 auto t_w = getFTensor0IntegrationWeight();
833 int row_nb_base_functions = row_data.getN().size2();
834 auto t_row_base_fun = row_data.getFTensor0N();
835 for (int gg = 0; gg != nb_integration_pts; ++gg) {
836 double a = v * t_w;
837 int rr = 0;
838 for (; rr != row_nb_dofs / 3; ++rr) {
839 auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
840 auto t_m = get_ftensor1(K, 3 * rr, 0);
841 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
842 double div_col_base = t_col_diff_base_fun(i, i);
843 t_m(i) += a * t_row_base_fun * div_col_base;
844 ++t_m;
845 ++t_col_diff_base_fun;
846 }
847 ++t_row_base_fun;
848 }
849 for (; rr != row_nb_base_functions; ++rr)
850 ++t_row_base_fun;
851 ++t_w;
852 }
854}
855
857 EntData &col_data) {
859
860 if (alphaW < std::numeric_limits<double>::epsilon() &&
861 alphaRho < std::numeric_limits<double>::epsilon())
863
864 const int nb_integration_pts = row_data.getN().size1();
865 const int row_nb_dofs = row_data.getIndices().size();
866 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
868 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
869
870 );
871 };
873
874 auto v = getVolume();
875 auto t_w = getFTensor0IntegrationWeight();
876
877 int row_nb_base_functions = row_data.getN().size2();
878 auto t_row_base_fun = row_data.getFTensor0N();
879
880 double ts_scale = alphaW * getTSa();
881 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
882 ts_scale += alphaRho * getTSaa();
883
884 for (int gg = 0; gg != nb_integration_pts; ++gg) {
885 double a = v * t_w * ts_scale;
886
887 int rr = 0;
888 for (; rr != row_nb_dofs / 3; ++rr) {
889
890 auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
891 auto t_m = get_ftensor1(K, 3 * rr, 0);
892 for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
893 const double b = a * t_row_base_fun * t_col_base_fun;
894 t_m(i) -= b;
895 ++t_m;
896 ++t_col_base_fun;
897 }
898
899 ++t_row_base_fun;
900 }
901
902 for (; rr != row_nb_base_functions; ++rr)
903 ++t_row_base_fun;
904
905 ++t_w;
906 }
907
909}
910
912 EntData &col_data) {
914 int nb_integration_pts = row_data.getN().size1();
915 int row_nb_dofs = row_data.getIndices().size();
916 int col_nb_dofs = col_data.getIndices().size();
917 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
919
920 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
921
922 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
923
924 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
925
926 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
927
928 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
929
930 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
931 };
932
939
940 auto v = getVolume();
941 auto t_w = getFTensor0IntegrationWeight();
942 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
943
944 auto t_diff_u =
945 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
946
947 int row_nb_base_functions = row_data.getN().size2();
948 auto t_row_base_fun = row_data.getFTensor0N();
949 for (int gg = 0; gg != nb_integration_pts; ++gg) {
950 double a = v * t_w;
951
952 int rr = 0;
953 for (; rr != row_nb_dofs / 6; ++rr) {
954
955 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
956 auto t_m = get_ftensor3(K, 6 * rr, 0);
957
958 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
959
961 t_RTP_dP(i, j, k) = t_R(k, i) * t_col_base_fun(j) * t_row_base_fun;
962
963 for (int kk = 0; kk != 3; ++kk)
964 for (int ii = 0; ii != 3; ++ii)
965 for (int jj = ii; jj != 3; ++jj) {
966 for (int mm = 0; mm != 3; ++mm)
967 for (int nn = 0; nn != 3; ++nn)
968 t_m(ii, jj, kk) +=
969 a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
970 }
971
972 ++t_col_base_fun;
973 ++t_m;
974 }
975
976 ++t_row_base_fun;
977 }
978 for (; rr != row_nb_base_functions; ++rr)
979 ++t_row_base_fun;
980 ++t_w;
981 ++t_R;
982 ++t_diff_u;
983 }
985}
986
988 EntData &col_data) {
990 int nb_integration_pts = row_data.getN().size1();
991 int row_nb_dofs = row_data.getIndices().size();
992 int col_nb_dofs = col_data.getIndices().size();
993 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
995 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
996 &m(r + 5, c));
997 };
998
1004
1005 auto v = getVolume();
1006 auto t_w = getFTensor0IntegrationWeight();
1007 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1008 auto t_diff_u =
1009 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1010
1011 int row_nb_base_functions = row_data.getN().size2();
1012 auto t_row_base_fun = row_data.getFTensor0N();
1013 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1014 double a = v * t_w;
1015
1016 int rr = 0;
1017 for (; rr != row_nb_dofs / 6; ++rr) {
1018
1019 auto t_m = get_ftensor2(K, 6 * rr, 0);
1020
1021 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1022 for (int cc = 0; cc != col_nb_dofs; ++cc) {
1023
1024 FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
1025 t_RTP_dBubble(i, j) =
1026 a * t_row_base_fun * t_R(k, i) * t_col_base_fun(k, j);
1027 t_m(i, j) += t_diff_u(m, n, i, j) * t_RTP_dBubble(m, n);
1028
1029 ++t_m;
1030 ++t_col_base_fun;
1031 }
1032
1033 ++t_row_base_fun;
1034 }
1035 for (; rr != row_nb_base_functions; ++rr)
1036 ++t_row_base_fun;
1037 ++t_w;
1038 ++t_R;
1039 ++t_diff_u;
1040 }
1042}
1043
1045 EntData &col_data) {
1047
1048 int nb_integration_pts = row_data.getN().size1();
1049 int row_nb_dofs = row_data.getIndices().size();
1050 int col_nb_dofs = col_data.getIndices().size();
1051 auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
1053
1054 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1055 &m(r + 0, c + 4), &m(r + 0, c + 5),
1056
1057 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1058 &m(r + 1, c + 4), &m(r + 1, c + 5),
1059
1060 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1061 &m(r + 2, c + 4), &m(r + 2, c + 5),
1062
1063 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1064 &m(r + 3, c + 4), &m(r + 3, c + 5),
1065
1066 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1067 &m(r + 4, c + 4), &m(r + 4, c + 5),
1068
1069 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1070 &m(r + 5, c + 4), &m(r + 5, c + 5)
1071
1072 );
1073 };
1080
1084
1085 auto v = getVolume();
1086 auto t_w = getFTensor0IntegrationWeight();
1087 auto t_P_dh0 = getFTensor3FromMat(dataAtPts->P_dh0);
1088 auto t_P_dh1 = getFTensor3FromMat(dataAtPts->P_dh1);
1089 auto t_P_dh2 = getFTensor3FromMat(dataAtPts->P_dh2);
1090 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1091 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1092 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1093 auto t_dot_log_u =
1094 getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
1095 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1096 auto t_diff_u =
1097 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1098 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1099 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1100 auto &nbUniq = dataAtPts->nbUniq;
1101
1102 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1104 t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
1105
1106 FTensor::Tensor4<double, 3, 3, 3, 3> t_one4_symmetric;
1107 t_one4_symmetric(i, j, k, l) = 0;
1108 for (auto ii : {0, 1, 2})
1109 for (auto jj : {0, 1, 2})
1110 for (auto kk : {0, 1, 2})
1111 for (auto ll : {0, 1, 2}) {
1112
1113 if (ll != kk)
1114 t_one4_symmetric(ii, jj, kk, ll) =
1115 t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1116 else
1117 t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1118 }
1119
1120 int row_nb_base_functions = row_data.getN().size2();
1121 auto t_row_base_fun = row_data.getFTensor0N();
1122
1124 t_one(i, j) = 0;
1125 for (auto ii : {0, 1, 2})
1126 t_one(ii, ii) = 1;
1127
1128 const double ts_a = getTSa();
1129 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1130 double a = v * t_w;
1131
1133 t_P_dh(i, j, N0, k) = t_P_dh0(i, j, k);
1134 t_P_dh(i, j, N1, k) = t_P_dh1(i, j, k);
1135 t_P_dh(i, j, N2, k) = t_P_dh2(i, j, k);
1136
1138 t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
1140 t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
1141
1143 t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
1144
1146 t_stress_diff(i, j, k, l) = 0;
1147 for (int ii = 0; ii != 3; ++ii)
1148 for (int jj = 0; jj != 3; ++jj)
1149 for (int kk = 0; kk != 3; ++kk)
1150 for (int ll = 0; ll != 3; ++ll)
1151 for (int mm = 0; mm != 3; ++mm)
1152 for (int nn = 0; nn != 3; ++nn)
1153 t_stress_diff(ii, jj, mm, nn) -=
1154 t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1155
1157 t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
1158 FTensor::Tensor4<double, 3, 3, 3, 3> t_viscous_stress_diff;
1159 t_viscous_stress_diff(i, j, k, l) = t_one4_symmetric(i, j, k, l);
1160 t_viscous_stress_diff(i, j, k, l) *= (ts_a * alphaU);
1161
1163 t_deltaP2(i, j) = t_deltaP(i, j) - t_viscous_stress(i, j);
1164
1165 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1166 t_eigen_vals, t_eigen_vecs, f, d_f, dd_f, t_deltaP2, nbUniq[gg]);
1167 for (int ii = 0; ii != 3; ++ii) {
1168 for (int jj = 0; jj < ii; ++jj) {
1169 for (int kk = 0; kk != 3; ++kk) {
1170 for (int ll = 0; ll != kk; ++ll) {
1171 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1172 }
1173 }
1174 }
1175 }
1176 for (int ii = 0; ii != 3; ++ii) {
1177 for (int jj = ii; jj != 3; ++jj) {
1178 for (int kk = 0; kk != 3; ++kk) {
1179 for (int ll = 0; ll < kk; ++ll) {
1180 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1181 }
1182 }
1183 }
1184 }
1185
1186 int rr = 0;
1187 for (; rr != row_nb_dofs / 6; ++rr) {
1188
1189 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1190 auto t_m = get_ftensor4(K, 6 * rr, 0);
1191
1192 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1193
1194 const double b = a * t_row_base_fun * t_col_base_fun;
1195
1196 for (int ii = 0; ii != 3; ++ii)
1197 for (int jj = ii; jj != 3; ++jj)
1198 for (int kk = 0; kk != 3; ++kk)
1199 for (int ll = kk; ll != 3; ++ll)
1200 for (int mm = 0; mm != 3; ++mm)
1201 for (int nn = 0; nn != 3; ++nn)
1202 t_m(ii, jj, kk, ll) +=
1203 b * t_diff_u(mm, nn, ii, jj) *
1204 (t_stress_diff(mm, nn, kk, ll) -
1205 t_viscous_stress_diff(mm, nn, kk, ll));
1206
1207 t_m(i, j, k, l) += b * t_diff2_uP2(i, j, k, l);
1208
1209 ++t_m;
1210 ++t_col_base_fun;
1211 }
1212
1213 ++t_row_base_fun;
1214 }
1215 for (; rr != row_nb_base_functions; ++rr)
1216 ++t_row_base_fun;
1217
1218 ++t_w;
1219 ++t_P_dh0;
1220 ++t_P_dh1;
1221 ++t_P_dh2;
1222 ++t_R;
1223 ++t_approx_P;
1224 ++t_P;
1225 ++t_dot_log_u;
1226 ++t_u;
1227 ++t_diff_u;
1228 ++t_eigen_vals;
1229 ++t_eigen_vecs;
1230 }
1232}
1233
1235 EntData &col_data) {
1237 int row_nb_dofs = row_data.getIndices().size();
1238 int col_nb_dofs = col_data.getIndices().size();
1239 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1241
1242 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1243
1244 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1245
1246 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1247
1248 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1249
1250 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1251
1252 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1253
1254 );
1255 };
1260
1261 auto v = getVolume();
1262 auto t_w = getFTensor0IntegrationWeight();
1263 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1264 auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1265 auto t_diff_u =
1266 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1267
1268 int row_nb_base_functions = row_data.getN().size2();
1269 auto t_row_base_fun = row_data.getFTensor0N();
1270
1271 int nb_integration_pts = row_data.getN().size1();
1272
1273 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1274 double a = v * t_w;
1275
1277 t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
1279 t(i, j, k) = 0;
1280 for (int ii = 0; ii != 3; ++ii)
1281 for (int jj = ii; jj != 3; ++jj)
1282 for (int mm = 0; mm != 3; ++mm)
1283 for (int nn = 0; nn != 3; ++nn)
1284 for (int kk = 0; kk != 3; ++kk)
1285 t(ii, jj, kk) +=
1286 t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
1287
1288 int rr = 0;
1289 for (; rr != row_nb_dofs / 6; ++rr) {
1290 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1291 auto t_m = get_ftensor3(K, 6 * rr, 0);
1292 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1293
1294 double v = a * t_row_base_fun * t_col_base_fun;
1295 t_m(i, j, k) += v * t(i, j, k);
1296
1297 ++t_m;
1298 ++t_col_base_fun;
1299 }
1300
1301 ++t_row_base_fun;
1302 }
1303
1304 for (; rr != row_nb_base_functions; ++rr)
1305 ++t_row_base_fun;
1306
1307 ++t_w;
1308 ++t_P;
1309 ++t_diff_R;
1310 ++t_diff_u;
1311 }
1312
1314}
1315
1317 EntData &col_data) {
1319 int nb_integration_pts = row_data.getN().size1();
1320 int row_nb_dofs = row_data.getIndices().size();
1321 int col_nb_dofs = col_data.getIndices().size();
1322 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1324
1325 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1326
1327 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1328
1329 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1330
1331 );
1332 };
1339
1340 auto v = getVolume();
1341 auto t_w = getFTensor0IntegrationWeight();
1342 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1343
1344 int row_nb_base_functions = row_data.getN().size2();
1345 auto t_row_base_fun = row_data.getFTensor0N();
1346
1347 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1348 double a = v * t_w;
1349
1350 int rr = 0;
1351 for (; rr != row_nb_dofs / 3; ++rr) {
1352 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1353 auto t_m = get_ftensor2(K, 3 * rr, 0);
1355 t_levi_R(i, j, k) =
1356 (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1357 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1358 FTensor::Tensor1<double, 3> t_pushed_base;
1359 t_m(k, i) += t_levi_R(i, j, k) * t_col_base_fun(j);
1360 ++t_m;
1361 ++t_col_base_fun;
1362 }
1363
1364 ++t_row_base_fun;
1365 }
1366
1367 for (; rr != row_nb_base_functions; ++rr)
1368 ++t_row_base_fun;
1369 ++t_w;
1370 ++t_R;
1371 }
1373}
1374
1376 EntData &col_data) {
1378 int nb_integration_pts = row_data.getN().size1();
1379 int row_nb_dofs = row_data.getIndices().size();
1380 int col_nb_dofs = col_data.getIndices().size();
1381 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1383 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1384 };
1390
1391 auto v = getVolume();
1392 auto t_w = getFTensor0IntegrationWeight();
1393 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1394
1395 int row_nb_base_functions = row_data.getN().size2();
1396 auto t_row_base_fun = row_data.getFTensor0N();
1397
1398 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1399 double a = v * t_w;
1400
1401 int rr = 0;
1402 for (; rr != row_nb_dofs / 3; ++rr) {
1403
1404 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1405 auto t_m = get_ftensor1(K, 3 * rr, 0);
1406
1408 t_levi_r(i, j, k) =
1409 (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1410
1411 for (int cc = 0; cc != col_nb_dofs; ++cc) {
1412 FTensor::Tensor2<double, 3, 3> t_pushed_base;
1413 t_m(k) += t_levi_r(i, j, k) * t_col_base_fun(i, j);
1414 ++t_m;
1415 ++t_col_base_fun;
1416 }
1417
1418 ++t_row_base_fun;
1419 }
1420
1421 for (; rr != row_nb_base_functions; ++rr)
1422 ++t_row_base_fun;
1423 ++t_w;
1424 ++t_R;
1425 }
1427}
1428
1430 EntData &col_data) {
1432 int nb_integration_pts = row_data.getN().size1();
1433 int row_nb_dofs = row_data.getIndices().size();
1434 int col_nb_dofs = col_data.getIndices().size();
1435 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1437
1438 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1439
1440 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1441
1442 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1443
1444 );
1445 };
1452
1453 auto v = getVolume();
1454 auto t_w = getFTensor0IntegrationWeight();
1455 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1456 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1457 auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1458
1459 int row_nb_base_functions = row_data.getN().size2();
1460 auto t_row_base_fun = row_data.getFTensor0N();
1461
1462 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1463 double a = v * t_w;
1464
1466 t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
1467
1468 FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
1469 t_leviPRT_domega(n, l) =
1470 levi_civita(i, j, n) * (t_approx_P(i, k) * t_diff_R(j, k, l));
1471
1472 int rr = 0;
1473 for (; rr != row_nb_dofs / 3; ++rr) {
1474
1475 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1476 auto t_m = get_ftensor2(K, 3 * rr, 0);
1477 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1478
1479 const double b = a * t_row_base_fun * t_col_base_fun;
1480
1481 t_m(k, m) += b * t_leviPRT_domega(k, m);
1482
1483 ++t_m;
1484 ++t_col_base_fun;
1485 }
1486
1487 ++t_row_base_fun;
1488 }
1489
1490 for (; rr != row_nb_base_functions; ++rr)
1491 ++t_row_base_fun;
1492 ++t_w;
1493 ++t_approx_P;
1494 ++t_R;
1495 ++t_diff_R;
1496 }
1498}
1499
1501 EntData &col_data) {
1503
1504 int nb_integration_pts = row_data.getN().size1();
1505 int row_nb_dofs = row_data.getIndices().size();
1506 int col_nb_dofs = col_data.getIndices().size();
1507
1508 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1510
1511 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1512
1513 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1514
1515 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1516
1517 );
1518 };
1519
1526
1527 auto v = getVolume();
1528 auto t_w = getFTensor0IntegrationWeight();
1529 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1530 auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1531 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1532 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
1533 int row_nb_base_functions = row_data.getN().size2() / 3;
1534 auto t_row_base_fun = row_data.getFTensor1N<3>();
1535
1536 const double ts_a = getTSa();
1537
1538 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1539 double a = v * t_w;
1540
1541 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1542
1543 int rr = 0;
1544 for (; rr != row_nb_dofs / 3; ++rr) {
1545
1547 t_PRT(i, k) = t_row_base_fun(j) * (t_diff_R(i, l, k) * t_u(l, j));
1548
1549 FTensor::Tensor2<double, 3, 3> t_levi1, t_levi2;
1550 t_levi1(i, k) = (levi_civita(i, m, n) * t_omega_dot(n)) *
1551 (t_row_base_fun(j) * t_diff_R(m, j, k));
1552 t_levi2(i, k) = levi_civita(i, m, k) * (t_row_base_fun(j) * t_R(m, j));
1553
1554 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1555 auto t_m = get_ftensor2(K, 3 * rr, 0);
1556 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1557 t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1558 t_m(i, j) += (a * t_col_base_fun) * t_levi1(i, j);
1559 t_m(i, j) += ((a * ts_a) * t_col_base_fun) * t_levi2(i, j);
1560 ++t_m;
1561 ++t_col_base_fun;
1562 }
1563
1564 ++t_row_base_fun;
1565 }
1566
1567 for (; rr != row_nb_base_functions; ++rr)
1568 ++t_row_base_fun;
1569 ++t_w;
1570 ++t_R;
1571 ++t_diff_R;
1572 ++t_u;
1573 ++t_omega_dot;
1574 }
1576}
1577
1580 EntData &col_data) {
1582
1583 int nb_integration_pts = row_data.getN().size1();
1584 int row_nb_dofs = row_data.getIndices().size();
1585 int col_nb_dofs = col_data.getIndices().size();
1586
1587 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1589 &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1590 };
1591
1598
1599 auto v = getVolume();
1600 auto t_w = getFTensor0IntegrationWeight();
1601 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1602 auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1603 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1604 auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
1605 int row_nb_base_functions = row_data.getN().size2() / 9;
1606 auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1607
1608 const double ts_a = getTSa();
1609
1610 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1611 double a = v * t_w;
1612
1613 int rr = 0;
1614 for (; rr != row_nb_dofs; ++rr) {
1615
1617 t_PRT(k) = t_row_base_fun(i, j) * (t_diff_R(i, l, k) * t_u(l, j));
1618
1619 FTensor::Tensor2<double, 3, 3> t_levi_omega;
1620 t_levi_omega(i, m) = levi_civita(i, m, n) * t_omega_dot(n);
1621 FTensor::Tensor3<double, 3, 3, 3> t_levi_omegaDiffR;
1622 t_levi_omegaDiffR(i, j, k) = t_levi_omega(i, m) * t_diff_R(m, j, k);
1624 t_levi_omegaR(i, j, k) = levi_civita(i, m, k) * t_R(m, j);
1625
1626 FTensor::Tensor1<double, 3> t_levi1, t_levi2;
1627 t_levi1(k) = t_row_base_fun(i, j) * t_levi_omegaDiffR(i, j, k);
1628 t_levi2(k) = t_row_base_fun(i, j) * t_levi_omegaR(i, j, k);
1629
1630 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1631 auto t_m = get_ftensor2(K, rr, 0);
1632 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1633 t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1634 t_m(j) += (a * t_col_base_fun) * t_levi1(j);
1635 t_m(j) += ((a * ts_a) * t_col_base_fun) * t_levi2(j);
1636 ++t_m;
1637 ++t_col_base_fun;
1638 }
1639
1640 ++t_row_base_fun;
1641 }
1642
1643 for (; rr != row_nb_base_functions; ++rr)
1644 ++t_row_base_fun;
1645
1646 ++t_w;
1647 ++t_R;
1648 ++t_diff_R;
1649 ++t_u;
1650 ++t_omega_dot;
1651 }
1653}
1654
1656 EntData &data) {
1658 if (type != MBVERTEX)
1660
1661 auto create_tag = [this](const std::string tag_name, const int size) {
1662 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1663 Tag th;
1664 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
1665 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1666 def_VAL);
1667 return th;
1668 };
1669
1670 Tag th_w = create_tag("SpatialDisplacement", 3);
1671 Tag th_omega = create_tag("Omega", 3);
1672 Tag th_approxP = create_tag("Piola1Stress", 9);
1673 Tag th_sigma = create_tag("CauchyStress", 9);
1674 Tag th_log_u = create_tag("LogSpatialStreach", 9);
1675 Tag th_u = create_tag("SpatialStreach", 9);
1676 Tag th_h = create_tag("h", 9);
1677 Tag th_X = create_tag("X", 3);
1678 Tag th_detF = create_tag("detF", 1);
1679 Tag th_angular_momentum = create_tag("AngularMomentum", 3);
1680
1681 Tag th_u_eig_vec = create_tag("SpatialStreachEigenVec", 9);
1682 Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
1683
1684 int nb_gauss_pts = data.getN().size1();
1685 if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
1686 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1687 "Nb. of integration points is not equal to number points on "
1688 "post-processing mesh");
1689 }
1690
1691 auto t_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
1692 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
1693 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1694 auto t_log_u =
1695 getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
1696 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1697 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1698 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1699 auto t_coords = getFTensor1CoordsAtGaussPts();
1700
1705
1706 // vectors
1707 VectorDouble3 v(3);
1708 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
1709 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
1711 t_v(i) = t_d(i);
1712 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1713 &*v.data().begin());
1715 };
1716
1717 MatrixDouble3by3 m(3, 3);
1719 &m(0, 0), &m(0, 1), &m(0, 2),
1720
1721 &m(1, 0), &m(1, 1), &m(1, 2),
1722
1723 &m(2, 0), &m(2, 1), &m(2, 2));
1724
1725 auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
1727 t_m(i, j) = t_d(i, j);
1728 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1729 &*m.data().begin());
1731 };
1732
1733 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1734
1735 // vetors
1736 CHKERR save_vec_tag(th_w, t_w, gg);
1737 CHKERR save_vec_tag(th_X, t_coords, gg);
1738 CHKERR save_vec_tag(th_omega, t_omega, gg);
1739
1740 // tensors
1741 CHKERR save_mat_tag(th_h, t_h, gg);
1742
1744 for (int ii = 0; ii != 3; ++ii)
1745 for (int jj = 0; jj != 3; ++jj)
1746 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
1747
1748 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
1749
1751 for (int ii = 0; ii != 3; ++ii)
1752 for (int jj = 0; jj != 3; ++jj)
1753 t_u_tmp(ii, jj) = t_u(ii, jj);
1754
1755 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
1756 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
1757
1758 const double jac = determinantTensor3by3(t_h);
1760 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1761 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
1762 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
1763
1765 t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
1767 t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
1768
1769 CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
1770 &t_leviPRT(0));
1771
1772 auto get_eiegn_vector_symmetric = [&](auto &t_u) {
1774
1775 for (int ii = 0; ii != 3; ++ii)
1776 for (int jj = 0; jj != 3; ++jj)
1777 t_m(ii, jj) = t_u(ii, jj);
1778
1779 VectorDouble3 eigen_values(3);
1780
1781 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
1782 // memory space
1783 int n = 3, lda = 3, info, lwork = -1;
1784 double wkopt;
1785 info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
1786 &(eigen_values.data()[0]), &wkopt, lwork);
1787 if (info != 0)
1788 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1789 "is something wrong with lapack_dsyev info = %d", info);
1790 lwork = (int)wkopt;
1791 double work[lwork];
1792 info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
1793 &(eigen_values.data()[0]), work, lwork);
1794 if (info != 0)
1795 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1796 "is something wrong with lapack_dsyev info = %d", info);
1797
1798 CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
1799 &*m.data().begin());
1800 CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
1801 &*eigen_values.data().begin());
1802
1804 };
1805
1806 CHKERR get_eiegn_vector_symmetric(t_u);
1807
1808 ++t_w;
1809 ++t_h;
1810 ++t_log_u;
1811 ++t_u;
1812 ++t_omega;
1813 ++t_R;
1814 ++t_approx_P;
1815 ++t_coords;
1816 }
1817
1819}
1820
1825 if (type == MBTET) {
1826 int nb_integration_pts = data.getN().size1();
1827 auto v = getVolume();
1828 auto t_w = getFTensor0IntegrationWeight();
1829 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1830 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1831
1834
1835 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1836 const double a = t_w * v;
1837 (*energy) += a * t_P(i, J) * t_h(i, J);
1838 ++t_w;
1839 ++t_P;
1840 ++t_h;
1841 }
1842 }
1844}
1845
1846} // 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.
@ 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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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)
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
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:261
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
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)
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
MoFEMErrorCode get_eigen_val_and_proj_lapack(T1 &X, T2 &eig, T3 &eig_vec)
bool is_eq(const double &a, const double &b)
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
double f(const double v)
void sort_eigen_vals(T1 &eig, T2 &eigen_vec)
double dd_f(const double v)
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
double d_f(const double v)
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
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
constexpr double t
plate stiffness
Definition: plate.cpp:59
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 integrate(EntData &row_data, EntData &col_data)
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 &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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