#include <boost/math/constants/constants.hpp>
int nb_integration_pts = getGaussPts().size2();
auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
false);
auto t_eshelby_stress =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
++t_energy;
++t_P;
++t_F;
++t_eshelby_stress;
}
}
int nb_integration_pts = getGaussPts().size2();
dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts, false);
dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 *
size_symm,
nb_integration_pts, false);
dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
false);
dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
dataAtPts->nbUniq.resize(nb_integration_pts, false);
dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
auto t_h_dlog_u =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
auto t_levi_kirchhoff_domega =
getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
dataAtPts->leviKirchhoffdLogStreatchAtPts);
auto t_levi_kirchhoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
auto t_approx_P_adjoint_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_approx_P_adjoint_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
auto t_approx_P_adjoint_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
auto t_approx_P_adjoint_log_du_domega =
getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
auto t_log_stretch_total =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
auto t_log_u2_h1 =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
auto &nbUniq = dataAtPts->nbUniq;
auto t_nb_uniq =
auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
auto t_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
auto next = [&]() {
++t_h;
++t_h_domega;
++t_h_dlog_u;
++t_levi_kirchhoff;
++t_levi_kirchhoff_domega;
++t_levi_kirchhoff_dstreach;
++t_levi_kirchhoff_dP;
++t_approx_P_adjoint_dstretch;
++t_approx_P_adjoint_log_du;
++t_approx_P_adjoint_log_du_dP;
++t_approx_P_adjoint_log_du_domega;
++t_R;
++t_u;
++t_diff_u;
++t_eigen_vals;
++t_eigen_vecs;
++t_nb_uniq;
++t_log_u2_h1;
++t_log_stretch_total;
++t_omega;
++t_grad_h1;
++t_approx_P;
++t_log_u;
};
auto bound_eig = [&](auto &eig) {
const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
const auto large = std::exp(std::numeric_limits<double>::max_exponent);
for (int ii = 0; ii != 3; ++ii) {
eig(ii) = std::min(std::max(zero, eig(ii)), large);
}
};
auto calculate_log_stretch = [&]() {
eigen_vec(
i,
j) = t_log_u(
i,
j);
MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
}
t_nb_uniq = get_uniq_nb<3>(&eig(0));
if (t_nb_uniq < 3) {
}
t_eigen_vals(
i) = eig(
i);
t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
auto get_t_diff_u = [&]() {
t_nb_uniq);
};
t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
return t_Ldiff_u;
};
auto calculate_total_stretch = [&](auto &t_h1) {
t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
return std::make_pair(t_R_h1, t_inv_u_h1);
} else {
t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
eigen_vec(
i,
j) = t_C_h1(
i,
j);
MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
}
}
t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
auto f_inv_sqrt = [](auto e) { return 1. / std::sqrt(e); };
t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
return std::make_pair(t_R_h1, t_inv_u_h1);
}
};
auto no_h1_loop = [&]() {
break;
break;
default:
"rotSelector should be large");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_Ldiff_u = calculate_log_stretch();
auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
break;
default:
"rotationSelector not handled");
}
constexpr double half_r = 1 / 2.;
constexpr double half_l = 1 / 2.;
t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
t_approx_P_adjoint_dstretch(
l,
k) =
(t_R(
i,
l) * t_approx_P(
i,
k) + t_R(
i,
k) * t_approx_P(
i,
l));
t_approx_P_adjoint_dstretch(
l,
k) /= 2.;
t_approx_P_adjoint_log_du(
L) =
t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
half_r * (t_diff_Rr(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
t_diff_Rr(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)))
+
half_l * (t_diff_Rl(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
t_diff_Rl(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)));
t_levi_kirchhoff(
m) /= 2.;
t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
+
half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
} else {
"symmetrySelector not handled");
}
t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
t_A(
k,
l,
m) = t_diff_Rr(
i,
l,
m) * t_approx_P(
i,
k) +
t_diff_Rr(
i,
k,
m) * t_approx_P(
i,
l);
t_B(
k,
l,
m) = t_diff_Rl(
i,
l,
m) * t_approx_P(
i,
k) +
t_diff_Rl(
i,
k,
m) * t_approx_P(
i,
l);
t_approx_P_adjoint_log_du_domega(
m,
L) =
half_r * (t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L)) +
half_l * (t_B(
k,
l,
m) * t_Ldiff_u(
k,
l,
L));
t_levi_kirchhoff_domega(
m,
n) =
half_r *
(t_diff_diff_Rr(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
t_diff_diff_Rr(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)))
+
half_l *
(t_diff_diff_Rl(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
t_diff_diff_Rl(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)));
t_levi_kirchhoff_domega(
m,
n) /= 2.;
}
next();
}
};
auto large_loop = [&]() {
break;
break;
default:
"rotSelector should be large");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
break;
break;
default:
"gradApproximator not handled");
};
auto t_Ldiff_u = calculate_log_stretch();
auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
t_h_u(
l,
k) = t_u(
l, o) * t_h1(o,
k);
t_Ldiff_h_u(
l,
k,
L) = t_Ldiff_u(
l, o,
L) * t_h1(o,
k);
break;
t_diff_diff_Rr(
i,
j,
l,
m) = 0;
t_diff_diff_Rl(
i,
j,
l,
m) = 0;
break;
default:
"rotationSelector not handled");
}
constexpr double half_r = 1 / 2.;
constexpr double half_l = 1 / 2.;
t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
t_approx_P_adjoint_dstretch(
l, o) =
(t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
t_approx_P_adjoint_log_du(
L) =
t_approx_P_adjoint_dstretch(
l, o) * t_Ldiff_u(
l, o,
L);
half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)))
+
half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)));
t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
+
half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
} else {
t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
}
t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
t_approx_P_adjoint_log_du_domega(
m,
L) =
half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
+
half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
} else {
t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
t_approx_P_adjoint_log_du_domega(
m,
L) =
t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
}
t_levi_kirchhoff_dstreach(
m,
L) =
half_r *
(t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
+
half_l * (t_diff_Rl(
i,
l,
m) *
(t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
t_levi_kirchhoff_dP(
m,
i,
k) =
half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
+
half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
t_levi_kirchhoff_domega(
m,
n) =
half_r *
(t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
+
half_l *
(t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
}
next();
}
};
auto moderate_loop = [&]() {
break;
break;
default:
"rotSelector should be large");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
break;
default:
"gradApproximator not handled");
};
auto t_Ldiff_u = calculate_log_stretch();
auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
t_h_u(
l,
k) = (t_symm_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
t_L_h(
l,
k,
L) = t_L(
l, o,
L) * t_h1(o,
k);
break;
t_diff_diff_Rr(
i,
j,
l,
m) = 0;
t_diff_diff_Rl(
i,
j,
l,
m) = 0;
break;
default:
"rotationSelector not handled");
}
constexpr double half_r = 1 / 2.;
constexpr double half_l = 1 / 2.;
t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
t_approx_P_adjoint_dstretch(
l, o) =
(t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
t_approx_P_adjoint_log_du(
L) =
t_approx_P_adjoint_dstretch(
l, o) * t_L(
l, o,
L);
half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)))
+
half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)));
t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
+
half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
} else {
t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
}
t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_L_h(
l,
k,
L);
t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_L_h(
l,
k,
L);
t_approx_P_adjoint_log_du_domega(
m,
L) =
half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
+
half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
} else {
t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_L_h(
l,
k,
L);
t_approx_P_adjoint_log_du_domega(
m,
L) =
t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
}
t_levi_kirchhoff_dstreach(
m,
L) =
half_r * (t_diff_Rr(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)))
+
half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
t_levi_kirchhoff_dP(
m,
i,
k) =
half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
+
half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
t_levi_kirchhoff_domega(
m,
n) =
half_r *
(t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
+
half_l *
(t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
}
next();
}
};
auto small_loop = [&]() {
break;
default:
"rotSelector should be small");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
break;
default:
"gradApproximator not handled");
};
t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
} else {
t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
}
t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
t_approx_P_adjoint_dstretch(
i,
j) = t_approx_P(
i,
j);
t_approx_P_adjoint_log_du(
L) =
t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
t_levi_kirchhoff_dstreach(
m,
L) = 0;
t_levi_kirchhoff_domega(
m,
n) = 0;
next();
}
};
break;
break;
break;
break;
default:
"gradApproximator not handled");
break;
};
}
auto n_total_side_eles = getLoopSize();
auto N_InLoop = getNinTheLoop();
auto sensee = getSkeletonSense();
auto nb_gauss_pts = getGaussPts().size2();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_sigma_ptr =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
if (N_InLoop == 0) {
dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
dataAtPts->tractionAtPts.clear();
}
auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
++t_traction;
++t_sigma_ptr;
++t_normal;
}
}
if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
};
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
loc_reaction_forces(
i) += (t_traction(
i))*t_w * getMeasure();
t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
(FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
t_traction(
k) * t_w * getMeasure();
++t_coords;
++t_spatial_disp;
++t_w;
++t_traction;
}
reactionVec[0] += loc_reaction_forces(0);
reactionVec[1] += loc_reaction_forces(1);
reactionVec[2] += loc_reaction_forces(2);
reactionVec[3] += loc_moment_forces(0);
reactionVec[4] += loc_moment_forces(1);
reactionVec[5] += loc_moment_forces(2);
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
dataAtPts->wL2DotDotAtPts.clear();
}
auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
auto piola_scale = dataAtPts->piolaScale;
auto alpha_w = alphaW / piola_scale;
auto alpha_rho = alphaRho / piola_scale;
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_div_P;
++t_s_dot_w;
++t_s_dot_dot_w;
}
}
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
auto t_omega_grad_dot =
getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
(
a * alphaOmega) * (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
++t_nf;
++t_row_base_fun;
++t_row_grad_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
++t_row_grad_fun;
}
++t_w;
++t_levi_kirchhoff;
++t_omega_grad_dot;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_R;
++t_u;
}
} else {
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_h;
}
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 9;
auto get_ftensor0 = [](
auto &
v) {
};
auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor0(nF);
int bb = 0;
for (; bb != nb_dofs; ++bb) {
t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_R;
++t_u;
}
} else {
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor0(nF);
t_residuum(
i,
j) = t_h(
i,
j);
int bb = 0;
for (; bb != nb_dofs; ++bb) {
t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_h;
}
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
double div_row_base = t_row_diff_base_fun(
i,
i);
t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
++t_nf;
++t_row_diff_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_diff_base_fun;
}
++t_w;
++t_w_l2;
}
}
template <AssemblyType A>
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 3;
for (auto &sm : scalingMethodsVec) {
} else {
scale *= sm->getScale(OP::getFEMethod()->ts_t);
}
}
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
int bb = 0;
t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_normal;
}
}
}
}
return OP::iNtegrate(data);
}
template <AssemblyType A>
for (auto &bc : (*bcRotPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
double theta = bc.theta;
} else {
theta *= OP::getFEMethod()->ts_t;
}
const double a = sqrt(t_normal(
i) * t_normal(
i));
t_omega(
i) = t_normal(
i) * (theta /
a);
? 0.
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_delta(
i) = t_center(
i) - t_coords(
i);
t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_normal;
++t_coords;
}
}
}
}
return OP::iNtegrate(data);
}
int nb_integration_pts = getGaussPts().size2();
int nb_base_functions = data.
getN().size2();
double ts_t = getFEMethod()->ts_t;
}
double time_scale = 1;
for (auto &sm : scalingMethodsVec) {
time_scale *= sm->getScale(ts_t);
}
#ifndef NDEBUG
if (this->locF.size() != nb_dofs)
"Size of locF %d != nb_dofs %d", this->locF.size(), nb_dofs);
#endif // NDEBUG
auto integrate_rhs = [&](auto &bc, auto calc_tau) {
auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
int rr = 0;
t_f(
i) -= time_scale * t_w * t_row_base * tau * (t_val(
i) *
scale);
++t_row_base;
++t_f;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_coords;
}
this->locF *= 2. * getMeasure();
};
for (auto &bc : *(bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
if (nb_dofs) {
if (std::regex_match(bc.blockName, std::regex(".*COOK.*"))) {
y -= 44;
y /= (60 - 44);
return -y * (y - 1) / 0.25;
};
CHKERR integrate_rhs(bc, calc_tau);
} else {
CHKERR integrate_rhs(bc, [](
double,
double,
double) {
return 1; });
}
}
}
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
double div_col_base = t_col_diff_base_fun(
i,
i);
t_m(
i) -=
a * t_row_base_fun * div_col_base;
++t_m;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
if (alphaW < std::numeric_limits<double>::epsilon() &&
alphaRho < std::numeric_limits<double>::epsilon())
const int nb_integration_pts = row_data.
getN().size1();
const int row_nb_dofs = row_data.
getIndices().size();
&
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto piola_scale = dataAtPts->piolaScale;
auto alpha_w = alphaW / piola_scale;
auto alpha_rho = alphaRho / piola_scale;
int row_nb_base_functions = row_data.
getN().size2();
double ts_scale = alpha_w * getTSa();
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
ts_scale += alpha_rho * getTSaa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a =
v * t_w * ts_scale;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
K, 3 * rr, 0);
for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
const double b =
a * t_row_base_fun * t_col_base_fun;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
&
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
&
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
&
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
auto t_approx_P_adjoint_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor3(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
t_row_base_fun;
++t_col_base_fun;
++t_m;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_approx_P_adjoint_log_du_dP;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c), &
m(
r + 3,
c), &
m(
r + 4,
c),
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
auto t_approx_P_adjoint_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor2(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
t_row_base_fun;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_approx_P_adjoint_log_du_dP;
}
}
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
&
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
&
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
&
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjoint_log_du_domega =
getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
int row_nb_base_functions = row_data.
getN().size2();
int nb_integration_pts = row_data.
getN().size1();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor3(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
double v =
a * t_row_base_fun * t_col_base_fun;
t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_approx_P_adjoint_log_du_domega;
}
}
int nb_integration_pts = getGaussPts().size2();
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
&
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
&
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
&
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5)
};
};
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
dataAtPts->leviKirchhoffdLogStreatchAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(
K, 3 * rr, 0);
const double b =
a * t_row_base_fun;
for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr) {
++t_row_base_fun;
}
++t_w;
++t_levi_kirchhoff_du;
}
}
int nb_integration_pts = getGaussPts().size2();
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
auto t_levi_kirchhoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
double b =
a * t_row_base_fun;
auto t_m = get_ftensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr) {
++t_row_base_fun;
}
++t_w;
++t_levi_kirchhoff_dP;
}
}
int nb_integration_pts = getGaussPts().size2();
};
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
double b =
a * t_row_base_fun;
auto t_m = get_ftensor1(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr) {
++t_row_base_fun;
}
++t_w;
++t_levi_kirchoff_dP;
}
}
int nb_integration_pts = getGaussPts().size2();
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
);
};
auto ts_a = getTSa();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchhoff_domega =
getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double c =
a * alphaOmega * ts_a;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(
K, 3 * rr, 0);
const double b =
a * t_row_base_fun;
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
++t_m;
++t_col_base_fun;
++t_col_grad_fun;
}
++t_row_base_fun;
++t_row_grad_fun;
}
for (; rr != row_nb_base_functions; ++rr) {
++t_row_base_fun;
++t_row_grad_fun;
}
++t_w;
++t_levi_kirchhoff_domega;
}
}
if (dataAtPts->matInvD.size1() ==
size_symm &&
return integrateImpl<0>(row_data, col_data);
} else {
return integrateImpl<size_symm>(row_data, col_data);
}
};
template <int S>
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
);
};
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 3;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
&*dataAtPts->matInvD.data().begin());
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
++t_m;
++t_col_base;
}
++t_row_base;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_inv_D;
}
}
if (dataAtPts->matInvD.size1() ==
size_symm &&
return integrateImpl<0>(row_data, col_data);
} else {
return integrateImpl<size_symm>(row_data, col_data);
}
};
template <int S>
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 9;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
&*dataAtPts->matInvD.data().begin());
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
for (int cc = 0; cc != col_nb_dofs; ++cc) {
a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
++t_col_base;
}
++t_row_base;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_inv_D;
}
}
if (dataAtPts->matInvD.size1() ==
size_symm &&
return integrateImpl<0>(row_data, col_data);
} else {
return integrateImpl<size_symm>(row_data, col_data);
}
};
template <int S>
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2)
);
};
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 9;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
&*dataAtPts->matInvD.data().begin());
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_m = get_ftensor1(
K, 0, 0);
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
++t_col_base;
++t_m;
}
++t_row_base;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_inv_D;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 3;
const double ts_a = getTSa();
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
auto t_m = get_ftensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_h_domega;
}
}
int nb_integration_pts = row_data.
getN().size1();
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 9;
const double ts_a = getTSa();
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
auto t_m = get_ftensor2(
K, rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_h_domega;
}
}
if (tagSense != getSkeletonSense())
auto get_tag = [&](auto name) {
auto &mob = getPtrFE()->mField.get_moab();
Tag tag;
return tag;
};
auto get_tag_value = [&](auto &&tag, int dim) {
auto &mob = getPtrFE()->mField.get_moab();
auto face = getSidePtrFE()->getFEEntityHandle();
std::vector<double> value(dim);
CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
return value;
};
auto create_tag = [this](const std::string tag_name, const int size) {
double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
th, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
};
Tag th_cauchy_streess = create_tag("CauchyStress", 9);
Tag th_detF = create_tag("detF", 1);
Tag th_traction = create_tag("traction", 3);
Tag th_disp_error = create_tag("DisplacementError", 1);
Tag th_energy = create_tag("Energy", 1);
auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
auto next = [&]() {
++t_w;
++t_h;
++t_approx_P;
++t_normal;
++t_disp;
};
if (dataAtPts->energyAtPts.size() == 0) {
dataAtPts->energyAtPts.resize(getGaussPts().size2());
dataAtPts->energyAtPts.clear();
}
auto set_float_precision = [](const double x) {
if (std::abs(x) < std::numeric_limits<float>::epsilon())
return 0.;
else
return x;
};
auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
v = set_float_precision(
v);
CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
};
auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
a = set_float_precision(
a);
CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
};
&
m(0, 0), &
m(0, 1), &
m(0, 2),
&
m(1, 0), &
m(1, 1), &
m(1, 2),
&
m(2, 0), &
m(2, 1), &
m(2, 2));
auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
v = set_float_precision(
v);
CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
};
const auto nb_gauss_pts = getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
CHKERR save_vec_tag(th_traction, t_traction, gg);
double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
CHKERR save_scal_tag(th_disp_error, u_error, gg);
CHKERR save_scal_tag(th_energy, t_energy, gg);
t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
next();
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::vector<FieldSpace> spaces, std::string geom_field_name,
boost::shared_ptr<Range> crack_front_edges_ptr) {
constexpr bool scale_l2 = false;
if (scale_l2) {
boost::shared_ptr<Range> edges_ptr)
: 0;
if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
return 0;
} else {
return OP::doWork(side,
type, data);
}
};
private:
boost::shared_ptr<Range> edgesPtr;
};
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
? crack_front_edges_ptr
: boost::make_shared<Range>()));
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::vector<FieldSpace> spaces, std::string geom_field_name,
boost::shared_ptr<Range> crack_front_edges_ptr) {
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
geom_field_name, jac));
pipeline.push_back(
}
constexpr bool scale_l2_ainsworth_legendre_base = false;
if (scale_l2_ainsworth_legendre_base) {
boost::shared_ptr<MatrixDouble> jac,
boost::shared_ptr<Range> edges_ptr)
: 0;
if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
return 0;
} else {
return OP::doWork(side,
type, data);
}
};
private:
boost::shared_ptr<Range> edgesPtr;
};
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
geom_field_name, jac,
: boost::make_shared<Range>()));
}
nullptr, nullptr, nullptr);
}
const auto nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
auto t_var_omega = getFTensor1FromMat<3>(dataAtPts->varRotAxis);
auto t_var_log_u = getFTensor2SymmetricFromMat<3>(dataAtPts->varLogStreach);
auto t_var_P = getFTensor2FromMat<3, 3>(dataAtPts->varPiola);
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
auto t_h_dlog_u =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
auto t_approx_P_adjoint_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
double var_element_energy = 0.;
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
t_var_L_u(
L) = t_L(
J,
K,
L) * t_var_log_u(
J,
K);
auto var_energy = t_P(
i,
J) * (t_h_domega(
i,
J,
j) * t_var_omega(
j) +
t_h_dlog_u(
i,
J,
L) * t_var_L_u(
L));
var_element_energy += t_w * var_energy;
auto var_complementary = t_var_P(
i,
J) * t_h(
i,
J);
var_element_energy += t_w * var_complementary;
++t_w;
++t_h;
++t_P;
++t_var_omega;
++t_var_log_u;
++t_var_P;
++t_h_domega;
++t_h_dlog_u;
++t_approx_P_adjoint_log_du;
}
var_element_energy *= getMeasure();
auto get_tag = [&]() {
auto &mob = getPtrFE()->mField.get_moab();
Tag tag;
double def_val[] = {0.};
CHK_MOAB_THROW(mob.tag_get_handle(
"ReleaseEnergy", 1, MB_TYPE_DOUBLE, tag,
MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
"create tag");
return tag;
};
auto set_tag = [&](auto &&tag, auto &energy) {
auto &mob = getPtrFE()->mField.get_moab();
auto tet = getPtrFE()->getFEEntityHandle();
};
set_tag(get_tag(), var_element_energy);
*releaseEnergyPtr += var_element_energy;
}
}
