v0.15.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
13#include <boost/math/constants/constants.hpp>
14
15#include <EshelbianAux.hpp>
16
17#include <lapack_wrap.h>
18
19#include <Lie.hpp>
20#include <MatrixFunction.hpp>
21
22#ifdef ENABLE_PYTHON_BINDING
23#include <boost/python.hpp>
24#include <boost/python/def.hpp>
25#include <boost/python/numpy.hpp>
26namespace bp = boost::python;
27namespace np = boost::python::numpy;
28#endif
29
30namespace EshelbianPlasticity {
31
32inline MatrixDouble analytical_expr_function(double delta_t, double t,
33 int nb_gauss_pts,
34 MatrixDouble &m_ref_coords,
35 MatrixDouble &m_ref_normals,
36 const std::string block_name);
37
38template <typename OP_PTR>
39auto getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr,
40 const std::string block_name) {
41
42 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
43
44 auto ts_time = op_ptr->getTStime();
45 auto ts_time_step = op_ptr->getTStimeStep();
48 ts_time_step = EshelbianCore::dynamicStep;
49 }
50
51 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
52 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
53
54 auto v_analytical_expr = analytical_expr_function(
55 ts_time_step, ts_time, nb_gauss_pts, m_ref_coords, m_ref_normals, block_name);
56
57#ifndef NDEBUG
58 if (v_analytical_expr.size2() != nb_gauss_pts)
60 "Wrong number of integration pts");
61#endif // NDEBUG
62
63 return std::make_tuple(block_name, v_analytical_expr);
64};
65
67 EntData &data) {
72
73 int nb_integration_pts = getGaussPts().size2();
74
77 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
78
79 dataAtPts->SigmaAtPts.resize(SPACE_DIM * SPACE_DIM, nb_integration_pts,
80 false);
81 auto t_eshelby_stress =
83
85
86 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
87 t_eshelby_stress(i, j) = t_energy * t_kd(i, j) - t_F(m, i) * t_P(m, j);
88 ++t_energy;
89 ++t_P;
90 ++t_F;
91 ++t_eshelby_stress;
92 }
93
95}
96
98 EntityType type,
99 EntData &data) {
101
102 auto ts_ctx = getTSCtx();
103 int nb_integration_pts = getGaussPts().size2();
104
105 // space size indices
113
114 // sym size indices
116
117 auto t_L = symm_L_tensor();
118
119 dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
120 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
121 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
122
123 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
124 dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts, false);
125 dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 * size_symm,
126 nb_integration_pts, false);
127 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
128
129 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
130 dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
131 dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
132 dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
133 false);
134 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
135 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
136 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
137 dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
138 dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
139 dataAtPts->nbUniq.resize(nb_integration_pts, false);
140 dataAtPts->eigenValsC.resize(3, nb_integration_pts, false);
141 dataAtPts->eigenVecsC.resize(9, nb_integration_pts, false);
142 dataAtPts->nbUniqC.resize(nb_integration_pts, false);
143
144 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
145 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
146
147 dataAtPts->internalStressAtPts.resize(9, nb_integration_pts, false);
148 dataAtPts->internalStressAtPts.clear();
149
150 // Calculated values
151 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
152 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
153 auto t_h_dlog_u =
155 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
156 auto t_levi_kirchhoff_domega =
157 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
158 auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
159 dataAtPts->leviKirchhoffdLogStreatchAtPts);
160 auto t_levi_kirchhoff_dP =
161 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
162 auto t_approx_P_adjoint__dstretch =
163 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
164 auto t_approx_P_adjoint_log_du =
166 auto t_approx_P_adjoint_log_du_dP =
168 auto t_approx_P_adjoint_log_du_domega =
169 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
170 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
171 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
172 auto t_diff_u =
173 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
174 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
175 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
176 auto &nbUniq = dataAtPts->nbUniq;
177 auto t_nb_uniq =
178 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
179 auto t_eigen_vals_C = getFTensor1FromMat<3>(dataAtPts->eigenValsC);
180 auto t_eigen_vecs_C = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecsC);
181 auto &nbUniqC = dataAtPts->nbUniqC;
182 auto t_nb_uniq_C =
183 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniqC.data().data());
184
185 auto t_log_stretch_total =
186 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
187 auto t_log_u2_h1 =
188 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
189
190 // Field values
191 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
192 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
193 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
194 auto t_log_u =
195 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
196
197 auto next = [&]() {
198 // calculated values
199 ++t_h;
200 ++t_h_domega;
201 ++t_h_dlog_u;
202 ++t_levi_kirchhoff;
203 ++t_levi_kirchhoff_domega;
204 ++t_levi_kirchhoff_dstreach;
205 ++t_levi_kirchhoff_dP;
206 ++t_approx_P_adjoint__dstretch;
207 ++t_approx_P_adjoint_log_du;
208 ++t_approx_P_adjoint_log_du_dP;
209 ++t_approx_P_adjoint_log_du_domega;
210 ++t_R;
211 ++t_u;
212 ++t_diff_u;
213 ++t_eigen_vals;
214 ++t_eigen_vecs;
215 ++t_nb_uniq;
216 ++t_eigen_vals_C;
217 ++t_eigen_vecs_C;
218 ++t_nb_uniq_C;
219 ++t_log_u2_h1;
220 ++t_log_stretch_total;
221 // field values
222 ++t_omega;
223 ++t_grad_h1;
224 ++t_approx_P;
225 ++t_log_u;
226 };
227
230
231 auto bound_eig = [&](auto &eig) {
233 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
234 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
235 for (int ii = 0; ii != 3; ++ii) {
236 eig(ii) = std::min(std::max(zero, eig(ii)), large);
237 }
239 };
240
241 auto calculate_log_stretch = [&]() {
244 eigen_vec(i, j) = t_log_u(i, j);
245 if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
246 MOFEM_LOG("SELF", Sev::error) << "Failed to compute eigen values";
247 }
248 // CHKERR bound_eig(eig);
249 // rare case when two eigen values are equal
250 t_nb_uniq = get_uniq_nb<3>(&eig(0));
251 if (t_nb_uniq < 3) {
252 sort_eigen_vals(eig, eigen_vec);
253 }
254 t_eigen_vals(i) = eig(i);
255 t_eigen_vecs(i, j) = eigen_vec(i, j);
256 t_u(i, j) =
257 EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f)(i, j);
258 auto get_t_diff_u = [&]() {
259 return EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs,
261 t_nb_uniq);
262 };
263 t_diff_u(i, j, k, l) = get_t_diff_u()(i, j, k, l);
265 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
266 return t_Ldiff_u;
267 };
268
269 auto calculate_total_stretch = [&](auto &t_h1) {
271
272 t_log_u2_h1(i, j) = 0;
273 t_log_stretch_total(i, j) = t_log_u(i, j);
274
277 t_R_h1(i, j) = t_kd(i, j);
278 t_inv_u_h1(i, j) = t_symm_kd(i, j);
279
280 return std::make_pair(t_R_h1, t_inv_u_h1);
281
282 } else {
283
286
288 t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
289 eigen_vec(i, j) = t_C_h1(i, j);
290 if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
291 MOFEM_LOG("SELF", Sev::error) << "Failed to compute eigen values";
292 }
293 // rare case when two eigen values are equal
294 t_nb_uniq_C = get_uniq_nb<3>(&eig(0));
295 if (t_nb_uniq_C < 3) {
296 sort_eigen_vals(eig, eigen_vec);
297 }
299 CHKERR bound_eig(eig);
300 }
301 t_eigen_vals_C(i) = eig(i);
302 t_eigen_vecs_C(i, j) = eigen_vec(i, j);
303
304 t_log_u2_h1(i, j) =
305 EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f)(i, j);
306 t_log_stretch_total(i, j) = t_log_u2_h1(i, j) / 2 + t_log_u(i, j);
307
308 auto f_inv_sqrt = [](auto e) { return 1. / std::sqrt(e); };
310 t_inv_u_h1(i, j) = EigenMatrix::getMat(eig, eigen_vec, f_inv_sqrt)(i, j);
312 t_R_h1(i, j) = t_h1(i, o) * t_inv_u_h1(o, j);
313
314 return std::make_pair(t_R_h1, t_inv_u_h1);
315 }
316 };
317
318 auto no_h1_loop = [&]() {
320
322 case LARGE_ROT:
323 break;
324 case SMALL_ROT:
325 break;
326 default:
327 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
328 "rotSelector should be large");
329 };
330
331 for (int gg = 0; gg != nb_integration_pts; ++gg) {
332
334
336 t_h1(i, j) = t_kd(i, j);
337
338 // calculate streach
339 auto t_Ldiff_u = calculate_log_stretch();
340
341 // calculate total stretch
342 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
343
348
350
351 // rotation
353 case LARGE_ROT:
354 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
355 t_diff_R(i, j, k) =
356 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
357 // left
358 t_diff_Rr(i, j, l) = t_R(i, k) * levi_civita(k, j, l);
359 t_diff_diff_Rr(i, j, l, m) = t_diff_R(i, k, m) * levi_civita(k, j, l);
360 // right
361 t_diff_Rl(i, j, l) = levi_civita(i, k, l) * t_R(k, j);
362 t_diff_diff_Rl(i, j, l, m) = levi_civita(i, k, l) * t_diff_R(k, j, m);
363 break;
364
365 default:
366 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
367 "rotationSelector not handled");
368 }
369
370 constexpr double half_r = 1 / 2.;
371 constexpr double half_l = 1 / 2.;
372
373 // calculate gradient
374 t_h(i, k) = t_R(i, l) * t_u(l, k);
375
376 // Adjoint stress
377 t_approx_P_adjoint__dstretch(l, k) =
378 (t_R(i, l) * t_approx_P(i, k) + t_R(i, k) * t_approx_P(i, l));
379 t_approx_P_adjoint__dstretch(l, k) /= 2.;
380
381 t_approx_P_adjoint_log_du(L) =
382 (t_approx_P_adjoint__dstretch(l, k) * t_Ldiff_u(l, k, L) +
383 t_approx_P_adjoint__dstretch(k, l) * t_Ldiff_u(l, k, L) +
384 t_approx_P_adjoint__dstretch(l, k) * t_Ldiff_u(k, l, L) +
385 t_approx_P_adjoint__dstretch(k, l) * t_Ldiff_u(k, l, L)) /
386 4.;
387
388 // Kirchhoff stress
389 t_levi_kirchhoff(m) =
390
391 half_r * (t_diff_Rr(i, l, m) * (t_u(l, k) * t_approx_P(i, k)) +
392 t_diff_Rr(i, k, m) * (t_u(l, k) * t_approx_P(i, l)))
393
394 +
395
396 half_l * (t_diff_Rl(i, l, m) * (t_u(l, k) * t_approx_P(i, k)) +
397 t_diff_Rl(i, k, m) * (t_u(l, k) * t_approx_P(i, l)));
398 t_levi_kirchhoff(m) /= 2.;
399
401
403 t_h_domega(i, k, m) = half_r * (t_diff_Rr(i, l, m) * t_u(l, k))
404
405 +
406
407 half_l * (t_diff_Rl(i, l, m) * t_u(l, k));
408 } else {
409 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u(l, k);
410 }
411
412 t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u(l, k, L);
413
414 t_approx_P_adjoint_log_du_dP(i, k, L) =
415 t_R(i, l) * (t_Ldiff_u(l, k, L) + t_Ldiff_u(k, l, L)) / 2.;
416
419 t_A(k, l, m) = t_diff_Rr(i, l, m) * t_approx_P(i, k) +
420 t_diff_Rr(i, k, m) * t_approx_P(i, l);
421 t_A(k, l, m) /= 2.;
423 t_B(k, l, m) = t_diff_Rl(i, l, m) * t_approx_P(i, k) +
424 t_diff_Rl(i, k, m) * t_approx_P(i, l);
425 t_B(k, l, m) /= 2.;
426
427 t_approx_P_adjoint_log_du_domega(m, L) =
428 half_r * (t_A(k, l, m) *
429 (t_Ldiff_u(k, l, L) + t_Ldiff_u(k, l, L)) / 2) +
430 half_l * (t_B(k, l, m) *
431 (t_Ldiff_u(k, l, L) + t_Ldiff_u(k, l, L)) / 2);
432 } else {
434 t_A(k, l, m) = t_diff_R(i, l, m) * t_approx_P(i, k);
435 t_approx_P_adjoint_log_du_domega(m, L) =
436 t_A(k, l, m) * t_Ldiff_u(k, l, L);
437 }
438
439 t_levi_kirchhoff_dstreach(m, L) =
440 half_r *
441 (t_diff_Rr(i, l, m) * (t_Ldiff_u(l, k, L) * t_approx_P(i, k)))
442
443 +
444
445 half_l *
446 (t_diff_Rl(i, l, m) * (t_Ldiff_u(l, k, L) * t_approx_P(i, k)));
447
448 t_levi_kirchhoff_dP(m, i, k) =
449 half_r * t_diff_Rr(i, l, m) * (t_u(l, k))
450
451 +
452
453 half_l * t_diff_Rl(i, l, m) * (t_u(l, k));
454
455 t_levi_kirchhoff_domega(m, n) =
456 half_r *
457 (t_diff_diff_Rr(i, l, m, n) * (t_u(l, k) * t_approx_P(i, k)) +
458 t_diff_diff_Rr(i, k, m, n) * (t_u(l, k) * t_approx_P(i, l)))
459
460 +
461
462 half_l *
463 (t_diff_diff_Rl(i, l, m, n) * (t_u(l, k) * t_approx_P(i, k)) +
464 t_diff_diff_Rl(i, k, m, n) * (t_u(l, k) * t_approx_P(i, l)));
465 t_levi_kirchhoff_domega(m, n) /= 2.;
466 }
467
468 next();
469 }
470
472 };
473
474 auto large_loop = [&]() {
476
478 case LARGE_ROT:
479 break;
480 case SMALL_ROT:
481 break;
482 default:
483 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
484 "rotSelector should be large");
485 };
486
487 for (int gg = 0; gg != nb_integration_pts; ++gg) {
488
490
494 t_h1(i, j) = t_kd(i, j);
495 break;
496 case LARGE_ROT:
497 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
498 break;
499 default:
500 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
501 "gradApproximator not handled");
502 };
503
504 // calculate streach
505 auto t_Ldiff_u = calculate_log_stretch();
506 // calculate total stretch
507 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
508
510 t_h_u(l, k) = t_u(l, o) * t_h1(o, k);
512 t_Ldiff_h_u(l, k, L) = t_Ldiff_u(l, o, L) * t_h1(o, k);
513
518
520
521 // rotation
523 case LARGE_ROT:
524 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
525 t_diff_R(i, j, k) =
526 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
527 // left
528 t_diff_Rr(i, j, l) = t_R(i, k) * levi_civita(k, j, l);
529 t_diff_diff_Rr(i, j, l, m) = t_diff_R(i, k, m) * levi_civita(k, j, l);
530 // right
531 t_diff_Rl(i, j, l) = levi_civita(i, k, l) * t_R(k, j);
532 t_diff_diff_Rl(i, j, l, m) = levi_civita(i, k, l) * t_diff_R(k, j, m);
533 break;
534 case SMALL_ROT:
535 t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
536 t_diff_R(i, j, k) = levi_civita(i, j, k);
537 // left
538 t_diff_Rr(i, j, l) = levi_civita(i, j, l);
539 t_diff_diff_Rr(i, j, l, m) = 0;
540 // right
541 t_diff_Rl(i, j, l) = levi_civita(i, j, l);
542 t_diff_diff_Rl(i, j, l, m) = 0;
543 break;
544
545 default:
546 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
547 "rotationSelector not handled");
548 }
549
550 constexpr double half_r = 1 / 2.;
551 constexpr double half_l = 1 / 2.;
552
553 // calculate gradient
554 t_h(i, k) = t_R(i, l) * t_h_u(l, k);
555
556 // Adjoint stress
557 t_approx_P_adjoint__dstretch(l, o) =
558 (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
559 t_approx_P_adjoint_log_du(L) =
560 t_approx_P_adjoint__dstretch(l, o) * t_Ldiff_u(l, o, L);
561
562 // Kirchhoff stress
563 t_levi_kirchhoff(m) =
564
565 half_r * ((t_diff_Rr(i, l, m) * (t_h_u(l, k)) * t_approx_P(i, k)))
566
567 +
568
569 half_l * ((t_diff_Rl(i, l, m) * (t_h_u(l, k)) * t_approx_P(i, k)));
570
572
574 t_h_domega(i, k, m) = half_r * (t_diff_Rr(i, l, m) * t_h_u(l, k))
575
576 +
577
578 half_l * (t_diff_Rl(i, l, m) * t_h_u(l, k));
579 } else {
580 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_h_u(l, k);
581 }
582
583 t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_h_u(l, k, L);
584
585 t_approx_P_adjoint_log_du_dP(i, k, L) =
586 t_R(i, l) * t_Ldiff_h_u(l, k, L);
587
590 t_A(m, L, i, k) = t_diff_Rr(i, l, m) * t_Ldiff_h_u(l, k, L);
592 t_B(m, L, i, k) = t_diff_Rl(i, l, m) * t_Ldiff_h_u(l, k, L);
593
594 t_approx_P_adjoint_log_du_domega(m, L) =
595 half_r * (t_A(m, L, i, k) * t_approx_P(i, k))
596
597 +
598
599 half_l * (t_B(m, L, i, k) * t_approx_P(i, k));
600 } else {
602 t_A(m, L, i, k) = t_diff_R(i, l, m) * t_Ldiff_h_u(l, k, L);
603 t_approx_P_adjoint_log_du_domega(m, L) =
604 t_A(m, L, i, k) * t_approx_P(i, k);
605 }
606
607 t_levi_kirchhoff_dstreach(m, L) =
608 half_r *
609 (t_diff_Rr(i, l, m) * (t_Ldiff_h_u(l, k, L) * t_approx_P(i, k)))
610
611 +
612
613 half_l * (t_diff_Rl(i, l, m) *
614 (t_Ldiff_h_u(l, k, L) * t_approx_P(i, k)));
615
616 t_levi_kirchhoff_dP(m, i, k) =
617
618 half_r * (t_diff_Rr(i, l, m) * t_h_u(l, k))
619
620 +
621
622 half_l * (t_diff_Rl(i, l, m) * t_h_u(l, k));
623
624 t_levi_kirchhoff_domega(m, n) =
625 half_r *
626 (t_diff_diff_Rr(i, l, m, n) * (t_h_u(l, k) * t_approx_P(i, k)))
627
628 +
629
630 half_l *
631 (t_diff_diff_Rl(i, l, m, n) * (t_h_u(l, k) * t_approx_P(i, k)));
632 }
633
634 next();
635 }
636
638 };
639
640 auto moderate_loop = [&]() {
642
644 case LARGE_ROT:
645 break;
646 case SMALL_ROT:
647 break;
648 default:
649 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
650 "rotSelector should be large");
651 };
652
653 for (int gg = 0; gg != nb_integration_pts; ++gg) {
654
656
659 case MODERATE_ROT:
660 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
661 break;
662 default:
663 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
664 "gradApproximator not handled");
665 };
666
667 // calculate streach
668 auto t_Ldiff_u = calculate_log_stretch();
669 // calculate total stretch
670 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
671
673 t_h_u(l, k) = (t_symm_kd(l, o) + t_log_u(l, o)) * t_h1(o, k);
675 t_L_h(l, k, L) = t_L(l, o, L) * t_h1(o, k);
676
681
683
684 // rotation
686 case LARGE_ROT:
687 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
688 t_diff_R(i, j, k) =
689 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
690 // left
691 t_diff_Rr(i, j, l) = t_R(i, k) * levi_civita(k, j, l);
692 t_diff_diff_Rr(i, j, l, m) = t_diff_R(i, k, m) * levi_civita(k, j, l);
693 // right
694 t_diff_Rl(i, j, l) = levi_civita(i, k, l) * t_R(k, j);
695 t_diff_diff_Rl(i, j, l, m) = levi_civita(i, k, l) * t_diff_R(k, j, m);
696 break;
697 case SMALL_ROT:
698 t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
699 t_diff_R(i, j, l) = levi_civita(i, j, l);
700 // left
701 t_diff_Rr(i, j, l) = levi_civita(i, j, l);
702 t_diff_diff_Rr(i, j, l, m) = 0;
703 // right
704 t_diff_Rl(i, j, l) = levi_civita(i, j, l);
705 t_diff_diff_Rl(i, j, l, m) = 0;
706 break;
707
708 default:
709 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
710 "rotationSelector not handled");
711 }
712
713 constexpr double half_r = 1 / 2.;
714 constexpr double half_l = 1 / 2.;
715
716 // calculate gradient
717 t_h(i, k) = t_R(i, l) * t_h_u(l, k);
718
719 // Adjoint stress
720 t_approx_P_adjoint__dstretch(l, o) =
721 (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
722
723 t_approx_P_adjoint_log_du(L) =
724 t_approx_P_adjoint__dstretch(l, o) * t_L(l, o, L);
725
726 // Kirchhoff stress
727 t_levi_kirchhoff(m) =
728
729 half_r * ((t_diff_Rr(i, l, m) * (t_h_u(l, k)) * t_approx_P(i, k)))
730
731 +
732
733 half_l * ((t_diff_Rl(i, l, m) * (t_h_u(l, k)) * t_approx_P(i, k)));
734
736
738 t_h_domega(i, k, m) = half_r * (t_diff_Rr(i, l, m) * t_h_u(l, k))
739
740 +
741
742 half_l * (t_diff_Rl(i, l, m) * t_h_u(l, k));
743 } else {
744 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_h_u(l, k);
745 }
746
747 t_h_dlog_u(i, k, L) = t_R(i, l) * t_L_h(l, k, L);
748
749 t_approx_P_adjoint_log_du_dP(i, k, L) = t_R(i, l) * t_L_h(l, k, L);
750
753 t_A(m, L, i, k) = t_diff_Rr(i, l, m) * t_L_h(l, k, L);
755 t_B(m, L, i, k) = t_diff_Rl(i, l, m) * t_L_h(l, k, L);
756 t_approx_P_adjoint_log_du_domega(m, L) =
757 half_r * (t_A(m, L, i, k) * t_approx_P(i, k))
758
759 +
760
761 half_l * (t_B(m, L, i, k) * t_approx_P(i, k));
762 } else {
764 t_A(m, L, i, k) = t_diff_R(i, l, m) * t_L_h(l, k, L);
765 t_approx_P_adjoint_log_du_domega(m, L) =
766 t_A(m, L, i, k) * t_approx_P(i, k);
767 }
768
769 t_levi_kirchhoff_dstreach(m, L) =
770 half_r * (t_diff_Rr(i, l, m) * (t_L_h(l, k, L) * t_approx_P(i, k)))
771
772 +
773
774 half_l * (t_diff_Rl(i, l, m) * (t_L_h(l, k, L) * t_approx_P(i, k)));
775
776 t_levi_kirchhoff_dP(m, i, k) =
777
778 half_r * (t_diff_Rr(i, l, m) * t_h_u(l, k))
779
780 +
781
782 half_l * (t_diff_Rl(i, l, m) * t_h_u(l, k));
783
784 t_levi_kirchhoff_domega(m, n) =
785 half_r *
786 (t_diff_diff_Rr(i, l, m, n) * (t_h_u(l, k) * t_approx_P(i, k)))
787
788 +
789
790 half_l *
791 (t_diff_diff_Rl(i, l, m, n) * (t_h_u(l, k) * t_approx_P(i, k)));
792 }
793
794 next();
795 }
796
798 };
799
800 auto small_loop = [&]() {
803 case SMALL_ROT:
804 break;
805 default:
806 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
807 "rotSelector should be small");
808 };
809
810 for (int gg = 0; gg != nb_integration_pts; ++gg) {
811
814 case SMALL_ROT:
815 t_h1(i, j) = t_kd(i, j);
816 break;
817 default:
818 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
819 "gradApproximator not handled");
820 };
821
822 // calculate streach
824
826 t_Ldiff_u(i, j, L) = calculate_log_stretch()(i, j, L);
827 } else {
828 t_u(i, j) = t_symm_kd(i, j) + t_log_u(i, j);
829 t_Ldiff_u(i, j, L) = t_L(i, j, L);
830 }
831 t_log_u2_h1(i, j) = 0;
832 t_log_stretch_total(i, j) = t_log_u(i, j);
833
835 t_h(i, j) = levi_civita(i, j, k) * t_omega(k) + t_u(i, j);
836
837 t_h_domega(i, j, k) = levi_civita(i, j, k);
838 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
839
840 // Adjoint stress
841 t_approx_P_adjoint__dstretch(i, j) = t_approx_P(i, j);
842 t_approx_P_adjoint_log_du(L) =
843 t_approx_P_adjoint__dstretch(i, j) * t_Ldiff_u(i, j, L);
844 t_approx_P_adjoint_log_du_dP(i, j, L) = t_Ldiff_u(i, j, L);
845 t_approx_P_adjoint_log_du_domega(m, L) = 0;
846
847 // Kirchhoff stress
848 t_levi_kirchhoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
849 t_levi_kirchhoff_dstreach(m, L) = 0;
850 t_levi_kirchhoff_dP(k, i, j) = levi_civita(i, j, k);
851 t_levi_kirchhoff_domega(m, n) = 0;
852
853 next();
854 }
855
857 };
858
861 CHKERR no_h1_loop();
862 break;
863 case LARGE_ROT:
864 CHKERR large_loop();
866 break;
867 case MODERATE_ROT:
868 CHKERR moderate_loop();
870 break;
871 case SMALL_ROT:
872 CHKERR small_loop();
874 break;
875 default:
876 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
877 "gradApproximator not handled");
878 break;
879 };
880
882}
883
885 EntData &data) {
889
890 auto N_InLoop = getNinTheLoop();
891 auto sensee = getSkeletonSense();
892 auto nb_gauss_pts = getGaussPts().size2();
893 auto t_normal = getFTensor1NormalsAtGaussPts();
894
895 auto t_sigma_ptr =
897 if (N_InLoop == 0) {
898 dataAtPts->tractionAtPts.resize(SPACE_DIM, nb_gauss_pts, false);
899 dataAtPts->tractionAtPts.clear();
900 }
901 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
902
903 for (int gg = 0; gg != nb_gauss_pts; gg++) {
904 t_traction(i) = t_sigma_ptr(i, j) * sensee * (t_normal(j) / t_normal.l2());
905 ++t_traction;
906 ++t_sigma_ptr;
907 ++t_normal;
908 }
909
911}
912
914 EntData &data) {
916 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
918 };
922 int nb_integration_pts = getGaussPts().size2();
923 auto t_w = getFTensor0IntegrationWeight();
924 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
925 auto t_coords = getFTensor1CoordsAtGaussPts();
926 auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
927
928 FTensor::Tensor1<double, 3> t_coords_spatial{0., 0., 0.};
929 // Offset for center of mass. Can be added in the future.
930 FTensor::Tensor1<double, 3> t_off{0.0, 0.0, 0.0};
931 FTensor::Tensor1<double, 3> loc_reaction_forces{0., 0., 0.};
932 FTensor::Tensor1<double, 3> loc_moment_forces{0., 0., 0.};
933
934 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
935 loc_reaction_forces(i) += (t_traction(i)) * t_w * getMeasure();
936 t_coords_spatial(i) = t_coords(i) + t_spatial_disp(i);
937 // t_coords_spatial(i) -= t_off(i);
938 loc_moment_forces(i) +=
939 (FTensor::levi_civita<double>(i, j, k) * t_coords_spatial(j)) *
940 t_traction(k) * t_w * getMeasure();
941 ++t_coords;
942 ++t_spatial_disp;
943 ++t_w;
944 ++t_traction;
945 }
946
947 reactionVec[0] += loc_reaction_forces(0);
948 reactionVec[1] += loc_reaction_forces(1);
949 reactionVec[2] += loc_reaction_forces(2);
950 reactionVec[3] += loc_moment_forces(0);
951 reactionVec[4] += loc_moment_forces(1);
952 reactionVec[5] += loc_moment_forces(2);
953
955}
956
959 int nb_dofs = data.getIndices().size();
960 int nb_integration_pts = data.getN().size1();
961 auto v = getVolume();
962 auto t_w = getFTensor0IntegrationWeight();
963 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
964 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
965 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
966 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
967 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
968 dataAtPts->wL2DotDotAtPts.clear();
969 }
970 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
971
972 auto piola_scale = dataAtPts->piolaScale;
973 auto alpha_w = alphaW / piola_scale;
974 auto alpha_rho = alphaRho / piola_scale;
975
976 int nb_base_functions = data.getN().size2();
977 auto t_row_base_fun = data.getFTensor0N();
978
979 FTensor::Index<'i', 3> i;
980 auto get_ftensor1 = [](auto &v) {
982 &v[2]);
983 };
984
985 for (int gg = 0; gg != nb_integration_pts; ++gg) {
986 double a = v * t_w;
987 auto t_nf = get_ftensor1(nF);
988 int bb = 0;
989 for (; bb != nb_dofs / 3; ++bb) {
990 t_nf(i) -= a * t_row_base_fun * t_div_P(i);
991 t_nf(i) += a * t_row_base_fun * alpha_w * t_s_dot_w(i);
992 t_nf(i) += a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(i);
993 ++t_nf;
994 ++t_row_base_fun;
995 }
996 for (; bb != nb_base_functions; ++bb)
997 ++t_row_base_fun;
998 ++t_w;
999 ++t_div_P;
1000 ++t_s_dot_w;
1001 ++t_s_dot_dot_w;
1002 }
1003
1005}
1006
1009 int nb_dofs = data.getIndices().size();
1010 int nb_integration_pts = getGaussPts().size2();
1011 auto v = getVolume();
1012 auto t_w = getFTensor0IntegrationWeight();
1013 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
1014 auto t_omega_grad_dot =
1015 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
1016 int nb_base_functions = data.getN().size2();
1017 auto t_row_base_fun = data.getFTensor0N();
1018 auto t_row_grad_fun = data.getFTensor1DiffN<3>();
1019 FTensor::Index<'i', 3> i;
1020 FTensor::Index<'j', 3> j;
1021 FTensor::Index<'k', 3> k;
1022 auto get_ftensor1 = [](auto &v) {
1024 &v[2]);
1025 };
1026 auto time_step = getTStimeStep();
1027
1028 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1029
1030 double a = v * t_w;
1031 auto t_nf = get_ftensor1(nF);
1032 int bb = 0;
1033 for (; bb != nb_dofs / 3; ++bb) {
1034 t_nf(k) -= (a * t_row_base_fun) * t_levi_kirchhoff(k);
1035 t_nf(k) += (a * alphaOmega /*/ time_step*/) *
1036 (t_row_grad_fun(i) * t_omega_grad_dot(k, i));
1037 ++t_nf;
1038 ++t_row_base_fun;
1039 ++t_row_grad_fun;
1040 }
1041 for (; bb != nb_base_functions; ++bb) {
1042 ++t_row_base_fun;
1043 ++t_row_grad_fun;
1044 }
1045 ++t_w;
1046 ++t_levi_kirchhoff;
1047 ++t_omega_grad_dot;
1048 }
1050}
1051
1054 int nb_dofs = data.getIndices().size();
1055 int nb_integration_pts = data.getN().size1();
1056 auto v = getVolume();
1057 auto t_w = getFTensor0IntegrationWeight();
1058
1059 int nb_base_functions = data.getN().size2() / 3;
1060 auto t_row_base_fun = data.getFTensor1N<3>();
1061 FTENSOR_INDEX(3, i);
1062 FTENSOR_INDEX(3, j);
1063 FTENSOR_INDEX(3, k);
1064 FTENSOR_INDEX(3, m);
1065 FTENSOR_INDEX(3, l);
1066
1067 auto get_ftensor1 = [](auto &v) {
1069 &v[2]);
1070 };
1071
1073
1074 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1075 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1076
1077 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1078 double a = v * t_w;
1079 auto t_nf = get_ftensor1(nF);
1080
1081 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1082
1083 int bb = 0;
1084 for (; bb != nb_dofs / 3; ++bb) {
1085 t_nf(i) -= a * (t_row_base_fun(k) * (t_R(i, l) * t_u(l, k)) / 2);
1086 t_nf(i) -= a * (t_row_base_fun(l) * (t_R(i, k) * t_u(l, k)) / 2);
1087 t_nf(i) += a * (t_row_base_fun(j) * t_kd(i, j));
1088 ++t_nf;
1089 ++t_row_base_fun;
1090 }
1091
1092 for (; bb != nb_base_functions; ++bb)
1093 ++t_row_base_fun;
1094
1095 ++t_w;
1096 ++t_R;
1097 ++t_u;
1098 }
1099
1100 } else {
1101
1102 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1103
1104 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1105 double a = v * t_w;
1106 auto t_nf = get_ftensor1(nF);
1107
1108 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1110
1111 t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
1112
1113 int bb = 0;
1114 for (; bb != nb_dofs / 3; ++bb) {
1115 t_nf(i) -= a * t_row_base_fun(j) * t_residuum(i, j);
1116 ++t_nf;
1117 ++t_row_base_fun;
1118 }
1119
1120 for (; bb != nb_base_functions; ++bb)
1121 ++t_row_base_fun;
1122
1123 ++t_w;
1124 ++t_h;
1125 }
1126 }
1127
1129}
1130
1133 int nb_dofs = data.getIndices().size();
1134 int nb_integration_pts = data.getN().size1();
1135 auto v = getVolume();
1136 auto t_w = getFTensor0IntegrationWeight();
1137
1138 int nb_base_functions = data.getN().size2() / 9;
1139 auto t_row_base_fun = data.getFTensor2N<3, 3>();
1140 FTENSOR_INDEX(3, i);
1141 FTENSOR_INDEX(3, j);
1142 FTENSOR_INDEX(3, k);
1143 FTENSOR_INDEX(3, m);
1144 FTENSOR_INDEX(3, l);
1145
1146 auto get_ftensor0 = [](auto &v) {
1148 };
1149
1151
1152 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1153 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1154
1155 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1156 double a = v * t_w;
1157 auto t_nf = get_ftensor0(nF);
1158
1159 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1160
1161 int bb = 0;
1162 for (; bb != nb_dofs; ++bb) {
1163 t_nf -= a * t_row_base_fun(i, k) * (t_R(i, l) * t_u(l, k)) / 2;
1164 t_nf -= a * t_row_base_fun(i, l) * (t_R(i, k) * t_u(l, k)) / 2;
1165 ++t_nf;
1166 ++t_row_base_fun;
1167 }
1168 for (; bb != nb_base_functions; ++bb) {
1169 ++t_row_base_fun;
1170 }
1171 ++t_w;
1172 ++t_R;
1173 ++t_u;
1174 }
1175
1176 } else {
1177
1178 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1179
1180 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1181 double a = v * t_w;
1182 auto t_nf = get_ftensor0(nF);
1183
1184 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1186 t_residuum(i, j) = t_h(i, j);
1187
1188 int bb = 0;
1189 for (; bb != nb_dofs; ++bb) {
1190 t_nf -= a * t_row_base_fun(i, j) * t_residuum(i, j);
1191 ++t_nf;
1192 ++t_row_base_fun;
1193 }
1194 for (; bb != nb_base_functions; ++bb) {
1195 ++t_row_base_fun;
1196 }
1197 ++t_w;
1198 ++t_h;
1199 }
1200 }
1201
1203}
1204
1207 int nb_dofs = data.getIndices().size();
1208 int nb_integration_pts = data.getN().size1();
1209 auto v = getVolume();
1210 auto t_w = getFTensor0IntegrationWeight();
1211 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1212 int nb_base_functions = data.getN().size2() / 3;
1213 auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
1214 FTensor::Index<'i', 3> i;
1215 auto get_ftensor1 = [](auto &v) {
1217 &v[2]);
1218 };
1219
1220 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1221 double a = v * t_w;
1222 auto t_nf = get_ftensor1(nF);
1223 int bb = 0;
1224 for (; bb != nb_dofs / 3; ++bb) {
1225 double div_row_base = t_row_diff_base_fun(i, i);
1226 t_nf(i) -= a * div_row_base * t_w_l2(i);
1227 ++t_nf;
1228 ++t_row_diff_base_fun;
1229 }
1230 for (; bb != nb_base_functions; ++bb) {
1231 ++t_row_diff_base_fun;
1232 }
1233 ++t_w;
1234 ++t_w_l2;
1235 }
1236
1238}
1239
1240template <>
1242 EntData &data) {
1244
1245 int nb_integration_pts = getGaussPts().size2();
1246
1247 Tag tag;
1248 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1249 int tag_length;
1250 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1251 if (tag_length != 9) {
1252 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1253 "Number of internal stress components should be 9 but is %d",
1254 tag_length);
1255 }
1256
1257 VectorDouble const_stress_vec(9);
1258 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1259 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1260 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1261 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1262
1263 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts, false);
1264 dataAtPts->internalStressAtPts.clear();
1265 auto t_internal_stress =
1266 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1267
1268 FTensor::Index<'L', 9> L;
1269 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1270 t_internal_stress(L) = t_const_stress(L);
1271 ++t_internal_stress;
1272 }
1273
1275}
1276
1277template <>
1279 EntityType type,
1280 EntData &data) {
1282
1283 int nb_integration_pts = getGaussPts().size2();
1284
1285 Tag tag;
1286 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1287 int tag_length;
1288 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1289 if (tag_length != 9) {
1290 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1291 "Number of internal stress components should be 9 but is %d",
1292 tag_length);
1293 }
1294
1295 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1296 const EntityHandle *vert_conn;
1297 int vert_num;
1298 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1299 vert_num, true);
1300 VectorDouble vert_data(vert_num * tag_length);
1301 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1302 &vert_data[0]);
1303
1304 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts, false);
1305 dataAtPts->internalStressAtPts.clear();
1306 auto t_internal_stress =
1307 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1308
1309 auto t_shape_n = data.getFTensor0N();
1310 int nb_shape_fn = data.getN(NOBASE).size2();
1311 FTensor::Index<'L', 9> L;
1312 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1313 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1314 for (int bb = 0; bb != nb_shape_fn; ++bb) {
1315 t_internal_stress(L) += t_vert_data(L) * t_shape_n;
1316 ++t_vert_data;
1317 ++t_shape_n;
1318 }
1319 ++t_internal_stress;
1320 }
1321
1323}
1324
1325template <>
1328
1329 int nb_dofs = data.getIndices().size();
1330 int nb_integration_pts = data.getN().size1();
1331 auto v = getVolume();
1332 auto t_w = getFTensor0IntegrationWeight();
1333
1334 FTensor::Index<'i', 3> i;
1335 FTensor::Index<'j', 3> j;
1336
1337 auto get_ftensor2 = [](auto &v) {
1339 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
1340 };
1341
1342 auto t_internal_stress =
1343 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1344
1346 auto t_L = symm_L_tensor();
1347
1348 int nb_base_functions = data.getN().size2();
1349 auto t_row_base_fun = data.getFTensor0N();
1350 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1351 double a = v * t_w;
1352 auto t_nf = get_ftensor2(nF);
1353
1354 FTensor::Tensor2<double, 3, 3> t_symm_stress;
1355 t_symm_stress(i, j) =
1356 (t_internal_stress(i, j) + t_internal_stress(j, i)) / 2;
1357
1359 t_residual(L) = t_L(i, j, L) * t_symm_stress(i, j);
1360
1361 int bb = 0;
1362 for (; bb != nb_dofs / 6; ++bb) {
1363 t_nf(L) += a * t_row_base_fun * t_residual(L);
1364 ++t_nf;
1365 ++t_row_base_fun;
1366 }
1367 for (; bb != nb_base_functions; ++bb)
1368 ++t_row_base_fun;
1369
1370 ++t_w;
1371 ++t_internal_stress;
1372 }
1374}
1375
1376template <>
1379
1380 int nb_dofs = data.getIndices().size();
1381 int nb_integration_pts = data.getN().size1();
1382 auto v = getVolume();
1383 auto t_w = getFTensor0IntegrationWeight();
1384
1385 auto get_ftensor2 = [](auto &v) {
1387 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
1388 };
1389
1390 auto t_internal_stress =
1391 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1392
1394 FTensor::Index<'M', size_symm> M;
1396 t_L = voigt_to_symm();
1397
1398 int nb_base_functions = data.getN().size2();
1399 auto t_row_base_fun = data.getFTensor0N();
1400 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1401 double a = v * t_w;
1402 auto t_nf = get_ftensor2(nF);
1403
1405 t_residual(L) = t_L(M, L) * t_internal_stress(M);
1406
1407 int bb = 0;
1408 for (; bb != nb_dofs / 6; ++bb) {
1409 t_nf(L) += a * t_row_base_fun * t_residual(L);
1410 ++t_nf;
1411 ++t_row_base_fun;
1412 }
1413 for (; bb != nb_base_functions; ++bb)
1414 ++t_row_base_fun;
1415
1416 ++t_w;
1417 ++t_internal_stress;
1418 }
1420}
1421
1422template <AssemblyType A>
1425 // get entity of face
1426 EntityHandle fe_ent = OP::getFEEntityHandle();
1427 // iterate over all boundary data
1428 for (auto &bc : (*bcDispPtr)) {
1429 // check if finite element entity is part of boundary condition
1430 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1431 int nb_dofs = data.getIndices().size();
1432
1433 int nb_integration_pts = OP::getGaussPts().size2();
1434 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1435 auto t_w = OP::getFTensor0IntegrationWeight();
1436 int nb_base_functions = data.getN().size2() / 3;
1437 auto t_row_base_fun = data.getFTensor1N<3>();
1438
1441
1442 double scale = 1;
1443 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1445 scale *= scalingMethodsMap.at(bc.blockName)
1446 ->getScale(EshelbianCore::dynamicTime);
1447 } else {
1448 scale *= scalingMethodsMap.at(bc.blockName)
1449 ->getScale(OP::getFEMethod()->ts_t);
1450 }
1451 } else {
1452 MOFEM_LOG("SELF", Sev::warning)
1453 << "No scaling method found for " << bc.blockName;
1454 }
1455
1456 // get bc data
1457 FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
1458 t_bc_disp(i) *= scale;
1459
1460 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1461 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1462 int bb = 0;
1463 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
1464 t_nf(i) +=
1465 t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
1466 ++t_nf;
1467 ++t_row_base_fun;
1468 }
1469 for (; bb != nb_base_functions; ++bb)
1470 ++t_row_base_fun;
1471
1472 ++t_w;
1473 ++t_normal;
1474 }
1475 }
1476 }
1478}
1479
1481 return OP::iNtegrate(data);
1482}
1483
1484template <AssemblyType A>
1487
1488 FTENSOR_INDEX(3, i);
1489 FTENSOR_INDEX(3, j);
1490 FTENSOR_INDEX(3, k);
1491
1492 double time = OP::getFEMethod()->ts_t;
1495 }
1496
1497 // get entity of face
1498 EntityHandle fe_ent = OP::getFEEntityHandle();
1499 // interate over all boundary data
1500 for (auto &bc : (*bcRotPtr)) {
1501 // check if finite element entity is part of boundary condition
1502 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1503 int nb_dofs = data.getIndices().size();
1504 int nb_integration_pts = OP::getGaussPts().size2();
1505 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1506 auto t_w = OP::getFTensor0IntegrationWeight();
1507
1508 int nb_base_functions = data.getN().size2() / 3;
1509 auto t_row_base_fun = data.getFTensor1N<3>();
1510
1511 auto get_ftensor1 = [](auto &v) {
1513 &v[2]);
1514 };
1515
1516 // Note: First three values of bc.vals are the center of rotation
1517 // 4th is rotation angle in radians, and remaining values are axis of
1518 // rotation. Also, if rotation axis is not provided, it defaults to the
1519 // normal vector of the face.
1520
1521 // get bc data
1522 FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
1523
1524 auto get_rotation_angle = [&]() {
1525 double theta = bc.theta;
1526 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1527 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1528 }
1529 return theta;
1530 };
1531
1532 auto get_rotation = [&](auto theta) {
1534 if (bc.vals.size() == 7) {
1535 t_omega(0) = bc.vals[4];
1536 t_omega(1) = bc.vals[5];
1537 t_omega(2) = bc.vals[6];
1538 } else {
1539 // Use gemetric face normal as rotation axis
1540 t_omega(i) = OP::getFTensor1Normal()(i);
1541 }
1542 t_omega.normalize();
1543 t_omega(i) *= theta;
1546 ? 0.
1547 : t_omega.l2());
1548 };
1549
1550 auto t_R = get_rotation(get_rotation_angle());
1551 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1552
1553 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1555 t_delta(i) = t_center(i) - t_coords(i);
1557 t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
1558
1559 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1560 int bb = 0;
1561 for (; bb != nb_dofs / 3; ++bb) {
1562 t_nf(i) += t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
1563 ++t_nf;
1564 ++t_row_base_fun;
1565 }
1566 for (; bb != nb_base_functions; ++bb)
1567 ++t_row_base_fun;
1568
1569 ++t_w;
1570 ++t_normal;
1571 ++t_coords;
1572 }
1573 }
1574 }
1576}
1577
1579 return OP::iNtegrate(data);
1580}
1581
1582template <AssemblyType A>
1585
1586 double time = OP::getFEMethod()->ts_t;
1589 }
1590
1591 // get entity of face
1592 EntityHandle fe_ent = OP::getFEEntityHandle();
1593 // iterate over all boundary data
1594 for (auto &bc : (*bcDispPtr)) {
1595 // check if finite element entity is part of boundary condition
1596 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1597
1598 auto t_approx_P = getFTensor2FromMat<3, 3>(*piolaStressPtr);
1599 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1600 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1601 auto t_w = OP::getFTensor0IntegrationWeight();
1602
1605
1607
1608 double scale = 1;
1609 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1610 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1611 } else {
1612 MOFEM_LOG("SELF", Sev::warning)
1613 << "No scaling method found for " << bc.blockName;
1614 }
1615
1616 // get bc data
1617 double val = scale * bc.val;
1618
1619 int nb_dofs = data.getIndices().size();
1620 int nb_integration_pts = OP::getGaussPts().size2();
1621 int nb_base_functions = data.getN().size2();
1622 auto t_row_base = data.getFTensor0N();
1623 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1624
1626 t_N(i) = t_normal(i);
1627 t_N.normalize();
1628
1630 t_P(i, j) = t_N(i) * t_N(j);
1632 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
1633
1635 t_tracton(i) = t_approx_P(i, j) * t_N(j);
1636
1638 t_res(i) = t_Q(i, j) * t_tracton(j) + t_P(i, j) * 2* t_u(j) - t_N(i) * val;
1639
1640 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1641 int bb = 0;
1642 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
1643 t_nf(i) += (t_w * t_row_base) * t_res(i);
1644 ++t_nf;
1645 ++t_row_base;
1646 }
1647 for (; bb != nb_base_functions; ++bb)
1648 ++t_row_base;
1649
1650 ++t_w;
1651 ++t_normal;
1652 ++t_u;
1653 ++t_approx_P;
1654 }
1655
1656 OP::locF *= OP::getMeasure();
1657 }
1658 }
1660}
1661
1662template <AssemblyType A>
1665 EntData &col_data) {
1667
1668 double time = OP::getFEMethod()->ts_t;
1671 }
1672
1673 int row_nb_dofs = row_data.getIndices().size();
1674 int col_nb_dofs = col_data.getIndices().size();
1675 auto &locMat = OP::locMat;
1676 locMat.resize(row_nb_dofs, col_nb_dofs, false);
1677 locMat.clear();
1678
1679 // get entity of face
1680 EntityHandle fe_ent = OP::getFEEntityHandle();
1681 // iterate over all boundary data
1682 for (auto &bc : (*bcDispPtr)) {
1683 // check if finite element entity is part of boundary condition
1684 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1685
1686 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1687 auto t_w = OP::getFTensor0IntegrationWeight();
1688
1691
1692 double scale = 1;
1693 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1694 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1695 } else {
1696 MOFEM_LOG("SELF", Sev::warning)
1697 << "No scaling method found for " << bc.blockName;
1698 }
1699
1700 int nb_integration_pts = OP::getGaussPts().size2();
1701 int row_nb_dofs = row_data.getIndices().size();
1702 int col_nb_dofs = col_data.getIndices().size();
1703 int nb_base_functions = row_data.getN().size2();
1704 auto t_row_base = row_data.getFTensor0N();
1705
1707
1708 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1709
1711 t_N(i) = t_normal(i);
1712 t_N.normalize();
1713
1715 t_P(i, j) = t_N(i) * t_N(j);
1716
1718 t_d_res(i, j) = 2.0 * t_P(i, j);
1719
1720 int rr = 0;
1721 for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
1723 locMat, SPACE_DIM * rr);
1724 auto t_col_base = col_data.getFTensor0N(gg, 0);
1725 for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
1726 t_mat(i, j) += (t_w * t_row_base * t_col_base) * t_d_res(i, j);
1727 ++t_mat;
1728 ++t_col_base;
1729 }
1730 ++t_row_base;
1731 }
1732
1733 for (; rr != nb_base_functions; ++rr)
1734 ++t_row_base;
1735
1736 ++t_w;
1737 ++t_normal;
1738 }
1739
1740 locMat *= OP::getMeasure();
1741
1742 }
1743 }
1745}
1746
1747template <AssemblyType A>
1750 EntData &col_data) {
1752
1753 double time = OP::getFEMethod()->ts_t;
1756 }
1757
1758 int row_nb_dofs = row_data.getIndices().size();
1759 int col_nb_dofs = col_data.getIndices().size();
1760 auto &locMat = OP::locMat;
1761 locMat.resize(row_nb_dofs, col_nb_dofs, false);
1762 locMat.clear();
1763
1764 // get entity of face
1765 EntityHandle fe_ent = OP::getFEEntityHandle();
1766 // iterate over all boundary data
1767 for (auto &bc : (*bcDispPtr)) {
1768 // check if finite element entity is part of boundary condition
1769 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1770
1771 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1772 auto t_w = OP::getFTensor0IntegrationWeight();
1773
1777
1779
1780 double scale = 1;
1781 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1782 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1783 } else {
1784 MOFEM_LOG("SELF", Sev::warning)
1785 << "No scaling method found for " << bc.blockName;
1786 }
1787
1788 int nb_integration_pts = OP::getGaussPts().size2();
1789 int nb_base_functions = row_data.getN().size2();
1790 auto t_row_base = row_data.getFTensor0N();
1791
1792 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1793
1795 t_N(i) = t_normal(i);
1796 t_N.normalize();
1797
1799 t_P(i, j) = t_N(i) * t_N(j);
1801 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
1802
1804 t_d_res(i, j) = t_Q(i, j);
1805
1806 int rr = 0;
1807 for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
1809 OP::locMat, SPACE_DIM * rr);
1810 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1811 for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
1812 t_mat(i, j) +=
1813 ((t_w * t_row_base) * (t_N(k) * t_col_base(k))) * t_d_res(i, j);
1814 ++t_mat;
1815 ++t_col_base;
1816 }
1817 ++t_row_base;
1818 }
1819
1820 for (; rr != nb_base_functions; ++rr)
1821 ++t_row_base;
1822
1823 ++t_w;
1824 ++t_normal;
1825 }
1826
1827 locMat *= OP::getMeasure();
1828 }
1829 }
1831}
1832
1834 return OP::iNtegrate(data);
1835}
1836
1838 EntData &col_data) {
1839 return OP::iNtegrate(row_data, col_data);
1840}
1841
1843 EntData &col_data) {
1844 return OP::iNtegrate(row_data, col_data);
1845}
1846
1847template <AssemblyType A>
1850
1851 double time = OP::getFEMethod()->ts_t;
1854 }
1855
1856 // get entity of face
1857 EntityHandle fe_ent = OP::getFEEntityHandle();
1858 // iterate over all boundary data
1859 for (auto &bc : (*bcDispPtr)) {
1860 // check if finite element entity is part of boundary condition
1861 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1862
1863 MatrixDouble analytical_expr;
1864 // placeholder to pass boundary block id to python
1865
1866 auto [block_name, v_analytical_expr] =
1867 getAnalyticalExpr(this, analytical_expr, bc.blockName);
1868
1869 int nb_dofs = data.getIndices().size();
1870
1871 int nb_integration_pts = OP::getGaussPts().size2();
1872 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1873 auto t_w = OP::getFTensor0IntegrationWeight();
1874 int nb_base_functions = data.getN().size2() / 3;
1875 auto t_row_base_fun = data.getFTensor1N<3>();
1876 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1877
1880
1881 // get bc data
1882 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1883
1884 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1885 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1886
1887 int bb = 0;
1888 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
1889 t_nf(i) +=
1890 t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
1891 ++t_nf;
1892 ++t_row_base_fun;
1893 }
1894 for (; bb != nb_base_functions; ++bb)
1895 ++t_row_base_fun;
1896
1897 ++t_bc_disp;
1898 ++t_coord;
1899 ++t_w;
1900 ++t_normal;
1901 }
1902 }
1903 }
1905}
1906
1908 return OP::iNtegrate(data);
1909}
1910
1913
1914 FTENSOR_INDEX(3, i);
1915
1916 int nb_dofs = data.getFieldData().size();
1917 int nb_integration_pts = getGaussPts().size2();
1918 int nb_base_functions = data.getN().size2();
1919
1920 double time = getFEMethod()->ts_t;
1923 }
1924
1925#ifndef NDEBUG
1926 if (this->locF.size() != nb_dofs)
1927 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1928 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1929#endif // NDEBUG
1930
1931 auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
1933
1934 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1935 auto t_row_base = data.getFTensor0N();
1936 auto t_w = getFTensor0IntegrationWeight();
1937 auto t_coords = getFTensor1CoordsAtGaussPts();
1938
1939 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1940
1941 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1942
1943 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1944 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1945 int rr = 0;
1946 for (; rr != nb_dofs / SPACE_DIM; ++rr) {
1947 t_f(i) -= (time_scale * t_w * t_row_base * tau) * (t_val(i) * scale);
1948 ++t_row_base;
1949 ++t_f;
1950 }
1951
1952 for (; rr != nb_base_functions; ++rr)
1953 ++t_row_base;
1954 ++t_w;
1955 ++t_coords;
1956 }
1957 this->locF *= getMeasure();
1959 };
1960
1961 // get entity of face
1962 EntityHandle fe_ent = getFEEntityHandle();
1963 for (auto &bc : *(bcData)) {
1964 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1965
1966 double time_scale = 1;
1967 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1968 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1969 }
1970
1971 int nb_dofs = data.getFieldData().size();
1972 if (nb_dofs) {
1973
1974 if (std::regex_match(bc.blockName, std::regex(".*COOK.*"))) {
1975 auto calc_tau = [](double, double y, double) {
1976 y -= 44;
1977 y /= (60 - 44);
1978 return -y * (y - 1) / 0.25;
1979 };
1980 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1981 } else {
1982 CHKERR integrate_rhs(
1983 bc, [](double, double, double) { return 1; }, time_scale);
1984 }
1985 }
1986 }
1987 }
1989}
1990
1993
1994 FTENSOR_INDEX(3, i);
1995
1996 int nb_dofs = data.getFieldData().size();
1997 int nb_integration_pts = getGaussPts().size2();
1998 int nb_base_functions = data.getN().size2();
1999
2000 double time = getFEMethod()->ts_t;
2003 }
2004
2005#ifndef NDEBUG
2006 if (this->locF.size() != nb_dofs)
2007 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2008 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2009#endif // NDEBUG
2010
2011 auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
2013
2014 auto val = bc.val;
2015 auto t_row_base = data.getFTensor0N();
2016 auto t_w = getFTensor0IntegrationWeight();
2017 auto t_coords = getFTensor1CoordsAtGaussPts();
2018 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2019 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2020
2021 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2022
2023 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2024
2025 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2026
2032
2036
2037 t_normal(i) = (FTensor::levi_civita<double>(i, j, k) * t_tangent1(j)) *
2038 t_tangent2(k);
2039 } else {
2040 t_normal(i) = (FTensor::levi_civita<double>(i, j, k) *
2041 (t_tangent1(j) + t_grad_gamma_u(j, N0))) *
2042 (t_tangent2(k) + t_grad_gamma_u(k, N1));
2043 }
2044 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2045 auto t_val = FTensor::Tensor1<double, 3>();
2046 t_val(i) = (time_scale * t_w * tau * scale * val) * t_normal(i);
2047
2048 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2049 int rr = 0;
2050 for (; rr != nb_dofs / SPACE_DIM; ++rr) {
2051 t_f(i) += t_row_base * t_val(i);
2052 ++t_row_base;
2053 ++t_f;
2054 }
2055
2056 for (; rr != nb_base_functions; ++rr)
2057 ++t_row_base;
2058 ++t_w;
2059 ++t_coords;
2060 ++t_tangent1;
2061 ++t_tangent2;
2062 ++t_grad_gamma_u;
2063 }
2064 this->locF /= 2.;
2065
2067 };
2068
2069 // get entity of face
2070 EntityHandle fe_ent = getFEEntityHandle();
2071 for (auto &bc : *(bcData)) {
2072 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2073
2074 double time_scale = 1;
2075 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2076 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2077 }
2078
2079 int nb_dofs = data.getFieldData().size();
2080 if (nb_dofs) {
2081 CHKERR integrate_rhs(
2082 bc, [](double, double, double) { return 1; }, time_scale);
2083 }
2084 }
2085 }
2087}
2088
2089template <AssemblyType A>
2092 EntData &col_data) {
2094
2098 }
2099
2100 double time = OP::getFEMethod()->ts_t;
2103 }
2104
2105 int nb_base_functions = row_data.getN().size2();
2106 int row_nb_dofs = row_data.getIndices().size();
2107 int col_nb_dofs = col_data.getIndices().size();
2108 int nb_integration_pts = OP::getGaussPts().size2();
2109 auto &locMat = OP::locMat;
2110 locMat.resize(row_nb_dofs, col_nb_dofs, false);
2111 locMat.clear();
2112
2113auto integrate_lhs = [&](auto &bc, auto calc_tau, double time_scale) {
2115
2116 auto val = bc.val;
2117 auto t_row_base = row_data.getFTensor0N();
2118 auto t_w = OP::getFTensor0IntegrationWeight();
2119 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2120 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2121 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2122
2123 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2125
2126 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2127
2132
2135
2136 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2137 auto t_val = time_scale * t_w * tau * val;
2138
2139 int rr = 0;
2140 for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
2142 locMat, SPACE_DIM * rr);
2143 auto t_diff_col_base = col_data.getFTensor1DiffN<2>(gg, 0);
2144 for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
2146 t_normal_du(i, l) = (FTensor::levi_civita<double>(i, j, k) *
2147 (t_tangent2(k) + t_grad_gamma_u(k, N1))) *
2148 t_kd(j, l) * t_diff_col_base(N0)
2149
2150 +
2151
2153 (t_tangent1(j) + t_grad_gamma_u(j, N0))) *
2154 t_kd(k, l) * t_diff_col_base(N1);
2155
2156 t_mat(i, j) += (t_w * t_row_base) * t_val * t_normal_du(i, j);
2157 ++t_mat;
2158 ++t_diff_col_base;
2159 }
2160 ++t_row_base;
2161 }
2162
2163 for (; rr != nb_base_functions; ++rr)
2164 ++t_row_base;
2165 ++t_w;
2166 ++t_coords;
2167 ++t_tangent1;
2168 ++t_tangent2;
2169 ++t_grad_gamma_u;
2170 }
2171
2172 OP::locMat /= 2.;
2173
2175 };
2176
2177 // get entity of face
2178 EntityHandle fe_ent = OP::getFEEntityHandle();
2179 for (auto &bc : *(bcData)) {
2180 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2181
2182 double time_scale = 1;
2183 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2184 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2185 }
2186
2187 int nb_dofs = row_data.getFieldData().size();
2188 if (nb_dofs) {
2189 CHKERR integrate_lhs(
2190 bc, [](double, double, double) { return 1; }, time_scale);
2191 }
2192 }
2193 }
2194
2196
2197}
2198
2200 EntData &col_data) {
2201 return OP::iNtegrate(row_data, col_data);
2202}
2203
2204
2207
2208 FTENSOR_INDEX(3, i);
2209
2210 int nb_dofs = data.getFieldData().size();
2211 int nb_integration_pts = getGaussPts().size2();
2212 int nb_base_functions = data.getN().size2();
2213
2214 double time = getFEMethod()->ts_t;
2217 }
2218
2219#ifndef NDEBUG
2220 if (this->locF.size() != nb_dofs)
2221 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2222 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2223#endif // NDEBUG
2224
2225 // get entity of face
2226 EntityHandle fe_ent = getFEEntityHandle();
2227 for (auto &bc : *(bcData)) {
2228 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2229
2230 MatrixDouble analytical_expr;
2231 // placeholder to pass boundary block id to python
2232 auto [block_name, v_analytical_expr] =
2233 getAnalyticalExpr(this, analytical_expr, bc.blockName);
2234 auto t_val =
2235 getFTensor1FromMat<3>(v_analytical_expr);
2236 auto t_row_base = data.getFTensor0N();
2237 auto t_w = getFTensor0IntegrationWeight();
2238 auto t_coords = getFTensor1CoordsAtGaussPts();
2239
2240 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2241
2242 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2243
2244 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2245 int rr = 0;
2246 for (; rr != nb_dofs / SPACE_DIM; ++rr) {
2247 t_f(i) -= t_w * t_row_base * (t_val(i) * scale);
2248 ++t_row_base;
2249 ++t_f;
2250 }
2251
2252 for (; rr != nb_base_functions; ++rr)
2253 ++t_row_base;
2254 ++t_w;
2255 ++t_coords;
2256 ++t_val;
2257 }
2258 this->locF *= getMeasure();
2259 }
2260 }
2262}
2263
2265 EntData &col_data) {
2267 int nb_integration_pts = row_data.getN().size1();
2268 int row_nb_dofs = row_data.getIndices().size();
2269 int col_nb_dofs = col_data.getIndices().size();
2270 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2272 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
2273 };
2274 FTensor::Index<'i', 3> i;
2275 auto v = getVolume();
2276 auto t_w = getFTensor0IntegrationWeight();
2277 int row_nb_base_functions = row_data.getN().size2();
2278 auto t_row_base_fun = row_data.getFTensor0N();
2279 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2280 double a = v * t_w;
2281 int rr = 0;
2282 for (; rr != row_nb_dofs / 3; ++rr) {
2283 auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
2284 auto t_m = get_ftensor1(K, 3 * rr, 0);
2285 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2286 double div_col_base = t_col_diff_base_fun(i, i);
2287 t_m(i) -= a * t_row_base_fun * div_col_base;
2288 ++t_m;
2289 ++t_col_diff_base_fun;
2290 }
2291 ++t_row_base_fun;
2292 }
2293 for (; rr != row_nb_base_functions; ++rr)
2294 ++t_row_base_fun;
2295 ++t_w;
2296 }
2298}
2299
2301 EntData &col_data) {
2303
2304 if (alphaW < std::numeric_limits<double>::epsilon() &&
2305 alphaRho < std::numeric_limits<double>::epsilon())
2307
2308 const int nb_integration_pts = row_data.getN().size1();
2309 const int row_nb_dofs = row_data.getIndices().size();
2310 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2312 &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
2313
2314 );
2315 };
2316 FTensor::Index<'i', 3> i;
2317
2318 auto v = getVolume();
2319 auto t_w = getFTensor0IntegrationWeight();
2320
2321 auto piola_scale = dataAtPts->piolaScale;
2322 auto alpha_w = alphaW / piola_scale;
2323 auto alpha_rho = alphaRho / piola_scale;
2324
2325 int row_nb_base_functions = row_data.getN().size2();
2326 auto t_row_base_fun = row_data.getFTensor0N();
2327
2328 double ts_scale = alpha_w * getTSa();
2329 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2330 ts_scale += alpha_rho * getTSaa();
2331
2332 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2333 double a = v * t_w * ts_scale;
2334
2335 int rr = 0;
2336 for (; rr != row_nb_dofs / 3; ++rr) {
2337
2338 auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
2339 auto t_m = get_ftensor1(K, 3 * rr, 0);
2340 for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2341 const double b = a * t_row_base_fun * t_col_base_fun;
2342 t_m(i) += b;
2343 ++t_m;
2344 ++t_col_base_fun;
2345 }
2346
2347 ++t_row_base_fun;
2348 }
2349
2350 for (; rr != row_nb_base_functions; ++rr)
2351 ++t_row_base_fun;
2352
2353 ++t_w;
2354 }
2355
2357}
2358
2360 EntData &col_data) {
2362
2368
2369 int nb_integration_pts = row_data.getN().size1();
2370 int row_nb_dofs = row_data.getIndices().size();
2371 int col_nb_dofs = col_data.getIndices().size();
2372 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
2374
2375 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2376
2377 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2378
2379 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
2380
2381 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
2382
2383 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
2384
2385 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
2386 };
2387
2388 auto v = getVolume();
2389 auto t_w = getFTensor0IntegrationWeight();
2390
2391 int row_nb_base_functions = row_data.getN().size2();
2392 auto t_row_base_fun = row_data.getFTensor0N();
2393
2394 auto t_approx_P_adjoint_log_du_dP =
2396
2397 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2398 double a = v * t_w;
2399 int rr = 0;
2400 for (; rr != row_nb_dofs / 6; ++rr) {
2401
2402 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
2403 auto t_m = get_ftensor3(K, 6 * rr, 0);
2404
2405 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2406 t_m(L, i) -=
2407 a * (t_approx_P_adjoint_log_du_dP(i, j, L) * t_col_base_fun(j)) *
2408 t_row_base_fun;
2409 ++t_col_base_fun;
2410 ++t_m;
2411 }
2412
2413 ++t_row_base_fun;
2414 }
2415 for (; rr != row_nb_base_functions; ++rr)
2416 ++t_row_base_fun;
2417 ++t_w;
2418 ++t_approx_P_adjoint_log_du_dP;
2419 }
2420
2422}
2423
2425 EntData &col_data) {
2427
2433
2434 int nb_integration_pts = row_data.getN().size1();
2435 int row_nb_dofs = row_data.getIndices().size();
2436 int col_nb_dofs = col_data.getIndices().size();
2437 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2439 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
2440 &m(r + 5, c));
2441 };
2442
2443 auto v = getVolume();
2444 auto t_w = getFTensor0IntegrationWeight();
2445 auto t_row_base_fun = row_data.getFTensor0N();
2446
2447 int row_nb_base_functions = row_data.getN().size2();
2448
2449 auto t_approx_P_adjoint_log_du_dP =
2451
2452 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2453 double a = v * t_w;
2454 int rr = 0;
2455 for (; rr != row_nb_dofs / 6; ++rr) {
2456 auto t_m = get_ftensor2(K, 6 * rr, 0);
2457 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
2458 for (int cc = 0; cc != col_nb_dofs; ++cc) {
2459 t_m(L) -=
2460 a * (t_approx_P_adjoint_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
2461 t_row_base_fun;
2462 ++t_m;
2463 ++t_col_base_fun;
2464 }
2465 ++t_row_base_fun;
2466 }
2467 for (; rr != row_nb_base_functions; ++rr)
2468 ++t_row_base_fun;
2469 ++t_w;
2470 ++t_approx_P_adjoint_log_du_dP;
2471 }
2473}
2474
2476 EntData &col_data) {
2478
2480 auto t_L = symm_L_tensor();
2481
2482 int row_nb_dofs = row_data.getIndices().size();
2483 int col_nb_dofs = col_data.getIndices().size();
2484 auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
2486
2487 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2488
2489 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2490
2491 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
2492
2493 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
2494
2495 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
2496
2497 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
2498
2499 );
2500 };
2501 FTensor::Index<'i', 3> i;
2502 FTensor::Index<'j', 3> j;
2503 FTensor::Index<'k', 3> k;
2504 FTensor::Index<'m', 3> m;
2505 FTensor::Index<'n', 3> n;
2506
2507 auto v = getVolume();
2508 auto t_w = getFTensor0IntegrationWeight();
2509 auto t_approx_P_adjoint_log_du_domega =
2510 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2511
2512 int row_nb_base_functions = row_data.getN().size2();
2513 auto t_row_base_fun = row_data.getFTensor0N();
2514
2515 int nb_integration_pts = row_data.getN().size1();
2516
2517 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2518 double a = v * t_w;
2519
2520 int rr = 0;
2521 for (; rr != row_nb_dofs / 6; ++rr) {
2522 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2523 auto t_m = get_ftensor3(K, 6 * rr, 0);
2524 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2525 double v = a * t_row_base_fun * t_col_base_fun;
2526 t_m(L, k) -= v * t_approx_P_adjoint_log_du_domega(k, L);
2527 ++t_m;
2528 ++t_col_base_fun;
2529 }
2530 ++t_row_base_fun;
2531 }
2532
2533 for (; rr != row_nb_base_functions; ++rr)
2534 ++t_row_base_fun;
2535
2536 ++t_w;
2537 ++t_approx_P_adjoint_log_du_domega;
2538 }
2539
2541}
2542
2544 EntData &col_data) {
2546 int nb_integration_pts = getGaussPts().size2();
2547 int row_nb_dofs = row_data.getIndices().size();
2548 int col_nb_dofs = col_data.getIndices().size();
2549 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2551 size_symm>{
2552
2553 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
2554 &m(r + 0, c + 4), &m(r + 0, c + 5),
2555
2556 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
2557 &m(r + 1, c + 4), &m(r + 1, c + 5),
2558
2559 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
2560 &m(r + 2, c + 4), &m(r + 2, c + 5)
2561
2562 };
2563 };
2564
2567
2568 auto v = getVolume();
2569 auto t_w = getFTensor0IntegrationWeight();
2570 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2571 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2572 int row_nb_base_functions = row_data.getN().size2();
2573 auto t_row_base_fun = row_data.getFTensor0N();
2574 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2575 double a = v * t_w;
2576 int rr = 0;
2577 for (; rr != row_nb_dofs / 3; ++rr) {
2578 auto t_m = get_ftensor2(K, 3 * rr, 0);
2579 const double b = a * t_row_base_fun;
2580 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2581 for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
2582 t_m(k, L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(k, L);
2583 ++t_m;
2584 ++t_col_base_fun;
2585 }
2586 ++t_row_base_fun;
2587 }
2588 for (; rr != row_nb_base_functions; ++rr) {
2589 ++t_row_base_fun;
2590 }
2591 ++t_w;
2592 ++t_levi_kirchhoff_du;
2593 }
2595}
2596
2598 EntData &col_data) {
2600
2607
2608 int nb_integration_pts = getGaussPts().size2();
2609 int row_nb_dofs = row_data.getIndices().size();
2610 int col_nb_dofs = col_data.getIndices().size();
2611 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2613
2614 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2615
2616 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2617
2618 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2619
2620 );
2621 };
2622
2623 auto v = getVolume();
2624 auto t_w = getFTensor0IntegrationWeight();
2625
2626 int row_nb_base_functions = row_data.getN().size2();
2627 auto t_row_base_fun = row_data.getFTensor0N();
2628 auto t_levi_kirchhoff_dP =
2629 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2630
2631 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2632 double a = v * t_w;
2633 int rr = 0;
2634 for (; rr != row_nb_dofs / 3; ++rr) {
2635 double b = a * t_row_base_fun;
2636 auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
2637 auto t_m = get_ftensor2(K, 3 * rr, 0);
2638 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2639 t_m(m, i) -= b * (t_levi_kirchhoff_dP(m, i, k) * t_col_base_fun(k));
2640 ++t_m;
2641 ++t_col_base_fun;
2642 }
2643 ++t_row_base_fun;
2644 }
2645 for (; rr != row_nb_base_functions; ++rr) {
2646 ++t_row_base_fun;
2647 }
2648
2649 ++t_w;
2650 ++t_levi_kirchhoff_dP;
2651 }
2653}
2654
2656 EntData &col_data) {
2658 int nb_integration_pts = getGaussPts().size2();
2659 int row_nb_dofs = row_data.getIndices().size();
2660 int col_nb_dofs = col_data.getIndices().size();
2661
2662 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2664 &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
2665 };
2666
2667 FTENSOR_INDEX(3, i);
2668 FTENSOR_INDEX(3, k);
2669 FTENSOR_INDEX(3, m);
2670
2671 auto v = getVolume();
2672 auto t_w = getFTensor0IntegrationWeight();
2673 auto t_levi_kirchoff_dP =
2674 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2675
2676 int row_nb_base_functions = row_data.getN().size2();
2677 auto t_row_base_fun = row_data.getFTensor0N();
2678
2679 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2680 double a = v * t_w;
2681 int rr = 0;
2682 for (; rr != row_nb_dofs / 3; ++rr) {
2683 double b = a * t_row_base_fun;
2684 auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
2685 auto t_m = get_ftensor1(K, 3 * rr, 0);
2686 for (int cc = 0; cc != col_nb_dofs; ++cc) {
2687 t_m(m) -= b * (t_levi_kirchoff_dP(m, i, k) * t_col_base_fun(i, k));
2688 ++t_m;
2689 ++t_col_base_fun;
2690 }
2691 ++t_row_base_fun;
2692 }
2693
2694 for (; rr != row_nb_base_functions; ++rr) {
2695 ++t_row_base_fun;
2696 }
2697 ++t_w;
2698 ++t_levi_kirchoff_dP;
2699 }
2701}
2702
2704 EntData &col_data) {
2706 int nb_integration_pts = getGaussPts().size2();
2707 int row_nb_dofs = row_data.getIndices().size();
2708 int col_nb_dofs = col_data.getIndices().size();
2709 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2711
2712 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2713
2714 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2715
2716 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2717
2718 );
2719 };
2720 FTensor::Index<'i', 3> i;
2721 FTensor::Index<'j', 3> j;
2722 FTensor::Index<'k', 3> k;
2723 FTensor::Index<'l', 3> l;
2724 FTensor::Index<'m', 3> m;
2725 FTensor::Index<'n', 3> n;
2726
2728
2729 auto v = getVolume();
2730 auto ts_a = getTSa();
2731 auto t_w = getFTensor0IntegrationWeight();
2732 auto t_levi_kirchhoff_domega =
2733 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2734 int row_nb_base_functions = row_data.getN().size2();
2735 auto t_row_base_fun = row_data.getFTensor0N();
2736 auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
2737
2738 auto time_step = getTStimeStep();
2739 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2740 double a = v * t_w;
2741 double c = (a * alphaOmega) * (ts_a /*/ time_step*/);
2742
2743 int rr = 0;
2744 for (; rr != row_nb_dofs / 3; ++rr) {
2745 auto t_m = get_ftensor2(K, 3 * rr, 0);
2746 const double b = a * t_row_base_fun;
2747 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2748 auto t_col_grad_fun = col_data.getFTensor1DiffN<3>(gg, 0);
2749 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2750 t_m(k, l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(k, l);
2751 t_m(k, l) += t_kd(k, l) * (c * (t_row_grad_fun(i) * t_col_grad_fun(i)));
2752 ++t_m;
2753 ++t_col_base_fun;
2754 ++t_col_grad_fun;
2755 }
2756 ++t_row_base_fun;
2757 ++t_row_grad_fun;
2758 }
2759 for (; rr != row_nb_base_functions; ++rr) {
2760 ++t_row_base_fun;
2761 ++t_row_grad_fun;
2762 }
2763 ++t_w;
2764 ++t_levi_kirchhoff_domega;
2765 }
2767}
2768
2770 EntData &col_data) {
2772 if (dataAtPts->matInvD.size1() == size_symm &&
2773 dataAtPts->matInvD.size2() == size_symm) {
2774 MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
2775 } else {
2777 }
2779};
2780
2781template <int S>
2783 EntData &col_data) {
2785
2786 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2788
2789 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2790
2791 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2792
2793 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2794
2795 );
2796 };
2797
2798 int nb_integration_pts = getGaussPts().size2();
2799 int row_nb_dofs = row_data.getIndices().size();
2800 int col_nb_dofs = col_data.getIndices().size();
2801
2802 auto v = getVolume();
2803 auto t_w = getFTensor0IntegrationWeight();
2804 int row_nb_base_functions = row_data.getN().size2() / 3;
2805
2810
2812 &*dataAtPts->matInvD.data().begin());
2813
2814 auto t_row_base = row_data.getFTensor1N<3>();
2815 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2816 double a = v * t_w;
2817
2818 int rr = 0;
2819 for (; rr != row_nb_dofs / 3; ++rr) {
2820 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
2821 auto t_m = get_ftensor2(K, 3 * rr, 0);
2822 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2823 t_m(i, k) -= a * t_row_base(j) * (t_inv_D(i, j, k, l) * t_col_base(l));
2824 ++t_m;
2825 ++t_col_base;
2826 }
2827
2828 ++t_row_base;
2829 }
2830
2831 for (; rr != row_nb_base_functions; ++rr)
2832 ++t_row_base;
2833
2834 ++t_w;
2835 ++t_inv_D;
2836 }
2838}
2839
2842 EntData &col_data) {
2844 if (dataAtPts->matInvD.size1() == size_symm &&
2845 dataAtPts->matInvD.size2() == size_symm) {
2846 MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
2847 } else {
2849 }
2851};
2852
2853template <int S>
2856 EntData &col_data) {
2858
2859 int nb_integration_pts = getGaussPts().size2();
2860 int row_nb_dofs = row_data.getIndices().size();
2861 int col_nb_dofs = col_data.getIndices().size();
2862
2863 auto v = getVolume();
2864 auto t_w = getFTensor0IntegrationWeight();
2865 int row_nb_base_functions = row_data.getN().size2() / 9;
2866
2871
2873 &*dataAtPts->matInvD.data().begin());
2874
2875 auto t_row_base = row_data.getFTensor2N<3, 3>();
2876 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2877 double a = v * t_w;
2878
2879 int rr = 0;
2880 for (; rr != row_nb_dofs; ++rr) {
2881 auto t_col_base = col_data.getFTensor2N<3, 3>(gg, 0);
2882 for (int cc = 0; cc != col_nb_dofs; ++cc) {
2883 K(rr, cc) -=
2884 a * (t_row_base(i, j) * (t_inv_D(i, j, k, l) * t_col_base(k, l)));
2885 ++t_col_base;
2886 }
2887
2888 ++t_row_base;
2889 }
2890
2891 for (; rr != row_nb_base_functions; ++rr)
2892 ++t_row_base;
2893 ++t_w;
2894 ++t_inv_D;
2895 }
2897}
2898
2900 EntData &col_data) {
2902 if (dataAtPts->matInvD.size1() == size_symm &&
2903 dataAtPts->matInvD.size2() == size_symm) {
2904 MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
2905 } else {
2907 }
2909};
2910
2911template <int S>
2914 EntData &col_data) {
2916
2917 auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2919
2920 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2)
2921
2922 );
2923 };
2924
2925 int nb_integration_pts = getGaussPts().size2();
2926 int row_nb_dofs = row_data.getIndices().size();
2927 int col_nb_dofs = col_data.getIndices().size();
2928
2929 auto v = getVolume();
2930 auto t_w = getFTensor0IntegrationWeight();
2931 int row_nb_base_functions = row_data.getN().size2() / 9;
2932
2939
2941 &*dataAtPts->matInvD.data().begin());
2942
2943 auto t_row_base = row_data.getFTensor2N<3, 3>();
2944 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2945 double a = v * t_w;
2946
2947 auto t_m = get_ftensor1(K, 0, 0);
2948
2949 int rr = 0;
2950 for (; rr != row_nb_dofs; ++rr) {
2951 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
2952 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2953 t_m(k) -= a * (t_row_base(i, j) * t_inv_D(i, j, k, l)) * t_col_base(l);
2954 ++t_col_base;
2955 ++t_m;
2956 }
2957
2958 ++t_row_base;
2959 }
2960
2961 for (; rr != row_nb_base_functions; ++rr)
2962 ++t_row_base;
2963 ++t_w;
2964 ++t_inv_D;
2965 }
2967}
2968
2970 EntData &col_data) {
2972
2979
2980 int nb_integration_pts = row_data.getN().size1();
2981 int row_nb_dofs = row_data.getIndices().size();
2982 int col_nb_dofs = col_data.getIndices().size();
2983
2984 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2986
2987 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2988
2989 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2990
2991 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2992
2993 );
2994 };
2995
2996 auto v = getVolume();
2997 auto t_w = getFTensor0IntegrationWeight();
2998 int row_nb_base_functions = row_data.getN().size2() / 3;
2999 auto t_row_base_fun = row_data.getFTensor1N<3>();
3000
3001 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3002
3003 for (int gg = 0; gg != nb_integration_pts; ++gg) {
3004 double a = v * t_w;
3005
3006 int rr = 0;
3007 for (; rr != row_nb_dofs / 3; ++rr) {
3008
3010 t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
3011
3012 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
3013 auto t_m = get_ftensor2(K, 3 * rr, 0);
3014 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3015 t_m(i, j) -= (a * t_col_base_fun) * t_PRT(i, j);
3016 ++t_m;
3017 ++t_col_base_fun;
3018 }
3019
3020 ++t_row_base_fun;
3021 }
3022
3023 for (; rr != row_nb_base_functions; ++rr)
3024 ++t_row_base_fun;
3025 ++t_w;
3026 ++t_h_domega;
3027 }
3029}
3030
3033 EntData &col_data) {
3035
3042
3043 int nb_integration_pts = row_data.getN().size1();
3044 int row_nb_dofs = row_data.getIndices().size();
3045 int col_nb_dofs = col_data.getIndices().size();
3046
3047 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
3049 &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
3050 };
3051
3052 auto v = getVolume();
3053 auto t_w = getFTensor0IntegrationWeight();
3054 int row_nb_base_functions = row_data.getN().size2() / 9;
3055 auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
3056
3057 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3058 for (int gg = 0; gg != nb_integration_pts; ++gg) {
3059 double a = v * t_w;
3060
3061 int rr = 0;
3062 for (; rr != row_nb_dofs; ++rr) {
3063
3065 t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
3066
3067 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
3068 auto t_m = get_ftensor2(K, rr, 0);
3069 for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3070 t_m(j) -= (a * t_col_base_fun) * t_PRT(j);
3071 ++t_m;
3072 ++t_col_base_fun;
3073 }
3074
3075 ++t_row_base_fun;
3076 }
3077
3078 for (; rr != row_nb_base_functions; ++rr)
3079 ++t_row_base_fun;
3080
3081 ++t_w;
3082 ++t_h_domega;
3083 }
3085}
3086
3088 EntData &data) {
3090
3091 if (tagSense != getSkeletonSense())
3093
3094 auto get_tag = [&](auto name) {
3095 auto &mob = getPtrFE()->mField.get_moab();
3096 Tag tag;
3097 CHK_MOAB_THROW(mob.tag_get_handle(name, tag), "get tag");
3098 return tag;
3099 };
3100
3101 auto get_tag_value = [&](auto &&tag, int dim) {
3102 auto &mob = getPtrFE()->mField.get_moab();
3103 auto face = getSidePtrFE()->getFEEntityHandle();
3104 std::vector<double> value(dim);
3105 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()), "set tag");
3106 return value;
3107 };
3108
3109 auto create_tag = [this](const std::string tag_name, const int size) {
3110 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3111 Tag th;
3112 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3113 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3114 def_VAL);
3115 return th;
3116 };
3117
3118 Tag th_cauchy_streess = create_tag("CauchyStress", 9);
3119 Tag th_detF = create_tag("detF", 1);
3120 Tag th_traction = create_tag("traction", 3);
3121 Tag th_disp_error = create_tag("DisplacementError", 1);
3122
3123 Tag th_energy = create_tag("Energy", 1);
3124
3125 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3126 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3127 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3128
3129 auto t_normal = getFTensor1NormalsAtGaussPts();
3130 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3131
3132 auto next = [&]() {
3133 ++t_w;
3134 ++t_h;
3135 ++t_approx_P;
3136 ++t_normal;
3137 ++t_disp;
3138 };
3139
3140 if (dataAtPts->energyAtPts.size() == 0) {
3141 // that is for case that energy is not calculated
3142 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3143 dataAtPts->energyAtPts.clear();
3144 }
3145 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
3146
3147 FTensor::Index<'i', 3> i;
3148 FTensor::Index<'j', 3> j;
3149 FTensor::Index<'k', 3> k;
3150 FTensor::Index<'l', 3> l;
3151
3152 auto set_float_precision = [](const double x) {
3153 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3154 return 0.;
3155 else
3156 return x;
3157 };
3158
3159 // scalars
3160 auto save_scal_tag = [&](auto &th, auto v, const int gg) {
3162 v = set_float_precision(v);
3163 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
3165 };
3166
3167 // vectors
3168 VectorDouble3 v(3);
3169 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
3170 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
3172 t_v(i) = t_d(i);
3173 for (auto &a : v.data())
3174 a = set_float_precision(a);
3175 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
3176 &*v.data().begin());
3178 };
3179
3180 // tensors
3181
3182 MatrixDouble3by3 m(3, 3);
3184 &m(0, 0), &m(0, 1), &m(0, 2),
3185
3186 &m(1, 0), &m(1, 1), &m(1, 2),
3187
3188 &m(2, 0), &m(2, 1), &m(2, 2));
3189
3190 auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
3192 t_m(i, j) = t_d(i, j);
3193 for (auto &v : m.data())
3194 v = set_float_precision(v);
3195 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
3196 &*m.data().begin());
3198 };
3199
3200 const auto nb_gauss_pts = getGaussPts().size2();
3201 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
3202
3203 FTensor::Tensor1<double, 3> t_traction;
3204 t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
3205 // vectors
3206 CHKERR save_vec_tag(th_traction, t_traction, gg);
3207
3208 double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
3209 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3210 CHKERR save_scal_tag(th_energy, t_energy, gg);
3211
3212 const double jac = determinantTensor3by3(t_h);
3214 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
3215 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3216 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3217
3218 next();
3219 }
3220
3222}
3223
3225 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3226 std::vector<FieldSpace> spaces, std::string geom_field_name,
3227 boost::shared_ptr<Range> crack_front_edges_ptr) {
3229
3230 constexpr bool scale_l2 = false;
3231
3232 if (scale_l2) {
3233
3235 : public MoFEM::OpGetHONormalsOnFace<SPACE_DIM> {
3236
3238
3239 OpGetHONormalsOnFace(std::string &field_name,
3240 boost::shared_ptr<Range> edges_ptr)
3241 : OP(field_name), edgesPtr(edges_ptr) {}
3242
3243 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
3244
3245 auto ent = data.getFieldEntities().size()
3246 ? data.getFieldEntities()[0]->getEnt()
3247 : 0;
3248
3249 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3250 return 0;
3251 } else {
3252 return OP::doWork(side, type, data);
3253 }
3254 };
3255
3256 private:
3257 boost::shared_ptr<Range> edgesPtr;
3258 };
3259
3260 auto jac = boost::make_shared<MatrixDouble>();
3261 auto det = boost::make_shared<VectorDouble>();
3262 pipeline.push_back(new OpGetHONormalsOnFace(
3263 geom_field_name, EshelbianCore::setSingularity
3264 ? crack_front_edges_ptr
3265 : boost::make_shared<Range>()));
3266 pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac));
3267 pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
3268 pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
3270 }
3271
3272 CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
3273
3275}
3276
3278 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3279 std::vector<FieldSpace> spaces, std::string geom_field_name,
3280 boost::shared_ptr<Range> crack_front_edges_ptr) {
3282
3283 constexpr bool scale_l2 = false;
3284
3285 if (scale_l2) {
3286 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
3287 "Scale L2 Ainsworth Legendre base is not implemented");
3288 }
3289
3290 CHKERR MoFEM::AddHOOps<2, 2, 3>::add(pipeline, spaces, geom_field_name);
3291
3293}
3294
3296 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3297 std::vector<FieldSpace> spaces, std::string geom_field_name,
3298 boost::shared_ptr<Range> crack_front_edges_ptr) {
3300
3302 auto jac = boost::make_shared<MatrixDouble>();
3303 auto det = boost::make_shared<VectorDouble>();
3305 geom_field_name, jac));
3306 pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
3307 pipeline.push_back(
3309 }
3310
3311 constexpr bool scale_l2_ainsworth_legendre_base = false;
3312
3313 if (scale_l2_ainsworth_legendre_base) {
3314
3316 : public MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM> {
3317
3319
3320 OpCalculateVectorFieldGradient(const std::string &field_name,
3321 boost::shared_ptr<MatrixDouble> jac,
3322 boost::shared_ptr<Range> edges_ptr)
3323 : OP(field_name, jac), edgesPtr(edges_ptr) {}
3324
3325 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
3326
3327 auto ent = data.getFieldEntities().size()
3328 ? data.getFieldEntities()[0]->getEnt()
3329 : 0;
3330
3331 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3332 return 0;
3333 } else {
3334 return OP::doWork(side, type, data);
3335 }
3336 };
3337
3338 private:
3339 boost::shared_ptr<Range> edgesPtr;
3340 };
3341
3342 auto jac = boost::make_shared<MatrixDouble>();
3343 auto det = boost::make_shared<VectorDouble>();
3344 pipeline.push_back(new OpCalculateVectorFieldGradient(
3345 geom_field_name, jac,
3346 EshelbianCore::setSingularity ? crack_front_edges_ptr
3347 : boost::make_shared<Range>()));
3348 pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
3349 pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
3351 }
3352
3353 CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name,
3354 nullptr, nullptr, nullptr);
3355
3357}
3358
3360 EntData &data) {
3362
3368
3369 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(), false);
3370 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(), false);
3371 if (getNinTheLoop() == 0) {
3372 dataAtPts->faceMaterialForceAtPts.clear();
3373 dataAtPts->normalPressureAtPts.clear();
3374 }
3375 auto loop_size = getLoopSize();
3376 if (loop_size == 1) {
3377 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3378 auto pstatus = numebered_fe_ptr->getPStatus();
3379 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3380 loop_size = 2;
3381 }
3382 }
3383
3385
3386 auto t_normal = getFTensor1NormalsAtGaussPts();
3388 dataAtPts->faceMaterialForceAtPts); //< face material force
3389 auto t_p =
3390 getFTensor0FromVec(dataAtPts->normalPressureAtPts); //< normal pressure
3391 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3392 auto t_grad_u_gamma =
3394 auto t_strain =
3395 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3396
3398 case GRIFFITH_FORCE:
3399 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3401 t_N(I) = t_normal(I);
3402 t_N.normalize();
3403
3405 t_grad_u(i, j) =
3406 (t_grad_u_gamma(i, j) - t_grad_u_gamma(j, i)) / 2. + t_strain(i, j);
3407 t_T(I) += t_N(J) * (t_grad_u(i, I) * t_P(i, J)) / loop_size;
3408 // note that works only for Hooke material,
3409 t_T(I) -= t_N(I) * ((t_strain(i, K) * t_P(i, K)) / 2) / loop_size;
3410
3411 t_p += t_N(I) *
3412 (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
3413 loop_size;
3414
3415 ++t_normal;
3416 ++t_P;
3417 ++t_grad_u_gamma;
3418 ++t_strain;
3419 ++t_T;
3420 ++t_p;
3421 }
3422 break;
3423 case GRIFFITH_SKELETON:
3424 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3426 t_N(I) = t_normal(I);
3427 t_N.normalize();
3428
3429 t_T(I) += t_N(J) * (t_grad_u_gamma(i, I) * t_P(i, J)) / loop_size;
3430
3431 t_p += t_N(I) *
3432 (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
3433 loop_size;
3434
3435 ++t_normal;
3436 ++t_P;
3437 ++t_grad_u_gamma;
3438 ++t_T;
3439 ++t_p;
3440 }
3441 break;
3442
3443 default:
3444 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
3445 "Grffith energy release "
3446 "selector not implemented");
3447 };
3448
3449#ifndef NDEBUG
3450 auto side_fe_ptr = getSidePtrFE();
3451 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3452 auto pstatus = side_fe_mi_ptr->getPStatus();
3453 if (pstatus) {
3454 auto owner = side_fe_mi_ptr->getOwnerProc();
3455 MOFEM_LOG("SELF", Sev::noisy)
3456 << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3457 << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
3458 << getNinTheLoop() << " loop size " << getLoopSize();
3459 }
3460#endif // NDEBUG
3461
3463}
3464
3466 EntData &data) {
3468
3469#ifndef NDEBUG
3470 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3471 auto pstatus = fe_mi_ptr->getPStatus();
3472 if (pstatus) {
3473 auto owner = fe_mi_ptr->getOwnerProc();
3474 MOFEM_LOG("SELF", Sev::noisy)
3475 << "OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3476 << " " << getPtrFE()->mField.get_comm_rank();
3477 }
3478#endif // NDEBUG
3479
3481
3483 t_face_T(I) = 0.;
3484 double face_pressure = 0.;
3486 dataAtPts->faceMaterialForceAtPts); //< face material force
3487 auto t_p =
3488 getFTensor0FromVec(dataAtPts->normalPressureAtPts); //< normal pressure
3489 auto t_w = getFTensor0IntegrationWeight();
3490 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3491 t_face_T(I) += t_w * t_T(I);
3492 face_pressure += t_w * t_p;
3493 ++t_w;
3494 ++t_T;
3495 ++t_p;
3496 }
3497 t_face_T(I) *= getMeasure();
3498 face_pressure *= getMeasure();
3499
3500 auto get_tag = [&](auto name, auto dim) {
3501 auto &moab = getPtrFE()->mField.get_moab();
3502 Tag tag;
3503 double def_val[] = {0., 0., 0.};
3504 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3505 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3506 "create tag");
3507 return tag;
3508 };
3509
3510 auto set_tag = [&](auto &&tag, auto ptr) {
3511 auto &moab = getPtrFE()->mField.get_moab();
3512 auto face = getPtrFE()->getFEEntityHandle();
3513 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr), "set tag");
3514 };
3515
3516 set_tag(get_tag("MaterialForce", 3), &t_face_T(0));
3517 set_tag(get_tag("FacePressure", 1), &face_pressure);
3518
3520}
3521
3522//! [Analytical displacement function from python]
3523#ifdef ENABLE_PYTHON_BINDING
3524
3525 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(const std::string py_file) {
3527 try {
3528
3529 // create main module
3530 auto main_module = bp::import("__main__");
3531 mainNamespace = main_module.attr("__dict__");
3532 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3533 // create a reference to python function
3534 analyticalDispFun = mainNamespace["analytical_disp"];
3535 analyticalTractionFun = mainNamespace["analytical_traction"];
3536
3537 } catch (bp::error_already_set const &) {
3538 // print all other errors to stderr
3539 PyErr_Print();
3541 }
3543 }
3544
3545 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalDisp(
3546
3547 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
3548 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3549 const std::string &block_name, np::ndarray &analytical_expr
3550
3551 ) {
3553 try {
3554
3555 // call python function
3556 analytical_expr = bp::extract<np::ndarray>(
3557 analyticalDispFun(delta_t, t, x, y, z, nx, ny, nz, block_name));
3558
3559 } catch (bp::error_already_set const &) {
3560 // print all other errors to stderr
3561 PyErr_Print();
3563 }
3565 }
3566
3567 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalTraction(
3568
3569 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
3570 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3571 const std::string& block_name, np::ndarray &analytical_expr
3572
3573 ) {
3575 try {
3576
3577 // call python function
3578 analytical_expr = bp::extract<np::ndarray>(
3579 analyticalTractionFun(delta_t, t, x, y, z, nx, ny, nz, block_name));
3580
3581 } catch (bp::error_already_set const &) {
3582 // print all other errors to stderr
3583 PyErr_Print();
3585 }
3587 }
3588
3589boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3590
3591inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
3592 int id) {
3593 auto dtype = np::dtype::get_builtin<double>();
3594 auto size = bp::make_tuple(nb_gauss_pts);
3595 auto stride = bp::make_tuple(3 * sizeof(double));
3596 return (np::from_data(&data[id], dtype, size, stride, bp::object()));
3597};
3598#endif
3599//! Analytical displacement function from python]
3600
3601
3602inline MatrixDouble analytical_expr_function(double delta_t, double t,
3603 int nb_gauss_pts,
3604 MatrixDouble &m_ref_coords,
3605 MatrixDouble &m_ref_normals,
3606 const std::string block_name) {
3607
3608#ifdef ENABLE_PYTHON_BINDING
3609 if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3610
3611 VectorDouble v_ref_coords = m_ref_coords.data();
3612 VectorDouble v_ref_normals = m_ref_normals.data();
3613
3614 bp::list python_coords;
3615 bp::list python_normals;
3616
3617 for (int idx = 0; idx < 3; ++idx) {
3618 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3619 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3620 }
3621
3622 np::ndarray np_analytical_expr = np::empty(bp::make_tuple(nb_gauss_pts, 3),
3623 np::dtype::get_builtin<double>());
3624
3625 auto disp_block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
3626 auto traction_block_name = "(.*)ANALYTICAL_TRACTION(.*)";
3627
3628 std::regex reg_disp_name(disp_block_name);
3629 std::regex reg_traction_name(traction_block_name);
3630 if (std::regex_match(block_name, reg_disp_name)) {
3631 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalDisp(
3632 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
3633 bp::extract<np::ndarray>(python_coords[1]),
3634 bp::extract<np::ndarray>(python_coords[2]),
3635 bp::extract<np::ndarray>(python_normals[0]),
3636 bp::extract<np::ndarray>(python_normals[1]),
3637 bp::extract<np::ndarray>(python_normals[2]),
3638 block_name, np_analytical_expr),
3639 "Failed analytical_disp() python call");
3640 } else if (std::regex_match(block_name, reg_traction_name)) {
3641 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalTraction(
3642 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
3643 bp::extract<np::ndarray>(python_coords[1]),
3644 bp::extract<np::ndarray>(python_coords[2]),
3645 bp::extract<np::ndarray>(python_normals[0]),
3646 bp::extract<np::ndarray>(python_normals[1]),
3647 bp::extract<np::ndarray>(python_normals[2]),
3648 block_name, np_analytical_expr),
3649 "Failed analytical_traction() python call");
3650 } else {
3652 "Unknown analytical block");
3653 }
3654
3655 double *analytical_expr_val_ptr =
3656 reinterpret_cast<double *>(np_analytical_expr.get_data());
3657
3658 MatrixDouble v_analytical_expr;
3659 v_analytical_expr.resize(3, nb_gauss_pts, false);
3660 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3661 for (int idx = 0; idx < 3; ++idx)
3662 v_analytical_expr(idx, gg) = *(analytical_expr_val_ptr + (3 * gg + idx));
3663 }
3664 return v_analytical_expr;
3665 }
3666#endif
3667 return MatrixDouble();
3668}
3669
3671 EntData &data) {
3674
3675 // auto calc_eigen_values = [&](auto &m) {
3676 // int n = m.size1();
3677 // const int lda = n;
3678 // const int size = (n + 2) * n;
3679 // int lwork = size;
3680 // double *work = new double[size];
3681
3682 // VectorDouble eig(n, 0.0);
3683 // MatrixDouble eigen_vec = m;
3684 // if (lapack_dsyev('V', 'U', n, eigen_vec.data().data(), lda,
3685 // eig.data().data(), work, lwork) > 0)
3686 // CHK_THROW_MESSAGE(MOFEM_INVALID_DATA,
3687 // "The algorithm failed to compute eigenvalues.");
3688 // delete[] work;
3689
3690 // for (auto e : eig) {
3691 // if (e < 0) {
3692 // cerr << "Negative eigenvalue: " << e << endl;
3693 // // CHK_THROW_MESSAGE(MOFEM_INVALID_DATA,
3694 // // "Negative eigenvalue found in matrix.");
3695 // }
3696 // // cerr << e << endl;
3697 // }
3698 // // cerr << endl;
3699 // };
3700
3701 // auto &m_uu = mapMatrix[{"u", "u"}];
3702 // CHKERR computeMatrixInverse(m_uu);
3703 // m_uu *= -1;
3704
3705 // auto &m_PP = mapMatrix[{"P", "P"}];
3706
3707 // auto &m_Pu = mapMatrix[{"P", "u"}];
3708 // auto &m_uP = mapMatrix[{"u", "P"}];
3709 // auto &m_Bu = mapMatrix[{"bubble", "u"}];
3710 // auto &m_uB = mapMatrix[{"u", "bubble"}];
3711
3712 // auto &m_uo = mapMatrix[{"u", "omega"}];
3713 // auto &m_ou = mapMatrix[{"omega", "u"}];
3714 // auto &m_Po = mapMatrix[{"P", "omega"}];
3715 // auto &m_oP = mapMatrix[{"omega", "P"}];
3716 // auto &m_Bo = mapMatrix[{"bubble", "omega"}];
3717 // auto &m_oB = mapMatrix[{"omega", "bubble"}];
3718 // auto &m_oo = mapMatrix[{"omega", "omega"}];
3719
3720 // auto &m_ww = mapMatrix[{"wL2", "wL2"}];
3721 // auto &m_wP = mapMatrix[{"wL2", "P"}];
3722 // auto &m_Pw = mapMatrix[{"P", "wL2"}];
3723
3724 // MatrixDouble m_uuuP;
3725 // m_uuuP = prod(m_uu, m_uP);
3726 // MatrixDouble m_uuuB;
3727 // m_uuuB = prod(m_uu, m_uB);
3728 // m_PP += prod(m_Pu, m_uuuP);
3729 // MatrixDouble m_BB;
3730 // m_BB = prod(m_Bu, m_uuuB);
3731 // MatrixDouble m_BP;
3732 // m_BP = prod(m_Bu, m_uuuP);
3733 // MatrixDouble m_PB;
3734 // m_PB = prod(m_Pu, m_uuuB);
3735 // MatrixDouble m_uuuo;
3736 // m_uuuo = prod(m_uu, m_uo);
3737 // m_oP += prod(m_ou, m_uuuP);
3738 // m_oB += prod(m_ou, m_uuuB);
3739 // m_Po += prod(m_Pu, m_uuuo);
3740 // m_Bo += prod(m_Bu, m_uuuo);
3741 // m_oo += prod(m_ou, m_uuuo);
3742 // CHKERR computeMatrixInverse(m_BB);
3743 // m_BB *= -1;
3744 // MatrixDouble m_BBBP;
3745 // m_BBBP = prod(m_BB, m_BP);
3746 // m_PP += prod(m_PB, m_BBBP);
3747 // MatrixDouble m_BBBo;
3748 // m_BBBo = prod(m_BB, m_Bo);
3749 // m_oo += prod(m_oB, m_BBBo);
3750
3751 // calc_eigen_values(m_PP);
3752
3753 // m_PP += m_PP_org;
3754 // calc_eigen_values(m_PP);
3755
3756 // CHKERR computeMatrixInverse(m_PP);
3757 // m_PP *= -1;
3758
3759
3760 // MatrixDouble m_PPPo;
3761 // m_PPPo = prod(m_PP, m_Po);
3762 // m_oo += prod(m_oP, m_PPPo);
3763
3764 // calc_eigen_values(m_oo);
3765
3766 // auto nbp = m_PP.size1();
3767 // auto nbo = m_oo.size1();
3768 // auto nb = nbp + nbo;
3769
3770 // MatrixDouble m(nb, nb);
3771 // noalias(subrange(m, 0, nbp, 0, nbp)) = m_PP;
3772 // noalias(subrange(m, nbp, nb, nbp, nb)) = m_oo;
3773 // noalias(subrange(m, 0, nbp, nbp, nb)) = m_Po;
3774 // noalias(subrange(m, nbp, nb, 0, nbp)) = m_oP;
3775
3776 // MatrixDouble inv_m_PP = m;
3777 // CHKERR computeMatrixInverse(inv_m_PP);
3778 // inv_m_PP *= -1;
3779 // MatrixDouble m_PPPo;
3780 // m_PPPo = prod(subrange(inv_m_PP, 0, nbp, 0, nbp), m_Po);
3781 // m_oo += prod(m_oP, m_PPPo);
3782
3783 // MatrixDouble m_inv_oo = m_oo;
3784 // CHKERR computeMatrixInverse(m_inv_oo);
3785 // MatrixDouble m_oooP;
3786 // m_oooP = prod(m_inv_oo, m_oP);
3787 // MatrixDouble m_PooooP;
3788 // m_PooooP = m_PP - prod(m_Po, m_oooP);
3789
3790
3791
3792 // MatrixDouble test;
3793 // test = m_oo - trans(m_oo);
3794 // cerr << test << endl;
3795
3796 // cerr << "AAA" << endl;
3797 // calc_eigen_values(m);
3798 // calc_eigen_values(m_PP);
3799 // calc_eigen_values(m_BB);
3800 // calc_eigen_values(m_PP);
3801 // calc_eigen_values(m_uu);
3802 // calc_eigen_values(m_oo);
3803
3805}
3806
3807} // namespace EshelbianPlasticity
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
constexpr double a
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ USER_BASE
user implemented approximation base
Definition definitions.h:68
@ NOBASE
Definition definitions.h:59
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
Definition definitions.h:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr auto t_kd
#define MOFEM_LOG(channel, severity)
Log.
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)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
MoFEM::TsCtx * ts_ctx
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
auto get_uniq_nb(double *ptr)
auto sort_eigen_vals(A &eig, B &eigen_vec)
static constexpr auto size_symm
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, MatrixDouble &m_ref_normals, const std::string block_name)
[Analytical displacement function from python]
auto getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static FTensor::Tensor3< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 3 (non symmetries) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static auto getFTensor4DdgFromPtr(T *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
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.
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
constexpr IntegrationType I
double scale
Definition plastic.cpp:123
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static enum SymmetrySelector symmetrySelector
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> f
static boost::function< double(const double)> d_f
static boost::function< double(const double)> inv_f
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< AnalyticalTractionBcVec > bcData
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< PressureBcVec > bcData
boost::shared_ptr< MatrixDouble > hybridGradDispPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< TractionBcVec > bcData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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, 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 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 integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(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 integrateImpl(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)
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)
static auto diffExp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:79
static auto exp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:48
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::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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.
const VectorFieldEntities & getFieldEntities() const
get field entities
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 VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Scale base functions by inverses of measure of element.