v0.14.0
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
// #undef ADD_CONTACT
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
#include <IntegrationRules.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps;
template <> struct ElementsAndOps<2> {
using DomainEle = PipelineManager::FaceEle;
using BoundaryEle = PipelineManager::EdgeEle;
using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
template <> struct ElementsAndOps<3> {
using DomainEle = PipelineManager::VolEle;
using BoundaryEle = PipelineManager::FaceEle;
using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
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::BLOCK_SCHUR
: AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
IntegrationType::GAUSS; //< selected integration type
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
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
kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,
double C1_k) {
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Tensor2_symmetric<double, DIM> t_alpha;
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>
inline auto kinematic_hardening_dplastic_strain(double C1_k) {
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
FTensor::Ddg<double, DIM, DIM> t_diff;
constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
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>
#include <PlasticNaturalBCs.hpp>
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
#include <ContactOps.hpp>
#endif // ADD_CONTACT
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
using OpDomainRhsBCs =
DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
using OpBoundaryRhsBCs =
BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
using OpBoundaryLhsBCs =
BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
using namespace PlasticOps;
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField;
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode bC();
MoFEMErrorCode OPs();
MoFEMErrorCode tsSolve();
MoFEMErrorCode testOperators();
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 {
using MoFEM::TimeScale::TimeScale;
double getScale(const double time) {
return scale * MoFEM::TimeScale::getScale(time);
};
};
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif // ADD_CONTACT
};
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
MoFEMFunctionBegin;
CHKERR createCommonData();
CHKERR setupProblem();
CHKERR bC();
CHKERR OPs();
PetscBool test_ops = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_operators", &test_ops,
PETSC_NULL);
if (test_ops == PETSC_FALSE) {
CHKERR tsSolve();
} else {
CHKERR testOperators();
}
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
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())
CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
const auto type = type_from_handle(domain_ents[0]);
switch (type) {
case MBQUAD:
return DEMKOWICZ_JACOBI_BASE;
case MBHEX:
return DEMKOWICZ_JACOBI_BASE;
case MBTRI:
return AINSWORTH_LEGENDRE_BASE;
case MBTET:
return AINSWORTH_LEGENDRE_BASE;
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,
SPACE_DIM);
CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
SPACE_DIM);
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 = true;
Range contact_range;
for (auto m :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true; ///< blocs interation is collective, so that is set irrespective
///< if there are entities 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();
MOFEM_LOG_SYNCHRONISE(mField.get_comm());
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);
};
PetscBool project_geometry = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
&project_geometry, PETSC_NULL);
if (project_geometry) {
CHKERR project_ho_geometry();
}
MoFEMFunctionReturn(0);
}
//! [Set up problem]
//! [Create common data]
MoFEMErrorCode Example::createCommonData() {
MoFEMFunctionBegin;
auto get_command_line_parameters = [&]() {
MoFEMFunctionBegin;
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)
tau_order = order - 2;
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) {
scale /= young_modulus;
}
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");
MoFEMFunctionReturn(0);
};
CHKERR get_command_line_parameters();
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdfPythonPtr->sdfInit("sdf.py");
ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif // ADD_CONTACT
MoFEMFunctionReturn(0);
}
//! [Create common data]
//! [Boundary condition]
MoFEMErrorCode Example::bC() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
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;
MoFEMFunctionReturn(0);
}
//! [Boundary condition]
//! [Push operators to pipeline]
MoFEMErrorCode Example::OPs() {
MoFEMFunctionBegin;
auto pip_mng = mField.getInterface<PipelineManager>();
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) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
"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
auto simple = mField.getInterface<Simple>();
CHKERR
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U");
CHKERR
ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
vol_rule);
#endif // ADD_CONTACT
MoFEMFunctionReturn(0);
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
"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
MoFEMFunctionReturn(0);
};
auto add_domain_ops_lhs = [this](auto &pip) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
"GEOMETRY");
if (is_quasi_static == PETSC_FALSE) {
//! [Only used for dynamics]
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
GAUSS>::OpMass<1, SPACE_DIM>;
//! [Only used for dynamics]
auto get_inertia_and_mass_damping = [this](const double, const double,
const double) {
auto *pip = mField.getInterface<PipelineManager>();
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");
MoFEMFunctionReturn(0);
};
auto add_domain_ops_rhs = [this](auto &pip) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
"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]
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Only used for dynamics]
auto mat_acceleration = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
"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(
new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
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
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
"GEOMETRY");
CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
MoFEMFunctionReturn(0);
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = vol_rule;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
reactionFe->postProcessHook =
EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
MoFEMFunctionReturn(0);
}
//! [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(
MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
);
virtual MoFEMErrorCode setUp(TS solver) = 0;
protected:
SetUpSchur() = default;
};
MoFEMErrorCode Example::tsSolve() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
auto set_section_monitor = [&](auto solver) {
MoFEMFunctionBegin;
SNES snes;
CHKERR TSGetSNES(solver, &snes);
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
MoFEMFunctionReturn(0);
};
auto create_post_process_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(mField);
auto push_vol_ops = [this](auto &pip) {
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
"GEOMETRY");
auto [common_plastic_ptr, common_hencky_ptr] =
PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
if (common_hencky_ptr) {
if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
}
return std::make_pair(common_plastic_ptr, common_hencky_ptr);
};
auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
MoFEMFunctionBegin;
auto &pip = pp_fe->getOpPtrVector();
auto [common_plastic_ptr, common_hencky_ptr] = p;
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
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));
if (is_large_strains) {
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_hencky_ptr->matGradPtr},
{"FIRST_PIOLA", common_hencky_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()}}
)
);
}
MoFEMFunctionReturn(0);
};
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);
CHK_MOAB_THROW(
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) {
SmartPetscObj<IS> is;
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),
SmartPetscObj<VecScatter>(scatter));
};
auto set_time_monitor = [&](auto dm, auto solver) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
auto set_schur_pc = [&](auto solver,
boost::shared_ptr<SetUpSchur> &schur_ptr) {
MoFEMFunctionBeginHot;
auto name_prb = simple->getProblemName();
// create sub dm for Schur complement
auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
SmartPetscObj<DM> &dm_sub) {
MoFEMFunctionBegin;
dm_sub = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
for (auto f : {"U"}) {
CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
}
CHKERR DMSetUp(dm_sub);
MoFEMFunctionReturn(0);
};
auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
SmartPetscObj<DM> &dm_sub) {
MoFEMFunctionBegin;
dm_sub = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
#ifdef ADD_CONTACT
for (auto f : {"SIGMA", "EP", "TAU"}) {
CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
}
#else
for (auto f : {"EP", "TAU"}) {
CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
}
#endif
CHKERR DMSetUp(dm_sub);
MoFEMFunctionReturn(0);
};
// Create nested (sub BC) Schur DM
if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
SmartPetscObj<DM> dm_schur;
CHKERR create_schur_dm(simple->getDM(), dm_schur);
SmartPetscObj<DM> dm_block;
CHKERR create_block_dm(simple->getDM(), dm_block);
#ifdef ADD_CONTACT
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
auto block_mat_data = createBlockMatStructure(
simple->getDM(),
{
{simple->getDomainFEName(),
{{"U", "U"},
{"SIGMA", "SIGMA"},
{"U", "SIGMA"},
{"SIGMA", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "EP"},
{"TAU", "U"}
}},
{simple->getBoundaryFEName(),
{{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
}}
}
);
return createSchurNestedMatrixStruture(
{dm_schur, dm_block}, block_mat_data,
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
);
};
#else
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
auto block_mat_data =
createBlockMatStructure(simple->getDM(),
{{simple->getDomainFEName(),
{{"U", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "U"},
{"TAU", "EP"}
}}}
);
return createSchurNestedMatrixStruture(
{dm_schur, dm_block}, block_mat_data,
{"EP", "TAU"}, {nullptr, nullptr}, false
);
};
#endif
auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
auto block_is = getDMSubData(dm_block)->getSmartRowIs();
auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
// Indices has to be map fro very to level, while assembling Schur
// complement.
schur_ptr =
SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
CHKERR schur_ptr->setUp(solver);
}
MoFEMFunctionReturnHot(0);
};
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 = []() {
MoFEMFunctionBegin;
std::fill(PlasticOps::CommonData::activityData.begin(),
PlasticOps::CommonData::activityData.end(), 0);
MoFEMFunctionReturn(0);
};
auto active_post_lhs = [&]() {
MoFEMFunctionBegin;
auto get_iter = [&]() {
SNES snes;
CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
int iter;
CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
"Can not get iter");
return iter;
};
auto iter = get_iter();
if (iter >= 0) {
std::array<int, 5> activity_data;
std::fill(activity_data.begin(), activity_data.end(), 0);
MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
activity_data.data(), activity_data.size(), MPI_INT,
MPI_SUM, mField.get_comm());
int &active_points = activity_data[0];
int &avtive_full_elems = activity_data[1];
int &avtive_elems = activity_data[2];
int &nb_points = activity_data[3];
int &nb_elements = activity_data[4];
if (nb_points) {
double proc_nb_points =
100 * static_cast<double>(active_points) / nb_points;
double proc_nb_active =
100 * static_cast<double>(avtive_elems) / nb_elements;
double proc_nb_full_active = 100;
if (avtive_elems)
proc_nb_full_active =
100 * static_cast<double>(avtive_full_elems) / avtive_elems;
MOFEM_LOG_C("PLASTICITY", Sev::inform,
"Iter %d nb pts %d nb active 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);
}
}
MoFEMFunctionReturn(0);
};
auto add_active_dofs_elem = [&](auto dm) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
auto set_essential_bc = [&](auto dm, auto solver) {
MoFEMFunctionBegin;
// 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]() {
MoFEMFunctionBegin;
CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
mField, post_proc_rhs_ptr, 1.)();
MoFEMFunctionReturn(0);
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
return EssentialPostProcLhs<DisplacementCubitBcData>(
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);
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
MoFEMFunctionReturn(0);
};
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_CHANNEL("TIMER");
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";
MoFEMFunctionReturn(0);
}
//! [Solve]
//! [TestOperators]
MoFEMErrorCode Example::testOperators() {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto opt = mField.getInterface<OperatorsTester>(); // get interface to
// OperatorsTester
auto pip = mField.getInterface<PipelineManager>(); // get interface to
// pipeline manager
constexpr double eps = 1e-9;
auto x = opt->setRandomFields(simple->getDM(), {
{"U", {-1e-6, 1e-6}},
{"EP", {-1e-6, 1e-6}},
{"TAU", {0, 1e-4}}
});
auto dot_x_plastic_active =
opt->setRandomFields(simple->getDM(), {
{"U", {-1, 1}},
{"EP", {-1, 1}},
{"TAU", {0.1, 1}}
});
auto diff_x_plastic_active =
opt->setRandomFields(simple->getDM(), {
{"U", {-1, 1}},
{"EP", {-1, 1}},
{"TAU", {-1, 1}}
});
auto dot_x_elastic =
opt->setRandomFields(simple->getDM(), {
{"U", {-1, 1}},
{"EP", {-1, 1}},
{"TAU", {-1, -0.1}}
});
auto diff_x_elastic =
opt->setRandomFields(simple->getDM(), {
{"U", {-1, 1}},
{"EP", {-1, 1}},
{"TAU", {-1, 1}}
});
auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
auto dot_x, auto diff_x) {
MoFEMFunctionBegin;
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
// Calculate norm of difference between directional derivative calculated
// from finite difference, and tangent matrix.
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C("PLASTICITY", Sev::inform,
"Test consistency of tangent matrix %3.4e", fnorm);
constexpr double err = 1e-5;
if (fnorm > err)
SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Norm of directional derivative too large err = %3.4e", fnorm);
MoFEMFunctionReturn(0);
};
MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
pip->getDomainRhsFE(), dot_x_plastic_active,
diff_x_plastic_active);
MoFEMFunctionReturn(0);
};
//! [TestOperators]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
Py_Initialize();
np::initialize();
#endif
#endif // ADD_CONTACT
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// 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";
CHKERR DMRegister_MoFEM(dm_name);
//! [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]
Simple *simple = m_field.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif // ADD_CONTACT
return 0;
}
struct SetUpSchurImpl : public SetUpSchur {
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
: SetUpSchur(), mField(m_field), subDM(sub_dm),
fieldSplitIS(field_split_is), aoSchur(ao_up) {
if (S) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Is expected that schur matrix is not "
"allocated. This is "
"possible only is PC is set up twice");
}
}
virtual ~SetUpSchurImpl() { S.reset(); }
MoFEMErrorCode setUp(TS solver);
MoFEMErrorCode preProc();
MoFEMErrorCode postProc();
private:
SmartPetscObj<Mat> S;
MoFEM::Interface &mField;
SmartPetscObj<DM> subDM; ///< field split sub dm
SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
SmartPetscObj<AO> aoSchur; ///> main DM to subDM
};
MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto pip_mng = mField.getInterface<PipelineManager>();
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
if (S) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Is expected that schur matrix is not "
"allocated. This is "
"possible only is PC is set up twice");
}
S = createDMMatrix(subDM);
CHKERR MatSetBlockSize(S, SPACE_DIM);
// Set DM to use shell block matrix
DM solver_dm;
CHKERR TSGetDM(solver, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
auto ts_ctx_ptr = getDMTsCtx(solver_dm);
auto A = createDMBlockMat(simple->getDM());
auto P = createDMNestSchurMat(simple->getDM());
if (is_quasi_static == PETSC_TRUE) {
auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
Mat A, Mat B, void *ctx) {
return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
};
CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
} else {
auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
PetscReal a, PetscReal aa, Mat A, Mat B,
void *ctx) {
return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
};
CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
}
CHKERR KSPSetOperators(ksp, A, P);
auto set_ops = [&]() {
MoFEMFunctionBegin;
auto pip_mng = mField.getInterface<PipelineManager>();
#ifndef ADD_CONTACT
// Boundary
pip_mng->getOpBoundaryLhsPipeline().push_front(
createOpSchurAssembleBegin());
pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
{"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
));
// Domain
pip_mng->getOpDomainLhsPipeline().push_front(
createOpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
{"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
));
#else
double eps_stab = 1e-4;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
PETSC_NULL);
using B = FormsIntegrators<BoundaryEleOp>::Assembly<
BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
// Boundary
pip_mng->getOpBoundaryLhsPipeline().push_front(
createOpSchurAssembleBegin());
pip_mng->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
return eps_stab;
}));
pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
false, false
));
// Domain
pip_mng->getOpDomainLhsPipeline().push_front(
createOpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
false, false
));
#endif // ADD_CONTACT
MoFEMFunctionReturn(0);
};
auto set_assemble_elems = [&]() {
MoFEMFunctionBegin;
auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
schur_asmb_pre_proc->preProcessHook = [this]() {
MoFEMFunctionBegin;
CHKERR MatZeroEntries(S);
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
MoFEMFunctionReturn(0);
};
auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
MoFEMFunctionBegin;
MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
// Apply essential constrains to Schur complement
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
mField, schur_asmb_post_proc, 1, S, aoSchur)();
MoFEMFunctionReturn(0);
};
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);
MoFEMFunctionReturn(0);
};
auto set_pc = [&]() {
MoFEMFunctionBegin;
CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
MoFEMFunctionReturn(0);
};
auto set_diagonal_pc = [&]() {
MoFEMFunctionBegin;
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
return SmartPetscObj<PC>(pc_raw,
true); // bump reference
};
CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
CHKERR PetscFree(subksp);
MoFEMFunctionReturn(0);
};
CHKERR set_ops();
CHKERR set_pc();
CHKERR set_assemble_elems();
CHKERR TSSetUp(solver);
CHKERR KSPSetUp(ksp);
CHKERR set_diagonal_pc();
} else {
pip_mng->getOpBoundaryLhsPipeline().push_front(
createOpSchurAssembleBegin());
pip_mng->getOpBoundaryLhsPipeline().push_back(
createOpSchurAssembleEnd({}, {}));
pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
pip_mng->getOpDomainLhsPipeline().push_back(
createOpSchurAssembleEnd({}, {}));
}
// fieldSplitIS.reset();
// aoSchur.reset();
MoFEMFunctionReturn(0);
}
boost::shared_ptr<SetUpSchur>
SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
SmartPetscObj<AO> ao_up) {
return boost::shared_ptr<SetUpSchur>(
new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
}
/** \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> {
enum ParamsIndexes {
YOUNG_MODULUS,
POISSON_RATIO,
SIGMA_Y,
H,
VIS_H,
QINF,
BISO,
C1_k,
LAST_PARAM
};
using BlockParams = std::array<double, LAST_PARAM>;
BlockParams blockParams;
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]
VectorDouble plasticSurface;
MatrixDouble plasticFlow;
VectorDouble plasticTau;
VectorDouble plasticTauDot;
MatrixDouble plasticStrain;
MatrixDouble plasticStrainDot;
VectorDouble resC;
VectorDouble resCdTau;
MatrixDouble resCdStrain;
MatrixDouble resCdPlasticStrain;
MatrixDouble resFlow;
MatrixDouble resFlowDtau;
MatrixDouble resFlowDstrain;
MatrixDouble resFlowDstrainDot;
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(),
&plasticStrainDot);
}
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]
FTensor::Index<'I', 3> I;
FTensor::Index<'J', 3> J;
FTensor::Index<'M', 3> M;
FTensor::Index<'N', 3> N;
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>
using OpCalculatePlasticSurface =
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 {
using AssemblyDomainEleOp =
typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowRhs =
OpCalculatePlasticFlowRhsImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
using OpCalculateConstraintsRhs =
OpCalculateConstraintsRhsImpl<I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticInternalForceLhs_dEP =
OpCalculatePlasticInternalForceLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticInternalForceLhs_LogStrain_dEP =
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl<
DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowLhs_dU =
OpCalculatePlasticFlowLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowLhs_LogStrain_dU =
OpCalculatePlasticFlowLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowLhs_dEP =
OpCalculatePlasticFlowLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowLhs_dTAU =
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculateConstraintsLhs_dU =
OpCalculateConstraintsLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculateConstraintsLhs_LogStrain_dU =
OpCalculateConstraintsLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculateConstraintsLhs_dEP =
OpCalculateConstraintsLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
using OpCalculateConstraintsLhs_dTAU =
OpCalculateConstraintsLhs_dTAUImpl<I, AssemblyDomainEleOp>;
};
};
}; // namespace PlasticOps
#include <PlasticOpsGeneric.hpp>
#include <PlasticOpsSmallStrains.hpp>
#include <PlasticOpsLargeStrains.hpp>
#include <PlasticOpsMonitor.hpp>
namespace PlasticOps {
using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
using CommonPlasticPtr = boost::shared_ptr<PlasticOps::CommonData>;
using CommonHenckyPtr = boost::shared_ptr<HenckyOps::CommonData>;
template <int DIM>
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) {
MoFEMFunctionBegin;
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) {
MoFEMFunctionBegin;
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));
MoFEMFunctionReturnHot(0);
}
}
(*matParamsPtr) = {young_modulus, poisson_ratio, sigmaY, H,
visH, Qinf, b_iso, C1_k};
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
MoFEMFunctionReturn(0);
}
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
*/
MoFEMErrorCode
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
MoFEMFunctionBegin;
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()});
}
MOFEM_LOG_CHANNEL("WORLD");
MoFEMFunctionReturn(0);
}
/**
* @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) {
MoFEMFunctionBegin;
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
double A = (DIM == 2)
? 2 * shear_modulus_G /
(bulk_modulus_K + (4. / 3.) * shear_modulus_G)
: 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();
MoFEMFunctionReturn(0);
}
};
// 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()
))
));
MoFEMFunctionReturn(0);
}
template <int DIM, IntegrationType I, typename DomainEleOp>
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) {
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()));
pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
u, common_plastic_ptr->mGradPtr));
CommonHenckyPtr common_hencky_ptr;
if (is_large_strains) {
common_hencky_ptr = boost::make_shared<HenckyOps::CommonData>();
common_hencky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
common_hencky_ptr->matDPtr = common_plastic_ptr->mDPtr;
common_hencky_ptr->matLogCPlastic =
common_plastic_ptr->getPlasticStrainPtr();
common_plastic_ptr->mStrainPtr = common_hencky_ptr->getMatLogC();
common_plastic_ptr->mStressPtr = common_hencky_ptr->getMatHenckyStress();
using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
u, common_hencky_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(u, common_hencky_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
u, common_hencky_ptr));
pip.push_back(
new typename H::template OpCalculateHenckyPlasticStress<DIM, I, 0>(
u, common_hencky_ptr, m_D_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
u, common_hencky_ptr));
} else {
pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
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_hencky_ptr);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode
opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
MoFEMFunctionBegin;
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
using OpInternalForceCauchy =
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
using OpInternalForcePiola =
typename B::template OpGradTimesTensor<1, DIM, DIM>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_hencky_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()));
pip.push_back(new OpCalculateScalarFieldValuesDot(
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_hencky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_hencky_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));
MoFEMFunctionReturn(0);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode
opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
std::string u, std::string ep, std::string tau) {
MoFEMFunctionBegin;
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_hencky_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()));
pip.push_back(new OpCalculateScalarFieldValuesDot(
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
if (common_hencky_ptr) {
using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
u, common_hencky_ptr, m_D_ptr));
pip.push_back(new OpKPiola(u, u, common_hencky_ptr->getMatTangent()));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
u, ep, common_plastic_ptr, common_hencky_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_hencky_ptr) {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
tau, u, common_plastic_ptr, common_hencky_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_hencky_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));
MoFEMFunctionReturn(0);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field,
std::string block_name, Pip &pip,
std::string u, std::string ep,
std::string tau) {
MoFEMFunctionBegin;
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
using OpInternalForceCauchy =
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
using OpInternalForcePiola =
typename B::template OpGradTimesTensor<1, DIM, DIM>;
auto [common_plastic_ptr, common_hencky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, 1., Sev::inform);
// Calculate internal forces
if (common_hencky_ptr) {
pip.push_back(new OpInternalForcePiola(
u, common_hencky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
}
MoFEMFunctionReturn(0);
}
} // namespace PlasticOps
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
PlasticOps::CommonData::BISO
@ BISO
Definition: PlasticOps.hpp:52
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
PlasticOps::CommonHenckyPtr
boost::shared_ptr< HenckyOps::CommonData > CommonHenckyPtr
Definition: PlasticOps.hpp:244
Example::tsSolve
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:740
SetUpSchurImpl
Definition: test_broken_space.cpp:515
PlasticOps::CommonData::resCdStrain
MatrixDouble resCdStrain
Definition: PlasticOps.hpp:80
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
set_timer
PetscBool set_timer
Set timer.
Definition: plastic.cpp:117
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:157
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Example::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:556
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticInternalForceLhs_LogStrain_dEP
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_LogStrain_dEP
Definition: PlasticOps.hpp:196
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
PlasticOps::J
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:116
SetUpSchurImpl::aoSchur
SmartPetscObj< AO > aoSchur
Definition: plastic.cpp:1416
PlasticOps::CommonData::BlockParams
std::array< double, LAST_PARAM > BlockParams
Definition: PlasticOps.hpp:57
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
ContactOps
Definition: contact.cpp:99
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:528
rho
double rho
Definition: plastic.cpp:140
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:123
PlasticOps::PlasticityIntegrators::Assembly::OpCalculateConstraintsLhs_dU
OpCalculateConstraintsLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dU
Definition: PlasticOps.hpp:216
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
PlasticOps::CommonData::LAST_PARAM
@ LAST_PARAM
Definition: PlasticOps.hpp:54
PlasticOps::opFactoryDomainRhs
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
MoFEM::BLOCK_PRECONDITIONER_SCHUR
@ BLOCK_PRECONDITIONER_SCHUR
Definition: FormsIntegrators.hpp:109
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
NOBASE
@ NOBASE
Definition: definitions.h:59
PlasticOps::CommonData::mStressPtr
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: PlasticOps.hpp:68
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
PlasticOps::CommonData::C1_k
@ C1_k
Definition: PlasticOps.hpp:53
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
PlasticOps::PlasticityIntegrators::Assembly::OpCalculateConstraintsLhs_LogStrain_dU
OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_LogStrain_dU
Definition: PlasticOps.hpp:220
PlasticOps::CommonData::resC
VectorDouble resC
Definition: PlasticOps.hpp:78
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PlasticOps::CommonData::getPlasticStrainPtr
auto getPlasticStrainPtr()
Definition: PlasticOps.hpp:96
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::TsSetI2Jacobian
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition: TsCtx.cpp:511
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
iso_hardening_dtau
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:77
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
PlasticOps::CommonData::plasticTau
VectorDouble plasticTau
Definition: PlasticOps.hpp:73
HenckyOps
Definition: HenckyOps.hpp:11
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
PlasticOps::CommonPlasticPtr
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
Definition: PlasticOps.hpp:243
scale
double scale
Definition: plastic.cpp:119
PlasticOps::createCommonPlasticOps
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
HenckyOps.hpp
PlasticOps::opFactoryDomainLhs
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
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:137
PlasticOps::CommonData::getPlasticFlowPtr
auto getPlasticFlowPtr()
Definition: PlasticOps.hpp:103
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1414
PlasticOps::CommonData::getParamsPtr
auto getParamsPtr()
Definition: PlasticOps.hpp:60
sdf.r
int r
Definition: sdf.py:8
PlasticOpsGeneric.hpp
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition: DMMoFEM.hpp:1076
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
PlasticOps::N
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:118
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
iso_hardening_exp
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:63
PlasticOps::CommonData::getPlasticSurfacePtr
auto getPlasticSurfacePtr()
Definition: PlasticOps.hpp:87
MoFEM::PipelineManager::EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: PipelineManager.hpp:36
PlasticOps::CommonData::mGradPtr
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: PlasticOps.hpp:66
PlasticOps::Pip
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
Definition: PlasticOps.hpp:242
PlasticOps::PlasticityIntegrators::OpCalculatePlasticity
OpCalculatePlasticityImpl< DIM, I, DomainEleOp > OpCalculatePlasticity
Definition: PlasticOps.hpp:171
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
PlasticOps.hpp
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
PlasticOps::CommonData::H
@ H
Definition: PlasticOps.hpp:49
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2186
AT
constexpr AssemblyType AT
Definition: test_broken_space.cpp:20
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
PlasticOps::CommonData::getPlasticStrainDotPtr
auto getPlasticStrainDotPtr()
Definition: PlasticOps.hpp:99
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:525
PlasticOps::CommonData::blockParams
BlockParams blockParams
Definition: PlasticOps.hpp:58
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticFlowLhs_dU
OpCalculatePlasticFlowLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dU
Definition: PlasticOps.hpp:200
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticInternalForceLhs_dEP
OpCalculatePlasticInternalForceLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_dEP
Definition: PlasticOps.hpp:191
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
BlockData
Definition: ElasticityMixedFormulation.hpp:12
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1415
BoundaryEleOp
PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
kinematic_hardening
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:91
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticFlowLhs_LogStrain_dU
OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_LogStrain_dU
Definition: PlasticOps.hpp:204
is_quasi_static
PetscBool is_quasi_static
Definition: plastic.cpp:139
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:177
MoFEM::DMMoFEMCreateSubDM
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:215
double
convert.type
type
Definition: convert.py:64
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticFlowLhs_dEP
OpCalculatePlasticFlowLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dEP
Definition: PlasticOps.hpp:208
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:524
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
cn1
double cn1
Definition: plastic.cpp:132
PlasticOps::CommonData::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
Definition: PlasticOps.hpp:62
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
C1_k
double C1_k
Kinematic hardening.
Definition: plastic.cpp:129
PlasticOps::CommonData::ParamsIndexes
ParamsIndexes
Definition: PlasticOps.hpp:45
kinematic_hardening_dplastic_strain
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:105
PlasticOps::opFactoryDomainReactions
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
PlasticOps::CommonData::mStrainPtr
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: PlasticOps.hpp:67
PlasticOps::CommonData::plasticTauDot
VectorDouble plasticTauDot
Definition: PlasticOps.hpp:74
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2181
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
cn0
double cn0
Definition: plastic.cpp:131
MatrixFunction.hpp
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:36
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
PlasticOps::CommonData::getPlasticTauPtr
auto getPlasticTauPtr()
Definition: PlasticOps.hpp:90
PlasticOps::CommonData::plasticStrain
MatrixDouble plasticStrain
Definition: PlasticOps.hpp:75
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1898
visH
double visH
Viscous hardening.
Definition: plastic.cpp:125
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
PlasticOps::CommonData::YOUNG_MODULUS
@ YOUNG_MODULUS
Definition: PlasticOps.hpp:46
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
PlasticOps::PlasticityIntegrators::Assembly::OpCalculateConstraintsLhs_dEP
OpCalculateConstraintsLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dEP
Definition: PlasticOps.hpp:224
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:40
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
ContactOps.hpp
FTensor::Index< 'i', DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
PlasticOps::CommonData::resFlowDstrainDot
MatrixDouble resFlowDstrainDot
Definition: PlasticOps.hpp:85
PlasticOps::CommonData::QINF
@ QINF
Definition: PlasticOps.hpp:51
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2223
b_iso
double b_iso
Saturation exponent.
Definition: plastic.cpp:128
PlasticOps::CommonData::resCdPlasticStrain
MatrixDouble resCdPlasticStrain
Definition: PlasticOps.hpp:81
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
MoFEM::PipelineManager::VolEle
MoFEM::VolumeElementForcesAndSourcesCore VolEle
Definition: PipelineManager.hpp:34
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()=default
PlasticOps
Definition: PlasticNaturalBCs.hpp:13
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
Example::testOperators
MoFEMErrorCode testOperators()
[Solve]
Definition: plastic.cpp:1214
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
Example::bC
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:512
PlasticOpsLargeStrains.hpp
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticFlowLhs_dTAU
OpCalculatePlasticFlowLhs_dTAUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dTAU
Definition: PlasticOps.hpp:212
IntegrationRules.hpp
iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:72
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
is_large_strains
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:116
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:233
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
PlasticOps::PlasticityIntegrators::OpCalculatePlasticSurface
OpCalculatePlasticSurfaceImpl< DIM, I, DomainEleOp > OpCalculatePlasticSurface
Definition: PlasticOps.hpp:168
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
FTensor::Ddg
Definition: Ddg_value.hpp:7
MoFEM::createSchurNestedMatrixStruture
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:1944
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::DMMoFEMTSSetMonitor
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:1056
IT
constexpr IntegrationType IT
Definition: test_broken_space.cpp:24
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ep_order
int ep_order
Order of ep files.
Definition: plastic.cpp:136
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
PlasticOps::PlasticityIntegrators::Assembly::OpCalculateConstraintsLhs_dTAU
OpCalculateConstraintsLhs_dTAUImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dTAU
Definition: PlasticOps.hpp:228
Qinf
double Qinf
Saturation yield stress.
Definition: plastic.cpp:127
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
PlasticOps::CommonData::resFlowDtau
MatrixDouble resFlowDtau
Definition: PlasticOps.hpp:83
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:214
PlasticOps::CommonData::resFlow
MatrixDouble resFlow
Definition: PlasticOps.hpp:82
PlasticOpsSmallStrains.hpp
MoFEM::FormsIntegrators
Integrator forms.
Definition: FormsIntegrators.hpp:306
PlasticOps::CommonData::plasticSurface
VectorDouble plasticSurface
[Common data set externally]
Definition: PlasticOps.hpp:71
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:517
PlasticOps::CommonData::plasticStrainDot
MatrixDouble plasticStrainDot
Definition: PlasticOps.hpp:76
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
PlasticOps::CommonData::VIS_H
@ VIS_H
Definition: PlasticOps.hpp:50
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:184
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:40
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
PlasticOps::CommonData::getPlasticTauDotPtr
auto getPlasticTauDotPtr()
Definition: PlasticOps.hpp:93
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1009
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
PlasticOps::I
FTensor::Index< 'I', 3 > I
[Common data]
Definition: PlasticOps.hpp:115
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:404
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
MoFEM::SmartPetscObj< VecScatter >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:768
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
PlasticOps::CommonData::resFlowDstrain
MatrixDouble resFlowDstrain
Definition: PlasticOps.hpp:84
PlasticOps::CommonData::POISSON_RATIO
@ POISSON_RATIO
Definition: PlasticOps.hpp:47
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
PlasticOps::PlasticityIntegrators::Assembly::OpCalculateConstraintsRhs
OpCalculateConstraintsRhsImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsRhs
Definition: PlasticOps.hpp:187
PlasticNaturalBCs.hpp
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HenckyOps::HenckyIntegrators
Definition: HenckyOps.hpp:569
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
H
double H
Hardening.
Definition: plastic.cpp:124
tau_order
int tau_order
Order of tau files.
Definition: plastic.cpp:135
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
alpha_damping
double alpha_damping
Definition: plastic.cpp:141
PlasticOps::CommonData::plasticFlow
MatrixDouble plasticFlow
Definition: PlasticOps.hpp:72
PlasticOpsMonitor.hpp
PlasticOps::PlasticityIntegrators::Assembly::OpCalculatePlasticFlowRhs
OpCalculatePlasticFlowRhsImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowRhs
Definition: PlasticOps.hpp:183
PlasticOps::CommonData::resCdTau
VectorDouble resCdTau
Definition: PlasticOps.hpp:79
PlasticOps::CommonData::SIGMA_Y
@ SIGMA_Y
Definition: PlasticOps.hpp:48
PlasticOps::PlasticityIntegrators::OpPlasticStress
OpPlasticStressImpl< DIM, I, DomainEleOp > OpPlasticStress
Definition: PlasticOps.hpp:174
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
PlasticOps::addMatBlockOps
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow