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

Plasticity in 2d and 3d

/**
* \file plastic.cpp
* \example plastic.cpp
*
* Plasticity in 2d and 3d
*
*/
/* The above code is a preprocessor directive in C++ that checks if the macro
"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it replaces
the " */
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
// #define ADD_CONTACT
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps;
template <> struct ElementsAndOps<2> {
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
template <> struct ElementsAndOps<3> {
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
constexpr AssemblyType AT =
(SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
: AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
IntegrationType::GAUSS; //< selected integration type
#ifdef ADD_CONTACT
//! [Specialisation for assembly]
// Assemble to A matrix, by default, however, some terms are assembled only to
// preconditioning.
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
};
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
};
/**
* @brief Element used to specialise assembly
*
*/
using BoundaryEleOp::BoundaryEleOp;
};
/**
* @brief Specialise assembly for Stabilised matrix
*
* @tparam
*/
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
};
//! [Specialisation for assembly]
#endif // ADD_CONTACT
inline double iso_hardening_exp(double tau, double b_iso) {
return std::exp(
std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
-b_iso * tau));
}
/**
* Isotropic hardening
*/
inline double iso_hardening(double tau, double H, double Qinf, double b_iso,
double sigmaY) {
return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;
}
inline double iso_hardening_dtau(double tau, double H, double Qinf,
double b_iso) {
auto r = [&](auto tau) {
return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);
};
constexpr double eps = 1e-12;
return std::max(r(tau), eps * r(0));
}
/**
* Kinematic hardening
*/
template <typename T, int DIM>
inline auto
double C1_k) {
t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
return t_alpha;
}
template <int DIM>
t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
return t_diff;
}
PetscBool is_large_strains = PETSC_TRUE; ///< Large strains
PetscBool set_timer = PETSC_FALSE; ///< Set timer
double scale = 1.;
double young_modulus = 206913; ///< Young modulus
double poisson_ratio = 0.29; ///< Poisson ratio
double sigmaY = 450; ///< Yield stress
double H = 129; ///< Hardening
double visH = 0; ///< Viscous hardening
double zeta = 5e-2; ///< Viscous hardening
double Qinf = 265; ///< Saturation yield stress
double b_iso = 16.93; ///< Saturation exponent
double C1_k = 0; ///< Kinematic hardening
double cn0 = 1;
double cn1 = 1;
int order = 2; ///< Order if fixed.
int geom_order = 2; ///< Order if fixed.
PetscBool is_quasi_static = PETSC_TRUE;
double rho = 0.0;
double alpha_damping = 0;
#include <HenckyOps.hpp>
#include <PlasticOps.hpp>
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
#endif
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
#include <ContactOps.hpp>
#endif // ADD_CONTACT
DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
using namespace PlasticOps;
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
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;
struct ScaledTimeScale : public MoFEM::TimeScale {
double getScale(const double time) {
};
};
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif // ADD_CONTACT
};
//! [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);
CHKERR simple->setFieldOrder("EP", order - 1);
CHKERR simple->setFieldOrder("TAU", order - 2);
CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
#ifdef ADD_CONTACT
CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
auto get_skin = [&]() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_blocks = [&](auto skin) {
Range contact_range;
for (auto m :
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
MOFEM_LOG("CONTACT", Sev::inform)
<< "Find contact block set: " << m->getName();
auto meshset = m->getMeshset();
CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
contact_range, true);
MOFEM_LOG("SYNC", Sev::inform)
<< "Nb entities in contact surface: " << contact_range.size();
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
contact_range);
skin = intersect(skin, contact_range);
}
return skin;
};
auto filter_true_skin = [&](auto skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
CHKERR simple->setFieldOrder("SIGMA", 0);
CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
#endif
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, "", "-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 PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, 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);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
&alpha_damping, PETSC_NULL);
MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
MOFEM_LOG("PLASTICITY", Sev::inform)
<< "Geometry approximation order " << geom_order;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
PetscBool is_scale = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
PETSC_NULL);
if (is_scale) {
}
MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
#ifdef ADD_CONTACT
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
&ContactOps::cn_contact, PETSC_NULL);
MOFEM_LOG("CONTACT", Sev::inform)
<< "cn_contact " << ContactOps::cn_contact;
#endif // ADD_CONTACT
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
&is_quasi_static, PETSC_NULL);
MOFEM_LOG("PLASTICITY", Sev::inform)
<< "Is quasi static: " << (is_quasi_static ? "true" : "false");
};
CHKERR get_command_line_parameters();
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdfPythonPtr->sdfInit("sdf.py");
ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif // ADD_CONTACT
}
//! [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);
#ifdef ADD_CONTACT
for (auto b : {"FIX_X", "REMOVE_X"})
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
"SIGMA", 0, 0, false, true);
for (auto b : {"FIX_Y", "REMOVE_Y"})
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
"SIGMA", 1, 1, false, true);
for (auto b : {"FIX_Z", "REMOVE_Z"})
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
"SIGMA", 2, 2, false, true);
for (auto b : {"FIX_ALL", "REMOVE_ALL"})
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
"SIGMA", 0, 3, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
#endif
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map)
MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto pip_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
"GEOMETRY");
pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
// Add Natural BCs to LHS
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip, mField, "U", Sev::inform);
#ifdef ADD_CONTACT
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U");
ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
vol_rule);
#endif // ADD_CONTACT
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
"GEOMETRY");
pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
// Add Natural BCs to RHS
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
pip, "SIGMA", "U");
#endif // ADD_CONTACT
};
auto add_domain_ops_lhs = [this](auto &pip) {
"GEOMETRY");
if (is_quasi_static == PETSC_FALSE) {
//! [Only used for dynamics]
//! [Only used for dynamics]
auto get_inertia_and_mass_damping = [this](const double, const double,
const double) {
auto &fe_domain_lhs = pip->getDomainLhsFE();
return (rho / scale) * fe_domain_lhs->ts_aa +
(alpha_damping / scale) * fe_domain_lhs->ts_a;
};
pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
}
CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
};
auto add_domain_ops_rhs = [this](auto &pip) {
"GEOMETRY");
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip, mField, "U",
{boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
Sev::inform);
// only in case of dynamics
if (is_quasi_static == PETSC_FALSE) {
//! [Only used for dynamics]
//! [Only used for dynamics]
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(
new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
return rho / scale;
}));
if (alpha_damping > 0) {
auto mat_velocity = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
new OpInertiaForce("U", mat_velocity, [](double, double, double) {
return alpha_damping / scale;
}));
}
}
CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
pip, "SIGMA", "U");
#endif // ADD_CONTACT
};
CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
// Boundary
CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto create_reaction_pipeline = [&](auto &pip) {
"GEOMETRY");
CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = vol_rule;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
reactionFe->postProcessHook =
}
//! [Push operators to pipeline]
//! [Solve]
struct SetUpSchur {
/**
* @brief Create data structure for handling Schur complement
*
* @param m_field
* @param sub_dm Schur complement sub dm
* @param field_split_it IS of Schur block
* @param ao_map AO map from sub dm to main problem
* @return boost::shared_ptr<SetUpSchur>
*/
static boost::shared_ptr<SetUpSchur> createSetUpSchur(
SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
);
virtual MoFEMErrorCode setUp(TS solver) = 0;
protected:
SetUpSchur() = default;
};
auto snes_ctx_ptr = getDMSnesCtx(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_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(mField);
auto &pip = pp_fe->getOpPtrVector();
auto push_vol_ops = [this](auto &pip) {
"GEOMETRY");
auto [common_plastic_ptr, common_henky_ptr] =
PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
if (common_henky_ptr) {
if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
}
return std::make_pair(common_plastic_ptr, common_henky_ptr);
};
auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
auto &pip = pp_fe->getOpPtrVector();
auto [common_plastic_ptr, common_henky_ptr] = p;
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
pip.push_back(
new OpPPMap(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()}},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{{"GRAD", common_plastic_ptr->mGradPtr},
{"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
{{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
} else {
pip.push_back(
new OpPPMap(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()}},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{},
{{"STRAIN", common_plastic_ptr->mStrainPtr},
{"STRESS", common_plastic_ptr->mStressPtr},
{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
}
};
auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
PetscBool post_proc_vol = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
&post_proc_vol, PETSC_NULL);
if (post_proc_vol == PETSC_FALSE)
return boost::shared_ptr<PostProcEle>();
auto pp_fe = boost::make_shared<PostProcEle>(mField);
push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
PetscBool post_proc_skin = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
&post_proc_skin, PETSC_NULL);
if (post_proc_skin == PETSC_FALSE)
return boost::shared_ptr<SkinPostProcEle>();
auto simple = mField.getInterface<Simple>();
auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
auto op_side = new OpLoopSide<DomainEle>(
mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
pp_fe->getOpPtrVector().push_back(op_side);
CHK_MOAB_THROW(push_vol_post_proc_ops(
pp_fe, push_vol_ops(op_side->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
return std::make_pair(vol_post_proc(), skin_post_proc());
};
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_elements(), reactionFe, uXScatter,
uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto set_schur_pc = [&](auto solver,
boost::shared_ptr<SetUpSchur> &schur_ptr) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
// create sub dm for Schur complement
auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,
SmartPetscObj<DM> &dm_sub) {
dm_sub = createDM(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);
};
// Create nested (sub BC) Schur DM
if constexpr (AT == AssemblyType::SCHUR) {
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
simple->getProblemName(), ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
simple->getProblemName(), 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);
#ifdef ADD_CONTACT
auto add_sigma_to_is = [&](auto is_union) {
SmartPetscObj<IS> is_union_sigma;
auto add_sigma_to_is_impl = [&]() {
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
simple->getProblemName(), ROW, "SIGMA", 0, MAX_DOFS_ON_ENTITY,
is_sigma);
IS is_union_raw_sigma;
CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);
is_union_sigma = SmartPetscObj<IS>(is_union_raw_sigma);
};
CHK_THROW_MESSAGE(add_sigma_to_is_impl(), "Can not add sigma to IS");
return is_union_sigma;
};
is_union = add_sigma_to_is(is_union);
#endif // ADD_CONTACT
CHKERR create_sub_u_dm(simple->getDM(), dm_u_sub);
// Indices has to be map fro very to level, while assembling Schur
// complement.
auto is_up = getDMSubData(dm_u_sub)->getSmartRowIs();
auto ao_up = createAOMappingIS(is_up, PETSC_NULL);
schur_ptr =
SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);
CHKERR schur_ptr->setUp(solver);
}
};
auto dm = simple->getDM();
auto D = createDMVector(dm);
auto DD = vectorDuplicate(D);
uXScatter = scatter_create(D, 0);
uYScatter = scatter_create(D, 1);
if constexpr (SPACE_DIM == 3)
uZScatter = scatter_create(D, 2);
auto create_solver = [pip_mng]() {
if (is_quasi_static == PETSC_TRUE)
return pip_mng->createTSIM();
else
return pip_mng->createTSIM2();
};
auto solver = create_solver();
auto active_pre_lhs = []() {
};
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);
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;
MOFEM_LOG_C("PLASTICITY", 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);
}
}
};
auto add_active_dofs_elem = [&](auto dm) {
auto fe_pre_proc = boost::make_shared<FEMethod>();
fe_pre_proc->preProcessHook = active_pre_lhs;
auto fe_post_proc = boost::make_shared<FEMethod>();
fe_post_proc->postProcessHook = active_post_lhs;
auto ts_ctx_ptr = getDMTsCtx(dm);
ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
};
auto set_essential_bc = [&](auto dm, auto solver) {
// This is low level pushing finite elements (pipelines) to solver
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
// Add boundary condition scaling
auto disp_time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
{disp_time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
auto ts_ctx_ptr = getDMTsCtx(dm);
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_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);
if (is_pcfs == PETSC_FALSE) {
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
}
};
if (is_quasi_static == PETSC_TRUE) {
CHKERR TSSetSolution(solver, D);
} else {
CHKERR TS2SetSolution(solver, D, DD);
}
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetUp(solver);
CHKERR add_active_dofs_elem(dm);
boost::shared_ptr<SetUpSchur> schur_ptr;
CHKERR set_schur_pc(solver, schur_ptr);
CHKERR set_essential_bc(dm, solver);
MOFEM_LOG_TAG("TIMER", "timer");
if (set_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";
}
//! [Solve]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
Py_Initialize();
#endif
#endif // ADD_CONTACT
// 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(), "PLASTICITY"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
LogManager::setLog("PLASTICITY");
MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
#ifdef ADD_CONTACT
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
MOFEM_LOG_TAG("CONTACT", "Contact");
#endif // ADD_CONTACT
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]
}
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif // ADD_CONTACT
return 0;
}
struct SetUpSchurImpl : public SetUpSchur {
SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
: SetUpSchur(), mField(m_field), subDM(sub_dm),
fieldSplitIS(field_split_is), aoUp(ao_up) {
if (S) {
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
}
virtual ~SetUpSchurImpl() {
#ifdef ADD_CONTACT
A.reset();
P.reset();
#endif // ADD_CONTACT
S.reset();
}
MoFEMErrorCode setUp(TS solver);
private:
#ifdef ADD_CONTACT
#endif // ADD_CONTACT
SmartPetscObj<DM> subDM; ///< field split sub dm
SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
SmartPetscObj<AO> aoUp; ///> main DM to subDM
};
auto pip_mng = mField.getInterface<PipelineManager>();
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPGetPC(ksp, &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");
}
#ifdef ADD_CONTACT
auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
A = createDMMatrix(simple->getDM());
P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
#endif // ADD_CONTACT
CHKERR MatSetBlockSize(S, SPACE_DIM);
auto set_ops = [&]() {
auto pip_mng = mField.getInterface<PipelineManager>();
#ifndef ADD_CONTACT
// Boundary
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
{"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
{SmartPetscObj<Mat>(), S}, {false, false}
));
// Domain
pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(
{"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
{SmartPetscObj<Mat>(), S}, {false, false}
));
#else
double eps_stab = 1e-4;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
PETSC_NULL);
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
// Boundary
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
return eps_stab;
}));
pip_mng->getOpBoundaryLhsPipeline().push_back(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
{false, false, false}
));
// Domain
pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
{false, false, false}
));
#endif // ADD_CONTACT
};
auto set_assemble_elems = [&]() {
auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
schur_asmb_pre_proc->preProcessHook = [this]() {
#ifdef ADD_CONTACT
CHKERR MatZeroEntries(A);
CHKERR MatZeroEntries(P);
#endif // ADD_CONTACT
CHKERR MatZeroEntries(S);
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
};
auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
#ifndef ADD_CONTACT
mField, schur_asmb_post_proc, 1)();
#else // ADD_CONTACT
CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
// Apply essential constrains to A matrix
mField, schur_asmb_post_proc, 1, A)();
CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
#endif // ADD_CONTACT
// Apply essential constrains to Schur complement
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
mField, schur_asmb_post_proc, 1, S, aoUp)();
};
auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
};
auto set_pc = [&]() {
CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
};
CHKERR set_ops();
CHKERR set_pc();
CHKERR set_assemble_elems();
} else {
pip_mng->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
pip_mng->getOpBoundaryLhsPipeline().push_back(
new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(
new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
}
// we do not those anymore
subDM.reset();
fieldSplitIS.reset();
// aoUp.reset();
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(
new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
}
static Index< 'p', 3 > p
std::string param_file
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#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
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ 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 int SPACE_DIM
Definition: elastic.cpp:17
const int dim
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
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:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
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.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
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:1761
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:1042
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
int r
Definition: sdf.py:8
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
Definition: plastic.cpp:165
constexpr AssemblyType AT
Definition: plastic.cpp:42
double C1_k
Kinematic hardening.
Definition: plastic.cpp:173
PostProcBrokenMeshInMoab< BoundaryEle > SkinPostProcEle
Definition: plastic.cpp:58
double Qinf
Saturation yield stress.
Definition: plastic.cpp:171
constexpr IntegrationType IT
Definition: plastic.cpp:45
NaturalBC< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT > BoundaryLhsBCs
Definition: plastic.cpp:211
NaturalBC< BoundaryEleOp >::Assembly< AT >::LinearForm< IT > BoundaryRhsBCs
Definition: plastic.cpp:208
double rho
Definition: plastic.cpp:182
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
PetscBool is_quasi_static
Definition: plastic.cpp:181
BoundaryLhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
Definition: plastic.cpp:213
double alpha_damping
Definition: plastic.cpp:183
double visH
Viscous hardening.
Definition: plastic.cpp:169
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:139
PetscBool set_timer
Set timer.
Definition: plastic.cpp:161
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:125
double scale
Definition: plastic.cpp:163
constexpr auto size_symm
Definition: plastic.cpp:40
NaturalBC< DomainEleOp >::Assembly< AT >::LinearForm< IT > DomainRhsBCs
Definition: plastic.cpp:205
double zeta
Viscous hardening.
Definition: plastic.cpp:170
double H
Hardening.
Definition: plastic.cpp:168
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:111
DomainRhsBCs::OpFlux< PlasticOps::DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
Definition: plastic.cpp:207
double cn0
Definition: plastic.cpp:175
double b_iso
Saturation exponent.
Definition: plastic.cpp:172
BoundaryRhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
Definition: plastic.cpp:210
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:160
int geom_order
Order if fixed.
Definition: plastic.cpp:179
double sigmaY
Yield stress.
Definition: plastic.cpp:167
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:120
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:149
double cn1
Definition: plastic.cpp:176
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:50
FTensor::Index< 'm', 3 > m
Element used to specialise assembly.
Definition: contact.cpp:86
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:240
[Example]
Definition: plastic.cpp:217
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:236
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:234
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:715
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:393
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:235
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:530
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:253
MoFEM::Interface & mField
Definition: plastic.cpp:224
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:265
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:485
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:232
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =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)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:47
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:55
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:63
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
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.
[Push operators to pipeline]
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
[Push operators to pipeline]
Definition: plastic.cpp:692
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1411
virtual MoFEMErrorCode setUp(TS solver)=0
SmartPetscObj< Mat > A
Definition: contact.cpp:844
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1233
SmartPetscObj< Mat > S
Definition: plastic.cpp:1230
MoFEMErrorCode setUp(TS solver)
Definition: plastic.cpp:1238
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1234
virtual ~SetUpSchurImpl()
Definition: plastic.cpp:1213
MoFEMErrorCode postProc()
SmartPetscObj< AO > aoUp
Definition: plastic.cpp:1235
SmartPetscObj< Mat > P
Definition: contact.cpp:845
MoFEMErrorCode preProc()
MoFEM::Interface & mField
Definition: plastic.cpp:1232