v0.13.2
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 = (SCHUR_ASSEMBLE)
? AssemblyType::SCHUR
: 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;
PetscBool set_timer = PETSC_FALSE;
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); };
constexpr double eps = 1e-12;
return std::max(r(tau), eps * 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 PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
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 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));
} else {
pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
"U", "EP", commonPlasticDataPtr, m_D_ptr));
}
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 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(pip->getOpDomainLhsPipeline());
CHKERR add_domain_stress_ops(pip->getOpDomainLhsPipeline(),
CHKERR add_domain_ops_lhs_mechanical(pip->getOpDomainLhsPipeline(),
CHKERR add_domain_ops_lhs_constrain(pip->getOpDomainLhsPipeline(),
add_boundary_ops_lhs_mechanical(pip->getOpBoundaryLhsPipeline());
CHKERR add_domain_base_ops(pip->getOpDomainRhsPipeline());
CHKERR add_domain_stress_ops(pip->getOpDomainRhsPipeline(),
CHKERR add_domain_ops_rhs_mechanical(pip->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_constrain(pip->getOpDomainRhsPipeline(),
// Boundary
add_boundary_ops_rhs_mechanical(pip->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 pip->setDomainRhsIntegrationRule(vol_rule);
CHKERR pip->setDomainLhsIntegrationRule(vol_rule);
CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pip->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]
struct SetUpSchur {
static boost::shared_ptr<SetUpSchur> createSetUpSchur(
);
virtual MoFEMErrorCode setUp(KSP solver) = 0;
virtual MoFEMErrorCode preProc() = 0;
virtual MoFEMErrorCode postProc() = 0;
protected:
SetUpSchur() = default;
};
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(
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(
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,
boost::shared_ptr<SetUpSchur>
&schur_ptr) {
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 create_sub_bc_dm = [&](SmartPetscObj<DM> base_dm,
SmartPetscObj<AO> &ao_sub) {
dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_BC");
CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
for (auto f : {"U", "EP", "TAU"}) {
}
CHKERR DMSetUp(dm_sub);
CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_X", "U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Y", "U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Z", "U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_ALL", "U", 0,
2);
auto *prb_ptr = getProblemPtr(dm_sub);
if (auto sub_data = prb_ptr->getSubData()) {
is_sub = sub_data->getSmartRowIs();
ao_sub = sub_data->getSmartRowMap();
int is_sub_size;
CHKERR ISGetSize(is_sub, &is_sub_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split second block size " << is_sub_size;
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
}
};
auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,
SmartPetscObj<AO> &ao_sub) {
dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_U");
CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
for (auto f : {"U"}) {
}
CHKERR DMSetUp(dm_sub);
auto *prb_ptr = getProblemPtr(dm_sub);
if (auto sub_data = prb_ptr->getSubData()) {
is_sub = sub_data->getSmartRowIs();
ao_sub = sub_data->getSmartRowMap();
// ISView(is_sub, PETSC_VIEWER_STDOUT_WORLD);
int is_sub_size;
CHKERR ISGetSize(is_sub, &is_sub_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split second block size " << is_sub_size;
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
}
CHKERR DMSetUp(dm_sub);
};
auto create_all_bc_is = [&](SmartPetscObj<IS> &is_all_bc) {
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 first block size " << is_all_bc_size;
};
SmartPetscObj<IS> is_all_bc;
SmartPetscObj<DM> dm_bc_sub;
SmartPetscObj<IS> is_bc_sub;
SmartPetscObj<AO> ao_bc_sub;
CHKERR create_all_bc_is(is_all_bc);
CHKERR create_sub_bc_dm(simple->getDM(), dm_bc_sub, is_bc_sub, ao_bc_sub);
CHKERR mField.getInterface<ISManager>()->isCreateProblem(
simple->getProblemName(), ROW, is_prb);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_bc_sub);
if constexpr (A == AssemblyType::SCHUR) {
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
"SUB_BC", ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
"SUB_BC", ROW, "TAU", 0, MAX_DOFS_ON_ENTITY, is_tau);
IS is_union_raw;
CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
SmartPetscObj<IS> is_union(is_union_raw);
CHKERR create_sub_u_dm(dm_bc_sub, dm_u_sub, is_u_sub, ao_u_sub);
// Indices has to be map fro very to level, while assembling Schur
// complement.
auto is_up = is_u_sub;
CHKERR AOPetscToApplicationIS(ao_bc_sub, is_up);
auto ao_up = createAOMappingIS(is_u_sub, PETSC_NULL);
schur_ptr =
SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);
PetscInt n;
KSP *ksps;
CHKERR PCFieldSplitGetSubKSP(pc, &n, &ksps);
CHKERR schur_ptr->setUp(ksps[1]);
}
}
};
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 = pip->createTSIM();
auto active_pre_lhs = [&]() {
std::fill(commonPlasticDataPtr->activityData.begin(),
commonPlasticDataPtr->activityData.end(), 0);
};
auto active_post_lhs = [&]() {
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 elements %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);
}
}
};
CHKERR TSSetSolution(solver, D);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
boost::shared_ptr<SetUpSchur> schur_ptr;
CHKERR set_fieldsplit_preconditioner(solver, schur_ptr);
mField.getInterface<PipelineManager>()->getDomainLhsFE()->preProcessHook =
[&]() {
if (schur_ptr)
CHKERR schur_ptr->preProc();
CHKERR active_pre_lhs();
};
mField.getInterface<PipelineManager>()->getDomainLhsFE()->postProcessHook =
[&]() {
// if (schur_ptr)
// CHKERR schur_ptr->postProc();
CHKERR active_post_lhs();
};
mField.getInterface<PipelineManager>()->getBoundaryLhsFE()->postProcessHook =
[&]() {
if (schur_ptr)
CHKERR schur_ptr->postProc();
};
MOFEM_LOG_TAG("TIMER", "timer");
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
CHKERR TSSetUp(solver);
MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
CHKERR TSSolve(solver, NULL);
MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
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"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
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 interface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
struct SetUpSchurImpl : public SetUpSchur {
: SetUpSchur(), mField(m_field), subDM(sub_dm), subIS(sub_is),
subAO(sub_ao) {
if (S) {
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
}
virtual ~SetUpSchurImpl() { S.reset(); }
MoFEMErrorCode setUp(KSP solver);
private:
};
PC pc;
CHKERR KSPSetFromOptions(solver);
CHKERR KSPGetPC(solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
if (S) {
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
} else {
pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpBoundaryLhsPipeline().push_back(
new OpSchurAssembleEnd({}, {}, {}, {}, {}));
pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(
new OpSchurAssembleEnd({}, {}, {}, {}, {}));
}
subDM.reset();
subIS.reset();
subAO.reset();
}
// Boundary
pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd(
{"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), subAO},
{SmartPetscObj<Mat>(), S}, {false, false}
));
// Domain
pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd(
{"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), subAO},
{SmartPetscObj<Mat>(), S}, {false, false}
));
}
CHKERR PCFieldSplitSetIS(pc, NULL, subIS);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
}
CHKERR MatZeroEntries(S);
}
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
}
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
if (S) {
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
}
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(
new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_sub));
}
std::string param_file
constexpr double third
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
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
static const double eps
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
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ 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:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
const int dim
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
constexpr auto t_kd
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:983
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:221
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:485
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:444
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:511
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:244
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:277
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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:1634
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: DMMoFEM.cpp:1044
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
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 createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:971
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1017
const double r
rate factor
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Definition: plastic.cpp:99
double Qinf
Definition: plastic.cpp:108
double hardening_dtau2(double tau)
Definition: plastic.cpp:132
constexpr size_t activ_history_sise
Definition: plastic.cpp:114
double rho
Definition: plastic.cpp:101
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:68
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:52
double visH
Definition: plastic.cpp:104
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:54
PetscBool set_timer
Definition: plastic.cpp:95
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:75
double scale
Definition: plastic.cpp:97
constexpr auto size_symm
Definition: plastic.cpp:33
double zeta
Definition: plastic.cpp:107
double H
Definition: plastic.cpp:103
double hardening(double tau)
Definition: plastic.cpp:122
double cn0
Definition: plastic.cpp:105
double b_iso
Definition: plastic.cpp:109
PetscBool is_large_strains
Definition: plastic.cpp:94
int geom_order
Order if fixed.
Definition: plastic.cpp:112
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:66
double sigmaY
Definition: plastic.cpp:102
double hardening_dtau(double tau)
Definition: plastic.cpp:126
constexpr IntegrationType G
Definition: plastic.cpp:38
double hardening_exp(double tau)
Definition: plastic.cpp:116
double cn1
Definition: plastic.cpp:106
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:77
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:169
[Example]
Definition: plastic.cpp:141
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:752
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:161
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:261
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:459
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:176
MoFEM::Interface & mField
Definition: plastic.cpp:148
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:162
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:164
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:397
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:165
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:158
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:156
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:157
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
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)
Essential 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.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:77
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:27
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.
[Push operators to pipeline]
Definition: plastic.cpp:736
SetUpSchur()=default
virtual MoFEMErrorCode setUp(KSP solver)=0
virtual MoFEMErrorCode postProc()=0
virtual MoFEMErrorCode preProc()=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM >, SmartPetscObj< IS >, SmartPetscObj< AO >)
Definition: plastic.cpp:1351
SmartPetscObj< PC > smartPC
Definition: plastic.cpp:1259
SmartPetscObj< DM > subDM
Definition: plastic.cpp:1262
MoFEMErrorCode setUp(KSP solver)
Definition: plastic.cpp:1267
SmartPetscObj< Mat > S
Definition: plastic.cpp:1258
SmartPetscObj< AO > subAO
Definition: plastic.cpp:1264
virtual ~SetUpSchurImpl()
Definition: plastic.cpp:1248
MoFEMErrorCode postProc()
Definition: plastic.cpp:1340
MoFEMErrorCode setPC(PC pc)
Definition: plastic.cpp:1323
SmartPetscObj< IS > subIS
Definition: plastic.cpp:1263
MoFEMErrorCode setOperator()
Definition: plastic.cpp:1301
MoFEMErrorCode preProc()
Definition: plastic.cpp:1331
MoFEM::Interface & mField
Definition: plastic.cpp:1261
/** \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 * std::sqrt(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 constraint(double eqiv, double dot_tau, double f, double sigma_y,
double abs_w) {
return visH * dot_tau + (sigmaY / 2) * ((cn0 * (dot_tau - eqiv) +
cn1 * (std::sqrt(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 * std::sqrt(eqiv) * (1 - sign));
};
inline double diff_constrain_eqiv(double sign, double eqiv, double dot_tau) {
return (sigmaY / 2) *
(-cn0 + cn1 * dot_tau * (0.5 / std::sqrt(eqiv)) * (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:550
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
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:559
double diff_constrain_eqiv(double sign, double eqiv, double dot_tau)
Definition: PlasticOps.hpp:540
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:545
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:547
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:576
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 constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w)
Definition: PlasticOps.hpp:528
double diff_constrain_ddot_tau(double sign, double eqiv)
Definition: PlasticOps.hpp:536
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_strain_dot)
Definition: PlasticOps.hpp:567
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:588
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