v0.13.1
Loading...
Searching...
No Matches
ADV-0: Plastic problem
Note
Prerequisites of this tutorial include VEC-0: Linear elasticity


Note
Intended learning outcome:
  • general structure of a program developed using MoFEM
  • multi-field problem
  • linearisation
  • tangent operators
/**
* \file plastic.cpp
* \example plastic.cpp
*
* Plasticity in 2d and 3d
*
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType G =
IntegrationType::GAUSS; //< selected integration type
//! [Only used with Hooke equation (linear material model)]
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
//! [Only used for dynamics]
//! [Only used with Hencky/nonlinear material]
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
//! [Only used with Hencky/nonlinear material]
//! [Essential boundary conditions]
//! [Essential boundary conditions]
using OpBodyForce =
using OpForce =
PetscBool is_large_strains = PETSC_TRUE;
double scale = 1.;
double young_modulus = 206913;
double poisson_ratio = 0.29;
double rho = 0;
double sigmaY = 450;
double H = 129;
double visH = 0;
double cn0 = 1;
double cn1 = 1;
double zeta = 5e-2;
double Qinf = 265;
double b_iso = 16.93;
int order = 2; ///< Order if fixed.
int geom_order = 2; ///< Order if fixed.
constexpr size_t activ_history_sise = 1;
inline double hardening_exp(double tau) {
return std::exp(
std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
-b_iso * tau));
}
inline double hardening(double tau) {
return H * tau + Qinf * (1. - hardening_exp(tau)) + sigmaY;
}
inline double hardening_dtau(double tau) {
auto r = [](auto tau) { return H + Qinf * b_iso * hardening_exp(tau); };
return std::max(r(tau), 1e-6 * r(0));
}
inline double hardening_dtau2(double tau) {
return -(Qinf * (b_iso * b_iso)) * hardening_exp(tau);
}
#include <HenckyOps.hpp>
#include <PlasticOps.hpp>
using namespace PlasticOps;
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<DomainEle> reactionFe;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
struct PlasticityTimeScale : public MoFEM::TimeScale {
double getScale(const double time) {
};
};
};
//! [Run problem]
}
//! [Run problem]
//! [Set up problem]
Range domain_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
true);
auto get_ents_by_dim = [&](const auto dim) {
if (dim == SPACE_DIM) {
return domain_ents;
} else {
Range ents;
if (dim == 0)
CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
else
CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(SPACE_DIM);
if (domain_ents.empty())
const auto type = type_from_handle(domain_ents[0]);
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
}
return NOBASE;
};
const auto base = get_base();
MOFEM_LOG("WORLD", Sev::inform) << "Base " << ApproximationBaseNames[base];
CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
CHKERR simple->addDomainField("EP", L2, base, size_symm);
CHKERR simple->addDomainField("TAU", L2, base, 1);
CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
auto ents = get_ents_by_dim(0);
ents.merge(get_ents_by_dim(1));
// ents.merge(get_ents_by_dim(2));
CHKERR simple->setFieldOrder("U", order, &ents);
CHKERR simple->setFieldOrder("EP", order - 1);
CHKERR simple->setFieldOrder("TAU", order - 2);
CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
CHKERR simple->setUp();
CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
}
//! [Set up problem]
//! [Create common data]
auto get_command_line_parameters = [&]() {
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
&young_modulus, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
&poisson_ratio, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
&is_large_strains, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
PETSC_NULL);
MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
MOFEM_LOG("EXAMPLE", Sev::inform) << "cn0 " << cn0;
MOFEM_LOG("EXAMPLE", Sev::inform) << "cn1 " << cn1;
MOFEM_LOG("EXAMPLE", Sev::inform) << "zeta " << zeta;
MOFEM_LOG("EXAMPLE", Sev::inform) << "order " << order;
MOFEM_LOG("EXAMPLE", Sev::inform) << "geom order " << geom_order;
PetscBool is_scale = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
PETSC_NULL);
if (is_scale) {
rho *= scale;
H *= scale;
Qinf *= scale;
visH *= scale;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Young modulus " << young_modulus;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Poisson ratio " << poisson_ratio;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Saturation yield stress " << Qinf;
}
};
auto set_matrial_stiffness = [&]() {
const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
// Plane stress or when 1, plane strain or 3d
const double A = (SPACE_DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Axiator);
auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Deviator);
constexpr double third = boost::math::constants::third<double>();
t_D_axiator(i, j, k, l) = A *
(bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
t_kd(i, j) * t_kd(k, l);
t_D_deviator(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
};
auto make_d_mat = []() {
return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
};
commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
commonPlasticDataPtr->mDPtr = make_d_mat();
commonPlasticDataPtr->mDPtr_Axiator = make_d_mat();
commonPlasticDataPtr->mDPtr_Deviator = make_d_mat();
commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
CHKERR get_command_line_parameters();
CHKERR set_matrial_stiffness();
commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
commonHenckyDataPtr->matLogCPlastic =
commonPlasticDataPtr->getPlasticStrainPtr();
commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
commonPlasticDataPtr->mStressPtr =
commonHenckyDataPtr->getMatHenckyStress();
}
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
auto &bc_map = bc_mng->getBcMapByBlockName();
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
"U", 0, 3);
for (auto bc : bc_map)
MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
// OK. We have problem with GMesh, it adding empty characters at the end of
// block. So first block is search by regexp. popMarkDOFsOnEntities should
// work with regexp.
std::string reaction_block_set;
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "REACTION")) {
reaction_block_set = bc.first;
break;
}
}
if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
reactionMarker = bc->getBcMarkersPtr();
// Only take reaction from nodes
Range nodes;
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
ProblemsManager::MarkOP::AND, nodes,
} else {
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto add_domain_base_ops = [&](auto &pipeline) {
"GEOMETRY");
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
pipeline.push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
"EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
"U", commonPlasticDataPtr->mGradPtr));
};
auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
"U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
} else {
pipeline.push_back(
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(
new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
}
if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator) {
pipeline.push_back(
pipeline.push_back(new OpCalculatePlasticity(
"U", mField, commonPlasticDataPtr, m_D_ptr));
}
};
auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(
new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
"U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
// pipeline.push_back(new OpCalculateArgLhs_LogStrain_dUdTau(
// "U", "TAU", commonPlasticDataPtr, commonHenckyDataPtr));
// pipeline.push_back(new OpCalculateArgLhs_LogStrain_dUdU(
// "U", "U", commonPlasticDataPtr, commonHenckyDataPtr));
} else {
pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
"U", "EP", commonPlasticDataPtr, m_D_ptr));
// pipeline.push_back(
// new OpCalculateArgLhs_dUdTau("U", "TAU", commonPlasticDataPtr));
// pipeline.push_back(
// new OpCalculateArgLhs_dUdU("U", "U", commonPlasticDataPtr));
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
"BODY_FORCE", Sev::inform);
// Calculate internal forces
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
"TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
"EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(new OpCalculateConstraintsLhs_dU(
"TAU", "U", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
"EP", "U", commonPlasticDataPtr, m_D_ptr));
}
pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
"EP", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculatePlasticFlowLhs_dTAU(
"EP", "TAU", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculateConstraintsLhs_dEP(
"TAU", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(
// pipeline.push_back(new OpCalculateArgRhs_LogStrain_dU(
// "U", commonPlasticDataPtr, commonHenckyDataPtr));
} else {
// pipeline.push_back(new OpCalculateArgRhs_dU("U", commonPlasticDataPtr));
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
mField, pipeline, simple->getProblemName(), "U");
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline, {NOSPACE}, "GEOMETRY");
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
"FORCE", Sev::inform);
pipeline.push_back(new OpUnSetBc("U"));
auto u_mat_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
mField, pipeline, simple->getProblemName(), "U", u_mat_ptr,
{boost::make_shared<TimeScale>()}); // note displacements have no
// scaling
};
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
add_boundary_ops_lhs_mechanical(pipeline_mng->getOpBoundaryLhsPipeline());
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline(),
// Boundary
add_boundary_ops_rhs_mechanical(pipeline_mng->getOpBoundaryRhsPipeline());
auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(vol_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(vol_rule);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto create_reaction_pipeline = [&](auto &pipeline) {
"GEOMETRY");
pipeline.push_back(
"U", commonPlasticDataPtr->mGradPtr));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(new OpPlasticStress(
}
pipeline.push_back(new OpSetBc("U", false, reactionMarker));
// Calculate internal force
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
}
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = vol_rule;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
}
//! [Push operators to pipeline]
//! [Solve]
auto snes_ctx_ptr = smartGetDMSnesCtx(simple->getDM());
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_element = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(mField);
pp_fe->getOpPtrVector(), {H1}, "GEOMETRY");
auto x_ptr = boost::make_shared<MatrixDouble>();
pp_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
pp_fe->getOpPtrVector().push_back(
pp_fe->getOpPtrVector().push_back(
"U", commonPlasticDataPtr->mGradPtr));
pp_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
pp_fe->getOpPtrVector().push_back(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
pp_fe->getOpPtrVector().push_back(
pp_fe->getOpPtrVector().push_back(
pp_fe->getOpPtrVector().push_back(
pp_fe->getOpPtrVector().push_back(
1 /*scale*/));
pp_fe->getOpPtrVector().push_back(
pp_fe->getOpPtrVector().push_back(
new OpPPMap(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{{"GRAD", commonPlasticDataPtr->mGradPtr},
{"FIRST_PIOLA", commonHenckyDataPtr->getMatFirstPiolaStress()}},
{}
)
);
} else {
pp_fe->getOpPtrVector().push_back(
commonPlasticDataPtr->mStrainPtr));
pp_fe->getOpPtrVector().push_back(new OpPlasticStress(
"U", commonPlasticDataPtr, commonPlasticDataPtr->mDPtr, 1 /*scale*/));
pp_fe->getOpPtrVector().push_back(
new OpPPMap(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{},
{{"STRAIN", commonPlasticDataPtr->mStrainPtr},
{"STRESS", commonPlasticDataPtr->mStressPtr}}
)
);
}
pp_fe->getOpPtrVector().push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
pp_fe->getOpPtrVector().push_back(
new OpPPMap(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()}},
{},
{},
{{"PLASTIC_STRAIN", commonPlasticDataPtr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", commonPlasticDataPtr->getPlasticFlowPtr()}}
)
);
return pp_fe;
};
auto scatter_create = [&](auto D, auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
ROW, "U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
Vec v;
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
VecScatter scatter;
CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(v),
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, create_post_process_element(), reactionFe, uXScatter,
uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Setup fieldsplit (block) solver - optional: yes/no
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
}
};
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
uXScatter = scatter_create(D, 0);
uYScatter = scatter_create(D, 1);
if constexpr (SPACE_DIM == 3)
uZScatter = scatter_create(D, 2);
auto solver = pipeline_mng->createTSIM();
auto pre_rhs = [&]() {
std::fill(commonPlasticDataPtr->activityData.begin(),
commonPlasticDataPtr->activityData.end(), 0);
};
auto post_rhs = [&]() {
auto get_iter = [&]() {
SNES snes;
CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
int iter;
CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
"Can not get iter");
return iter;
};
auto iter = get_iter();
if (iter >= 0) {
std::array<int, 5> activity_data;
std::fill(activity_data.begin(), activity_data.end(), 0);
MPI_Allreduce(commonPlasticDataPtr->activityData.data(),
activity_data.data(), activity_data.size(), MPI_INT,
MPI_SUM, mField.get_comm());
int &active_points = activity_data[0];
int &avtive_full_elems = activity_data[1];
int &avtive_elems = activity_data[2];
int &nb_points = activity_data[3];
int &nb_elements = activity_data[4];
if (nb_points) {
double proc_nb_points =
100 * static_cast<double>(active_points) / nb_points;
double proc_nb_active =
100 * static_cast<double>(avtive_elems) / nb_elements;
double proc_nb_full_active = 100;
if (avtive_elems)
proc_nb_full_active =
100 * static_cast<double>(avtive_full_elems) / avtive_elems;
"EXAMPLE", Sev::inform,
"Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active elemens %d "
"(%3.3f\%) nb full active elems %d (%3.3f\%)",
iter, nb_points, active_points, proc_nb_points, avtive_elems,
proc_nb_active, avtive_full_elems, proc_nb_full_active, iter);
}
}
};
mField.getInterface<PipelineManager>()->getDomainLhsFE()->preProcessHook =
pre_rhs;
mField.getInterface<PipelineManager>()->getDomainLhsFE()->postProcessHook =
post_rhs;
CHKERR TSSetSolution(solver, D);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
std::string param_file
constexpr double third
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr double shear_modulus_G
Definition: elastic.cpp:63
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
const int dim
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
constexpr auto t_kd
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1592
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:1003
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:987
const double r
rate factor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Definition: plastic.cpp:96
double Qinf
Definition: plastic.cpp:105
double hardening_dtau2(double tau)
Definition: plastic.cpp:128
constexpr size_t activ_history_sise
Definition: plastic.cpp:111
double rho
Definition: plastic.cpp:98
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:50
double visH
Definition: plastic.cpp:101
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:71
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:52
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:73
double scale
Definition: plastic.cpp:94
constexpr auto size_symm
Definition: plastic.cpp:33
double zeta
Definition: plastic.cpp:104
double H
Definition: plastic.cpp:100
double hardening(double tau)
Definition: plastic.cpp:119
constexpr AssemblyType A
Definition: plastic.cpp:35
double cn0
Definition: plastic.cpp:102
double b_iso
Definition: plastic.cpp:106
PetscBool is_large_strains
Definition: plastic.cpp:92
int geom_order
Order if fixed.
Definition: plastic.cpp:109
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:64
double sigmaY
Definition: plastic.cpp:99
double hardening_dtau(double tau)
Definition: plastic.cpp:123
constexpr IntegrationType G
Definition: plastic.cpp:36
double hardening_exp(double tau)
Definition: plastic.cpp:113
double cn1
Definition: plastic.cpp:103
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:75
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:165
[Example]
Definition: plastic.cpp:137
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:745
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:157
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:257
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:453
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:172
MoFEM::Interface & mField
Definition: plastic.cpp:144
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:158
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:160
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:184
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:391
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:161
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:154
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:152
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:156
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:153
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Natural boundary conditions.
Definition: Essential.hpp:101
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Assembly methods.
Definition: Natural.hpp:57
Get rate of scalar field at integration points.
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Definition: Essential.hpp:71
Enforce essential constrains on rhs.
Definition: Essential.hpp:55
Post post-proc data at points from hash maps.
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
PipelineManager interface.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
/** \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> {
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mDPtr_Axiator;
boost::shared_ptr<MatrixDouble> mDPtr_Deviator;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
VectorDouble plasticSurface;
MatrixDouble plasticFlow;
VectorDouble plasticTau;
VectorDouble plasticTauDot;
MatrixDouble plasticStrain;
MatrixDouble plasticStrainDot;
VectorDouble resC;
VectorDouble resCdTau;
MatrixDouble resCdStrain;
MatrixDouble resCdStrainDot;
MatrixDouble resFlow;
MatrixDouble resFlowDtau;
MatrixDouble resFlowDstrain;
MatrixDouble resFlowDstrainDot;
std::array<int, 5> activityData;
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);
}
};
//! [Common data]
//! [Operators definitions]
struct OpCalculatePlasticSurface : public DomainEleOp {
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculatePlasticity : public DomainEleOp {
OpCalculatePlasticity(const std::string field_name, MoFEM::Interface &m_field,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpPlasticStress : public DomainEleOp {
OpPlasticStress(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> mDPtr,
const double scale = 1);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<MatrixDouble> mDPtr;
const double scaleStress;
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculatePlasticFlowRhs : public AssemblyDomainEleOp {
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;
};
struct OpCalculateConstraintsRhs : public AssemblyDomainEleOp {
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;
};
struct OpCalculatePlasticInternalForceLhs_dEP : public AssemblyDomainEleOp {
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;
};
struct OpCalculatePlasticInternalForceLhs_LogStrain_dEP
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> common_henky_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<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_dU : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dU(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;
};
struct OpCalculatePlasticFlowLhs_LogStrain_dU : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> comman_henky_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<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_dEP : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dEP(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;
};
struct OpCalculatePlasticFlowLhs_dTAU : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dTAU(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;
};
struct OpCalculateConstraintsLhs_dU : public AssemblyDomainEleOp {
OpCalculateConstraintsLhs_dU(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;
};
struct OpCalculateConstraintsLhs_LogStrain_dU : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> comman_henky_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<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculateConstraintsLhs_dEP : public AssemblyDomainEleOp {
OpCalculateConstraintsLhs_dEP(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;
};
struct OpCalculateConstraintsLhs_dTAU : public AssemblyDomainEleOp {
OpCalculateConstraintsLhs_dTAU(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;
};
//! [Operators definitions]
//! [Lambda functions]
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;
if constexpr (SPACE_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;
}
inline auto diff_symmetrize() {
t_diff(i, j, k, l) = 0;
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 (SPACE_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>();
if constexpr (SPACE_DIM == 2)
return (t_stress(0, 0) + t_stress(1, 1)) * third;
else
return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) * third;
};
template <typename T>
double trace) {
t_dev(I, J) = 0;
for (int ii = 0; ii != SPACE_DIM; ++ii)
for (int jj = ii; jj != SPACE_DIM; ++jj)
t_dev(ii, jj) = t_stress(ii, jj);
t_dev(0, 0) -= trace;
t_dev(1, 1) -= trace;
t_dev(2, 2) -= trace;
return t_dev;
};
inline auto
t_diff_deviator(I, J, k, l) = 0;
for (int ii = 0; ii != SPACE_DIM; ++ii)
for (int jj = ii; jj != SPACE_DIM; ++jj)
for (int kk = 0; kk != SPACE_DIM; ++kk)
for (int ll = kk; ll != SPACE_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 (SPACE_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;
};
/**
*
\f[
\begin{split}
f&=\sqrt{s_{ij}s_{ij}}\\
A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
\frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
\frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
\sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
-A_{mn}A_{ij}
\right)\\
\frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
\\
\frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
\frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
\sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
\sigma_{mn}} D_{mnkl}
\end{split}
\f]
*/
inline double
return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
std::numeric_limits<double>::epsilon();
};
inline auto plastic_flow(long double f,
t_diff_f(k, l) =
(1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l))) / f;
return t_diff_f;
};
template <typename T>
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>
FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> &&t_diff_plastic_flow_dstress) {
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;
};
inline double constrain_diff_sign(double x) {
const auto y = x / zeta;
if (y > std::numeric_limits<float>::max_exponent10 ||
y < std::numeric_limits<float>::min_exponent10) {
return 0;
} else {
const auto e = std::exp(y);
const auto ep1 = e + 1;
return (2 / zeta) * (e / (ep1 * ep1));
}
};
inline double constrian_sign(double x) {
const auto y = x / zeta;
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);
}
};
inline double constrain_abs(double x) {
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) {
return (f - sigma_y) / sigmaY + cn1 * (dot_tau * eqiv);
};
/**
\f[
\dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
\| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0 \\
c_n \sigma_y \dot{\tau} - \frac{1}{2}\left\{c_n\sigma_y \dot{\tau} +
(f(\pmb\sigma) - \sigma_y) +
\| c_n \sigma_y \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
\f]
*/
inline double constrain(double eqiv, double dot_tau, double f, double sigma_y,
double abs_w) {
return visH * dot_tau +
(sigmaY / 2) * ((cn0 * (dot_tau - eqiv) + cn1 * (eqiv * dot_tau) -
(f - sigma_y) / sigmaY) -
abs_w);
};
inline double diff_constrain_ddot_tau(double sign, double eqiv) {
return visH + (sigmaY / 2) * (cn0 + cn1 * eqiv * (1 - sign));
};
inline double diff_constrain_equiv(double sign, double dot_tau) {
return (sigmaY / 2) * (-cn0 + cn1 * dot_tau * (1 - sign));
};
inline auto diff_constrain_df(double sign) { return (-1 - sign) / 2; };
inline auto diff_constrain_dsigma_y(double sign) { return (1 + sign) / 2; }
template <typename T>
t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
return t_diff_constrain_dstress;
};
template <typename T1, typename T2>
inline auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress) {
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>
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>
inline auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot,
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;
};
//! [Lambda functions]
//! [Auxiliary functions functions
static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
&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)};
}
static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
&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)};
}
//! [Auxiliary functions functions
}; // namespace PlasticOps
auto f
Definition: HenckyOps.hpp:5
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
double sign(double x)
auto diff_constrain_dstress(double diff_constrain_df, FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_flow)
Definition: PlasticOps.hpp:549
auto symm_L_tensor()
Definition: PlasticOps.hpp:315
auto deviator(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress, double trace)
Definition: PlasticOps.hpp:373
double constrian_sign(double x)
Definition: PlasticOps.hpp:488
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
Definition: PlasticOps.hpp:440
double diff_constrain_equiv(double sign, double dot_tau)
Definition: PlasticOps.hpp:540
FTensor::Index< 'p', SPACE_DIM > p
Definition: PlasticOps.hpp:99
FTensor::Index< 'j', SPACE_DIM > j
Definition: PlasticOps.hpp:93
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:103
auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
Definition: PlasticOps.hpp:558
double constrain_diff_sign(double x)
Definition: PlasticOps.hpp:476
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
Definition: PlasticOps.hpp:445
FTensor::Index< 'O', size_symm > O
Definition: PlasticOps.hpp:107
FTensor::Index< 'L', size_symm > L
Definition: PlasticOps.hpp:106
FTensor::Index< 'l', SPACE_DIM > l
Definition: PlasticOps.hpp:95
FTensor::Index< 'k', SPACE_DIM > k
Definition: PlasticOps.hpp:94
auto diff_tensor()
[Operators definitions]
Definition: PlasticOps.hpp:308
auto diff_constrain_df(double sign)
Definition: PlasticOps.hpp:544
double w(double eqiv, double dot_tau, double f, double sigma_y)
Definition: PlasticOps.hpp:513
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:104
auto diff_constrain_dsigma_y(double sign)
Definition: PlasticOps.hpp:546
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:101
auto diff_deviator(FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_stress)
Definition: PlasticOps.hpp:387
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: PlasticOps.hpp:92
FTensor::Index< 'o', SPACE_DIM > o
Definition: PlasticOps.hpp:98
FTensor::Index< 'm', SPACE_DIM > m
Definition: PlasticOps.hpp:96
double constrain_abs(double x)
Definition: PlasticOps.hpp:502
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:102
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain)
Definition: PlasticOps.hpp:575
double trace(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress)
Definition: PlasticOps.hpp:364
auto diff_symmetrize()
Definition: PlasticOps.hpp:333
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_flow, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
Definition: PlasticOps.hpp:455
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, SPACE_DIM, SPACE_DIM > &t_D, FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_plastic_flow_dstress)
Definition: PlasticOps.hpp:467
double diff_constrain_ddot_tau(double sign, double eqiv)
Definition: PlasticOps.hpp:536
double constrain(double eqiv, double dot_tau, double f, double sigma_y, double abs_w)
Definition: PlasticOps.hpp:528
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_strain_dot)
Definition: PlasticOps.hpp:566
FTensor::Index< 'n', SPACE_DIM > n
Definition: PlasticOps.hpp:97
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
Definition: PlasticOps.hpp:587
constexpr auto field_name
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: PlasticOps.hpp:48
MatrixDouble resFlowDstrainDot
Definition: PlasticOps.hpp:65
VectorDouble resCdTau
Definition: PlasticOps.hpp:59
MatrixDouble resCdStrain
Definition: PlasticOps.hpp:60
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: PlasticOps.hpp:47
MatrixDouble resFlow
Definition: PlasticOps.hpp:62
VectorDouble plasticTauDot
Definition: PlasticOps.hpp:54
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:44
VectorDouble plasticSurface
Definition: PlasticOps.hpp:51
boost::shared_ptr< MatrixDouble > mDPtr_Axiator
Definition: PlasticOps.hpp:45
MatrixDouble resFlowDtau
Definition: PlasticOps.hpp:63
MatrixDouble resCdStrainDot
Definition: PlasticOps.hpp:61
MatrixDouble plasticStrain
Definition: PlasticOps.hpp:55
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: PlasticOps.hpp:49
MatrixDouble plasticFlow
Definition: PlasticOps.hpp:52
MatrixDouble plasticStrainDot
Definition: PlasticOps.hpp:56
VectorDouble plasticTau
Definition: PlasticOps.hpp:53
MatrixDouble resFlowDstrain
Definition: PlasticOps.hpp:64
boost::shared_ptr< MatrixDouble > mDPtr_Deviator
Definition: PlasticOps.hpp:46
std::array< int, 5 > activityData
Definition: PlasticOps.hpp:67
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:277
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:275
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:276
OpCalculateConstraintsLhs_LogStrain_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > comman_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:290
OpCalculateConstraintsLhs_dEP(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)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:291
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:302
OpCalculateConstraintsLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpCalculateConstraintsLhs_dU(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)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:261
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:262
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculateConstraintsRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:162
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:163
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticFlowLhs_LogStrain_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > comman_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:222
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:221
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:220
OpCalculatePlasticFlowLhs_dEP(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)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:236
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:235
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:248
OpCalculatePlasticFlowLhs_dTAU(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)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:249
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticFlowLhs_dU(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)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:207
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:206
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:152
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
[Calculate stress]
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:151
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:191
OpCalculatePlasticInternalForceLhs_LogStrain_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > common_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticInternalForceLhs_dEP(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)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:175
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:176
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:116
OpCalculatePlasticity(const std::string field_name, MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:127
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:128
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:141
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mDPtr, const double scale=1)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:139