v0.14.0
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
*
*/
/* 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) {
if (C1_k < std::numeric_limits<double>::epsilon()) {
t_alpha(i, j) = 0;
return t_alpha;
}
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 displacement
int tau_order = order - 2; ///< Order of tau files
int ep_order = order - 1; ///< Order of ep files
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
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("PLASTICITY", 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);
PetscBool order_edge = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,
PETSC_NULL);
PetscBool order_face = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
PETSC_NULL);
PetscBool order_volume = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
PETSC_NULL);
if (order_edge || order_face || order_volume) {
MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
? "true"
: "false";
MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
? "true"
: "false";
MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
? "true"
: "false";
auto ents = get_ents_by_dim(0);
if (order_edge)
ents.merge(get_ents_by_dim(1));
if (order_face)
ents.merge(get_ents_by_dim(2));
if (order_volume)
ents.merge(get_ents_by_dim(3));
CHKERR simple->setFieldOrder("U", order, &ents);
} else {
CHKERR simple->setFieldOrder("U", order);
}
CHKERR simple->setFieldOrder("EP", ep_order);
CHKERR simple->setFieldOrder("TAU", tau_order);
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) {
bool is_contact_block = false;
Range contact_range;
for (auto m :
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true; ///< bloks interation is collectibe, so that is set irrespective
///< if there are enerities in given rank or not in the block
MOFEM_LOG("CONTACT", Sev::inform)
<< "Find contact block set: " << m->getName();
auto meshset = m->getMeshset();
Range contact_meshset_range;
CHKERR mField.get_moab().get_entities_by_dimension(
meshset, SPACE_DIM - 1, contact_meshset_range, true);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
contact_meshset_range);
contact_range.merge(contact_meshset_range);
}
if (is_contact_block) {
MOFEM_LOG("SYNC", Sev::inform)
<< "Nb entities in contact surface: " << contact_range.size();
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);
PetscBool tau_order_is_set; ///< true if tau order is set
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
&tau_order_is_set);
PetscBool ep_order_is_set; ///< true if tau order is set
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
&ep_order_is_set);
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;
if (tau_order_is_set == PETSC_FALSE)
if (ep_order_is_set == PETSC_FALSE)
ep_order = order - 1;
MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
MOFEM_LOG("PLASTICITY", Sev::inform)
<< "Ep approximation order " << ep_order;
MOFEM_LOG("PLASTICITY", Sev::inform)
<< "Tau approximation order " << tau_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
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
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");
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_henky_ptr->matGradPtr},
{"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<SideEle>(
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 order
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
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
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:1918
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:172
constexpr AssemblyType AT
Definition: plastic.cpp:44
double C1_k
Kinematic hardening.
Definition: plastic.cpp:180
double Qinf
Saturation yield stress.
Definition: plastic.cpp:178
constexpr IntegrationType IT
Definition: plastic.cpp:47
double rho
Definition: plastic.cpp:191
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
PetscBool is_quasi_static
Definition: plastic.cpp:190
double alpha_damping
Definition: plastic.cpp:192
double visH
Viscous hardening.
Definition: plastic.cpp:176
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:142
PetscBool set_timer
Set timer.
Definition: plastic.cpp:168
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:128
double scale
Definition: plastic.cpp:170
constexpr auto size_symm
Definition: plastic.cpp:42
double zeta
Viscous hardening.
Definition: plastic.cpp:177
double H
Hardening.
Definition: plastic.cpp:175
int tau_order
Order of tau files.
Definition: plastic.cpp:186
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:114
double cn0
Definition: plastic.cpp:182
double b_iso
Saturation exponent.
Definition: plastic.cpp:179
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:167
int geom_order
Order if fixed.
Definition: plastic.cpp:188
double sigmaY
Yield stress.
Definition: plastic.cpp:174
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:123
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:156
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
int ep_order
Order of ep files.
Definition: plastic.cpp:187
double cn1
Definition: plastic.cpp:183
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
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:249
[Example]
Definition: plastic.cpp:226
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:245
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:243
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:778
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:440
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:244
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:593
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:262
MoFEM::Interface & mField
Definition: plastic.cpp:233
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:274
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:548
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:241
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:755
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:1474
virtual MoFEMErrorCode setUp(TS solver)=0
SmartPetscObj< Mat > A
Definition: contact.cpp:854
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1296
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1297
virtual ~SetUpSchurImpl()
MoFEMErrorCode postProc()
SmartPetscObj< AO > aoUp
MoFEMErrorCode preProc()
MoFEM::Interface & mField
/** \file PlasticOps.hpp
* \example PlasticOps.hpp
\f[
\left\{
\begin{array}{ll}
\frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i} \right)\\
\sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
\dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
\sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
\sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
\Omega_u \cup \Omega_\sigma = \Omega \\
\Omega_u \cap \Omega_\sigma = \emptyset
\end{array}
\right.
\f]
\f[
\left\{
\begin{array}{ll}
\left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
\delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
\dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
& \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
\left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
(f(\pmb\sigma,\tau) - \sigma_y) +
\| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
\forall \delta\tau \in L^2(\Omega) \end{array} \right.
\f]
*/
namespace PlasticOps {
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
H,
};
using BlockParams = std::array<double, LAST_PARAM>;
inline auto getParamsPtr() {
return boost::shared_ptr<BlockParams>(shared_from_this(), &blockParams);
};
//! [Common data set externally]
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
//! [Common data set externally]
inline auto getPlasticSurfacePtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticSurface);
}
inline auto getPlasticTauPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTau);
}
inline auto getPlasticTauDotPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTauDot);
}
inline auto getPlasticStrainPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticStrain);
}
inline auto getPlasticStrainDotPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto getPlasticFlowPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticFlow);
}
static std::array<int, 5> activityData;
};
std::array<int, 5> CommonData::activityData = {0, 0, 0, 0, 0};
//! [Common data]
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticSurfaceImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticityImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpPlasticStressImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowRhsImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dTAUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dEPImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dTAUImpl;
template <typename DomainEleOp> struct PlasticityIntegrators {
template <int DIM, IntegrationType I>
OpCalculatePlasticSurfaceImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticity = OpCalculatePlasticityImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpPlasticStress = OpPlasticStressImpl<DIM, I, DomainEleOp>;
template <AssemblyType A> struct Assembly {
typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowRhsImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsRhsImpl<I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl<
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsLhs_dTAUImpl<I, AssemblyDomainEleOp>;
};
};
}; // namespace PlasticOps
namespace PlasticOps {
using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
using CommonPlasticPtr = boost::shared_ptr<PlasticOps::CommonData>;
using CommonHenkyPtr = boost::shared_ptr<HenckyOps::CommonData>;
template <int DIM>
addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
boost::shared_ptr<MatrixDouble> mat_D_Ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
double scale_value, Sev sev) {
struct OpMatBlocks : public DomainEleOp {
OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
double scale_value, MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m_D_ptr),
matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
auto getK = [](auto &p) {
auto young_modulus = p[CommonData::YOUNG_MODULUS];
auto poisson_ratio = p[CommonData::POISSON_RATIO];
return young_modulus / (3 * (1 - 2 * poisson_ratio));
};
auto getG = [](auto &p) {
auto young_modulus = p[CommonData::YOUNG_MODULUS];
auto poisson_ratio = p[CommonData::POISSON_RATIO];
return young_modulus / (2 * (1 + poisson_ratio));
};
auto scale_fun = [this](auto &p) {
p[CommonData::YOUNG_MODULUS] *= scaleVal;
p[CommonData::SIGMA_Y] *= scaleVal;
p[CommonData::H] *= scaleVal;
p[CommonData::VIS_H] *= scaleVal;
p[CommonData::QINF] *= scaleVal;
p[CommonData::C1_k] *= scaleVal;
};
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*matParamsPtr = b.bParams;
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
}
(*matParamsPtr) = {young_modulus, poisson_ratio, sigmaY, H,
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
const double scaleVal;
struct BlockData {
std::array<double, CommonData::LAST_PARAM> bParams;
Range blockEnts;
};
std::vector<BlockData> blockData;
/**
* @brief Extract block data from meshsets
*
* @param m_field
* @param meshset_vec_ptr
* @param sev
* @return MoFEMErrorCode
*/
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != CommonData::LAST_PARAM) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of block attribute");
}
auto get_block_ents = [&]() {
Range ents;
CHKERR m_field.get_moab().get_entities_by_handle(m->meshset, ents,
true);
return ents;
};
CommonData::BlockParams block_params;
for (auto i = 0; i != CommonData::LAST_PARAM; ++i) {
block_params[i] = block_data[i];
}
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< "E = " << block_params[CommonData::YOUNG_MODULUS]
<< " nu = " << block_params[CommonData::POISSON_RATIO];
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< std::endl
<< "sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
<< "h = " << block_params[CommonData::H] << std::endl
<< "vis_h = " << block_params[CommonData::VIS_H] << std::endl
<< "qinf = " << block_params[CommonData::QINF] << std::endl
<< "biso = " << block_params[CommonData::BISO] << std::endl
<< "C1_k = " << block_params[CommonData::C1_k] << std::endl;
blockData.push_back({block_params, get_block_ents()});
}
}
/**
* @brief Get elasticity tensor
*
* Calculate elasticity tensor for given material parameters
*
* @param mat_D_ptr
* @param bulk_modulus_K
* @param shear_modulus_G
* @return MoFEMErrorCode
*
*/
MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
double bulk_modulus_K, double shear_modulus_G) {
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
double A = (DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
t_D(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
t_kd(k, l);
};
//! [Calculate elasticity tensor]
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
mat_D_ptr->resize(size_symm * size_symm, 1);
set_material_stiffness();
}
};
// push operator to calculate material stiffness matrix for each block
pip.push_back(new OpMatBlocks(
mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
// Get blockset using regular expression
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
));
}
template <int DIM, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field, std::string block_name,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string u, std::string ep, std::string tau, double scale, Sev sev) {
using P = PlasticityIntegrators<DomainEleOp>;
auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
auto make_d_mat = []() {
return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
};
common_plastic_ptr->mDPtr = make_d_mat();
common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
auto m_D_ptr = common_plastic_ptr->mDPtr;
CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, block_name, pip, m_D_ptr,
common_plastic_ptr->getParamsPtr(),
scale, sev),
"add mat block ops");
pip.push_back(new OpCalculateScalarFieldValues(
tau, common_plastic_ptr->getPlasticTauPtr()));
pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
ep, common_plastic_ptr->getPlasticStrainPtr()));
u, common_plastic_ptr->mGradPtr));
CommonHenkyPtr common_henky_ptr;
common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
common_henky_ptr->matLogCPlastic =
common_plastic_ptr->getPlasticStrainPtr();
common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
u, common_henky_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
u, common_henky_ptr));
pip.push_back(new
typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
u, common_henky_ptr));
} else {
pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
u, common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
u, common_plastic_ptr));
return std::make_tuple(common_plastic_ptr, common_henky_ptr);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, scale, Sev::inform);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
// Calculate internal forces
if (common_henky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
DIM, I>(ep, common_plastic_ptr, m_D_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
using namespace HenckyOps;
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template BiLinearForm<I>;
using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, scale, Sev::verbose);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
if (common_henky_ptr) {
pip.push_back(new typename H::template OpHenckyTangent<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(new OpKCauchy(u, u, m_D_ptr));
pip.push_back(new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
u, ep, common_plastic_ptr, m_D_ptr));
}
if (common_henky_ptr) {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
DIM, I>(tau, u, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
DIM, I>(ep, u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
DIM, I>(ep, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
DIM, I>(ep, tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
DIM, I>(tau, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
I>(tau, tau, common_plastic_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string block_name, Pip &pip,
std::string u, std::string ep,
std::string tau) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, 1., Sev::inform);
// Calculate internal forces
if (common_henky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
}
}
} // namespace PlasticOps
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
@ NOSPACE
Definition: definitions.h:83
double bulk_modulus_K
double shear_modulus_G
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
Definition: PlasticOps.hpp:426
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:630
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:500
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
Definition: PlasticOps.hpp:242
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
Definition: PlasticOps.hpp:248
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
Definition: PlasticOps.hpp:243
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:118
FTensor::Index< 'I', 3 > I
[Common data]
Definition: PlasticOps.hpp:115
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:547
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:116
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
Definition: PlasticOps.hpp:244
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr AssemblyType A
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: PlasticOps.hpp:67
MatrixDouble resFlowDstrainDot
Definition: PlasticOps.hpp:85
VectorDouble resCdTau
Definition: PlasticOps.hpp:79
MatrixDouble resCdStrain
Definition: PlasticOps.hpp:80
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: PlasticOps.hpp:66
MatrixDouble resCdPlasticStrain
Definition: PlasticOps.hpp:81
MatrixDouble resFlow
Definition: PlasticOps.hpp:82
VectorDouble plasticTauDot
Definition: PlasticOps.hpp:74
BlockParams blockParams
Definition: PlasticOps.hpp:58
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
Definition: PlasticOps.hpp:65
VectorDouble plasticSurface
[Common data set externally]
Definition: PlasticOps.hpp:71
MatrixDouble resFlowDtau
Definition: PlasticOps.hpp:83
MatrixDouble plasticStrain
Definition: PlasticOps.hpp:75
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: PlasticOps.hpp:68
MatrixDouble plasticFlow
Definition: PlasticOps.hpp:72
MatrixDouble plasticStrainDot
Definition: PlasticOps.hpp:76
VectorDouble plasticTau
Definition: PlasticOps.hpp:73
MatrixDouble resFlowDstrain
Definition: PlasticOps.hpp:84
std::array< double, LAST_PARAM > BlockParams
Definition: PlasticOps.hpp:57
OpCalculatePlasticFlowRhsImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowRhs
Definition: PlasticOps.hpp:183
OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_LogStrain_dU
Definition: PlasticOps.hpp:204
OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_LogStrain_dU
Definition: PlasticOps.hpp:220
OpCalculatePlasticFlowLhs_dTAUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dTAU
Definition: PlasticOps.hpp:212
OpCalculateConstraintsLhs_dTAUImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dTAU
Definition: PlasticOps.hpp:228
OpCalculatePlasticInternalForceLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_dEP
Definition: PlasticOps.hpp:191
OpCalculatePlasticFlowLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dU
Definition: PlasticOps.hpp:200
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_LogStrain_dEP
Definition: PlasticOps.hpp:196
OpCalculatePlasticFlowLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dEP
Definition: PlasticOps.hpp:208
OpCalculateConstraintsRhsImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsRhs
Definition: PlasticOps.hpp:187
OpCalculateConstraintsLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dU
Definition: PlasticOps.hpp:216
OpCalculateConstraintsLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dEP
Definition: PlasticOps.hpp:224
OpPlasticStressImpl< DIM, I, DomainEleOp > OpPlasticStress
Definition: PlasticOps.hpp:174
OpCalculatePlasticSurfaceImpl< DIM, I, DomainEleOp > OpCalculatePlasticSurface
Definition: PlasticOps.hpp:168
OpCalculatePlasticityImpl< DIM, I, DomainEleOp > OpCalculatePlasticity
Definition: PlasticOps.hpp:171
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]