v0.13.2
Loading...
Searching...
No Matches
plastic.cpp

Plasticity in 2d and 3d

/**
* \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 =
DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
using OpForce =
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
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));
CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
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));
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
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: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
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:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
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
int r
Definition: sdf.py:34
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
BoundaryNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpForce
Definition: plastic.cpp:87
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
NaturalBC< DomainEleOp >::Assembly< A >::LinearForm< G > DomainNaturalBC
Definition: plastic.cpp:81
EssentialBC< BoundaryEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpEssentialLhs< DisplacementCubitBcData, 1, SPACE_DIM > OpEssentialLhs
Definition: plastic.cpp:90
double zeta
Definition: plastic.cpp:107
double H
Definition: plastic.cpp:103
double hardening(double tau)
Definition: plastic.cpp:122
NaturalBC< BoundaryEleOp >::Assembly< A >::LinearForm< G > BoundaryNaturalBC
Definition: plastic.cpp:85
EssentialBC< BoundaryEleOp >::Assembly< A >::LinearForm< GAUSS >::OpEssentialRhs< DisplacementCubitBcData, 1, SPACE_DIM > OpEssentialRhs
Definition: plastic.cpp:92
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
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:162
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:752
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:261
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:161
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:459
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:176
MoFEM::Interface & mField
Definition: plastic.cpp:148
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
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.
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:35
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