v0.13.1
plastic.cpp

Plasticity in 2d and 3d

/**
* \file plastic.cpp
* \example plastic.cpp
*
* Plasticity in 2d and 3d
*
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
//! [Only used with Hooke equation (linear material model)]
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
//! [Only used for dynamics]
//! [Only used with Hencky/nonlinear material]
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
//! [Only used with Hencky/nonlinear material]
//! [Essential boundary conditions]
//! [Essential boundary conditions]
using OpBodyForce =
DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
using OpForce =
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
PetscBool is_large_strains = PETSC_TRUE;
double scale = 1.;
double young_modulus = 206913;
double poisson_ratio = 0.29;
double rho = 0;
double sigmaY = 450;
double H = 129;
double visH = 0;
double cn = 1;
double Qinf = 265;
double b_iso = 16.93;
int order = 2;
inline long double hardening(long double tau) {
return H * tau + Qinf * (1. - std::exp(-b_iso * tau)) + sigmaY;
}
inline long double hardening_dtau(long double tau) {
return H + Qinf * b_iso * std::exp(-b_iso * tau);
}
#include <HenckyOps.hpp>
#include <PlasticOps.hpp>
using namespace PlasticOps;
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<PostProcEle> postProcFe;
boost::shared_ptr<DomainEle> reactionFe;
boost::shared_ptr<DomainEle> feAxiatorRhs;
boost::shared_ptr<DomainEle> feAxiatorLhs;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
std::vector<FTensor::Tensor1<double, 3>> bodyForces;
};
//! [Run problem]
}
//! [Run problem]
//! [Set up problem]
// Add field
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
CHKERR simple->addDomainField("TAU", L2, base, 1);
CHKERR simple->addDomainField("EP", L2, base, size_symm);
CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("TAU", order - 1);
CHKERR simple->setFieldOrder("EP", order - 1);
CHKERR simple->setUp();
CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
}
//! [Set up problem]
//! [Create common data]
auto get_command_line_parameters = [&]() {
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
&young_modulus, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
&poisson_ratio, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
&is_large_strains, PETSC_NULL);
MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
MOFEM_LOG("EXAMPLE", Sev::inform) << "cn " << cn;
PetscBool is_scale = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
PETSC_NULL);
if (is_scale) {
rho *= scale;
H *= scale;
Qinf *= scale;
visH *= scale;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Young modulus " << young_modulus;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Poisson ratio " << poisson_ratio;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled Saturation yield stress " << Qinf;
}
};
auto set_matrial_stiffness = [&]() {
const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
// Plane stress or when 1, plane strain or 3d
const double A = (SPACE_DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Axiator);
auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Deviator);
constexpr double third = boost::math::constants::third<double>();
t_D_axiator(i, j, k, l) = A *
(bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
t_kd(i, j) * t_kd(k, l);
t_D_deviator(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
};
commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
CHKERR get_command_line_parameters();
CHKERR set_matrial_stiffness();
commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
commonHenckyDataPtr->matLogCPlastic =
commonPlasticDataPtr->getPlasticStrainPtr();
commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
commonPlasticDataPtr->mStressPtr =
commonHenckyDataPtr->getMatHenckyStress();
}
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
auto &bc_map = bc_mng->getBcMapByBlockName();
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
"U", 0, 3);
for (auto bc : bc_map)
MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
// OK. We have problem with GMesh, it adding empty characters at the end of
// block. So first block is search by regexp. popMarkDOFsOnEntities should
// work with regexp.
std::string reaction_block_set;
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "REACTION")) {
reaction_block_set = bc.first;
break;
}
}
if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
reactionMarker = bc->getBcMarkersPtr();
// Only take reaction from nodes
Range nodes;
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
ProblemsManager::MarkOP::AND, nodes,
} else {
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
feAxiatorLhs = boost::make_shared<DomainEle>(mField);
feAxiatorRhs = boost::make_shared<DomainEle>(mField);
auto integration_rule_axiator = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
};
feAxiatorLhs->getRuleHook = integration_rule_axiator;
feAxiatorRhs->getRuleHook = integration_rule_axiator;
auto integration_rule_deviator = [](int o_row, int o_col, int approx_order) {
return 2 * (approx_order - 1);
};
auto integration_rule_bc = [](int, int, int approx_order) {
return 2 * approx_order;
};
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
if (SPACE_DIM == 2)
pipeline.push_back(new OpSetHOWeightsOnFace());
else
pipeline.push_back(new OpSetHOWeights(det_ptr));
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
"EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
"U", commonPlasticDataPtr->mGradPtr));
pipeline.push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
};
auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
"U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
} else {
pipeline.push_back(
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(
new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
}
if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
pipeline.push_back(
};
auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(
new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
"U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
"U", "EP", commonPlasticDataPtr, m_D_ptr));
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
pipeline, mField, "U", {boost::make_shared<TimeScale>()}, "BODY_FORCE",
Sev::inform);
// Calculate internal forces
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(new OpCalculateContrainsLhs_LogStrain_dU(
"TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
"EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
"EP", "U", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculateContrainsLhs_dU(
"TAU", "U", commonPlasticDataPtr, m_D_ptr));
}
pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
"EP", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
pipeline.push_back(new OpCalculateContrainsLhs_dEP(
"TAU", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
mField, pipeline, simple->getProblemName(), "U");
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
pipeline, mField, "U", {boost::make_shared<TimeScale>()}, "FORCE",
Sev::inform);
pipeline.push_back(new OpUnSetBc("U"));
auto u_mat_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
mField, pipeline, simple->getProblemName(), "U", u_mat_ptr,
{boost::make_shared<TimeScale>()});
};
// Axiator
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_boundary_ops_lhs_mechanical(
pipeline_mng->getOpBoundaryLhsPipeline());
CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
CHKERR add_domain_stress_ops(feAxiatorLhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_ops_rhs_mechanical(
pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto create_reaction_pipeline = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
if (SPACE_DIM == 2)
pipeline.push_back(new OpSetHOWeightsOnFace());
else
pipeline.push_back(new OpSetHOWeights(det_ptr));
pipeline.push_back(
"U", commonPlasticDataPtr->mGradPtr));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(new OpPlasticStress(
}
pipeline.push_back(new OpSetBc("U", false, reactionMarker));
// Calculate internal forece
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
}
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = integration_rule_deviator;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
}
//! [Push operators to pipeline]
//! [Solve]
auto snes_ctx_ptr = smartGetDMSnesCtx(simple->getDM());
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_element = [&]() {
postProcFe = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
"U", commonPlasticDataPtr->mGradPtr));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
postProcFe->getOpPtrVector().push_back(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
new OpPPMap(
postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
{},
{{"U", u_ptr}},
{{"GRAD", commonPlasticDataPtr->mGradPtr},
{"FIRST_PIOLA", commonHenckyDataPtr->getMatFirstPiolaStress()}},
{}
)
);
} else {
postProcFe->getOpPtrVector().push_back(
commonPlasticDataPtr->mStrainPtr));
postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
postProcFe->getOpPtrVector().push_back(
new OpPPMap(
postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
{},
{{"U", u_ptr}},
{},
{{"STRAIN", commonPlasticDataPtr->mStrainPtr},
{"STRESS", commonPlasticDataPtr->mStressPtr}}
)
);
}
postProcFe->getOpPtrVector().push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpPPMap(
postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
{{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()}},
{},
{},
{{"PLASTIC_STRAIN", commonPlasticDataPtr->getPlasticStrainPtr()}}
)
);
};
auto scatter_create = [&](auto D, auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
ROW, "U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
Vec v;
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
VecScatter scatter;
CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(v),
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Setup fieldsplit (block) solver - optional: yes/no
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
}
};
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
boost::shared_ptr<FEMethod> null;
CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feAxiatorLhs,
null, null);
CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feAxiatorRhs,
null, null);
CHKERR create_post_process_element();
uXScatter = scatter_create(D, 0);
uYScatter = scatter_create(D, 1);
if (SPACE_DIM == 3)
uZScatter = scatter_create(D, 2);
auto solver = pipeline_mng->createTSIM();
CHKERR TSSetSolution(solver, D);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
std::string param_file
constexpr double third
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr double shear_modulus_G
Definition: elastic.cpp:63
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr auto t_kd
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:747
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:800
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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: DMMMoFEM.cpp:1003
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:224
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:987
const double D
diffusivity
double A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpEssentialRhs< DisplacementCubitBcData, 1, SPACE_DIM > OpEssentialRhs
Definition: plastic.cpp:90
double young_modulus
Definition: plastic.cpp:96
int main(int argc, char *argv[])
Definition: plastic.cpp:915
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:63
EntitiesFieldData::EntData EntData
Definition: plastic.cpp:35
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
Definition: plastic.cpp:82
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:61
double Qinf
Definition: plastic.cpp:103
static char help[]
[Solve]
Definition: plastic.cpp:913
double rho
Definition: plastic.cpp:98
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:49
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
constexpr int SPACE_DIM
Definition: plastic.cpp:31
double visH
Definition: plastic.cpp:101
double poisson_ratio
Definition: plastic.cpp:97
BoundaryNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpForce
Definition: plastic.cpp:84
double scale
Definition: plastic.cpp:94
NaturalBC< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS > DomainNaturalBC
Definition: plastic.cpp:77
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpEssentialLhs< DisplacementCubitBcData, 1, SPACE_DIM > OpEssentialLhs
Definition: plastic.cpp:88
double H
Definition: plastic.cpp:100
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:70
long double hardening(long double tau)
Definition: plastic.cpp:107
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:47
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:72
int order
Definition: plastic.cpp:105
double b_iso
Definition: plastic.cpp:104
PetscBool is_large_strains
Definition: plastic.cpp:92
double sigmaY
Definition: plastic.cpp:99
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:68
long double hardening_dtau(long double tau)
Definition: plastic.cpp:111
double cn
Definition: plastic.cpp:102
constexpr EntityType boundary_ent
Definition: plastic.cpp:34
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
[Example]
Definition: plastic.cpp:120
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:137
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:683
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:143
boost::shared_ptr< DomainEle > feAxiatorLhs
Definition: plastic.cpp:140
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:186
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:122
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:371
boost::shared_ptr< DomainEle > feAxiatorRhs
Definition: plastic.cpp:139
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Definition: plastic.cpp:149
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:153
MoFEM::Interface & mField
Definition: plastic.cpp:127
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:144
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:146
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:165
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:309
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:147
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:138
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:135
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:142
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:136
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Natural boundary conditions.
Definition: Essential.hpp:101
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Assembly methods.
Definition: Natural.hpp:57
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
VolumeElementForcesAndSourcesCore VolEle