Implementation of operators.
#include <boost/math/constants/constants.hpp>
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
int nb_gauss_pts,
const std::string block_name);
template <typename OP_PTR>
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();
ts_time_step, ts_time, nb_gauss_pts, m_ref_coords, m_ref_normals, block_name);
#ifndef NDEBUG
if (v_analytical_expr.size2() != nb_gauss_pts)
"Wrong number of integration pts");
#endif
return std::make_tuple(block_name, v_analytical_expr);
};
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);
auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
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;
}
}
EntityType type,
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->eigenValsC.resize(3, nb_integration_pts, false);
dataAtPts->eigenVecsC.resize(9, nb_integration_pts, false);
dataAtPts->nbUniqC.resize(nb_integration_pts, false);
dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
dataAtPts->internalStressAtPts.resize(9, nb_integration_pts, false);
dataAtPts->internalStressAtPts.clear();
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 &nbUniq = dataAtPts->nbUniq;
auto t_nb_uniq =
auto t_eigen_vals_C = getFTensor1FromMat<3>(dataAtPts->eigenValsC);
auto t_eigen_vecs_C = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecsC);
auto &nbUniqC = dataAtPts->nbUniqC;
auto t_nb_uniq_C =
auto t_log_stretch_total =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
auto t_log_u2_h1 =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
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_eigen_vals_C;
++t_eigen_vecs_C;
++t_nb_uniq_C;
++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);
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);
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);
if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
}
t_nb_uniq_C = get_uniq_nb<3>(&eig(0));
if (t_nb_uniq_C < 3) {
}
}
t_eigen_vals_C(
i) = eig(
i);
t_eigen_vecs_C(
i,
j) = eigen_vec(
i,
j);
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);
t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
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) +
t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
l,
k, L) +
t_approx_P_adjoint__dstretch(
l,
k) * t_Ldiff_u(
k,
l, L) +
t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
k,
l, L)) /
4.;
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.;
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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 {
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_Ldiff_u(
k,
l, L)) / 2.;
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) =
(t_Ldiff_u(
k,
l, L) + t_Ldiff_u(
k,
l, L)) / 2) +
(t_Ldiff_u(
k,
l, L) + t_Ldiff_u(
k,
l, L)) / 2);
} else {
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) =
half_r *
(t_diff_Rr(
i,
l,
m) * (t_Ldiff_u(
l,
k, L) * t_approx_P(
i,
k)))
+
half_l *
(t_diff_Rl(
i,
l,
m) * (t_Ldiff_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_u(
l,
k))
+
half_l * t_diff_Rl(
i,
l,
m) * (t_u(
l,
k));
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);
t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
break;
t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
t_diff_diff_Rr(
i,
j,
l,
m) = 0;
t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
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)));
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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);
t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
break;
t_diff_R(
i,
j,
l) = levi_civita(
i,
j,
l);
t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
t_diff_diff_Rr(
i,
j,
l,
m) = 0;
t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
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)));
if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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(
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 =
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) {
};
auto time_step = getTStimeStep();
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_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 <>
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);
dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts, false);
dataAtPts->internalStressAtPts.clear();
auto t_internal_stress =
getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_internal_stress(L) = t_const_stress(L);
++t_internal_stress;
}
}
template <>
EntityType type,
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]);
dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts, false);
dataAtPts->internalStressAtPts.clear();
auto t_internal_stress =
getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
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 =
getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
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) * 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 =
getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
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) * 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);
}
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);
}
? 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;
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);
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
auto t_approx_P = getFTensor2FromMat<3, 3>(*piolaStressPtr);
auto t_u = getFTensor1FromMat<3>(*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_tracton(
i) = t_approx_P(
i,
j) * t_N(
j);
t_res(
i) = t_Q(
i,
j) * t_tracton(
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) * 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;
}
OP::locF *= 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) = 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>(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();
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 *= getMeasure();
};
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>(*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>(*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>(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 =
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) {
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;
}
}
&
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();
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
};
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();
auto time_step = getTStimeStep();
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 &&
} 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 = 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 &&
} 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 = 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) {
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 (dataAtPts->matInvD.size1() ==
size_symm &&
} 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 = 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;
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();
&
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 = 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();
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);
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 t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
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);
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) {
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>()));
}
jac, det, inv_jac);
}
dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(), false);
dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(), 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 = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->faceMaterialForceAtPts);
auto t_p =
getFTensor0FromVec(dataAtPts->normalPressureAtPts);
auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
auto t_grad_P = getFTensor3FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
dataAtPts->gradPAtPts);
auto t_u_gamma =
getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
auto t_grad_u_gamma =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
auto t_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
auto next = [&]() {
++t_normal;
++t_P;
++t_grad_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_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
next();
}
break;
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
(t_grad_u_gamma(
i,
j) - t_grad_u_gamma(
j,
i)) / 2. +
t_T(
I) -= t_N(
J) * (t_u_gamma(
i) * t_grad_P(
i,
J,
I)) / loop_size;
t_T(
I) -= t_N(
I) * ((t_strain(
i, K) * t_P(
i, K)) / 2.) / loop_size;
t_p += t_N(
I) * (t_N(
J) * (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 = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->faceMaterialForceAtPts);
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);
}
#ifdef ENABLE_PYTHON_BINDING
MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
try {
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
analyticalDispFun = mainNamespace["analytical_disp"];
analyticalTractionFun = mainNamespace["analytical_traction"];
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray nx, np::ndarray ny, np::ndarray nz,
const std::string &block_name, np::ndarray &analytical_expr
) {
try {
analytical_expr = bp::extract<np::ndarray>(
analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray nx, np::ndarray ny, np::ndarray nz,
const std::string& block_name, np::ndarray &analytical_expr
) {
try {
analytical_expr = bp::extract<np::ndarray>(
analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
int id) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(nb_gauss_pts);
auto stride = bp::make_tuple(3 * sizeof(double));
return (np::from_data(&data[id], dtype, size, stride, bp::object()));
};
#endif
int nb_gauss_pts,
const std::string block_name) {
#ifdef ENABLE_PYTHON_BINDING
if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
}
np::ndarray np_analytical_expr = np::empty(bp::make_tuple(nb_gauss_pts, 3),
np::dtype::get_builtin<double>());
auto disp_block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
auto traction_block_name = "(.*)ANALYTICAL_TRACTION(.*)";
std::regex reg_disp_name(disp_block_name);
std::regex reg_traction_name(traction_block_name);
if (std::regex_match(block_name, reg_disp_name)) {
delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]),
block_name, np_analytical_expr),
"Failed analytical_disp() python call");
} else if (std::regex_match(block_name, reg_traction_name)) {
delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]),
block_name, np_analytical_expr),
"Failed analytical_traction() python call");
} else {
"Unknown analytical block");
}
double *analytical_expr_val_ptr =
reinterpret_cast<double *>(np_analytical_expr.get_data());
v_analytical_expr.resize(3, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 3; ++idx)
v_analytical_expr(idx, gg) = *(analytical_expr_val_ptr + (3 * gg + idx));
}
return v_analytical_expr;
}
#endif
}
}
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_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
auto sort_eigen_vals(A &eig, B &eigen_vec)
static constexpr auto size_symm
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, MatrixDouble &m_ref_normals, const std::string block_name)
[Analytical displacement function from python]
auto getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
constexpr IntegrationType I
constexpr double t
plate stiffness
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static double dynamicTime
static PetscBool l2UserBaseScale
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static enum SymmetrySelector symmetrySelector
static boost::function< double(const double)> f
static PetscBool setSingularity
static boost::function< double(const double)> d_f
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_f
static auto diffExp(A &&t_w_vee, B &&theta)
static auto exp(A &&t_w_vee, B &&theta)
Add operators pushing bases from local to physical configuration.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
Get DOF values on entity.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Operator for inverting matrices at integration points.
Scale base functions by inverses of measure of element.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
