v0.14.0
EshelbianOperators.cpp

Implementation of operators

/**
* \file EshelbianOperators.cpp
* \example EshelbianOperators.cpp
*
* \brief Implementation of operators
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <boost/math/constants/constants.hpp>
#include <lapack_wrap.h>
inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
FTensor::Ddg<double, 3, 3> t_diff_deviator;
t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
constexpr double third = boost::math::constants::third<double>();
t_diff_deviator(0, 0, 0, 0) -= third;
t_diff_deviator(0, 0, 1, 1) -= third;
t_diff_deviator(1, 1, 0, 0) -= third;
t_diff_deviator(1, 1, 1, 1) -= third;
t_diff_deviator(2, 2, 0, 0) -= third;
t_diff_deviator(2, 2, 1, 1) -= third;
t_diff_deviator(0, 0, 2, 2) -= third;
t_diff_deviator(1, 1, 2, 2) -= third;
t_diff_deviator(2, 2, 2, 2) -= third;
return t_diff_deviator;
}
constexpr static auto size_symm = (3 * (3 + 1)) / 2;
inline auto diff_tensor() {
t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
return t_diff;
};
inline auto symm_L_tensor() {
t_L(i, j, L) = 0;
t_L(0, 0, 0) = 1;
t_L(1, 1, 3) = 1;
t_L(2, 2, 5) = 1;
t_L(1, 0, 1) = 1;
t_L(2, 0, 2) = 1;
t_L(2, 1, 4) = 1;
return t_L;
}
template <class T> struct TensorTypeExtractor {
typedef typename std::remove_pointer<T>::type Type;
};
template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
typedef typename std::remove_pointer<T>::type Type;
};
template <typename T>
RotSelector rotSelector = LARGE_ROT) {
using D = typename TensorTypeExtractor<T>::Type;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
t_R(i, j) = t_kd(i, j);
t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
if (rotSelector == SMALL_ROT) {
t_R(i, j) += t_Omega(i, j);
return t_R;
}
const auto angle =
sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
const auto a = sin(angle) / angle;
t_R(i, j) += a * t_Omega(i, j);
if (rotSelector == MODERATE_ROT)
return t_R;
const auto ss_2 = sin(angle / 2.) / angle;
const auto b = 2. * ss_2 * ss_2;
t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
return t_R;
}
template <typename T>
RotSelector rotSelector = LARGE_ROT) {
using D = typename TensorTypeExtractor<T>::Type;
const auto angle =
sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
if (rotSelector == SMALL_ROT) {
t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
return t_diff_R;
}
const auto ss = sin(angle);
const auto a = ss / angle;
t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
const auto angle2 = angle * angle;
const auto cc = cos(angle);
const auto diff_a = (angle * cc - ss) / angle2;
t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
if (rotSelector == MODERATE_ROT)
return t_diff_R;
const auto ss_2 = sin(angle / 2.);
const auto cc_2 = cos(angle / 2.);
const auto b = 2. * ss_2 * ss_2 / angle2;
t_diff_R(i, j, k) +=
b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
const auto diff_b =
(2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
t_diff_R(i, j, k) +=
diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
return t_diff_R;
}
template <typename T>
RotSelector rotSelector = LARGE_ROT) {
using D = typename TensorTypeExtractor<T>::Type;
constexpr double eps = 1e-10;
for (int l = 0; l != 3; ++l) {
t_omega_c(i) = t_omega(i);
t_omega_c(l) += std::complex<double>(0, eps);
auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
for (int i = 0; i != 3; ++i) {
for (int j = 0; j != 3; ++j) {
for (int k = 0; k != 3; ++k) {
t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
}
}
}
}
return t_diff2_R;
}
struct isEq {
static inline auto check(const double &a, const double &b) {
constexpr double eps = std::numeric_limits<float>::epsilon();
return std::abs(a - b) / absMax < eps;
}
static double absMax;
};
double isEq::absMax = 1;
inline auto is_eq(const double &a, const double &b) {
return isEq::check(a, b);
};
template <int DIM> inline auto get_uniq_nb(double *ptr) {
std::array<double, DIM> tmp;
std::copy(ptr, &ptr[DIM], tmp.begin());
std::sort(tmp.begin(), tmp.end());
isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
constexpr double eps = std::numeric_limits<float>::epsilon();
isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
}
template <int DIM>
isEq::absMax =
std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
int i = 0, j = 1, k = 2;
if (is_eq(eig(0), eig(1))) {
i = 0;
j = 2;
k = 1;
} else if (is_eq(eig(0), eig(2))) {
i = 0;
j = 1;
k = 2;
} else if (is_eq(eig(1), eig(2))) {
i = 1;
j = 0;
k = 2;
}
eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
{
eig(i) = eig_c(i);
eigen_vec(i, j) = eigen_vec_c(i, j);
}
}
MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
EntityType type,
EntData &data) {
auto t_L = symm_L_tensor();
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->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
dataAtPts->leviKirchhoffOmegaAtPts.resize(9, 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->diffRotMatAtPts.resize(27, 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->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
auto t_h_dlog_u =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
auto t_levi_kirchoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
auto t_levi_kirchoff_domega =
getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
auto t_approx_P_adjont_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_approx_P_adjont_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
auto t_approx_P_adjont_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
auto t_approx_P_adjont_log_du_domega =
getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
auto t_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
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_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
auto t_log_stretch_total =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
auto calulate_rotation = [&]() {
auto t0_diff =
get_diff_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
auto t0 = get_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
t_diff_R(i, j, k) = t0_diff(i, j, k);
t_R(i, j) = t0(i, j);
};
auto calulate_streach = [&]() {
eigen_vec(i, j) = t_log_u(i, j);
// rare case when two eigen values are equal
nbUniq[gg] = get_uniq_nb<3>(&eig(0));
if (nbUniq[gg] < 3) {
sort_eigen_vals<3>(eig, eigen_vec);
}
t_eigen_vals(i) = eig(i);
t_eigen_vecs(i, j) = eigen_vec(i, j);
auto t_u_tmp =
EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
t_u(i, j) = t_u_tmp(i, j);
auto t_diff_u_tmp =
EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f,
EshelbianCore::d_f, nbUniq[gg]);
t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
};
calulate_rotation();
calulate_streach();
switch (EshelbianCore::gradApperoximator) {
case LARGE_ROT:
// rotation can be large, moderate or small
t_Ru(i, m) = t_R(i, l) * t_u(l, m);
t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
// conservation of angular momentum at current configuration
t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
t_approx_P_adjont_dstretch(l, m) =
t_approx_P_intermidiate(i, m) * t_R(i, l);
t_levi_kirchoff(k) =
levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
t_levi_kirchoff_dP(i, j, k) =
(levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
t_levi_kirchoff_domega(k, n) =
(t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
t_approx_P_adjont_log_du(L) =
t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
t_approx_P_adjont_log_du_domega(n, L) =
t_Ldiff_u(l, m, L) *
(t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
break;
{
// assume small current rotation
t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
// conservation of angular momentum at intermediate configuration
t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
t_levi_kirchoff(k) =
levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
t_levi_kirchoff_domega(k, n) = 0;
t_approx_P_adjont_log_du(L) =
t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
t_approx_P_adjont_log_du_domega(n, L) = 0;
}
break;
case SMALL_ROT:
{
// assume small current rotation
t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
t_h(i, j) = t_Omega(i, j) + t_u(i, j);
t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
// conservation of angular momentum at reference condition
t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
t_levi_kirchoff_domega(k, l) = 0;
t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
t_approx_P_adjont_log_du(L) =
t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
t_approx_P_adjont_log_du_domega(m, L) = 0;
}
break;
}
t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
eigen_vec(i, j) = t_C_h1(i, j);
switch (EshelbianCore::stretchSelector) {
break;
for (int ii = 0; ii != 3; ++ii) {
eig(ii) = std::max(eig(ii),
2 * EshelbianCore::inv_f(EshelbianCore::f(
std::numeric_limits<double>::min_exponent)));
}
break;
}
auto t_log_u2_h1_tmp =
EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
switch (EshelbianCore::gradApperoximator) {
case LARGE_ROT:
t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
break;
case SMALL_ROT:
t_log_stretch_total(i, j) = t_log_u(i, j);
break;
};
++t_h;
++t_h_domega;
++t_h_dlog_u;
++t_levi_kirchoff;
++t_levi_kirchoff_dP;
++t_levi_kirchoff_domega;
++t_approx_P_adjont_dstretch;
++t_approx_P_adjont_log_du;
++t_approx_P_adjont_log_du_dP;
++t_approx_P_adjont_log_du_domega;
++t_approx_P;
++t_R;
++t_diff_R;
++t_log_u;
++t_u;
++t_diff_u;
++t_eigen_vals;
++t_eigen_vecs;
++t_omega;
++t_grad_h1;
++t_log_stretch_total;
}
}
MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
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);
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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 * alphaW * t_s_dot_w(i);
t_nf(i) -= a * t_row_base_fun * alphaRho * 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;
}
}
MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = getGaussPts().size2();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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_kirchoff(k);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_levi_kirchoff;
}
}
MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
CHKERR integratePiola(data);
} else {
if (polyConvex) {
CHKERR integratePolyconvexHencky(data);
} else {
CHKERR integrateHencky(data);
}
}
}
MoFEMErrorCode OpSpatialPhysical::integrateHencky(EntData &data) {
auto t_L = symm_L_tensor();
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
auto t_log_streach_h1 =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
auto get_ftensor2 = [](auto &v) {
&v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
};
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor2(nF);
t_T(i, j) =
t_D(i, j, k, l) * (t_log_streach_h1(k, l) + alphaU * t_dot_log_u(k, l));
t_residual(L) =
a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(L) += 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_approx_P_adjont_log_du;
++t_dot_log_u;
++t_log_streach_h1;
}
}
MoFEMErrorCode OpSpatialPhysical::integratePolyconvexHencky(EntData &data) {
auto t_L = symm_L_tensor();
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
auto t_log_streach_h1 =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
auto get_ftensor2 = [](auto &v) {
&v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
};
constexpr double nohat_k = 1. / 4;
constexpr double hat_k = 1. / 8;
double mu = dataAtPts->mu;
double lambda = dataAtPts->lambda;
constexpr double third = boost::math::constants::third<double>();
auto t_diff_deviator = diff_deviator(diff_tensor());
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor2(nF);
double log_det = t_log_streach_h1(i, i);
double log_det2 = log_det * log_det;
t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
double dev_norm2 = t_dev(i, j) * t_dev(i, j);
auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
auto B = lambda * std::exp(hat_k * log_det2) * log_det;
t_T(i, j) =
A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
+
B * t_kd(i, j)
+
alphaU * t_D(i, j, k, l) * t_log_streach_h1(k, l);
t_residual(L) =
a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(L) += 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_approx_P_adjont_log_du;
++t_dot_log_u;
++t_log_streach_h1;
}
}
MoFEMErrorCode OpSpatialPhysical::integratePiola(EntData &data) {
auto t_L = symm_L_tensor();
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du =
getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
auto get_ftensor2 = [](auto &v) {
&v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
};
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor2(nF);
t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
t_viscous_P(i, j) =
alphaU *
t_dot_log_u(i,
j); // That is chit, should be split on axiator and deviator
t_residual(L) =
a * (t_approx_P_adjont_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
t_L(i, j, L) * t_viscous_P(i, j));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(L) += 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_approx_P_adjont_log_du;
++t_P;
++t_dot_log_u;
++t_diff_u;
}
}
MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor1(nF);
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
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;
++t_omega_dot;
}
}
MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
int nb_base_functions = data.getN().size2() / 9;
auto t_row_base_fun = data.getFTensor2N<3, 3>();
auto get_ftensor0 = [](auto &v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor0(nF);
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
t_residuum(i, j) = t_h(i, j) - t_kd(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;
++t_omega_dot;
}
}
MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
int nb_base_functions = data.getN().size2() / 3;
auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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;
}
}
MoFEMErrorCode OpDispBc::integrate(EntData &data) {
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcDispPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
double scale = 1;
for (auto &sm : scalingMethodsVec) {
scale *= sm->getScale(getFEMethod()->ts_t);
}
// get bc data
FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
t_bc_disp(i) *= scale;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_bc_res(i) = t_bc_disp(i);
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(i) -=
t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
}
}
}
}
MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcRotPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
// get bc data
FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
double theta = bc.theta;
theta *= getFEMethod()->ts_t;
const double a = sqrt(t_normal(i) * t_normal(i));
t_omega(i) = t_normal(i) * (theta / a);
auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
auto t_coords = 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 = get_ftensor1(nF);
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_coords;
}
}
}
}
MoFEMErrorCode FeTractionBc::preProcess() {
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const {
return 2 * (p_data + 1);
}
};
if (ts_ctx == CTX_TSSETIFUNCTION) {
// Loop boundary elements with traction boundary conditions
{HDIV});
fe.getOpPtrVector().push_back(
new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
fe.ts_t = ts_t;
CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
fe, nullptr, MF_ZERO,
this->getCacheWeakPtr());
}
}
MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
auto t_normal = getFTensor1Normal();
const double nrm2 = sqrt(t_normal(i) * t_normal(i));
t_unit_normal(i) = t_normal(i) / nrm2;
int nb_dofs = data.getFieldData().size();
int nb_integration_pts = data.getN().size1();
int nb_base_functions = data.getN().size2() / 3;
double ts_t = getFEMethod()->ts_t;
auto integrate_matrix = [&]() {
auto t_w = getFTensor0IntegrationWeight();
matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
matPP.clear();
auto t_row_base_fun = data.getFTensor1N<3>();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double a = t_row_base_fun(i) * t_unit_normal(i);
auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
for (int cc = 0; cc != nb_dofs / 3; ++cc) {
const double b = t_col_base_fun(i) * t_unit_normal(i);
matPP(rr, cc) += t_w * a * b;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
};
auto integrate_rhs = [&](auto &bc) {
auto t_w = getFTensor0IntegrationWeight();
vecPv.resize(3, nb_dofs / 3, false);
vecPv.clear();
auto t_row_base_fun = data.getFTensor1N<3>();
double ts_t = getFEMethod()->ts_t;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
vecPv(dd, rr) += t * bc.vals[dd];
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
};
auto integrate_rhs_cook = [&](auto &bc) {
vecPv.resize(3, nb_dofs / 3, false);
vecPv.clear();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base_fun = data.getFTensor1N<3>();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto calc_tau = [](double y) {
y -= 44;
y /= (60 - 44);
return -y * (y - 1) / 0.25;
};
const double tau = calc_tau(t_coords(1));
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
vecPv(dd, rr) += t * tau * bc.vals[dd];
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_coords;
}
};
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
for (auto &bc : *(bcFe.bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getFieldData().size();
if (nb_dofs) {
CHKERR integrate_matrix();
if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
CHKERR integrate_rhs_cook(bc);
else
CHKERR integrate_rhs(bc);
const auto info =
lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
if (info > 0)
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"The leading minor of order %i is "
"not positive; definite;\nthe "
"solution could not be computed",
info);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
for (int rr = 0; rr != nb_dofs / 3; ++rr)
data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
}
}
}
}
MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
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;
}
}
MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
EntData &col_data) {
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
double ts_scale = alphaW * getTSa();
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
ts_scale += alphaRho * 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_col_base_fun = row_data.getFTensor0N(gg, 0);
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(i) -= b;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
auto t_m = get_ftensor3(K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(L, i) +=
a * (t_approx_P_adjont_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_adjont_log_du_dP;
}
}
MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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),
&m(r + 5, c));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du_dP =
getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor2(K, 6 * rr, 0);
auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_m(L) +=
a * (t_approx_P_adjont_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_adjont_log_du_dP;
}
}
MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
EntData &col_data) {
if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
CHKERR integratePiola(row_data, col_data);
} else {
if (polyConvex) {
CHKERR integratePolyconvexHencky(row_data, col_data);
} else {
CHKERR integrateHencky(row_data, col_data);
}
}
}
MoFEMErrorCode OpSpatialPhysical_du_du::integrateHencky(EntData &row_data,
EntData &col_data) {
auto t_L = symm_L_tensor();
auto t_diff = diff_tensor();
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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),
&m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
&m(r + 3, c + 4), &m(r + 3, c + 5),
&m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
&m(r + 4, c + 4), &m(r + 4, c + 5),
&m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
&m(r + 5, c + 4), &m(r + 5, c + 5)
);
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
auto &nbUniq = dataAtPts->nbUniq;
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto get_dP = [&]() {
dP.resize(size_symm * size_symm, nb_integration_pts, false);
auto ts_a = getTSa();
auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
(t_L(i, j, L) *
((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
if (EshelbianCore::stretchSelector > LINEAR) {
auto t_approx_P_adjont_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
auto &nbUniq = dataAtPts->nbUniq;
auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
// Work of symmetric tensor on undefined tensor is equal to the work
// of the symmetric part of it
t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
t_approx_P_adjont_dstretch(j, i));
t_sym(i, j) /= 2.0;
auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
EshelbianCore::dd_f, t_sym, nbUniq[gg]);
t_dP(L, J) = t_L(i, j, L) *
((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
t_L(k, l, J)) /
2. +
t_dP_tmp(L, J);
++t_dP;
++t_approx_P_adjont_dstretch;
++t_eigen_vals;
++t_eigen_vecs;
}
} else {
auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
t_dP(L, J) = t_dP_tmp(L, J);
++t_dP;
}
}
return getFTensor2FromMat<size_symm, size_symm>(dP);
};
auto t_dP = get_dP();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
auto t_m = get_ftensor2(K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
const double b = a * t_row_base_fun * t_col_base_fun;
t_m(L, J) += b * t_dP(L, 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_dP;
}
}
OpSpatialPhysical_du_du::integratePolyconvexHencky(EntData &row_data,
EntData &col_data) {
auto t_L = symm_L_tensor();
auto t_diff = diff_tensor();
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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),
&m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
&m(r + 3, c + 4), &m(r + 3, c + 5),
&m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
&m(r + 4, c + 4), &m(r + 4, c + 5),
&m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
&m(r + 5, c + 4), &m(r + 5, c + 5)
);
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto get_dP = [&]() {
dP.resize(size_symm * size_symm, nb_integration_pts, false);
auto ts_a = getTSa();
auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
constexpr double nohat_k = 1. / 4;
constexpr double hat_k = 1. / 8;
double mu = dataAtPts->mu;
double lambda = dataAtPts->lambda;
constexpr double third = boost::math::constants::third<double>();
auto t_diff_deviator = diff_deviator(diff_tensor());
auto t_approx_P_adjont_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_log_streach_h1 =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
auto &nbUniq = dataAtPts->nbUniq;
auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
double log_det = t_log_streach_h1(i, i);
double log_det2 = log_det * log_det;
t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
double dev_norm2 = t_dev(i, j) * t_dev(i, j);
auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
auto B = lambda * std::exp(hat_k * log_det2) * log_det;
t_A_diff(i, j) =
(A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
t_dT(i, j, k, l) =
t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
+
A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
+
t_B_diff(i, j) * t_kd(k, l);
t_dP(L, J) = -t_L(i, j, L) *
((
t_dT(i, j, k, l)
+
(alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
)) *
t_L(k, l, J));
// Work of symmetric tensor on undefined tensor is equal to the work
// of the symmetric part of it
if (EshelbianCore::stretchSelector > LINEAR) {
t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
t_approx_P_adjont_dstretch(j, i));
t_sym(i, j) /= 2.0;
auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
EshelbianCore::dd_f, t_sym, nbUniq[gg]);
t_dP(L, J) += t_L(i, j, L) *
((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
t_L(k, l, J)) /
2.;
}
++t_dP;
++t_approx_P_adjont_dstretch;
++t_log_streach_h1;
++t_eigen_vals;
++t_eigen_vecs;
}
return getFTensor2FromMat<size_symm, size_symm>(dP);
};
auto t_dP = get_dP();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
auto t_m = get_ftensor2(K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
const double b = a * t_row_base_fun * t_col_base_fun;
t_m(L, J) += b * t_dP(L, 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_dP;
}
}
MoFEMErrorCode OpSpatialPhysical_du_du::integratePiola(EntData &row_data,
EntData &col_data) {
auto t_L = symm_L_tensor();
auto t_diff = diff_tensor();
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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),
&m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
&m(r + 3, c + 4), &m(r + 3, c + 5),
&m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
&m(r + 4, c + 4), &m(r + 4, c + 5),
&m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
&m(r + 5, c + 4), &m(r + 5, c + 5)
);
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto get_dP = [&]() {
dP.resize(size_symm * size_symm, nb_integration_pts, false);
auto ts_a = getTSa();
auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
auto t_approx_P_adjont_dstretch =
getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
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_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
t_dP(L, J) =
-t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
// viscous stress derivative
t_dP(L, J) -= (alphaU * ts_a) *
(t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
if (EshelbianCore::stretchSelector > LINEAR) {
t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
t_deltaP_sym(i, j) /= 2.0;
auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
}
++t_P;
++t_dP;
++r_P_du;
++t_approx_P_adjont_dstretch;
++t_dot_log_u;
++t_u;
++t_diff_u;
++t_eigen_vals;
++t_eigen_vecs;
}
return getFTensor2FromMat<size_symm, size_symm>(dP);
};
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto t_dP = get_dP();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
auto t_m = get_ftensor2(K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
const double b = a * t_row_base_fun * t_col_base_fun;
t_m(L, J) += b * t_dP(L, 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_dP;
}
}
MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
EntData &col_data) {
auto t_L = symm_L_tensor();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P_adjont_log_du_domega =
getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
int nb_integration_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
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_adjont_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_adjont_log_du_domega;
}
}
MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = getGaussPts().size2();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
double b = a * t_row_base_fun;
auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
auto t_m = get_ftensor2(K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(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_kirchoff_dP;
}
}
MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = getGaussPts().size2();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff_dP =
getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
double b = a * t_row_base_fun;
auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
auto t_m = get_ftensor1(K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(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_levi_kirchoff_dP;
}
}
MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = getGaussPts().size2();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_levi_kirchoff_domega =
getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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;
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_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_levi_kirchoff_domega;
}
}
MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
int row_nb_base_functions = row_data.getN().size2() / 3;
auto t_row_base_fun = row_data.getFTensor1N<3>();
const double ts_a = getTSa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
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_col_base_fun = col_data.getFTensor0N(gg, 0);
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;
}
}
OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
int row_nb_base_functions = row_data.getN().size2() / 9;
auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
const double ts_a = getTSa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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_col_base_fun = col_data.getFTensor0N(gg, 0);
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;
}
}
MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
EntData &data) {
auto create_tag = [this](const std::string tag_name, const int size) {
double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
Tag th;
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_w = create_tag("SpatialDisplacement", 3);
Tag th_omega = create_tag("Omega", 3);
Tag th_approxP = create_tag("Piola1Stress", 9);
Tag th_sigma = create_tag("CauchyStress", 9);
Tag th_log_u = create_tag("LogSpatialStretch", 9);
Tag th_u = create_tag("SpatialStretch", 9);
Tag th_h = create_tag("h", 9);
Tag th_X = create_tag("X", 3);
Tag th_detF = create_tag("detF", 1);
Tag th_angular_momentum = create_tag("AngularMomentum", 3);
Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
Tag th_traction = create_tag("traction", 3);
Tag th_disp = create_tag("H1Displacement", 3);
Tag th_disp_error = create_tag("DisplacementError", 1);
Tag th_lambda_disp = create_tag("ContactDisplacement", 3);
auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
auto t_log_u =
getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
auto set_float_precision = [](const double x) {
if (std::abs(x) < std::numeric_limits<float>::epsilon())
return 0.;
else
return x;
};
// scalars
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);
};
// vectors
auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
t_v(i) = t_d(i);
for (auto &a : v.data())
a = set_float_precision(a);
CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
&*v.data().begin());
};
// tensors
&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) {
t_m(i, j) = t_d(i, j);
for (auto &v : m.data())
v = set_float_precision(v);
CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
&*m.data().begin());
};
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();
// vectors
CHKERR save_vec_tag(th_w, t_w, gg);
CHKERR save_vec_tag(th_X, t_coords, gg);
CHKERR save_vec_tag(th_omega, t_omega, gg);
CHKERR save_vec_tag(th_traction, t_traction, gg);
// tensors
CHKERR save_mat_tag(th_h, t_h, gg);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_u_tmp(ii, jj) = t_u(ii, jj);
CHKERR save_mat_tag(th_u, t_u_tmp, gg);
CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
CHKERR save_vec_tag(th_disp, t_disp, gg);
CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, 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);
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_sigma, t_cauchy, gg);
CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
t_levi(k) = t_levi_kirchoff(k);
CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
&t_levi(0));
auto get_eiegn_vector_symmetric = [&](auto &t_u) {
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_m(ii, jj) = t_u(ii, jj);
VectorDouble3 eigen_values(3);
auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
CHKERR computeEigenValuesSymmetric(t_m, t_eigen_values);
CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
&*m.data().begin());
CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
&*eigen_values.data().begin());
};
CHKERR get_eiegn_vector_symmetric(t_u);
++t_w;
++t_h;
++t_log_u;
++t_u;
++t_omega;
++t_R;
++t_approx_P;
++t_levi_kirchoff;
++t_coords;
++t_normal;
++t_disp;
++t_lambda_disp;
}
}
OpCalculateStrainEnergy::doWork(int side, EntityType type,
if (type == MBTET) {
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = t_w * v;
// FIXME: this is wrong, energy should be calculated in material
(*energy) += a * t_P(i, J) * t_h(i, J);
++t_w;
++t_P;
++t_h;
}
}
}
} // namespace EshelbianPlasticity
EshelbianPlasticity::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: EshelbianOperators.cpp:218
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
EshelbianPlasticity::is_eq
auto is_eq(const double &a, const double &b)
Definition: EshelbianOperators.cpp:203
FTensor
JSON compatible output.
Definition: Christof_constructor.hpp:6
EshelbianPlasticity::get_rotation_form_vector
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:83
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::size_symm
constexpr static auto size_symm
Definition: EshelbianOperators.cpp:47
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::ForcesAndSourcesCore::getRuleHook
RuleHookFun getRuleHook
Hook to get rule.
Definition: ForcesAndSourcesCore.hpp:42
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::diff_tensor
auto diff_tensor()
Definition: EshelbianOperators.cpp:49
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FTensor::levi_civita
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
Definition: Levi_Civita.hpp:617
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
EshelbianPlasticity::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: EshelbianOperators.cpp:207
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
convert.type
type
Definition: convert.py:64
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
lapack_wrap.h
MatrixFunction.hpp
MoFEM::EntitiesFieldData::EntData::getFTensor2N
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:762
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
EshelbianPlasticity::diff_deviator
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
Definition: EshelbianOperators.cpp:20
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
EshelbianPlasticity::symm_L_tensor
auto symm_L_tensor()
Definition: EshelbianOperators.cpp:60
FTensor::Index< 'i', 3 >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:162
FTensor::Dg
Definition: Dg_value.hpp:9
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::Tensor0
Definition: Tensor0.hpp:16
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:21
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
EshelbianPlasticity::get_diff2_rotation_form_vector
auto get_diff2_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:167
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:20
FTensor::Ddg< double, 3, 3 >
lapack_dposv
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:211
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::getDiffDiffMat
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
Definition: MatrixFunction.cpp:78
EshelbianPlasticity::get_diff_rotation_form_vector
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:119
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1430