v0.13.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
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
//! [Body force]
using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, SPACE_DIM>;
//! [Body force]
//! [Only used with Hooke equation (linear material model)]
using OpKCauchy = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
using OpInternalForceCauchy = FormsIntegrators<DomainEleOp>::Assembly<
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
//! [Only used for dynamics]
//! [Only used with Hencky/nonlinear material]
using OpKPiola = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
using OpInternalForcePiola = FormsIntegrators<DomainEleOp>::Assembly<
//! [Only used with Hencky/nonlinear material]
//! [Essential boundary conditions]
using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryVec = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryInternal = FormsIntegrators<BoundaryEleOp>::Assembly<
//! [Essential boundary conditions]
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, double temp) {
return H * tau + Qinf * (1. - std::exp(-b_iso * tau)) + sigmaY;
}
inline long double hardening_dtau(long double tau, double temp) {
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) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode tsSolve();
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]
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR tsSolve();
}
//! [Run problem]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// 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();
}
//! [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>(
*commonPlasticDataPtr->mDPtr);
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->matDPtr = commonPlasticDataPtr->mDPtr;
commonHenckyDataPtr->matLogCPlastic =
commonPlasticDataPtr->getPlasticStrainPtr();
commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
commonPlasticDataPtr->mStressPtr =
commonHenckyDataPtr->getMatHenckyStress();
}
}
//! [Create common data]
//! [Boundary condition]
auto simple = mField.getInterface<Simple>();
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(simple->getProblemName(), "FIX_X", "U",
0, 0);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
1, 1);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
"U", 0, 3);
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,
*reactionMarker);
} else {
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
if (!reactionMarker) {
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
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) {
if (SPACE_DIM == 2) {
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 OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<SPACE_DIM>(
"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(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
} else {
pipeline.push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(
new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
}
if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
pipeline.push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
};
auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
pipeline.push_back(new OpCalculatePlasticInternalForceLhs_LogStrain_dEP(
"U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
"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));
const std::string block_name = "BODY_FORCE";
if (it->getName().compare(0, block_name.size(), block_name) == 0) {
std::vector<double> attr;
CHKERR it->getAttributes(attr);
if (attr.size() == 3) {
bodyForces.push_back(
FTensor::Tensor1<double, 3>{attr[0], attr[1], attr[2]});
} else {
SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"Should be three atributes in BODYFORCE blockset, but is %d",
attr.size());
}
}
}
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
t_source(i) = 0;
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
// hardcoded gravity load
for (auto &t_b : bodyForces)
t_source(i) += (scale * t_b(i)) * time;
return t_source;
};
pipeline.push_back(new OpBodyForce("U", get_body_force));
// Calculate internal forece
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
}
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(
new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateContrainsLhs_LogStrain_dU(
"TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
"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(
new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonPlasticDataPtr));
pipeline.push_back(new OpCalculateContrainsLhs_dEP(
"TAU", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
new OpCalculateContrainsLhs_dTAU("TAU", "TAU", commonPlasticDataPtr));
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(
new OpCalculatePlasticFlowRhs("EP", commonPlasticDataPtr));
pipeline.push_back(
new OpCalculateContrainsRhs("TAU", commonPlasticDataPtr));
};
auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "FIX_")){
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
pipeline.push_back(new OpBoundaryMass(
"U", "U", [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
auto get_time = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return fe_domain_rhs->ts_t;
};
auto get_time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return fe_domain_rhs->ts_t * scale;
};
auto get_minus_time = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return -fe_domain_rhs->ts_t;
};
auto time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return -fe_domain_rhs->ts_t;
};
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
if (it->getName().compare(0, 5, "FORCE") == 0) {
Range force_edges;
std::vector<double> attr_vec;
CHKERR it->getMeshsetIdEntitiesByDimension(
mField.get_moab(), SPACE_DIM - 1, force_edges, true);
it->getAttributes(attr_vec);
if (attr_vec.size() < SPACE_DIM)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong size of boundary attributes vector. Set right block "
"size attributes.");
auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
force_vec_ptr->data().begin());
pipeline.push_back(
new OpBoundaryVec("U", force_vec_ptr, get_time_scaled,
boost::make_shared<Range>(force_edges)));
}
}
pipeline.push_back(new OpUnSetBc("U"));
auto u_mat_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
if (bc_mng->checkBlock(bc, "FIX_")) {
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
attr_vec->clear();
if (bc.second->bcAttributes.size() < SPACE_DIM)
SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong size of boundary attributes vector. Set right block "
"size attributes. Size of attributes %d",
bc.second->bcAttributes.size());
std::copy(&bc.second->bcAttributes[0],
&bc.second->bcAttributes[SPACE_DIM],
attr_vec->data().begin());
pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpBoundaryInternal(
"U", u_mat_ptr, [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
};
// 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) {
if (reactionMarker) {
if (SPACE_DIM == 2) {
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 OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
pipeline.push_back(
"U", commonPlasticDataPtr->mGradPtr));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
pipeline.push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr));
pipeline.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(new OpPlasticStress("U", commonPlasticDataPtr,
commonPlasticDataPtr->mDPtr, 1));
}
pipeline.push_back(new OpSetBc("U", false, reactionMarker));
// Calculate internal forece
pipeline.push_back(new OpInternalForcePiola(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
}
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]
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
CHKERR SNESMonitorSet(
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
};
auto create_post_process_element = [&]() {
postProcFe = boost::make_shared<PostProcEle>(mField);
postProcFe->generateReferenceElementMesh();
if (SPACE_DIM == 2) {
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(
new OpCalculateHOJacForFace(jac_ptr));
postProcFe->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
postProcFe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
postProcFe->getOpPtrVector().push_back(
"U", commonPlasticDataPtr->mGradPtr));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
postProcFe->getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong pointer for grad");
postProcFe->getOpPtrVector().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr, scale));
postProcFe->getOpPtrVector().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(new OpPostProcHencky<SPACE_DIM>(
"U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
commonHenckyDataPtr));
} else {
postProcFe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
"U", commonPlasticDataPtr, commonPlasticDataPtr->mDPtr, scale));
postProcFe->getOpPtrVector().push_back(
"U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
commonPlasticDataPtr->mStrainPtr,
commonPlasticDataPtr->mStressPtr));
}
postProcFe->getOpPtrVector().push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpPostProcPlastic("U", postProcFe->postProcMesh,
postProcFe->mapGaussPts, commonPlasticDataPtr));
postProcFe->addFieldValuesPostProc("U");
};
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) {
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::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]
Simple *simple = m_field.getInterface<Simple>();
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
int main(int argc, char *argv[])
static char help[]
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
constexpr double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double young_modulus
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:758
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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:811
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
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
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:994
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
auto hardening_dtau(long double tau)
auto hardening(long double tau)
const double D
diffusivity
double A
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:88
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:86
double Qinf
Definition: plastic.cpp:112
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
double visH
Definition: plastic.cpp:110
double scale
Definition: plastic.cpp:103
double H
Definition: plastic.cpp:109
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
double b_iso
Definition: plastic.cpp:113
PetscBool is_large_strains
Definition: plastic.cpp:101
double sigmaY
Definition: plastic.cpp:108
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
double cn
Definition: plastic.cpp:111
constexpr EntityType boundary_ent
Definition: plastic.cpp:54
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
[Example]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:787
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Scale base functions by inverses of measure of element.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Postprocess on face.
Post processing.
[Class definition]
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
/** \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> {
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mDPtr_Axiator;
boost::shared_ptr<MatrixDouble> mDPtr_Deviator;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
inline auto getPlasticTauPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTau);
}
inline auto getPlasticTauDotPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTauDot);
}
inline auto getPlasticStrainPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticStrain);
}
inline auto getPlasticStrainDotPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto getTempValPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &tempVal);
}
};
//! [Common data]
FTensor::Index<'L', (SPACE_DIM * (SPACE_DIM + 1)) / 2> L;
//! [Operators definitions]
struct OpCalculatePlasticSurface : public DomainEleOp {
OpCalculatePlasticSurface(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpPlasticStress : public DomainEleOp {
OpPlasticStress(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> mDPtr,
const double scale = 1);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<MatrixDouble> mDPtr;
const double scaleStress;
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculatePlasticFlowRhs : public AssemblyDomainEleOp {
OpCalculatePlasticFlowRhs(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculateContrainsRhs : public AssemblyDomainEleOp {
OpCalculateContrainsRhs(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculatePlasticInternalForceLhs_dEP : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticInternalForceLhs_LogStrain_dEP
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> common_henky_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_dU : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_LogStrain_dU : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_dEP : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculatePlasticFlowLhs_dTAU : public AssemblyDomainEleOp {
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpCalculateContrainsLhs_dU : public AssemblyDomainEleOp {
OpCalculateContrainsLhs_dU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculateContrainsLhs_LogStrain_dU : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculateContrainsLhs_dEP : public AssemblyDomainEleOp {
OpCalculateContrainsLhs_dEP(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> mat_D_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
};
struct OpCalculateContrainsLhs_dTAU : public AssemblyDomainEleOp {
OpCalculateContrainsLhs_dTAU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpPostProcPlastic : public DomainEleOp {
OpPostProcPlastic(const std::string field_name,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
std::vector<EntityHandle> &mapGaussPts;
boost::shared_ptr<CommonData> commonDataPtr;
};
//! [Operators definitions]
//! [Lambda functions]
inline auto diff_tensor() {
t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
return t_diff;
};
inline auto symm_L_tensor() {
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
t_L(i, j, L) = 0;
if (SPACE_DIM == 2) {
t_L(0, 0, 0) = 1;
t_L(1, 0, 1) = 1;
t_L(1, 1, 2) = 1;
} else if (SPACE_DIM == 3) {
t_L(0, 0, 0) = 1;
t_L(1, 0, 1) = 1;
t_L(2, 0, 2) = 1;
t_L(1, 1, 3) = 1;
t_L(2, 1, 4) = 1;
t_L(2, 2, 5) = 1;
}
return t_L;
}
inline auto diff_symmetrize() {
t_diff(i, j, k, l) = 0;
t_diff(0, 0, 0, 0) = 1;
t_diff(1, 1, 1, 1) = 1;
t_diff(1, 0, 1, 0) = 0.5;
t_diff(1, 0, 0, 1) = 0.5;
t_diff(0, 1, 0, 1) = 0.5;
t_diff(0, 1, 1, 0) = 0.5;
if (SPACE_DIM == 3) {
t_diff(2, 2, 2, 2) = 1;
t_diff(2, 0, 2, 0) = 0.5;
t_diff(2, 0, 0, 2) = 0.5;
t_diff(0, 2, 0, 2) = 0.5;
t_diff(0, 2, 2, 0) = 0.5;
t_diff(2, 1, 2, 1) = 0.5;
t_diff(2, 1, 1, 2) = 0.5;
t_diff(1, 2, 1, 2) = 0.5;
t_diff(1, 2, 2, 1) = 0.5;
}
return t_diff;
};
template <typename T>
constexpr double third = boost::math::constants::third<double>();
if (SPACE_DIM == 2)
return (t_stress(0, 0) + t_stress(1, 1)) * third;
else
return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) * third;
};
template <typename T>
double trace) {
t_dev(I, J) = 0;
for (int ii = 0; ii != SPACE_DIM; ++ii)
for (int jj = ii; jj != SPACE_DIM; ++jj)
t_dev(ii, jj) = t_stress(ii, jj);
t_dev(0, 0) -= trace;
t_dev(1, 1) -= trace;
t_dev(2, 2) -= trace;
return t_dev;
};
inline auto
t_diff_deviator(I, J, k, l) = 0;
for (int ii = 0; ii != SPACE_DIM; ++ii)
for (int jj = ii; jj != SPACE_DIM; ++jj)
for (int kk = 0; kk != SPACE_DIM; ++kk)
for (int ll = kk; ll != SPACE_DIM; ++ll)
t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
constexpr double third = boost::math::constants::third<double>();
t_diff_deviator(0, 0, 0, 0) -= third;
t_diff_deviator(0, 0, 1, 1) -= third;
t_diff_deviator(1, 1, 0, 0) -= third;
t_diff_deviator(1, 1, 1, 1) -= third;
t_diff_deviator(2, 2, 0, 0) -= third;
t_diff_deviator(2, 2, 1, 1) -= third;
if (SPACE_DIM == 3) {
t_diff_deviator(0, 0, 2, 2) -= third;
t_diff_deviator(1, 1, 2, 2) -= third;
t_diff_deviator(2, 2, 2, 2) -= third;
}
return t_diff_deviator;
};
/**
*
\f[
\begin{split}
f&=\sqrt{s_{ij}s_{ij}}\\
A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
\frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
\frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
\sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
-A_{mn}A_{ij}
\right)\\
\frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
\\
\frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
\frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
\sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
\sigma_{mn}} D_{mnkl}
\end{split}
\f]
*/
inline double
return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
std::numeric_limits<double>::epsilon();
};
inline auto plastic_flow(long double f,
t_diff_f(k, l) =
(1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l))) / f;
return t_diff_f;
};
template <typename T>
t_diff_flow(i, j, k, l) =
(1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
(2. / 3.) * t_flow(i, j) * t_flow(k, l))) /
f;
return t_diff_flow;
};
template <typename T>
FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> &&t_diff_plastic_flow_dstress) {
t_diff_flow(i, j, k, l) =
t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
return t_diff_flow;
};
inline double constrain_abs(long double x) {
return std::abs(x);
};
inline double constrian_sign(long double x) {
if (x > 0)
return 1;
else if (x < 0)
return -1;
else
return 0;
};
inline double w(long double dot_tau, long double f, long double sigma_y) {
return (f - sigma_y) / sigmaY + cn * dot_tau;
};
/**
\f[
\dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
\| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0 \\
c_n \sigma_y \dot{\tau} - \frac{1}{2}\left\{c_n\sigma_y \dot{\tau} +
(f(\pmb\sigma) - \sigma_y) +
\| c_n \sigma_y \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
\f]
*/
inline double constrain(long double dot_tau, long double f,
long double sigma_y) {
return visH * dot_tau +
(sigmaY / 2) * ((cn * dot_tau - (f - sigma_y) / sigmaY) -
constrain_abs(w(dot_tau, f, sigma_y)));
};
inline double diff_constrain_ddot_tau(long double dot_tau, long double f,
long double sigma_y) {
return visH +
(sigmaY / 2) * (cn - cn * constrian_sign(w(dot_tau, f, sigma_y)));
};
inline auto diff_constrain_df(long double dot_tau, long double f,
long double sigma_y) {
return (-1 - constrian_sign(w(dot_tau, f, sigma_y))) / 2;
};
inline auto diff_constrain_dsigma_y(long double dot_tau, long double f,
long double sigma_y) {
return (1 + constrian_sign(w(dot_tau, f, sigma_y))) / 2;
}
template <typename T>
t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
return t_diff_constrain_dstress;
};
template <typename T1, typename T2>
inline auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress) {
t_diff_constrain_dstrain(k, l) =
t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
return t_diff_constrain_dstrain;
};
//! [Lambda functions]
//! [Auxiliary functions functions
&mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
&mat(3 * rr + 1, 1), &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1),
&mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
&mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
&mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
&mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
&mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
&mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
&mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
&mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
&mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
}
static inline auto get_mat_scalar_dvector(MatrixDouble &mat,
&mat(0, 1)};
}
static inline auto get_mat_scalar_dvector(MatrixDouble &mat,
&mat(0, 0), &mat(0, 1), &mat(0, 2)};
}
//! [Auxiliary functions functions
}; // namespace PlasticOps
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
auto symm_L_tensor()
Definition: PlasticOps.hpp:301
auto deviator(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress, double trace)
Definition: PlasticOps.hpp:360
auto diff_constrain_dstress(double &&diff_constrain_df, FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_flow)
Definition: PlasticOps.hpp:515
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
Definition: PlasticOps.hpp:427
FTensor::Index< 'j', SPACE_DIM > j
[Common data]
Definition: PlasticOps.hpp:91
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:100
auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
Definition: PlasticOps.hpp:524
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
Definition: PlasticOps.hpp:432
static auto get_mat_scalar_dvector(MatrixDouble &mat, FTensor::Number< 2 >)
Definition: PlasticOps.hpp:555
FTensor::Index< 'l', SPACE_DIM > l
Definition: PlasticOps.hpp:93
FTensor::Index< 'k', SPACE_DIM > k
Definition: PlasticOps.hpp:92
auto diff_tensor()
[Operators definitions]
Definition: PlasticOps.hpp:294
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:101
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:98
double constrain(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:491
auto diff_deviator(FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_stress)
Definition: PlasticOps.hpp:374
FTensor::Index< 'i', SPACE_DIM > i
Definition: PlasticOps.hpp:95
double w(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:476
FTensor::Index< 'm', SPACE_DIM > m
Definition: PlasticOps.hpp:94
auto diff_constrain_dsigma_y(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:509
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:99
double trace(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress)
Definition: PlasticOps.hpp:351
auto diff_symmetrize()
Definition: PlasticOps.hpp:320
double constrian_sign(long double x)
Definition: PlasticOps.hpp:467
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_flow, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
Definition: PlasticOps.hpp:442
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, SPACE_DIM, SPACE_DIM > &t_D, FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_plastic_flow_dstress)
Definition: PlasticOps.hpp:454
double diff_constrain_ddot_tau(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:498
double constrain_abs(long double x)
Definition: PlasticOps.hpp:463
auto diff_constrain_df(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:504
static FTensor::Tensor3< FTensor::PackPtr< double *, 2 >, 2, 2, 2 > get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
Definition: PlasticOps.hpp:534
FTensor::Index< 'L',(SPACE_DIM *(SPACE_DIM+1))/2 > L
Definition: PlasticOps.hpp:103
FTensor::Index< 'n', SPACE_DIM > n
Definition: PlasticOps.hpp:96
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: PlasticOps.hpp:60
VectorDouble tempVal
Definition: PlasticOps.hpp:69
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: PlasticOps.hpp:59
VectorDouble plasticTauDot
Definition: PlasticOps.hpp:66
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:56
VectorDouble plasticSurface
Definition: PlasticOps.hpp:63
boost::shared_ptr< MatrixDouble > mDPtr_Axiator
Definition: PlasticOps.hpp:57
MatrixDouble plasticStrain
Definition: PlasticOps.hpp:67
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: PlasticOps.hpp:61
MatrixDouble plasticFlow
Definition: PlasticOps.hpp:64
MatrixDouble plasticStrainDot
Definition: PlasticOps.hpp:68
VectorDouble plasticTau
Definition: PlasticOps.hpp:65
boost::shared_ptr< MatrixDouble > mDPtr_Deviator
Definition: PlasticOps.hpp:58
OpCalculateContrainsLhs_LogStrain_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > comman_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:251
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:252
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:253
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:266
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculateContrainsLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:265
OpCalculateContrainsLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:276
OpCalculateContrainsLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:238
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:237
OpCalculateContrainsRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:143
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticFlowLhs_LogStrain_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > comman_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:201
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:200
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:199
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:214
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:213
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:225
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticFlowLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:186
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:185
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:134
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOps.hpp:171
OpCalculatePlasticInternalForceLhs_LogStrain_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< HenckyOps::CommonData > common_henky_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpCalculatePlasticInternalForceLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:155
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:156
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:112
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:125
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mDPtr, const double scale=1)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:123
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:289
std::vector< EntityHandle > & mapGaussPts
Definition: PlasticOps.hpp:288
OpPostProcPlastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
moab::Interface & postProcMesh
Definition: PlasticOps.hpp:287
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]