13#include <boost/math/constants/constants.hpp>
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;
36 const std::string block_name);
38template <
typename OP_PTR>
40 const std::string block_name) {
42 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
44 auto ts_time = op_ptr->getTStime();
45 auto ts_time_step = op_ptr->getTStimeStep();
51 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
52 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
55 ts_time_step, ts_time, nb_gauss_pts, m_ref_coords, m_ref_normals, block_name);
58 if (v_analytical_expr.size2() != nb_gauss_pts)
60 "Wrong number of integration pts");
63 return std::make_tuple(block_name, v_analytical_expr);
73 int nb_integration_pts = getGaussPts().size2();
75 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
76 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
81 auto t_eshelby_stress =
82 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
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);
103 int nb_integration_pts = getGaussPts().size2();
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);
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);
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,
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);
144 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
145 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
147 dataAtPts->internalStressAtPts.resize(9, nb_integration_pts,
false);
148 dataAtPts->internalStressAtPts.clear();
151 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
152 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
154 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
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 =
165 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
166 auto t_approx_P_adjoint_log_du_dP =
167 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
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);
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;
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;
185 auto t_log_stretch_total =
186 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
188 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
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);
195 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
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;
220 ++t_log_stretch_total;
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);
241 auto calculate_log_stretch = [&]() {
244 eigen_vec(
i,
j) = t_log_u(
i,
j);
246 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
250 t_nb_uniq = get_uniq_nb<3>(&eig(0));
254 t_eigen_vals(
i) = eig(
i);
255 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
258 auto get_t_diff_u = [&]() {
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);
269 auto calculate_total_stretch = [&](
auto &t_h1) {
272 t_log_u2_h1(
i,
j) = 0;
273 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
278 t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
280 return std::make_pair(t_R_h1, t_inv_u_h1);
288 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
289 eigen_vec(
i,
j) = t_C_h1(
i,
j);
291 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
294 t_nb_uniq_C = get_uniq_nb<3>(&eig(0));
295 if (t_nb_uniq_C < 3) {
301 t_eigen_vals_C(
i) = eig(
i);
302 t_eigen_vecs_C(
i,
j) = eigen_vec(
i,
j);
306 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
308 auto f_inv_sqrt = [](
auto e) {
return 1. / std::sqrt(e); };
312 t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
314 return std::make_pair(t_R_h1, t_inv_u_h1);
318 auto no_h1_loop = [&]() {
328 "rotSelector should be large");
331 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
339 auto t_Ldiff_u = calculate_log_stretch();
342 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
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);
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);
367 "rotationSelector not handled");
370 constexpr double half_r = 1 / 2.;
371 constexpr double half_l = 1 / 2.;
374 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
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.;
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)) /
389 t_levi_kirchhoff(
m) =
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)))
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.;
403 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
407 half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
409 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
412 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
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.;
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);
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);
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);
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);
439 t_levi_kirchhoff_dstreach(
m,
L) =
441 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)))
446 (t_diff_Rl(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)));
448 t_levi_kirchhoff_dP(
m,
i,
k) =
449 half_r * t_diff_Rr(
i,
l,
m) * (t_u(
l,
k))
453 half_l * t_diff_Rl(
i,
l,
m) * (t_u(
l,
k));
455 t_levi_kirchhoff_domega(
m,
n) =
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)))
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.;
474 auto large_loop = [&]() {
484 "rotSelector should be large");
487 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
501 "gradApproximator not handled");
505 auto t_Ldiff_u = calculate_log_stretch();
507 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
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);
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);
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);
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);
538 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
539 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
541 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
542 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
547 "rotationSelector not handled");
550 constexpr double half_r = 1 / 2.;
551 constexpr double half_l = 1 / 2.;
554 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
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);
563 t_levi_kirchhoff(
m) =
565 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
569 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
574 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
578 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
580 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
583 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
585 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
586 t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
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);
594 t_approx_P_adjoint_log_du_domega(
m,
L) =
595 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
599 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
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);
607 t_levi_kirchhoff_dstreach(
m,
L) =
609 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
613 half_l * (t_diff_Rl(
i,
l,
m) *
614 (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
616 t_levi_kirchhoff_dP(
m,
i,
k) =
618 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
622 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
624 t_levi_kirchhoff_domega(
m,
n) =
626 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
631 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
640 auto moderate_loop = [&]() {
650 "rotSelector should be large");
653 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
664 "gradApproximator not handled");
668 auto t_Ldiff_u = calculate_log_stretch();
670 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
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);
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);
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);
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);
701 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
702 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
704 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
705 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
710 "rotationSelector not handled");
713 constexpr double half_r = 1 / 2.;
714 constexpr double half_l = 1 / 2.;
717 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
720 t_approx_P_adjoint__dstretch(
l, o) =
721 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
723 t_approx_P_adjoint_log_du(
L) =
724 t_approx_P_adjoint__dstretch(
l, o) * t_L(
l, o,
L);
727 t_levi_kirchhoff(
m) =
729 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
733 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
738 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
742 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
744 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
747 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
749 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
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))
761 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
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);
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)))
774 half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
776 t_levi_kirchhoff_dP(
m,
i,
k) =
778 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
782 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
784 t_levi_kirchhoff_domega(
m,
n) =
786 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
791 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
800 auto small_loop = [&]() {
807 "rotSelector should be small");
810 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
819 "gradApproximator not handled");
826 t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
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);
831 t_log_u2_h1(
i,
j) = 0;
832 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
835 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
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);
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;
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;
877 "gradApproximator not handled");
890 auto N_InLoop = getNinTheLoop();
891 auto sensee = getSkeletonSense();
892 auto nb_gauss_pts = getGaussPts().size2();
893 auto t_normal = getFTensor1NormalsAtGaussPts();
896 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
898 dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
899 dataAtPts->tractionAtPts.clear();
901 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
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());
916 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
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);
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);
938 loc_moment_forces(
i) +=
939 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
940 t_traction(
k) * t_w * getMeasure();
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);
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();
970 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
972 auto piola_scale = dataAtPts->piolaScale;
973 auto alpha_w = alphaW / piola_scale;
974 auto alpha_rho = alphaRho / piola_scale;
976 int nb_base_functions = data.
getN().size2();
980 auto get_ftensor1 = [](
auto &
v) {
985 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
987 auto t_nf = get_ftensor1(nF);
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);
996 for (; bb != nb_base_functions; ++bb)
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();
1022 auto get_ftensor1 = [](
auto &
v) {
1026 auto time_step = getTStimeStep();
1028 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1031 auto t_nf = get_ftensor1(nF);
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 ) *
1036 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
1041 for (; bb != nb_base_functions; ++bb) {
1055 int nb_integration_pts = data.
getN().size1();
1056 auto v = getVolume();
1057 auto t_w = getFTensor0IntegrationWeight();
1059 int nb_base_functions = data.
getN().size2() / 3;
1067 auto get_ftensor1 = [](
auto &
v) {
1074 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1075 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1077 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1079 auto t_nf = get_ftensor1(nF);
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));
1092 for (; bb != nb_base_functions; ++bb)
1102 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1104 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1106 auto t_nf = get_ftensor1(nF);
1114 for (; bb != nb_dofs / 3; ++bb) {
1115 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1120 for (; bb != nb_base_functions; ++bb)
1134 int nb_integration_pts = data.
getN().size1();
1135 auto v = getVolume();
1136 auto t_w = getFTensor0IntegrationWeight();
1138 int nb_base_functions = data.
getN().size2() / 9;
1146 auto get_ftensor0 = [](
auto &
v) {
1152 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1153 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1155 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1157 auto t_nf = get_ftensor0(nF);
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;
1168 for (; bb != nb_base_functions; ++bb) {
1178 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1180 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1182 auto t_nf = get_ftensor0(nF);
1186 t_residuum(
i,
j) = t_h(
i,
j);
1189 for (; bb != nb_dofs; ++bb) {
1190 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1194 for (; bb != nb_base_functions; ++bb) {
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;
1215 auto get_ftensor1 = [](
auto &
v) {
1220 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1222 auto t_nf = get_ftensor1(nF);
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);
1228 ++t_row_diff_base_fun;
1230 for (; bb != nb_base_functions; ++bb) {
1231 ++t_row_diff_base_fun;
1245 int nb_integration_pts = getGaussPts().size2();
1248 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1250 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1251 if (tag_length != 9) {
1253 "Number of internal stress components should be 9 but is %d",
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);
1263 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1264 dataAtPts->internalStressAtPts.clear();
1265 auto t_internal_stress =
1266 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1269 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1270 t_internal_stress(
L) = t_const_stress(
L);
1271 ++t_internal_stress;
1283 int nb_integration_pts = getGaussPts().size2();
1286 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1288 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1289 if (tag_length != 9) {
1291 "Number of internal stress components should be 9 but is %d",
1295 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1298 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1301 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1304 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1305 dataAtPts->internalStressAtPts.clear();
1306 auto t_internal_stress =
1307 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
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;
1319 ++t_internal_stress;
1330 int nb_integration_pts = data.
getN().size1();
1331 auto v = getVolume();
1332 auto t_w = getFTensor0IntegrationWeight();
1337 auto get_ftensor2 = [](
auto &
v) {
1339 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1342 auto t_internal_stress =
1343 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1348 int nb_base_functions = data.
getN().size2();
1350 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1352 auto t_nf = get_ftensor2(nF);
1355 t_symm_stress(
i,
j) =
1356 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1359 t_residual(
L) = t_L(
i,
j,
L) * t_symm_stress(
i,
j);
1362 for (; bb != nb_dofs / 6; ++bb) {
1363 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1367 for (; bb != nb_base_functions; ++bb)
1371 ++t_internal_stress;
1381 int nb_integration_pts = data.
getN().size1();
1382 auto v = getVolume();
1383 auto t_w = getFTensor0IntegrationWeight();
1385 auto get_ftensor2 = [](
auto &
v) {
1387 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1390 auto t_internal_stress =
1391 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1398 int nb_base_functions = data.
getN().size2();
1400 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1402 auto t_nf = get_ftensor2(nF);
1405 t_residual(
L) = t_L(M,
L) * t_internal_stress(M);
1408 for (; bb != nb_dofs / 6; ++bb) {
1409 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1413 for (; bb != nb_base_functions; ++bb)
1417 ++t_internal_stress;
1422template <AssemblyType A>
1428 for (
auto &bc : (*bcDispPtr)) {
1430 if (bc.faces.find(fe_ent) != bc.faces.end()) {
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;
1443 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1445 scale *= scalingMethodsMap.at(bc.blockName)
1448 scale *= scalingMethodsMap.at(bc.blockName)
1449 ->getScale(OP::getFEMethod()->ts_t);
1453 <<
"No scaling method found for " << bc.blockName;
1460 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1461 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1463 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1465 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1469 for (; bb != nb_base_functions; ++bb)
1481 return OP::iNtegrate(data);
1484template <AssemblyType A>
1492 double time = OP::getFEMethod()->ts_t;
1500 for (
auto &bc : (*bcRotPtr)) {
1502 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1504 int nb_integration_pts = OP::getGaussPts().size2();
1505 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1506 auto t_w = OP::getFTensor0IntegrationWeight();
1508 int nb_base_functions = data.
getN().size2() / 3;
1511 auto get_ftensor1 = [](
auto &
v) {
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);
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];
1540 t_omega(
i) = OP::getFTensor1Normal()(
i);
1543 t_omega(
i) *= theta;
1550 auto t_R = get_rotation(get_rotation_angle());
1551 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
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);
1559 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
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;
1566 for (; bb != nb_base_functions; ++bb)
1579 return OP::iNtegrate(data);
1582template <AssemblyType A>
1586 double time = OP::getFEMethod()->ts_t;
1594 for (
auto &bc : (*bcDispPtr)) {
1596 if (bc.faces.find(fe_ent) != bc.faces.end()) {
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();
1609 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1610 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1613 <<
"No scaling method found for " << bc.blockName;
1617 double val =
scale * bc.val;
1620 int nb_integration_pts = OP::getGaussPts().size2();
1621 int nb_base_functions = data.
getN().size2();
1623 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1626 t_N(
i) = t_normal(
i);
1630 t_P(
i,
j) = t_N(
i) * t_N(
j);
1635 t_tracton(
i) = t_approx_P(
i,
j) * t_N(
j);
1638 t_res(
i) = t_Q(
i,
j) * t_tracton(
j) + t_P(
i,
j) * 2* t_u(
j) - t_N(
i) * val;
1640 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1642 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1643 t_nf(
i) += (t_w * t_row_base) * t_res(
i);
1647 for (; bb != nb_base_functions; ++bb)
1656 OP::locF *= OP::getMeasure();
1662template <AssemblyType A>
1668 double time = OP::getFEMethod()->ts_t;
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);
1682 for (
auto &bc : (*bcDispPtr)) {
1684 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1686 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1687 auto t_w = OP::getFTensor0IntegrationWeight();
1693 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1694 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1697 <<
"No scaling method found for " << bc.blockName;
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();
1708 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1711 t_N(
i) = t_normal(
i);
1715 t_P(
i,
j) = t_N(
i) * t_N(
j);
1718 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1721 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1722 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
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);
1733 for (; rr != nb_base_functions; ++rr)
1740 locMat *= OP::getMeasure();
1747template <AssemblyType A>
1753 double time = OP::getFEMethod()->ts_t;
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);
1767 for (
auto &bc : (*bcDispPtr)) {
1769 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1771 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1772 auto t_w = OP::getFTensor0IntegrationWeight();
1781 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1782 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1785 <<
"No scaling method found for " << bc.blockName;
1788 int nb_integration_pts = OP::getGaussPts().size2();
1789 int nb_base_functions = row_data.
getN().size2();
1792 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1795 t_N(
i) = t_normal(
i);
1799 t_P(
i,
j) = t_N(
i) * t_N(
j);
1804 t_d_res(
i,
j) = t_Q(
i,
j);
1807 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1808 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1811 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1813 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1820 for (; rr != nb_base_functions; ++rr)
1827 locMat *= OP::getMeasure();
1834 return OP::iNtegrate(data);
1839 return OP::iNtegrate(row_data, col_data);
1844 return OP::iNtegrate(row_data, col_data);
1847template <AssemblyType A>
1851 double time = OP::getFEMethod()->ts_t;
1859 for (
auto &bc : (*bcDispPtr)) {
1861 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1866 auto [block_name, v_analytical_expr] =
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;
1876 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1882 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1884 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1885 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1888 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1890 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1894 for (; bb != nb_base_functions; ++bb)
1908 return OP::iNtegrate(data);
1917 int nb_integration_pts = getGaussPts().size2();
1918 int nb_base_functions = data.
getN().size2();
1920 double time = getFEMethod()->ts_t;
1926 if (this->locF.size() != nb_dofs)
1928 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1931 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1934 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1936 auto t_w = getFTensor0IntegrationWeight();
1937 auto t_coords = getFTensor1CoordsAtGaussPts();
1939 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1941 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1943 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1944 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1946 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1947 t_f(
i) -= (time_scale * t_w * t_row_base * tau) * (t_val(
i) *
scale);
1952 for (; rr != nb_base_functions; ++rr)
1957 this->locF *= getMeasure();
1963 for (
auto &bc : *(bcData)) {
1964 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1966 double time_scale = 1;
1967 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1968 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1974 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1978 return -y * (y - 1) / 0.25;
1980 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1983 bc, [](
double,
double,
double) {
return 1; }, time_scale);
1997 int nb_integration_pts = getGaussPts().size2();
1998 int nb_base_functions = data.
getN().size2();
2000 double time = getFEMethod()->ts_t;
2006 if (this->locF.size() != nb_dofs)
2008 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2011 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2016 auto t_w = getFTensor0IntegrationWeight();
2017 auto t_coords = getFTensor1CoordsAtGaussPts();
2018 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2019 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2021 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2023 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2025 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2037 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
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));
2044 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2046 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2048 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2050 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2051 t_f(
i) += t_row_base * t_val(
i);
2056 for (; rr != nb_base_functions; ++rr)
2071 for (
auto &bc : *(bcData)) {
2072 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2074 double time_scale = 1;
2075 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2076 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2082 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2089template <AssemblyType A>
2100 double time = OP::getFEMethod()->ts_t;
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);
2113auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2118 auto t_w = OP::getFTensor0IntegrationWeight();
2119 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2120 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2121 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2123 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2126 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2136 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2137 auto t_val = time_scale * t_w * tau * val;
2140 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2141 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
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)
2152 (FTensor::levi_civita<double>(
i,
j,
k) *
2153 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2154 t_kd(
k,
l) * t_diff_col_base(N1);
2156 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2163 for (; rr != nb_base_functions; ++rr)
2179 for (
auto &bc : *(bcData)) {
2180 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2182 double time_scale = 1;
2183 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2184 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2190 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2201 return OP::iNtegrate(row_data, col_data);
2211 int nb_integration_pts = getGaussPts().size2();
2212 int nb_base_functions = data.
getN().size2();
2214 double time = getFEMethod()->ts_t;
2220 if (this->locF.size() != nb_dofs)
2222 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2227 for (
auto &bc : *(bcData)) {
2228 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2232 auto [block_name, v_analytical_expr] =
2235 getFTensor1FromMat<3>(v_analytical_expr);
2237 auto t_w = getFTensor0IntegrationWeight();
2238 auto t_coords = getFTensor1CoordsAtGaussPts();
2240 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2242 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2244 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2246 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2247 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2252 for (; rr != nb_base_functions; ++rr)
2258 this->locF *= getMeasure();
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));
2275 auto v = getVolume();
2276 auto t_w = getFTensor0IntegrationWeight();
2277 int row_nb_base_functions = row_data.
getN().size2();
2279 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2282 for (; rr != row_nb_dofs / 3; ++rr) {
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;
2289 ++t_col_diff_base_fun;
2293 for (; rr != row_nb_base_functions; ++rr)
2304 if (alphaW < std::numeric_limits<double>::epsilon() &&
2305 alphaRho < std::numeric_limits<double>::epsilon())
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)
2318 auto v = getVolume();
2319 auto t_w = getFTensor0IntegrationWeight();
2321 auto piola_scale = dataAtPts->piolaScale;
2322 auto alpha_w = alphaW / piola_scale;
2323 auto alpha_rho = alphaRho / piola_scale;
2325 int row_nb_base_functions = row_data.
getN().size2();
2328 double ts_scale = alpha_w * getTSa();
2329 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2330 ts_scale += alpha_rho * getTSaa();
2332 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2333 double a =
v * t_w * ts_scale;
2336 for (; rr != row_nb_dofs / 3; ++rr) {
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;
2350 for (; rr != row_nb_base_functions; ++rr)
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) {
2375 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2377 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2379 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2381 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2383 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2385 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2388 auto v = getVolume();
2389 auto t_w = getFTensor0IntegrationWeight();
2391 int row_nb_base_functions = row_data.
getN().size2();
2394 auto t_approx_P_adjoint_log_du_dP =
2395 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2397 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2400 for (; rr != row_nb_dofs / 6; ++rr) {
2403 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2405 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2407 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2415 for (; rr != row_nb_base_functions; ++rr)
2418 ++t_approx_P_adjoint_log_du_dP;
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),
2443 auto v = getVolume();
2444 auto t_w = getFTensor0IntegrationWeight();
2447 int row_nb_base_functions = row_data.
getN().size2();
2449 auto t_approx_P_adjoint_log_du_dP =
2450 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2452 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
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) {
2460 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2467 for (; rr != row_nb_base_functions; ++rr)
2470 ++t_approx_P_adjoint_log_du_dP;
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) {
2487 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2489 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2491 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2493 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2495 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2497 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
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);
2512 int row_nb_base_functions = row_data.
getN().size2();
2515 int nb_integration_pts = row_data.
getN().size1();
2517 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2521 for (; rr != row_nb_dofs / 6; ++rr) {
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);
2533 for (; rr != row_nb_base_functions; ++rr)
2537 ++t_approx_P_adjoint_log_du_domega;
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) {
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),
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),
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)
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();
2574 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
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;
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);
2588 for (; rr != row_nb_base_functions; ++rr) {
2592 ++t_levi_kirchhoff_du;
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) {
2614 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2616 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2618 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2623 auto v = getVolume();
2624 auto t_w = getFTensor0IntegrationWeight();
2626 int row_nb_base_functions = row_data.
getN().size2();
2628 auto t_levi_kirchhoff_dP =
2629 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2631 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2634 for (; rr != row_nb_dofs / 3; ++rr) {
2635 double b =
a * t_row_base_fun;
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));
2645 for (; rr != row_nb_base_functions; ++rr) {
2650 ++t_levi_kirchhoff_dP;
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();
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));
2671 auto v = getVolume();
2672 auto t_w = getFTensor0IntegrationWeight();
2673 auto t_levi_kirchoff_dP =
2674 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2676 int row_nb_base_functions = row_data.
getN().size2();
2679 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
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));
2694 for (; rr != row_nb_base_functions; ++rr) {
2698 ++t_levi_kirchoff_dP;
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) {
2712 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2714 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2716 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
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();
2738 auto time_step = getTStimeStep();
2739 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2741 double c = (
a * alphaOmega) * (ts_a );
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;
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)));
2759 for (; rr != row_nb_base_functions; ++rr) {
2764 ++t_levi_kirchhoff_domega;
2772 if (dataAtPts->matInvD.size1() ==
size_symm &&
2773 dataAtPts->matInvD.size2() ==
size_symm) {
2786 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2789 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2791 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2793 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
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();
2802 auto v = getVolume();
2803 auto t_w = getFTensor0IntegrationWeight();
2804 int row_nb_base_functions = row_data.
getN().size2() / 3;
2811 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2812 &*dataAtPts->matInvD.data().begin());
2815 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2819 for (; rr != row_nb_dofs / 3; ++rr) {
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));
2831 for (; rr != row_nb_base_functions; ++rr)
2844 if (dataAtPts->matInvD.size1() ==
size_symm &&
2845 dataAtPts->matInvD.size2() ==
size_symm) {
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();
2863 auto v = getVolume();
2864 auto t_w = getFTensor0IntegrationWeight();
2865 int row_nb_base_functions = row_data.
getN().size2() / 9;
2872 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2873 &*dataAtPts->matInvD.data().begin());
2876 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2880 for (; rr != row_nb_dofs; ++rr) {
2882 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2884 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2891 for (; rr != row_nb_base_functions; ++rr)
2902 if (dataAtPts->matInvD.size1() ==
size_symm &&
2903 dataAtPts->matInvD.size2() ==
size_symm) {
2917 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2920 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
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();
2929 auto v = getVolume();
2930 auto t_w = getFTensor0IntegrationWeight();
2931 int row_nb_base_functions = row_data.
getN().size2() / 9;
2940 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2941 &*dataAtPts->matInvD.data().begin());
2944 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2947 auto t_m = get_ftensor1(
K, 0, 0);
2950 for (; rr != row_nb_dofs; ++rr) {
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);
2961 for (; rr != row_nb_base_functions; ++rr)
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();
2984 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2987 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2989 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2991 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2996 auto v = getVolume();
2997 auto t_w = getFTensor0IntegrationWeight();
2998 int row_nb_base_functions = row_data.
getN().size2() / 3;
3001 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3003 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3007 for (; rr != row_nb_dofs / 3; ++rr) {
3010 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
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);
3023 for (; rr != row_nb_base_functions; ++rr)
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();
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));
3052 auto v = getVolume();
3053 auto t_w = getFTensor0IntegrationWeight();
3054 int row_nb_base_functions = row_data.
getN().size2() / 9;
3057 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3058 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3062 for (; rr != row_nb_dofs; ++rr) {
3065 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
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);
3078 for (; rr != row_nb_base_functions; ++rr)
3091 if (tagSense != getSkeletonSense())
3094 auto get_tag = [&](
auto name) {
3095 auto &mob = getPtrFE()->mField.get_moab();
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");
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};
3112 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3113 th, MB_TAG_CREAT | MB_TAG_SPARSE,
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);
3123 Tag th_energy = create_tag(
"Energy", 1);
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);
3129 auto t_normal = getFTensor1NormalsAtGaussPts();
3130 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3140 if (dataAtPts->energyAtPts.size() == 0) {
3142 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3143 dataAtPts->energyAtPts.clear();
3152 auto set_float_precision = [](
const double x) {
3153 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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);
3170 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
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());
3184 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3186 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3188 &
m(2, 0), &
m(2, 1), &
m(2, 2));
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());
3200 const auto nb_gauss_pts = getGaussPts().size2();
3201 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3204 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3206 CHKERR save_vec_tag(th_traction, t_traction, gg);
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);
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);
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) {
3230 constexpr bool scale_l2 =
false;
3234 "Scale L2 Ainsworth Legendre base is not implemented");
3243 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3244 std::vector<FieldSpace> spaces, std::string geom_field_name,
3245 boost::shared_ptr<Range> crack_front_edges_ptr) {
3248 constexpr bool scale_l2 =
false;
3252 "Scale L2 Ainsworth Legendre base is not implemented");
3261 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3262 std::vector<FieldSpace> spaces, std::string geom_field_name,
3263 boost::shared_ptr<Range> crack_front_edges_ptr,
3264 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3265 boost::shared_ptr<MatrixDouble> inv_jac) {
3269 auto jac = boost::make_shared<MatrixDouble>();
3270 auto det = boost::make_shared<VectorDouble>();
3272 geom_field_name, jac));
3278 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3280 if (scale_l2_ainsworth_legendre_base) {
3288 boost::shared_ptr<MatrixDouble> jac,
3289 boost::shared_ptr<Range> edges_ptr)
3298 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3301 return OP::doWork(side, type, data);
3306 boost::shared_ptr<Range> edgesPtr;
3309 auto jac = boost::make_shared<MatrixDouble>();
3310 auto det = boost::make_shared<VectorDouble>();
3312 geom_field_name, jac,
3314 : boost::make_shared<
Range>()));
3384 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3385 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3386 if (getNinTheLoop() == 0) {
3387 dataAtPts->faceMaterialForceAtPts.clear();
3388 dataAtPts->normalPressureAtPts.clear();
3390 auto loop_size = getLoopSize();
3391 if (loop_size == 1) {
3392 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3393 auto pstatus = numebered_fe_ptr->getPStatus();
3394 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3401 auto t_normal = getFTensor1NormalsAtGaussPts();
3402 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3403 dataAtPts->faceMaterialForceAtPts);
3406 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3410 getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
3411 auto t_grad_u_gamma =
3412 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3414 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3415 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
3437 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3438 t_N(
I) = t_normal(
I);
3441 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3443 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3446 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3449 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3451 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3457 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3460 t_N(
I) = t_normal(
I);
3465 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3469 t_grad_u_gamma(
i,
J) +
3470 (2 * t_R(
i,
K) * t_N(
K) - (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3473 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3476 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3479 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3487 "Grffith energy release "
3488 "selector not implemented");
3492 auto side_fe_ptr = getSidePtrFE();
3493 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3494 auto pstatus = side_fe_mi_ptr->getPStatus();
3496 auto owner = side_fe_mi_ptr->getOwnerProc();
3498 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3499 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3500 << getNinTheLoop() <<
" loop size " << getLoopSize();
3512 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3513 auto pstatus = fe_mi_ptr->getPStatus();
3515 auto owner = fe_mi_ptr->getOwnerProc();
3517 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3518 <<
" " << getPtrFE()->mField.get_comm_rank();
3526 double face_pressure = 0.;
3527 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3528 dataAtPts->faceMaterialForceAtPts);
3531 auto t_w = getFTensor0IntegrationWeight();
3532 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3533 t_face_T(
I) += t_w * t_T(
I);
3534 face_pressure += t_w * t_p;
3539 t_face_T(
I) *= getMeasure();
3540 face_pressure *= getMeasure();
3542 auto get_tag = [&](
auto name,
auto dim) {
3543 auto &moab = getPtrFE()->mField.get_moab();
3545 double def_val[] = {0., 0., 0.};
3546 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3547 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3552 auto set_tag = [&](
auto &&tag,
auto ptr) {
3553 auto &moab = getPtrFE()->mField.get_moab();
3554 auto face = getPtrFE()->getFEEntityHandle();
3555 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3558 set_tag(get_tag(
"MaterialForce", 3), &t_face_T(0));
3559 set_tag(get_tag(
"FacePressure", 1), &face_pressure);
3565#ifdef ENABLE_PYTHON_BINDING
3567 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3572 auto main_module = bp::import(
"__main__");
3573 mainNamespace = main_module.attr(
"__dict__");
3574 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3576 analyticalDispFun = mainNamespace[
"analytical_disp"];
3577 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3579 }
catch (bp::error_already_set
const &) {
3589 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3590 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3591 const std::string &block_name, np::ndarray &analytical_expr
3598 analytical_expr = bp::extract<np::ndarray>(
3599 analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3601 }
catch (bp::error_already_set
const &) {
3611 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3612 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3613 const std::string& block_name, np::ndarray &analytical_expr
3620 analytical_expr = bp::extract<np::ndarray>(
3621 analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3623 }
catch (bp::error_already_set
const &) {
3631boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3633inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3635 auto dtype = np::dtype::get_builtin<double>();
3636 auto size = bp::make_tuple(nb_gauss_pts);
3637 auto stride = bp::make_tuple(3 *
sizeof(
double));
3638 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3648 const std::string block_name) {
3650#ifdef ENABLE_PYTHON_BINDING
3651 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3656 bp::list python_coords;
3657 bp::list python_normals;
3659 for (
int idx = 0; idx < 3; ++idx) {
3660 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3661 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3664 np::ndarray np_analytical_expr = np::empty(bp::make_tuple(nb_gauss_pts, 3),
3665 np::dtype::get_builtin<double>());
3667 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3668 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3670 std::regex reg_disp_name(disp_block_name);
3671 std::regex reg_traction_name(traction_block_name);
3672 if (std::regex_match(block_name, reg_disp_name)) {
3674 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3675 bp::extract<np::ndarray>(python_coords[1]),
3676 bp::extract<np::ndarray>(python_coords[2]),
3677 bp::extract<np::ndarray>(python_normals[0]),
3678 bp::extract<np::ndarray>(python_normals[1]),
3679 bp::extract<np::ndarray>(python_normals[2]),
3680 block_name, np_analytical_expr),
3681 "Failed analytical_disp() python call");
3682 }
else if (std::regex_match(block_name, reg_traction_name)) {
3684 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3685 bp::extract<np::ndarray>(python_coords[1]),
3686 bp::extract<np::ndarray>(python_coords[2]),
3687 bp::extract<np::ndarray>(python_normals[0]),
3688 bp::extract<np::ndarray>(python_normals[1]),
3689 bp::extract<np::ndarray>(python_normals[2]),
3690 block_name, np_analytical_expr),
3691 "Failed analytical_traction() python call");
3694 "Unknown analytical block");
3697 double *analytical_expr_val_ptr =
3698 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3701 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3702 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3703 for (
int idx = 0; idx < 3; ++idx)
3704 v_analytical_expr(idx, gg) = *(analytical_expr_val_ptr + (3 * gg + idx));
3706 return v_analytical_expr;
3724 "OpZeroNoSideBrokenDofs works only for discontinuous field");
3726 auto fe_type = OP::getFEType();
3729 auto face_side_number = getFaceSideNumber();
3733 auto it = side_dof_map.get<1>().begin();
3734 it != side_dof_map.get<1>().end(); ++it
3737 if ((3 * it->dof) >= data.
getIndices().size())
3739 if (it->type != MBTRI || it->side != face_side_number) {
3740 for (
auto dd = 0; dd < 3; ++dd) {
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
Tensor1< T, Tensor_Dim > normalize()
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#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
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
#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
@ MOFEM_DATA_INCONSISTENCY
#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 ...
#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
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 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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
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.
constexpr IntegrationType I
constexpr double t
plate stiffness
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static double dynamicTime
static PetscBool l2UserBaseScale
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static enum SymmetrySelector symmetrySelector
static boost::function< double(const double)> f
static PetscBool setSingularity
static boost::function< double(const double)> d_f
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_f
static auto diffExp(A &&t_w_vee, B &&theta)
static auto exp(A &&t_w_vee, B &&theta)
Add operators pushing bases from local to physical configuration.
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
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 (const version)
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 DOF values on entity.
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 degrees of freedom on entity.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Operator for inverting matrices at integration points.
Scale base functions by inverses of measure of element.
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Caluclate face material force and normal pressure at gauss points.
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 &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
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 &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)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
