return t_diff;
};
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
if constexpr (DIM == 2) {
t_L(0, 0, 0) = 1;
t_L(1, 0, 1) = 1;
t_L(1, 1, 2) = 1;
} else {
t_L(0, 0, 0) = 1;
t_L(1, 0, 1) = 1;
t_L(2, 0, 2) = 1;
t_L(1, 1, 3) = 1;
t_L(2, 1, 4) = 1;
t_L(2, 2, 5) = 1;
}
return t_L;
}
t_diff(0, 0, 0, 0) = 1;
t_diff(1, 1, 1, 1) = 1;
t_diff(1, 0, 1, 0) = 0.5;
t_diff(1, 0, 0, 1) = 0.5;
t_diff(0, 1, 0, 1) = 0.5;
t_diff(0, 1, 1, 0) = 0.5;
if constexpr (DIM == 3) {
t_diff(2, 2, 2, 2) = 1;
t_diff(2, 0, 2, 0) = 0.5;
t_diff(2, 0, 0, 2) = 0.5;
t_diff(0, 2, 0, 2) = 0.5;
t_diff(0, 2, 2, 0) = 0.5;
t_diff(2, 1, 2, 1) = 0.5;
t_diff(2, 1, 1, 2) = 0.5;
t_diff(1, 2, 1, 2) = 0.5;
t_diff(1, 2, 2, 1) = 0.5;
}
return t_diff;
};
template <typename T>
constexpr double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1)) *
third;
};
template <typename T>
constexpr double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) *
third;
};
template <typename T, int DIM>
) {
for (int ii = 0; ii != DIM; ++ii)
for (int jj = ii; jj != DIM; ++jj)
t_dev(ii, jj) = t_stress(ii, jj);
for (int ii = 0; ii != DIM; ++ii)
for (int jj = ii; jj != DIM; ++jj)
t_dev(ii, jj) -= t_alpha(ii, jj);
return t_dev;
};
template <typename T>
};
template <typename T>
};
template <int DIM>
t_diff_deviator(
I,
J,
k,
l) = 0;
for (int ii = 0; ii != DIM; ++ii)
for (int jj = ii; jj != DIM; ++jj)
for (int kk = 0; kk != DIM; ++kk)
for (int ll = kk; ll != DIM; ++ll)
t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
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;
if constexpr (DIM == 3) {
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;
};
}
}
inline double
return std::sqrt(1.5 * t_stress_deviator(
I,
J) * t_stress_deviator(
I,
J)) +
std::numeric_limits<double>::epsilon();
};
template <int DIM>
(1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l))) / f;
return t_diff_f;
};
template <typename T, int DIM>
inline auto
t_diff_flow(
i,
j,
k,
l) =
(1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l) -
(2. / 3.) * t_flow(
i,
j) * t_flow(
k,
l))) /
f;
return t_diff_flow;
};
template <typename T, int DIM>
t_diff_flow(
i,
j,
k,
l) =
t_diff_plastic_flow_dstress(
i,
j,
m,
n) * t_D(
m,
n,
k,
l);
return t_diff_flow;
};
if (y > std::numeric_limits<float>::max_exponent10 ||
y < std::numeric_limits<float>::min_exponent10) {
if (x > 0)
return 1.;
else
return -1.;
} else {
const auto e = std::exp(y);
return (e - 1) / (1 + e);
}
};
const auto y = -x /
zeta;
if (y > std::numeric_limits<float>::max_exponent10 ||
y < std::numeric_limits<float>::min_exponent10) {
return std::abs(x);
} else {
const double e = std::exp(y);
return x + 2 *
zeta * std::log1p(e);
}
};
inline double w(
double eqiv,
double dot_tau,
double f,
double sigma_y,
double sigma_Y) {
return (f - sigma_y) / sigma_Y +
cn1 * (eqiv);
};
inline double constraint(
double eqiv,
double dot_tau,
double f,
double sigma_y,
double abs_w, double vis_H, double sigma_Y) {
return vis_H * dot_tau +
(sigma_Y / 2) * ((
cn0 * (dot_tau - eqiv) +
cn1 * (eqiv) -
(f - sigma_y) / sigma_Y) -
abs_w);
};
double vis_H, double sigma_Y) {
return vis_H + (sigma_Y / 2) * (
cn0);
};
double sigma_Y) {
return (sigma_Y / 2) * (-
cn0 +
cn1 * (1 - sign));
};
template <typename T, int DIM>
inline auto
return t_diff_constrain_dstress;
};
template <typename T1, typename T2, int DIM>
t_diff_constrain_dstrain(
k,
l) =
t_diff_constrain_dstress(
i,
j) * t_D(
i,
j,
k,
l);
return t_diff_constrain_dstrain;
};
template <typename T, int DIM>
constexpr double A = 2. / 3;
return std::sqrt(
A * t_plastic_strain_dot(
i,
j) *
t_plastic_strain_dot(
i,
j)) +
std::numeric_limits<double>::epsilon();
};
template <typename T1, typename T2, typename T3, int DIM>
T3 &t_diff_plastic_strain,
constexpr double A = 2. / 3;
t_diff_eqiv(
i,
j) =
A * (t_plastic_strain_dot(
k,
l) / eqiv) *
t_diff_plastic_strain(
k,
l,
i,
j);
return t_diff_eqiv;
};
&mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
&mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
&mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
&mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
&mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
&mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
&mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
}
template <int DIM, typename DomainEleOp>
struct OpCalculatePlasticSurfaceImpl<DIM, GAUSS,
DomainEleOp>
OpCalculatePlasticSurfaceImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp>
OpCalculatePlasticSurfaceImpl<DIM, GAUSS, DomainEleOp>::
OpCalculatePlasticSurfaceImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
}
template <int DIM, typename DomainEleOp>
MoFEMErrorCode OpCalculatePlasticSurfaceImpl<DIM, GAUSS, DomainEleOp>::doWork(
const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
auto t_stress =
getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
auto t_plastic_strain =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
commonDataPtr->plasticSurface.resize(nb_gauss_pts, false);
commonDataPtr->plasticFlow.resize(
size_symm, nb_gauss_pts,
false);
auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
auto ¶ms = commonDataPtr->blockParams;
for (auto &f : commonDataPtr->plasticSurface) {
t_stress,
trace(t_stress),
);
auto t_flow_tmp =
t_flow(
i,
j) = t_flow_tmp(
i,
j);
++t_plastic_strain;
++t_flow;
++t_stress;
}
}
template <int DIM, typename DomainEleOp>
OpCalculatePlasticityImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename DomainEleOp>
OpCalculatePlasticityImpl<DIM, GAUSS, DomainEleOp>::OpCalculatePlasticityImpl(
const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
}
template <int DIM, typename DomainEleOp>
MoFEMErrorCode OpCalculatePlasticityImpl<DIM, GAUSS, DomainEleOp>::doWork(
auto ¶ms = commonDataPtr->blockParams;
auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
auto t_plastic_strain =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
auto t_stress =
getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
t_flow_dir_dstress(
i,
j,
k,
l) =
1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l));
t_flow_dir_dstrain(
i,
j,
k,
l) =
t_flow_dir_dstress(
i,
j,
m,
n) * t_D_Op(
m,
n,
k,
l);
auto t_alpha_dir =
commonDataPtr->resC.resize(nb_gauss_pts, false);
commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
commonDataPtr->resCdStrain.resize(
size_symm, nb_gauss_pts,
false);
commonDataPtr->resCdPlasticStrain.resize(
size_symm, nb_gauss_pts,
false);
commonDataPtr->resFlow.resize(
size_symm, nb_gauss_pts,
false);
commonDataPtr->resFlowDtau.resize(
size_symm, nb_gauss_pts,
false);
false);
false);
commonDataPtr->resC.clear();
commonDataPtr->resCdTau.clear();
commonDataPtr->resCdStrain.clear();
commonDataPtr->resCdPlasticStrain.clear();
commonDataPtr->resFlow.clear();
commonDataPtr->resFlowDtau.clear();
commonDataPtr->resFlowDstrain.clear();
commonDataPtr->resFlowDstrainDot.clear();
auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
auto t_res_c_dstrain =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
auto t_res_c_plastic_strain =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
auto t_res_flow_dtau =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
auto t_res_flow_dstrain =
getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
auto t_res_flow_dplastic_strain =
getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
auto next = [&]() {
++t_tau;
++t_tau_dot;
++t_f;
++t_flow;
++t_plastic_strain;
++t_plastic_strain_dot;
++t_stress;
++t_res_c;
++t_res_c_dtau;
++t_res_c_dstrain;
++t_res_c_plastic_strain;
++t_res_flow;
++t_res_flow_dtau;
++t_res_flow_dstrain;
++t_res_flow_dplastic_strain;
++t_w;
};
auto get_avtive_pts = [&]() {
int nb_points_avtive_on_elem = 0;
int nb_points_on_elem = 0;
auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
for (auto &f : commonDataPtr->plasticSurface) {
eqiv, t_tau_dot, t_f,
++nb_points_on_elem;
if (sign_ww > 0) {
++nb_points_avtive_on_elem;
}
++t_tau;
++t_tau_dot;
++t_f;
++t_plastic_strain_dot;
}
++nb_elements;
nb_points += nb_points_on_elem;
if (nb_points_avtive_on_elem > 0) {
++avtive_elems;
active_points += nb_points_avtive_on_elem;
if (nb_points_avtive_on_elem == nb_points_on_elem) {
++avtive_full_elems;
}
}
if (nb_points_avtive_on_elem != nb_points_on_elem)
return 1;
else
return 0;
};
get_avtive_pts();
}
for (auto &f : commonDataPtr->plasticSurface) {
t_diff_plastic_strain,
const auto sigma_y =
const auto d_sigma_y =
auto c =
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
t_stress,
trace(t_stress),
);
t_flow_dir(
k,
l) = 1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l));
t_flow_dstrain(
i,
j) = t_flow(
k,
l) * t_D_Op(
k,
l,
i,
j);
auto get_res_c = [&]() {
return c; };
auto get_res_c_dstrain = [&](auto &t_diff_res) {
t_diff_res(
i,
j) = c_f * t_flow_dstrain(
i,
j);
};
auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
t_diff_res(
i,
j) = (this->getTSa() * c_equiv) * t_diff_eqiv(
i,
j);
t_diff_res(
k,
l) -= c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
};
auto get_res_c_dtau = [&]() {
return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
};
auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
t_diff_res(
k,
l) = -c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
};
auto get_res_flow = [&](auto &t_res_flow) {
const auto b = t_tau_dot;
t_res_flow(
k,
l) =
a * t_plastic_strain_dot(
k,
l) - b * t_flow_dir(
k,
l);
};
auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
const auto da = d_sigma_y;
const auto db = this->getTSa();
da * t_plastic_strain_dot(
k,
l) - db * t_flow_dir(
k,
l);
};
auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
const auto b = t_tau_dot;
t_res_flow_dstrain(
m,
n,
k,
l) = -t_flow_dir_dstrain(
m,
n,
k,
l) * b;
};
auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
t_res_flow_dplastic_strain(
m,
n,
k,
l) =
(
a * this->getTSa()) * t_diff_plastic_strain(
m,
n,
k,
l);
const auto b = t_tau_dot;
t_res_flow_dplastic_strain(
m,
n,
i,
j) +=
(t_flow_dir_dstrain(
m,
n,
k,
l) * t_alpha_dir(
k,
l,
i,
j)) * b;
};
t_res_c = get_res_c();
get_res_flow(t_res_flow);
if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
t_res_c_dtau = get_res_c_dtau();
get_res_c_dstrain(t_res_c_dstrain);
get_res_c_dplastic_strain(t_res_c_plastic_strain);
get_res_flow_dtau(t_res_flow_dtau);
get_res_flow_dstrain(t_res_flow_dstrain);
get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
}
next();
}
}
template <int DIM, typename DomainEleOp>
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> mDPtr);
private:
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp>
OpPlasticStressImpl<DIM, GAUSS, DomainEleOp>::OpPlasticStressImpl(
const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
}
template <int DIM, typename DomainEleOp>
MoFEMErrorCode
OpPlasticStressImpl<DIM, GAUSS, DomainEleOp>::doWork(
int side,
EntityType type,
const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
auto t_strain =
getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
auto t_plastic_strain =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
auto t_stress =
getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_D(
i,
j,
k,
l) * (t_strain(
k,
l) - t_plastic_strain(
k,
l));
++t_strain;
++t_plastic_strain;
++t_stress;
}
}
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::
OpCalculatePlasticFlowRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculatePlasticFlowRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const auto nb_base_functions = data.getN().size2();
auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
auto next = [&]() { ++t_res_flow; };
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_base = data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
t_rhs(
L) = alpha * (t_res_flow(
i,
j) * t_L(
i,
j,
L));
next();
auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
size_t bb = 0;
t_nf(
L) += t_base * t_rhs(
L);
++t_base;
++t_nf;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
}
}
template <typename AssemblyDomainEleOp>
OpCalculateConstraintsRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <typename AssemblyDomainEleOp>
OpCalculateConstraintsRhsImpl<GAUSS, AssemblyDomainEleOp>::
OpCalculateConstraintsRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
template <typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateConstraintsRhsImpl<GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const size_t nb_base_functions = data.getN().size2();
auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
auto next = [&]() { ++t_res_c; };
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_base = data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
const auto res = alpha * t_res_c;
next();
size_t bb = 0;
nf[bb] += t_base * res;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
}
}
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowLhs_dEPImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowLhs_dEPImpl<DIM, GAUSS, AssemblyDomainEleOp>::
OpCalculatePlasticFlowLhs_dEPImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
AssemblyDomainEleOp::OPROWCOL),
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
AssemblyDomainEleOp::sYmm = false;
}
&mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
&mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
&mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
&mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
&mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
&mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
&mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
&mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
&mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
&mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
&mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
&mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
&mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
&mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
}
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculatePlasticFlowLhs_dEPImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_res_flow_dstrain =
getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
auto t_res_flow_dplastic_strain =
getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
auto next = [&]() {
++t_res_flow_dstrain;
++t_res_flow_dplastic_strain;
};
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
alpha * (t_L(
i,
j, O) * ((t_res_flow_dplastic_strain(
i,
j,
k,
l) -
t_res_flow_dstrain(
i,
j,
k,
l)) *
next();
size_t rr = 0;
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(O,
L) += ((t_row_base * t_col_base) * t_res_mat(O,
L));
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr)
++t_row_base;
}
}
template <int DIM, typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dUImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowLhs_dTAUImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, GAUSS, AssemblyDomainEleOp>::
OpCalculatePlasticFlowLhs_dTAUImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
AssemblyDomainEleOp::sYmm = false;
}
&mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
&mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
}
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const size_t nb_row_base_functions = row_data.getN().size2();
const auto type = AssemblyDomainEleOp::getFEType();
const auto nb_nodes = moab::CN::VerticesPerEntity(type);
auto t_res_flow_dtau =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
auto next = [&]() { ++t_res_flow_dtau; };
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
t_res_vec(
L) = alpha * (t_res_flow_dtau(
i,
j) * t_L(
i,
j,
L));
next();
size_t rr = 0;
auto t_mat =
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(
L) += t_row_base * t_col_base * t_res_vec(
L);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_base;
}
}
template <int DIM, typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dEPImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> mat_D_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
template <int DIM, typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dEPImpl<DIM, GAUSS, AssemblyDomainEleOp>::
OpCalculateConstraintsLhs_dEPImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
AssemblyDomainEleOp::sYmm = false;
}
&mat(0, 0), &mat(0, 1), &mat(0, 2)};
}
&mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
}
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateConstraintsLhs_dEPImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_c_dstrain =
getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
auto t_c_dplastic_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
auto next = [&]() {
++t_c_dstrain;
++t_c_dplastic_strain;
};
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
t_L(
i,
j,
L) * (t_c_dplastic_strain(
i,
j) - t_c_dstrain(
i,
j));
next();
size_t rr = 0;
const auto row_base = alpha * t_row_base;
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(
L) += (row_base * t_col_base) * t_res_vec(
L);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_base;
}
}
template <typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dTAUImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dTAUImpl<GAUSS, AssemblyDomainEleOp>::
OpCalculateConstraintsLhs_dTAUImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
AssemblyDomainEleOp::sYmm = false;
}
template <typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateConstraintsLhs_dTAUImpl<GAUSS, AssemblyDomainEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
auto next = [&]() { ++t_res_c_dtau; };
auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
++t_w;
const auto res = alpha * (t_res_c_dtau);
next();
size_t rr = 0;
auto t_col_base = col_data.getFTensor0N(gg, 0);
*mat_ptr += t_row_base * t_col_base * res;
++t_col_base;
++mat_ptr;
}
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr)
++t_row_base;
}
}
};
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, DIM > &t_flow, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
double constrian_sign(double x)
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
auto symm_L_tensor(FTensor::Number< DIM >)
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
auto diff_constrain_dstress(double diff_constrain_df, FTensor::Tensor2_symmetric< T, DIM > &t_plastic_flow)
FTensor::Index< 'M', 3 > M
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain, FTensor::Number< DIM >)
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w, double vis_H, double sigma_Y)
auto diff_constrain_df(double sign)
FTensor::Index< 'N', 3 > N
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
[Common data]
double constrain_abs(double x)
auto diff_symmetrize(FTensor::Number< DIM >)
FTensor::Index< 'J', 3 > J
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
auto diff_deviator(FTensor::Ddg< double, DIM, DIM > &&t_diff_stress, FTensor::Number< DIM >)
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, DIM, DIM > &t_D, FTensor::Ddg< double, DIM, DIM > &&t_diff_plastic_flow_dstress)
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
auto deviator(FTensor::Tensor2_symmetric< T, DIM > &t_stress, double trace, FTensor::Tensor2_symmetric< double, DIM > &t_alpha, FTensor::Number< DIM >)
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
double zeta
Viscous hardening.
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
const TSMethod::TSContext getTSCtx() const
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
static std::array< int, 5 > activityData
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp