Implementation of operators.
#include <boost/math/constants/constants.hpp>
int nb_integration_pts = getGaussPts().size2();
auto t_P = dataAtPts->getFTensorApproxP(getGaussPts().size2());
auto t_F = dataAtPts->getFTensorSmallH(getGaussPts().size2());
auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
auto get_eshelby_stress =
DL>::size(dataAtPts->SigmaAtPts, nb_integration_pts);
auto t_eshelby_stress = get_eshelby_stress();
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->getStretchTensorAtPts(), nb_integration_pts);
*dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
*dataAtPts->getStretchH1AtPts(), nb_integration_pts);
*dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
*dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
*dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
*dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
*dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
*dataAtPts->getDeformationGradient(), nb_integration_pts);
dataAtPts->hdOmegaAtPts, nb_integration_pts);
dataAtPts->hdLogStretchAtPts, nb_integration_pts);
dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
dataAtPts->rotMatAtPts, nb_integration_pts);
*dataAtPts->getEigenVals(), nb_integration_pts);
*dataAtPts->getEigenVecs(), nb_integration_pts);
dataAtPts->nbUniq.resize(nb_integration_pts, false);
dataAtPts->eigenValsC, nb_integration_pts);
dataAtPts->eigenVecsC, nb_integration_pts);
dataAtPts->nbUniqC.resize(nb_integration_pts, false);
dataAtPts->logStretch2H1AtPts, nb_integration_pts);
dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
dataAtPts->internalStressAtPts, nb_integration_pts);
dataAtPts->internalStressAtPts.clear();
auto t_h = dataAtPts->getFTensorSmallH(getGaussPts().size2());
auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(getGaussPts().size2());
auto t_h_dlog_u =
dataAtPts->getFTensorSmallHdLogStretch(getGaussPts().size2());
auto t_levi_kirchhoff =
dataAtPts->getFTensorLeviKirchhoff(getGaussPts().size2());
auto t_levi_kirchhoff_domega =
dataAtPts->getFTensorLeviKirchhoffdOmega(getGaussPts().size2());
auto t_levi_kirchhoff_dstreach =
dataAtPts->getFTensorLeviKirchhoffdLogStretch(getGaussPts().size2());
auto t_levi_kirchhoff_dP =
dataAtPts->getFTensorLeviKirchhoffP(getGaussPts().size2());
auto t_approx_P_adjoint_dstretch =
dataAtPts->getFTensorAdjointPdstretch(getGaussPts().size2());
auto t_approx_P_adjoint_log_du =
dataAtPts->getFTensorAdjointPdU(getGaussPts().size2());
auto t_approx_P_adjoint_log_du_dP =
dataAtPts->getFTensorAdjointPdUdP(getGaussPts().size2());
auto t_approx_P_adjoint_log_du_domega =
dataAtPts->getFTensorAdjointPdUdOmega(getGaussPts().size2());
auto t_R = dataAtPts->getFTensorRotMat(getGaussPts().size2());
auto t_u = dataAtPts->getFTensorStretch(getGaussPts().size2());
auto t_diff_u = dataAtPts->getFTensorDiffStretch(getGaussPts().size2());
auto t_eigen_vals = dataAtPts->getFTensorEigenVals(getGaussPts().size2());
auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(getGaussPts().size2());
auto &nbUniq = dataAtPts->nbUniq;
auto t_nb_uniq =
auto t_eigen_vals_C = dataAtPts->getFTensorEigenValsC(nb_integration_pts);
auto t_eigen_vecs_C = dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
auto &nbUniqC = dataAtPts->nbUniqC;
auto t_nb_uniq_C =
auto t_u_h1 = dataAtPts->getFTensorStretchH1(getGaussPts().size2());
auto t_diff_u_h1 = dataAtPts->getFTensorDiffStretchH1(getGaussPts().size2());
auto t_log_stretch_total =
dataAtPts->getFTensorLogStretchTotal(getGaussPts().size2());
auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(getGaussPts().size2());
auto t_grad_h1 = dataAtPts->getFTensorSmallWGradH1(getGaussPts().size2());
auto t_omega = dataAtPts->getFTensorRotAxis(getGaussPts().size2());
auto t_approx_P = dataAtPts->getFTensorApproxP(getGaussPts().size2());
auto t_log_u = dataAtPts->getFTensorLogStretch(getGaussPts().size2());
auto t_omega0 = dataAtPts->getFTensorRotAxis0(getGaussPts().size2());
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_eigen_vals_C;
++t_eigen_vecs_C;
++t_nb_uniq_C;
++t_u_h1;
++t_diff_u_h1;
++t_log_u2_h1;
++t_log_stretch_total;
++t_omega;
++t_omega0;
++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);
if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
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);
};
auto calculate_total_stretch = [&](auto &t_h1) {
t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
} else {
t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
t_eigen_vec(
i,
j) = t_C_h1(
i,
j);
if (computeEigenValuesSymmetric(t_eigen_vec, t_eig) != MB_SUCCESS) {
MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
}
t_nb_uniq_C = get_uniq_nb<3>(&t_eig(0));
if (t_nb_uniq_C < 3) {
}
}
t_eigen_vals_C(
i) = t_eig(
i);
t_eigen_vecs_C(
i,
j) = t_eigen_vec(
i,
j);
t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
}
};
auto no_h1_loop = [&]() {
break;
default:
"no_h1_loop is implemented only for LARGE_ROT");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
CHKERR calculate_log_stretch();
CHKERR calculate_total_stretch(t_h1);
t_u_h1(
i,
j) = t_u(
i,
j);
t_diff_u_h1(
i,
j,
k,
l) = 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);
auto large_rot = [&]() {
t_diff_diff_R(
i,
j,
k,
l) =
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_approx_P_adjoint_log_du(L) =
t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k, L);
t_diff_R(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k));
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
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_R(
i,
l,
m) * t_approx_P(
i,
k);
t_approx_P_adjoint_log_du_domega(
m, L) =
t_A(
k,
l,
m) * t_Ldiff_u(
k,
l, L);
t_levi_kirchhoff_dstreach(
m, L) =
t_diff_R(
i,
l,
m) * (t_Ldiff_u(
l,
k, L) * t_approx_P(
i,
k));
t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
t_levi_kirchhoff_domega(
m,
n) =
t_diff_diff_R(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k));
}
};
auto moderate_rot = [&]() {
t_delta_omega(
m) = t_omega(
m) - t_omega0(
m);
(t_u(
l,
k) + levi_civita(
l,
k,
m) * t_delta_omega(
m));
t_approx_P_adjoint_dstretch(
l,
k) = t_R0(
i,
l) * t_approx_P(
i,
k);
t_approx_P_adjoint_log_du(L) =
t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k, L);
t_levi_kirchhoff(
m) = t_R0(
i,
l) * levi_civita(
l,
k,
m) *
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
t_h_domega(
i,
k,
m) = t_R0(
i,
l) * levi_civita(
l,
k,
m);
t_h_dlog_u(
i,
k, L) = t_R0(
i,
l) * t_Ldiff_u(
l,
k, L);
t_approx_P_adjoint_log_du_dP(
i,
k, L) =
t_R0(
i,
l) * t_Ldiff_u(
l,
k, L);
t_approx_P_adjoint_log_du_domega(
m, L) = 0;
t_levi_kirchhoff_dstreach(
m, L) = 0;
t_levi_kirchhoff_dP(
m,
i,
k) =
t_R0(
i,
l) * levi_civita(
l,
k,
m);
t_levi_kirchhoff_domega(
m,
n) = 0;
}
};
large_rot();
break;
moderate_rot();
break;
default:
"rotationSelector not handled");
}
next();
}
};
auto large_loop = [&]() {
break;
break;
default:
"rotSelector should be large or small");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
break;
default:
"Selected grad approximator not handled");
};
CHKERR calculate_log_stretch();
CHKERR calculate_total_stretch(t_h1);
t_u_h1(
l,
k) = t_u(
l, o) * t_h1(o,
k);
t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i, o,
k,
l) * t_h1(o,
j);
t_Ldiff_u_h1(
l,
k, L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j, L);
t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
t_diff_diff_R(
i,
j,
l,
m) = 0;
break;
t_diff_diff_R(
i,
j,
k,
l) =
break;
default:
"rotationSelector not handled");
}
t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
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_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k, L);
t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
t_h_dlog_u(
i,
k, L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k, L);
t_approx_P_adjoint_log_du_dP(
i,
k, L) =
t_R(
i,
l) * t_Ldiff_u_h1(
l,
k, L);
t_A(
m, L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
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) =
t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k, L) * t_approx_P(
i,
k));
t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
t_levi_kirchhoff_domega(
m,
n) =
t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
}
next();
}
};
auto moderate_loop = [&]() {
break;
break;
default:
"rotSelector should be large or small");
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
break;
default:
"Selected grad approximator not handled");
};
CHKERR calculate_log_stretch();
CHKERR calculate_total_stretch(t_h1);
t_u_h1(
l,
k) = (
t_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
t_diff_u_h1(
i,
j,
k,
l) = t_diff(
i, o,
k,
l) * t_h1(o,
j);
t_Ldiff_u_h1(
l,
k, L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j, L);
t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
t_diff_diff_R(
i,
j,
l,
m) = 0;
break;
t_diff_diff_R(
i,
j,
k,
l) =
break;
default:
"rotationSelector not handled");
}
t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
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_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k, L);
t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
t_h_dlog_u(
i,
k, L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k, L);
t_approx_P_adjoint_log_du_dP(
i,
k, L) =
t_R(
i,
l) * t_Ldiff_u_h1(
l,
k, L);
t_A(
m, L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
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) =
t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k, L) * t_approx_P(
i,
k));
t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
t_levi_kirchhoff_domega(
m,
n) =
t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
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");
};
"stretchSelector should be linear for small loop");
} else {
t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
t_u_h1(
i,
j) = t_u(
i,
j);
t_diff_u_h1(
i,
j,
k,
l) =
t_diff_u_h1(
i,
j,
k,
l) /= 2.;
t_Ldiff_u(
i,
j, L) = t_L(
i,
j, L);
}
t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
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(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
t_levi_kirchhoff_dstreach(
m, L) = 0;
t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
t_levi_kirchhoff_domega(
m,
n) = 0;
next();
}
};
break;
break;
break;
break;
default:
"gradApproximator not handled");
break;
};
}
auto N_InLoop = getNinTheLoop();
auto sensee = getSkeletonSense();
auto nb_gauss_pts = getGaussPts().size2();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_sigma_ptr = dataAtPts->getFTensorApproxP(getGaussPts().size2());
auto get_tracion =
dataAtPts->tractionAtPts, nb_gauss_pts);
if (N_InLoop == 0) {
dataAtPts->tractionAtPts.clear();
}
auto t_traction = get_tracion();
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 = dataAtPts->getFTensorTraction(nb_integration_pts);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_spatial_disp = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
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 = dataAtPts->getFTensorDivP(nb_integration_pts);
auto t_s_dot_w = dataAtPts->getFTensorSmallWL2Dot(nb_integration_pts);
auto w_l2_dot_dot_at_pts = dataAtPts->getSmallWL2DotDotAtPts();
const bool reset_w_l2_dot_dot =
w_l2_dot_dot_at_pts->size1() != nb_integration_pts ||
w_l2_dot_dot_at_pts->size2() != 3;
*w_l2_dot_dot_at_pts, nb_integration_pts);
if (reset_w_l2_dot_dot) {
w_l2_dot_dot_at_pts->clear();
}
auto t_s_dot_dot_w = dataAtPts->getFTensorSmallWL2DotDot(nb_integration_pts);
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) {
};
auto next = [&]() {
++t_w;
++t_div_P;
++t_s_dot_w;
++t_s_dot_dot_w;
};
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;
next();
}
}
int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchhoff =
dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts);
auto t_omega_grad_dot =
dataAtPts->getFTensorRotAxisGradDot(nb_integration_pts);
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);
t_nf(
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_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
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_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
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 = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
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 <>
int nb_integration_pts = getGaussPts().size2();
CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
int tag_length;
CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
if (tag_length != 9) {
"Number of internal stress components should be 9 but is %d",
tag_length);
}
auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
tag, &fe_ent, 1, &*const_stress_vec.data().begin());
auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
auto get_internal_stress =
dataAtPts->internalStressAtPts, nb_integration_pts);
dataAtPts->internalStressAtPts.clear();
auto t_internal_stress = get_internal_stress();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_internal_stress(L) = t_const_stress(L);
++t_internal_stress;
}
}
template <>
int nb_integration_pts = getGaussPts().size2();
CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
int tag_length;
CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
if (tag_length != 9) {
"Number of internal stress components should be 9 but is %d",
tag_length);
}
auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
int vert_num;
CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
vert_num, true);
CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
&vert_data[0]);
auto get_internal_stress =
dataAtPts->internalStressAtPts, nb_integration_pts);
dataAtPts->internalStressAtPts.clear();
auto t_internal_stress = get_internal_stress();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
for (int bb = 0; bb != nb_shape_fn; ++bb) {
t_internal_stress(L) += t_vert_data(L) * t_shape_n;
++t_vert_data;
++t_shape_n;
}
++t_internal_stress;
}
}
template <>
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto get_ftensor2 = [](
auto &
v) {
&
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
};
auto t_internal_stress =
dataAtPts->getFTensorInternalStress(nb_integration_pts);
: getFEMethod()->ts_t;
double scale = scalingMethodPtr->getScale(time);
int nb_base_functions = data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor2(nF);
(t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
t_residual(L) = t_L(
i,
j, L) * (
scale * t_symm_stress(
i,
j));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(L) +=
a * t_row_base_fun * t_residual(L);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_internal_stress;
}
}
template <>
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto get_ftensor2 = [](
auto &
v) {
&
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
};
auto t_internal_stress =
dataAtPts->getFTensorInternalStressVec(nb_integration_pts);
: getFEMethod()->ts_t;
double scale = scalingMethodPtr->getScale(time);
int nb_base_functions = data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor2(nF);
t_residual(L) = t_L(
M, L) * (
scale * t_internal_stress(
M));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(L) +=
a * t_row_base_fun * t_residual(L);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_internal_stress;
}
}
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;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
scale *= scalingMethodsMap.at(bc.blockName)
} else {
scale *= scalingMethodsMap.at(bc.blockName)
->getScale(OP::getFEMethod()->ts_t);
}
} else {
<< "No scaling method found for " << bc.blockName;
}
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);
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2();
#ifndef NDEBUG
if (!this->sourceVec) {
"Source vector for OpTauStabilizationDispRhsBc is not set");
}
if (data.
getN().size1() != nb_integration_pts) {
"Number of integration points in data should be %d but is %d",
nb_integration_pts, (
int)data.
getN().size1());
}
if (nb_base_functions < nb_dofs /
SPACE_DIM) {
"Number of base functions in data should be %d but is %d",
}
#endif
auto t_disp_val =
*this->sourceVec, nb_integration_pts)();
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
scale *= scalingMethodsMap.at(bc.blockName)
} else {
scale *= scalingMethodsMap.at(bc.blockName)
->getScale(OP::getFEMethod()->ts_t);
}
} else {
<< "No scaling method found for " << bc.blockName;
}
auto area = getMeasure();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
int bb = 0;
if (bc.flags[ii]) {
t_nf(ii) += (tau_scale * t_row_base_fun) *
(t_disp_val(ii) - t_bc_disp(ii));
}
}
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_coords;
++t_disp_val;
}
}
}
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = row_data.
getN().size2();
auto get_t_vec = [&](const int rr) {
std::array<double *, SPACE_DIM> ptrs;
ptrs[
i] = &OP::locMat(rr +
i,
i);
};
auto area = getMeasure();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
int rr = 0;
for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
if (bc.flags[ii]) {
t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
}
}
++t_col_base_fun;
++t_mat;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_coords;
}
}
}
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto &v_analytical_expr = std::get<1>(analytical_data);
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2();
#ifndef NDEBUG
if (!this->sourceVec) {
"Source vector for OpTauStabilizationOpAnalyticalDispBc is not "
"set");
}
if (data.
getN().size1() != nb_integration_pts) {
"Number of integration points in data should be %d but is %d",
nb_integration_pts, (
int)data.
getN().size1());
}
if (nb_base_functions < nb_dofs /
SPACE_DIM) {
"Number of base functions in data should be at least %d but is "
"%d",
}
#endif
auto t_disp_val =
*this->sourceVec, nb_integration_pts)();
auto t_bc_disp = getFTensor1FromMat<3, -1,
DL>(v_analytical_expr);
auto area = getMeasure();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
int bb = 0;
if (bc.flags[ii]) {
t_nf(ii) +=
(tau_scale * t_row_base_fun) * (t_disp_val(ii) - t_bc_disp(ii));
}
}
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_coords;
++t_disp_val;
++t_bc_disp;
}
}
}
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = row_data.
getN().size2();
auto get_t_vec = [&](const int rr) {
std::array<double *, SPACE_DIM> ptrs;
ptrs[
i] = &OP::locMat(rr +
i,
i);
};
auto area = getMeasure();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
int rr = 0;
for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
if (bc.flags[ii]) {
t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
}
}
++t_col_base_fun;
++t_mat;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_coords;
}
}
}
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
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) {
};
auto get_rotation_angle = [&]() {
double theta = bc.theta;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
}
return theta;
};
auto get_rotation = [&](auto theta) {
if (bc.vals.size() == 7) {
t_omega(0) = bc.vals[4];
t_omega(1) = bc.vals[5];
t_omega(2) = bc.vals[6];
} else {
t_omega(
i) = OP::getFTensor1Normal()(
i);
}
if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
} else {
<< "Rotation axis is zero vector for block " << bc.blockName
<< ". This may lead to unexpected results.";
}
? 0.
};
auto t_R = get_rotation(get_rotation_angle());
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;
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);
}
double time = OP::getFEMethod()->ts_t;
}
for (auto &bc : (*bcRotPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
auto get_rotation_angle = [&]() {
double theta = bc.theta;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
}
return theta;
};
auto get_rotation = [&](auto theta) {
if (bc.vals.size() == 7) {
t_omega(0) = bc.vals[4];
t_omega(1) = bc.vals[5];
t_omega(2) = bc.vals[6];
} else {
t_omega(
i) = OP::getFTensor1Normal()(
i);
}
if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
}
? 0.
};
auto area = getMeasure();
auto t_R = get_rotation(get_rotation_angle());
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
auto t_disp_val =
*this->sourceVec, nb_integration_pts)();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
t_delta(
i) = t_center(
i) - t_coords(
i);
t_bc_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
int bb = 0;
(tau_scale * t_row_base_fun) * (t_disp_val(
i) - t_bc_disp(
i));
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_coords;
++t_disp_val;
}
}
}
}
for (auto &bc : (*bcRotPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = OP::getGaussPts().size2();
auto t_w = OP::getFTensor0IntegrationWeight();
int nb_base_functions = row_data.
getN().size2();
auto get_t_vec = [&](const int rr) {
std::array<double *, SPACE_DIM> ptrs;
ptrs[
i] = &OP::locMat(rr +
i,
i);
};
auto area = getMeasure();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau_scale =
area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
int rr = 0;
for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
}
++t_col_base_fun;
++t_mat;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_coords;
}
}
}
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
for (auto &bd : (*brokenBaseSideDataPtr)) {
auto t_approx_P = getFTensor2FromMat<3, 3, -1,
DL>(bd.getFlux());
auto t_u = getFTensor1FromMat<3, -1,
DL>(*hybridDispPtr);
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_w = OP::getFTensor0IntegrationWeight();
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
} else {
<< "No scaling method found for " << bc.blockName;
}
double val =
scale * bc.val;
int nb_integration_pts = OP::getGaussPts().size2();
int nb_base_functions = data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_P(
i,
j) = t_N(
i) * t_N(
j);
t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
int bb = 0;
t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
++t_nf;
++t_row_base;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base;
++t_w;
++t_normal;
++t_u;
++t_approx_P;
}
}
}
}
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
auto &locMat = OP::locMat;
locMat.resize(row_nb_dofs, col_nb_dofs, false);
locMat.clear();
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_w = OP::getFTensor0IntegrationWeight();
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
} else {
<< "No scaling method found for " << bc.blockName;
}
int nb_integration_pts = OP::getGaussPts().size2();
int nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_P(
i,
j) = t_N(
i) * t_N(
j);
t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
int rr = 0;
for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_normal;
}
locMat *= OP::getMeasure();
}
}
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
auto &locMat = OP::locMat;
locMat.resize(row_nb_dofs, col_nb_dofs, false);
locMat.clear();
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_w = OP::getFTensor0IntegrationWeight();
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
} else {
<< "No scaling method found for " << bc.blockName;
}
int nb_integration_pts = OP::getGaussPts().size2();
int nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_P(
i,
j) = t_N(
i) * t_N(
j);
t_d_res(
i,
j) = t_Q(
i,
j);
int rr = 0;
for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_normal;
}
locMat *= OP::getMeasure();
}
}
}
return OP::iNtegrate(data);
}
return OP::iNtegrate(row_data, col_data);
}
return OP::iNtegrate(row_data, col_data);
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto [block_name, v_analytical_expr] =
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 t_coord = OP::getFTensor1CoordsAtGaussPts();
auto t_bc_disp = getFTensor1FromMat<3, -1,
DL>(v_analytical_expr);
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_bc_disp;
++t_coord;
++t_w;
++t_normal;
}
}
}
}
return OP::iNtegrate(data);
}
int nb_integration_pts = getGaussPts().size2();
int nb_base_functions = data.
getN().size2();
double time = getFEMethod()->ts_t;
}
#ifndef NDEBUG
if (this->locF.size() != nb_dofs)
"Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif
auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = sqrt(t_normal(
i) * t_normal(
i));
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;
(time_scale *
a * 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;
++t_normal;
}
};
for (auto &bc : *(bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
double time_scale = 1;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
}
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, time_scale);
} else {
bc, [](double, double, double) { return 1; }, time_scale);
}
}
}
}
}
int nb_integration_pts = getGaussPts().size2();
int nb_base_functions = data.
getN().size2();
double time = getFEMethod()->ts_t;
}
#ifndef NDEBUG
if (this->locF.size() != nb_dofs)
"Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif
auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
auto val = bc.val;
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
auto t_grad_gamma_u = getFTensor2FromMat<3, 2, -1,
DL>(*hybridGradDispPtr);
double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
} else {
t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
(t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
(t_tangent2(
k) + t_grad_gamma_u(
k, N1));
}
auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
int rr = 0;
t_f(
i) += t_row_base * t_val(
i);
++t_row_base;
++t_f;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_coords;
++t_tangent1;
++t_tangent2;
++t_grad_gamma_u;
}
this->locF /= 2.;
};
for (auto &bc : *(bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
double time_scale = 1;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
}
if (nb_dofs) {
bc, [](double, double, double) { return 1; }, time_scale);
}
}
}
}
template <AssemblyType A>
}
double time = OP::getFEMethod()->ts_t;
}
int nb_base_functions = row_data.
getN().size2();
int nb_integration_pts = OP::getGaussPts().size2();
auto &locMat = OP::locMat;
locMat.resize(row_nb_dofs, col_nb_dofs, false);
locMat.clear();
auto integrate_lhs = [&](auto &bc, auto calc_tau, double time_scale) {
auto val = bc.val;
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_coords = OP::getFTensor1CoordsAtGaussPts();
auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
auto t_grad_gamma_u = getFTensor2FromMat<3, 2, -1,
DL>(*hybridGradDispPtr);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
auto t_val = time_scale * t_w * tau * val;
int rr = 0;
for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
(t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
t_kd(
j,
l) * t_diff_col_base(N0)
+
(FTensor::levi_civita<double>(
i,
j,
k) *
(t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
t_kd(
k,
l) * t_diff_col_base(N1);
t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
++t_mat;
++t_diff_col_base;
}
++t_row_base;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_coords;
++t_tangent1;
++t_tangent2;
++t_grad_gamma_u;
}
OP::locMat /= 2.;
};
for (auto &bc : *(bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
double time_scale = 1;
if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
}
if (nb_dofs) {
bc, [](double, double, double) { return 1; }, time_scale);
}
}
}
}
return OP::iNtegrate(row_data, col_data);
}
int nb_integration_pts = getGaussPts().size2();
int nb_base_functions = data.
getN().size2();
double time = getFEMethod()->ts_t;
}
#ifndef NDEBUG
if (this->locF.size() != nb_dofs)
"Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif
for (auto &bc : *(bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto [block_name, v_analytical_expr] =
auto t_val = getFTensor1FromMat<3, -1,
DL>(v_analytical_expr);
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
int rr = 0;
t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
++t_row_base;
++t_f;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_coords;
++t_val;
}
this->locF *= getMeasure();
}
}
}
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 =
dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
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 =
dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
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) {
t_m(L) -=
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;
}
}
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),
&
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 =
dataAtPts->getFTensorAdjointPdUdOmega(nb_integration_pts);
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 / 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 =
dataAtPts->getFTensorLeviKirchhoffdLogStretch(nb_integration_pts);
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 =
dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
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();
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
};
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff_dP =
dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
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 = std::abs(alphaOmega) > std::numeric_limits<double>::epsilon()
? getTSa()
: 0.0;
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchhoff_domega =
dataAtPts->getFTensorLeviKirchhoffdOmega(nb_integration_pts);
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 (hasNonhomogeneousMatBlock) {
integrateImpl<size_symm * size_symm>(row_data, col_data));
} else {
}
};
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 =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
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 (hasNonhomogeneousMatBlock) {
integrateImpl<size_symm * size_symm>(row_data, col_data));
} else {
}
};
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 =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
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) {
K(rr, 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 (hasNonhomogeneousMatBlock) {
integrateImpl<size_symm * size_symm>(row_data, col_data));
} else {
}
};
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 =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
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;
}
}
&
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;
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(
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;
}
}
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;
auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
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();
&
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2() / 9;
auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
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 &mob = getPtrFE()->mField.get_moab();
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);
return th;
};
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);
Tag th_young_modulus = create_tag(
"YoungModulus", 1);
const auto nb_gauss_pts = getGaussPts().size2();
auto t_w = dataAtPts->getFTensorSmallWL2(nb_gauss_pts);
auto t_h = dataAtPts->getFTensorSmallH(nb_gauss_pts);
auto t_approx_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_disp = dataAtPts->getFTensorSmallWH1(nb_gauss_pts);
if (dataAtPts->energyAtPts.size() == 0) {
dataAtPts->energyAtPts.resize(nb_gauss_pts);
dataAtPts->energyAtPts.clear();
}
auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
auto t_youngs_modulus = getFTensor0FromVec(dataAtPts->youngModulusAtPts);
auto next = [&]() {
++t_w;
++t_h;
++t_approx_P;
++t_normal;
++t_disp;
++t_youngs_modulus;
++t_energy;
};
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,
};
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
t_traction(
i) *= tagSense;
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)));
if (!std::isfinite(u_error))
u_error = -1.;
CHKERR save_scal_tag(th_disp_error, u_error, gg);
CHKERR save_scal_tag(th_energy, t_energy, gg);
CHKERR save_scal_tag(th_young_modulus, t_youngs_modulus, gg);
const double jac = determinantTensor3by3(t_h);
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) {
"Scale L2 Ainsworth Legendre base is not implemented");
}
}
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) {
"Scale L2 Ainsworth Legendre base is not implemented");
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::vector<FieldSpace> spaces, std::string geom_field_name,
boost::shared_ptr<Range> crack_front_edges_ptr,
boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
boost::shared_ptr<MatrixDouble> inv_jac) {
if (!geom_field_name.empty()) {
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
pipeline.push_back(
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;
};
if (!geom_field_name.empty()) {
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
geom_field_name, jac,
: boost::make_shared<Range>()));
}
}
det, inv_jac);
}
const auto nb_gauss_pts = getGaussPts().size2();
dataAtPts->faceMaterialForceAtPts, nb_gauss_pts);
dataAtPts->normalPressureAtPts.resize(nb_gauss_pts, false);
if (getNinTheLoop() == 0) {
dataAtPts->faceMaterialForceAtPts.clear();
dataAtPts->normalPressureAtPts.clear();
}
auto loop_size = getLoopSize();
if (loop_size == 1) {
auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
auto pstatus = numebered_fe_ptr->getPStatus();
if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
loop_size = 2;
}
}
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_T = dataAtPts->getFTensorFaceMaterialForce(
nb_gauss_pts);
auto t_p =
getFTensor0FromVec(dataAtPts->normalPressureAtPts);
auto t_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
auto t_u_gamma = dataAtPts->getFTensorSmallHybridDisp(nb_gauss_pts);
auto t_grad_u_gamma = dataAtPts->getFTensorGradHybridDisp(nb_gauss_pts);
auto t_strain = dataAtPts->getFTensorLogStretch(nb_gauss_pts);
auto t_omega = dataAtPts->getFTensorRotAxis(nb_gauss_pts);
auto next = [&]() {
++t_normal;
++t_P;
++t_omega;
++t_u_gamma;
++t_grad_u_gamma;
++t_strain;
++t_T;
++t_p;
};
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
t_T(
I) -= t_N(
I) * ((t_strain(
i, K) * t_P(
i, K)) / 2.) / loop_size;
(t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
loop_size;
next();
}
break;
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
(2 * t_R(
i, K) * t_N(K) - (t_R(
k, L) * t_N(
k) * t_N(L)) * t_N(
i)) *
t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
t_T(
I) -= t_N(
I) * ((t_strain(
i, K) * t_P(
i, K)) / 2.) / loop_size;
(t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
loop_size;
next();
}
break;
default:
"Grffith energy release "
"selector not implemented");
};
#ifndef NDEBUG
auto side_fe_ptr = getSidePtrFE();
auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
auto pstatus = side_fe_mi_ptr->getPStatus();
if (pstatus) {
auto owner = side_fe_mi_ptr->getOwnerProc();
<< "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
<< " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
<< getNinTheLoop() << " loop size " << getLoopSize();
}
#endif
}
#ifndef NDEBUG
auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
auto pstatus = fe_mi_ptr->getPStatus();
if (pstatus) {
auto owner = fe_mi_ptr->getOwnerProc();
<< "OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
<< " " << getPtrFE()->mField.get_comm_rank();
}
#endif
double face_pressure = 0.;
auto t_T = dataAtPts->getFTensorFaceMaterialForce(
getGaussPts().size2());
auto t_p =
getFTensor0FromVec(dataAtPts->normalPressureAtPts);
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_face_T(
I) += t_w * t_T(
I);
face_pressure += t_w * t_p;
++t_w;
++t_T;
++t_p;
}
t_face_T(
I) *= getMeasure();
face_pressure *= getMeasure();
auto get_tag = [&](
auto name,
auto dim) {
auto &moab = getPtrFE()->mField.get_moab();
double def_val[] = {0., 0., 0.};
MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
"create tag");
return tag;
};
auto set_tag = [&](auto &&tag, auto ptr) {
auto &moab = getPtrFE()->mField.get_moab();
auto face = getPtrFE()->getFEEntityHandle();
};
set_tag(
get_tag(
"MaterialForce", 3), &t_face_T(0));
set_tag(
get_tag(
"FacePressure", 1), &face_pressure);
}
template <typename OP_PTR>
std::tuple<std::string, MatrixDouble>
const std::string block_name) {
auto nb_gauss_pts = op_ptr->getGaussPts().size2();
auto ts_time = op_ptr->getTStime();
auto ts_time_step = op_ptr->getTStimeStep();
}
MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
auto v_analytical_expr =
m_ref_coords, m_ref_normals, block_name);
if (PetscUnlikely(!v_analytical_expr.size2())) {
"Analytical expression is empty or does not exist, "
"check python file");
}
return std::make_tuple(block_name, v_analytical_expr);
}
boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
boost::shared_ptr<Range> ents_ptr)
:
OP(broken_base_side_data, ents_ptr) {
this->sourceVec = vec;
this->betaCoeff = beta_coeff;
}
boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
:
OP(broken_base_side_data, ents_ptr) {
this->sourceVec = boost::shared_ptr<MatrixDouble>();
this->betaCoeff = beta_coeff;
}
if (OP::entsPtr) {
if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
}
#ifndef NDEBUG
if (!brokenBaseSideData) {
}
#endif
auto do_work_rhs = [this](int row_side, EntityType row_type,
if (!OP::nbRows)
OP::nbIntegrationPts = OP::getGaussPts().size2();
OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
OP::locF.resize(OP::nbRows, false);
OP::locF.clear();
CHKERR this->iNtegrate(row_data);
CHKERR this->aSsemble(row_data);
};
switch (OP::opType) {
case OP::OPSPACE:
for (auto &bd : *brokenBaseSideData) {
this->sourceVec =
boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
this->sourceVec.reset();
}
break;
default:
(std::string("wrong op type ") +
OpBaseDerivativesBase::OpTypeNames[OP::opType])
.c_str());
}
}
const std::string row_field,
boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
:
OP(row_field, boost::shared_ptr<MatrixDouble>(), beta_coeff, ents_ptr),
brokenBaseSideDataPtr(broken_base_side_data) {
this->betaCoeff = beta_coeff;
}
for (auto &bd : (*brokenBaseSideDataPtr)) {
this->sourceVec =
boost::shared_ptr<MatrixDouble>(brokenBaseSideDataPtr, &bd.getFlux());
#ifndef NDEBUG
"Inconsistent size of the source vector");
}
if (this->sourceVec->size1() != OP::getGaussPts().size2()) {
"Inconsistent size of the source vector");
}
#endif
this->sourceVec.reset();
}
}
std::string row_field,
boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
boost::shared_ptr<Range> ents_ptr)
:
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
ents_ptr) {
this->betaCoeff = beta;
this->sYmm = false;
}
boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
:
OP(broken_base_side_data, ents_ptr) {
this->sYmm = false;
this->betaCoeff = beta;
OP::assembleTranspose = false;
OP::onlyTranspose = false;
}
if (OP::entsPtr) {
if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
}
#ifndef NDEBUG
if (!brokenBaseSideData) {
}
#endif
auto do_work_lhs = [this](int row_side, int col_side, EntityType row_type,
EntityType col_type,
auto check_if_assemble_transpose = [&] {
if (this->sYmm) {
if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
return true;
else
return false;
} else if (OP::assembleTranspose) {
return true;
}
return false;
};
OP::rowSide = row_side;
OP::rowType = row_type;
OP::colSide = col_side;
OP::colType = col_type;
OP::locMat.resize(OP::nbRows, OP::nbCols, false);
OP::locMat.clear();
CHKERR this->iNtegrate(row_data, col_data);
CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
};
switch (OP::opType) {
case OP::OPSPACE:
for (auto &bd : *brokenBaseSideData) {
#ifndef NDEBUG
if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
"base functions not set");
}
#endif
OP::nbRows = bd.getData().getIndices().size();
if (!OP::nbRows)
OP::nbIntegrationPts = OP::getGaussPts().size2();
OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
if (!OP::nbRows)
bd.getSide(), bd.getSide(),
bd.getType(), bd.getType(),
bd.getData(), bd.getData()
);
}
break;
default:
(std::string("wrong op type ") +
OpBaseDerivativesBase::OpTypeNames[OP::opType])
.c_str());
}
}
}
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
#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_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 Tag get_tag(moab::Interface &moab, std::string tag_name, int size)
std::tuple< std::string, MatrixDouble > getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
decltype(GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
constexpr IntegrationType I
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static PetscBool l2UserBaseScale
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool physicalTimeFlg
static double currentPhysicalTime
static boost::function< double(const double)> f
static PetscBool setSingularity
static bool hasNonHomogeneousMaterialBlock
static boost::function< double(const double)> d_f
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_f
static auto diffDiffExp(A &&t_w_vee, B &&theta)
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.
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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
auto getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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.
auto 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.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseBrokenBase(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseTimesBrokenDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBrokenBaseTimesHybridDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
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)
OpHybridBaseTimesBrokenDisp(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
OpHyrbridBaseBrokenBase(std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
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 integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
