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,
MatrixDouble &m_ref_coords,
MatrixDouble &m_ref_normals,
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_coords = op_ptr->getCoordsAtGaussPts();
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);
nb_integration_pts, false);
dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
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);
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>(
auto t_levi_kirchhoff_dP =
getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchhoffPAtPts);
auto t_approx_P_adjoint__dstretch =
getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
auto t_approx_P_adjoint_log_du =
getFTensor1FromMat<size_symm>(
dataAtPts->adjointPdUAtPts);
auto t_approx_P_adjoint_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjointPdUdPAtPts);
auto t_approx_P_adjoint_log_du_domega =
getFTensor2FromMat<3, size_symm>(
dataAtPts->adjointPdUdOmegaAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
auto t_nb_uniq =
auto t_eigen_vals_C = getFTensor1FromMat<3>(
dataAtPts->eigenValsC);
auto t_eigen_vecs_C = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecsC);
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";
}
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";
}
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) {
}
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;
}
}
};
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);
t_traction(
k) * t_w * getMeasure();
++t_coords;
++t_spatial_disp;
++t_w;
++t_traction;
}
}
int nb_integration_pts = data.
getN().size1();
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);
}
auto t_s_dot_dot_w = getFTensor1FromMat<3>(
dataAtPts->wL2DotDotAtPts);
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;
}
}
auto t_levi_kirchhoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchhoffAtPts);
auto t_omega_grad_dot =
getFTensor2FromMat<3, 3>(
dataAtPts->rotAxisGradDotAtPts);
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(
nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
(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();
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();
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_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);
}
VectorDouble const_stress_vec(9);
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);
VectorDouble vert_data(vert_num * tag_length);
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);
}
}
}
template <AssemblyType A>
double time = OP::getFEMethod()->ts_t;
}
for (auto &bc : (*bcDispPtr)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
MatrixDouble analytical_expr;
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();
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();
};
if (bc.faces.find(fe_ent) != bc.faces.end()) {
double time_scale = 1;
}
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();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
} else {
(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.;
};
if (bc.faces.find(fe_ent) != bc.faces.end()) {
double time_scale = 1;
}
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_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
t_kd(
j,
l) * t_diff_col_base(N0)
+
(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);
}
}
}
}
}
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
if (bc.faces.find(fe_ent) != bc.faces.end()) {
MatrixDouble analytical_expr;
auto [block_name, v_analytical_expr] =
auto t_val =
getFTensor1FromMat<3>(v_analytical_expr);
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
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();
auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
};
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();
auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
);
};
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())
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();
auto get_ftensor3 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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));
};
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();
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
};
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;
}
}
auto get_ftensor3 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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_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;
}
}
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
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;
}
}
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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 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;
}
}
auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
};
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;
}
}
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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_levi_kirchhoff_domega =
getFTensor2FromMat<3, 3>(
dataAtPts->leviKirchhoffdOmegaAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
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;
}
}
} else {
}
};
template <int S>
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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 row_nb_base_functions = row_data.
getN().size2() / 3;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
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;
}
}
} else {
}
};
template <int S>
int row_nb_base_functions = row_data.
getN().size2() / 9;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
for (int cc = 0; cc != col_nb_dofs; ++cc) {
a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
++t_col_base;
}
++t_row_base;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base;
++t_w;
++t_inv_D;
}
}
} else {
}
};
template <int S>
auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
);
};
int row_nb_base_functions = row_data.
getN().size2() / 9;
auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
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();
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
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 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();
auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
&
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
};
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;
}
}
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};
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;
};
dataAtPts->energyAtPts.resize(getGaussPts().size2());
}
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);
};
auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
a = set_float_precision(
a);
};
&
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);
};
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);
next();
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::vector<FieldSpace> spaces, std::string geom_field_name,
boost::shared_ptr<Range> crack_front_edges_ptr) {
constexpr bool scale_l2 = false;
if (scale_l2) {
boost::shared_ptr<Range> edges_ptr)
: 0;
if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
return 0;
} else {
return OP::doWork(side, type, data);
}
};
private:
boost::shared_ptr<Range> edgesPtr;
};
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
? crack_front_edges_ptr
: boost::make_shared<Range>()));
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::vector<FieldSpace> spaces, std::string geom_field_name,
boost::shared_ptr<Range> crack_front_edges_ptr) {
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) {
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
geom_field_name, jac));
pipeline.push_back(
}
constexpr bool scale_l2_ainsworth_legendre_base = false;
if (scale_l2_ainsworth_legendre_base) {
boost::shared_ptr<MatrixDouble> jac,
boost::shared_ptr<Range> edges_ptr)
: 0;
if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
return 0;
} else {
return OP::doWork(side, type, data);
}
};
private:
boost::shared_ptr<Range> edgesPtr;
};
auto jac = boost::make_shared<MatrixDouble>();
auto det = boost::make_shared<VectorDouble>();
geom_field_name, jac,
: boost::make_shared<Range>()));
}
nullptr, nullptr, nullptr);
}
dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
if (getNinTheLoop() == 0) {
}
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>(
auto t_p =
getFTensor0FromVec(
dataAtPts->normalPressureAtPts);
auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->approxPAtPts);
auto t_grad_u_gamma =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->gradHybridDispAtPts);
auto t_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(
dataAtPts->logStretchTensorAtPts);
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
(t_grad_u_gamma(
i,
j) - t_grad_u_gamma(
j,
i)) / 2. + 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;
++t_normal;
++t_P;
++t_grad_u_gamma;
++t_strain;
++t_T;
++t_p;
}
break;
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_T(
I) += t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J)) / loop_size;
(t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
loop_size;
++t_normal;
++t_P;
++t_grad_u_gamma;
++t_T;
++t_p;
}
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>(
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,
MatrixDouble &m_ref_coords,
MatrixDouble &m_ref_normals,
const std::string block_name) {
#ifdef ENABLE_PYTHON_BINDING
if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
VectorDouble v_ref_coords = m_ref_coords.data();
VectorDouble v_ref_normals = m_ref_normals.data();
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());
MatrixDouble v_analytical_expr;
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
return MatrixDouble();
}
}
}
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 nme:nme847.
@ 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 get_uniq_nb(double *ptr)
auto sort_eigen_vals(A &eig, B &eigen_vec)
static constexpr auto size_symm
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, MatrixDouble &m_ref_normals, const std::string block_name)
[Analytical displacement function from python]
auto getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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 double dynamicTime
static enum StretchSelector stretchSelector
static PetscBool setSingularity
static enum SymmetrySelector symmetrySelector
static enum RotSelector rotSelector
static PetscBool l2UserBaseScale
static enum EnergyReleaseSelector energyReleaseSelector
static enum RotSelector gradApproximator
static boost::function< double(const double)> f
static boost::function< double(const double)> d_f
static PetscBool dynamicRelaxation
static boost::function< double(const double)> inv_f
MoFEMErrorCode iNtegrate(EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< double > piolaScalePtr
boost::shared_ptr< AnalyticalTractionBcVec > bcData
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< double > piolaScalePtr
MoFEMErrorCode iNtegrate(EntData &data)
boost::shared_ptr< PressureBcVec > bcData
boost::shared_ptr< MatrixDouble > hybridGradDispPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< double > piolaScalePtr
boost::shared_ptr< TractionBcVec > bcData
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
std::array< double, 6 > & reactionVec
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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)
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
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getTStimeStep() const
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Scale base functions by inverses of measure of element.
double getVolume() const
element volume (linear geometry)
