v0.14.0
Loading...
Searching...
No Matches
PlasticOps.hpp

\[ \left\{ \begin{array}{ll} \frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\ \varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)\\ \sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\ \dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial \sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\ f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\ u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\ \sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\ \Omega_u \cup \Omega_\sigma = \Omega \\ \Omega_u \cap \Omega_\sigma = \emptyset \end{array} \right. \]

\[ \left\{ \begin{array}{ll} \left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall \delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left( \dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0 & \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\ \left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) + \| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 & \forall \delta\tau \in L^2(\Omega) \end{array} \right. \]

/** \file PlasticOps.hpp
* \example PlasticOps.hpp
\f[
\left\{
\begin{array}{ll}
\frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i} \right)\\
\sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
\dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
\sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
\sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
\Omega_u \cup \Omega_\sigma = \Omega \\
\Omega_u \cap \Omega_\sigma = \emptyset
\end{array}
\right.
\f]
\f[
\left\{
\begin{array}{ll}
\left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
\delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
\dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
& \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
\left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
(f(\pmb\sigma,\tau) - \sigma_y) +
\| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
\forall \delta\tau \in L^2(\Omega) \end{array} \right.
\f]
*/
namespace PlasticOps {
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
H,
};
using BlockParams = std::array<double, LAST_PARAM>;
inline auto getParamsPtr() {
return boost::shared_ptr<BlockParams>(shared_from_this(), &blockParams);
};
//! [Common data set externally]
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
//! [Common data set externally]
VectorDouble plasticSurface;
MatrixDouble plasticFlow;
VectorDouble plasticTau;
VectorDouble plasticTauDot;
MatrixDouble plasticStrain;
MatrixDouble plasticStrainDot;
VectorDouble resC;
VectorDouble resCdTau;
MatrixDouble resCdStrain;
MatrixDouble resCdPlasticStrain;
MatrixDouble resFlow;
MatrixDouble resFlowDtau;
MatrixDouble resFlowDstrain;
MatrixDouble resFlowDstrainDot;
inline auto getPlasticSurfacePtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticSurface);
}
inline auto getPlasticTauPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTau);
}
inline auto getPlasticTauDotPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTauDot);
}
inline auto getPlasticStrainPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticStrain);
}
inline auto getPlasticStrainDotPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto getPlasticFlowPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticFlow);
}
static std::array<int, 5> activityData;
};
std::array<int, 5> CommonData::activityData = {0, 0, 0, 0, 0};
//! [Common data]
FTensor::Index<'I', 3> I;
FTensor::Index<'J', 3> J;
FTensor::Index<'M', 3> M;
FTensor::Index<'N', 3> N;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticSurfaceImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticityImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpPlasticStressImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowRhsImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dTAUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dEPImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dTAUImpl;
template <typename DomainEleOp> struct PlasticityIntegrators {
template <int DIM, IntegrationType I>
OpCalculatePlasticSurfaceImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticity = OpCalculatePlasticityImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpPlasticStress = OpPlasticStressImpl<DIM, I, DomainEleOp>;
template <AssemblyType A> struct Assembly {
typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowRhsImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsRhsImpl<I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl<
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsLhs_dTAUImpl<I, AssemblyDomainEleOp>;
};
};
}; // namespace PlasticOps
namespace PlasticOps {
using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
using CommonPlasticPtr = boost::shared_ptr<PlasticOps::CommonData>;
using CommonHenkyPtr = boost::shared_ptr<HenckyOps::CommonData>;
template <int DIM>
MoFEMErrorCode
addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
boost::shared_ptr<MatrixDouble> mat_D_Ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
double scale_value, Sev sev) {
struct OpMatBlocks : public DomainEleOp {
OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
double scale_value, MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m_D_ptr),
matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
auto getK = [](auto &p) {
auto young_modulus = p[CommonData::YOUNG_MODULUS];
auto poisson_ratio = p[CommonData::POISSON_RATIO];
return young_modulus / (3 * (1 - 2 * poisson_ratio));
};
auto getG = [](auto &p) {
auto young_modulus = p[CommonData::YOUNG_MODULUS];
auto poisson_ratio = p[CommonData::POISSON_RATIO];
return young_modulus / (2 * (1 + poisson_ratio));
};
auto scale_fun = [this](auto &p) {
p[CommonData::YOUNG_MODULUS] *= scaleVal;
p[CommonData::SIGMA_Y] *= scaleVal;
p[CommonData::H] *= scaleVal;
p[CommonData::VIS_H] *= scaleVal;
p[CommonData::QINF] *= scaleVal;
p[CommonData::C1_k] *= scaleVal;
};
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*matParamsPtr = b.bParams;
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
}
(*matParamsPtr) = {young_modulus, poisson_ratio, sigmaY, H,
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
const double scaleVal;
struct BlockData {
std::array<double, CommonData::LAST_PARAM> bParams;
Range blockEnts;
};
std::vector<BlockData> blockData;
/**
* @brief Extract block data from meshsets
*
* @param m_field
* @param meshset_vec_ptr
* @param sev
* @return MoFEMErrorCode
*/
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != CommonData::LAST_PARAM) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of block attribute");
}
auto get_block_ents = [&]() {
Range ents;
CHKERR m_field.get_moab().get_entities_by_handle(m->meshset, ents,
true);
return ents;
};
CommonData::BlockParams block_params;
for (auto i = 0; i != CommonData::LAST_PARAM; ++i) {
block_params[i] = block_data[i];
}
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< "E = " << block_params[CommonData::YOUNG_MODULUS]
<< " nu = " << block_params[CommonData::POISSON_RATIO];
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< std::endl
<< "sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
<< "h = " << block_params[CommonData::H] << std::endl
<< "vis_h = " << block_params[CommonData::VIS_H] << std::endl
<< "qinf = " << block_params[CommonData::QINF] << std::endl
<< "biso = " << block_params[CommonData::BISO] << std::endl
<< "C1_k = " << block_params[CommonData::C1_k] << std::endl;
blockData.push_back({block_params, get_block_ents()});
}
}
/**
* @brief Get elasticity tensor
*
* Calculate elasticity tensor for given material parameters
*
* @param mat_D_ptr
* @param bulk_modulus_K
* @param shear_modulus_G
* @return MoFEMErrorCode
*
*/
MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
double bulk_modulus_K, double shear_modulus_G) {
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
double A = (DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
t_D(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
t_kd(k, l);
};
//! [Calculate elasticity tensor]
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
mat_D_ptr->resize(size_symm * size_symm, 1);
set_material_stiffness();
}
};
// push operator to calculate material stiffness matrix for each block
pip.push_back(new OpMatBlocks(
mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
// Get blockset using regular expression
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
));
}
template <int DIM, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field, std::string block_name,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string u, std::string ep, std::string tau, double scale, Sev sev) {
using P = PlasticityIntegrators<DomainEleOp>;
auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
auto make_d_mat = []() {
return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
};
common_plastic_ptr->mDPtr = make_d_mat();
common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
auto m_D_ptr = common_plastic_ptr->mDPtr;
CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, block_name, pip, m_D_ptr,
common_plastic_ptr->getParamsPtr(),
scale, sev),
"add mat block ops");
pip.push_back(new OpCalculateScalarFieldValues(
tau, common_plastic_ptr->getPlasticTauPtr()));
pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
ep, common_plastic_ptr->getPlasticStrainPtr()));
u, common_plastic_ptr->mGradPtr));
CommonHenkyPtr common_henky_ptr;
common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
common_henky_ptr->matLogCPlastic =
common_plastic_ptr->getPlasticStrainPtr();
common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
u, common_henky_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
u, common_henky_ptr));
pip.push_back(new
typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
u, common_henky_ptr));
} else {
pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
common_plastic_ptr, m_D_ptr));
}
pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
u, common_plastic_ptr));
return std::make_tuple(common_plastic_ptr, common_henky_ptr);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, scale, Sev::inform);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
// Calculate internal forces
if (common_henky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
DIM, I>(ep, common_plastic_ptr, m_D_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
using namespace HenckyOps;
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template BiLinearForm<I>;
using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, scale, Sev::verbose);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
if (common_henky_ptr) {
pip.push_back(new typename H::template OpHenckyTangent<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(new OpKCauchy(u, u, m_D_ptr));
pip.push_back(new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
u, ep, common_plastic_ptr, m_D_ptr));
}
if (common_henky_ptr) {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
DIM, I>(tau, u, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
DIM, I>(ep, u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
DIM, I>(ep, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
DIM, I>(ep, tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
DIM, I>(tau, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
I>(tau, tau, common_plastic_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string block_name, Pip &pip,
std::string u, std::string ep,
std::string tau) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, 1., Sev::inform);
// Calculate internal forces
if (common_henky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
}
}
} // namespace PlasticOps
static Index< 'p', 3 > p
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double bulk_modulus_K
double shear_modulus_G
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'M', 3 > M
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
FTensor::Index< 'N', 3 > N
FTensor::Index< 'I', 3 > I
[Common data]
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'J', 3 > J
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
double young_modulus
Young modulus.
Definition plastic.cpp:172
double C1_k
Kinematic hardening.
Definition plastic.cpp:180
double Qinf
Saturation yield stress.
Definition plastic.cpp:178
double visH
Viscous hardening.
Definition plastic.cpp:176
double scale
Definition plastic.cpp:170
constexpr auto size_symm
Definition plastic.cpp:42
double H
Hardening.
Definition plastic.cpp:175
double b_iso
Saturation exponent.
Definition plastic.cpp:179
PetscBool is_large_strains
Large strains.
Definition plastic.cpp:167
double sigmaY
Yield stress.
Definition plastic.cpp:174
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition seepage.cpp:63
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition seepage.cpp:65
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< MatrixDouble > mStrainPtr
MatrixDouble resFlowDstrainDot
MatrixDouble resCdStrain
boost::shared_ptr< MatrixDouble > mGradPtr
static std::array< int, 5 > activityData
MatrixDouble resCdPlasticStrain
VectorDouble plasticTauDot
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
VectorDouble plasticSurface
[Common data set externally]
MatrixDouble resFlowDtau
MatrixDouble plasticStrain
boost::shared_ptr< MatrixDouble > mStressPtr
MatrixDouble plasticFlow
MatrixDouble plasticStrainDot
VectorDouble plasticTau
MatrixDouble resFlowDstrain
std::array< double, LAST_PARAM > BlockParams
OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_LogStrain_dU
OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_LogStrain_dU
OpCalculateConstraintsLhs_dTAUImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dTAU
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_LogStrain_dEP
OpCalculatePlasticFlowLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dEP
OpCalculateConstraintsLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dEP
OpCalculatePlasticFlowLhs_dTAUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dTAU
OpCalculatePlasticFlowLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dU
OpCalculatePlasticFlowRhsImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowRhs
OpCalculateConstraintsLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dU
OpCalculateConstraintsRhsImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsRhs
OpCalculatePlasticInternalForceLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_dEP
OpPlasticStressImpl< DIM, I, DomainEleOp > OpPlasticStress
OpCalculatePlasticityImpl< DIM, I, DomainEleOp > OpCalculatePlasticity
OpCalculatePlasticSurfaceImpl< DIM, I, DomainEleOp > OpCalculatePlasticSurface
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy