v0.13.0
Note
Prerequisites of this tutorial include VEC-0: Linear elasticity

Note
Intended learning outcome:
• multi-physics
• multi-field problem
/**
* \file thermo_plastic.cpp
* \example thermo_plastic.cpp
*
* Thermo plasticity
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* 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
*
* You should have received a copy of the GNU Lesser General Public
#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<
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<
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]
// Thermal operators
/**
* @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
*
*/
using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
/**
* @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
* and transpose of it, i.e. (T x FLAX)
*
*/
using OpHdivT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
/**
* @brief Integrate Lhs base of temerature times (heat capacity) times base of
* temperature (T x T)
*
*/
using OpCapacity = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
/**
* @brief Integrating Rhs flux base (1/k) flux (FLUX)
*/
using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpBaseTimesVector<3, 3, 1>;
/**
* @brief Integrate Rhs div flux base times temperature (T)
*
*/
using OpHDivTemp = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpMixDivTimesU<3, 1, 2>;
/**
* @brief Integrate Rhs base of temerature time heat capacity times heat rate
* (T)
*
*/
using OpBaseDotT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
/**
* @brief Integrate Rhs base of temerature times divergenc of flux (T)
*
*/
using OpHeatSource = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, 1>;
using OpTemperatureBC = FormsIntegrators<BoundaryEleOp>::Assembly<
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 = 0; // 265;
double b_iso = 16.93;
16.2; // Force / (time temerature ) or Power /
// (length temperature) // Time unit is hour. force unit kN
double heat_capacity =
5961.6; // length^2/(time^2 temerature) // length is millimeter time is hour
double omega_0 = 2e-3;
double omega_h = 2e-3;
double omega_inf = 0;
int order = 2;
double amplitude_cycle = 0.5;
double amplitude_shift = 0.5;
double phase_shift = 0.8;
inline long double hardening(long double tau, double temp) {
return H * tau * (1 - omega_h * temp) +
Qinf * (1. - std::exp(-b_iso * tau)) * (1 - omega_inf * temp) +
sigmaY * (1 - omega_0 * temp);
}
inline long double hardening_dtau(long double tau, double temp) {
return H * (1 - omega_h * temp) +
Qinf * b_iso * std::exp(-b_iso * tau) * (1 - omega_inf * temp);
}
inline long double hardening_dtemp(long double tau, double temp) {
return -H * tau * omega_h - Qinf * (1. - std::exp(-b_iso * tau)) * omega_inf -
}
#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<PlasticThermalOps::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;
boost::shared_ptr<DomainEle> feThermalRhs;
boost::shared_ptr<DomainEle> feThermalLhs;
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;
struct BcTempFun {
BcTempFun(double v, FEMethod &fe) : valTemp(v), fE(fe) {}
double operator()(const double, const double, const double) {
return -valTemp; // * fE.ts_t;
}
private:
double valTemp;
FEMethod &fE;
};
std::vector<BcTempFun> bcTemperatureFunctions;
};
//! [Run problem]
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR tsSolve();
}
//! [Run problem]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
// Mechanical fields
// Temerature
const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
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->setFieldOrder("FLUX", order);
CHKERR simple->setFieldOrder("T", 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 PetscOptionsGetScalar(PETSC_NULL, "", "-capacity", &heat_capacity,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
&heat_conductivity, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_0", &omega_0,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_h", &omega_0,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_inf", &omega_inf,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-fraction_of_dissipation",
&fraction_of_dissipation, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-number_of_cycles_in_total_time",
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_cycle",
&amplitude_cycle, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phase_shift", &phase_shift,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_shift",
&amplitude_shift, 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;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "fraction_of_dissipation " << fraction_of_dissipation;
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;
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Scaled heat capacity " << heat_capacity;
}
};
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<PlasticThermalOps::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->mStrainPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
CHKERR get_command_line_parameters();
CHKERR set_matrial_stiffness();
}
//! [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->removeBlockDOFsOnEntities(simple->getProblemName(),
"ZERO_FLUX", "FLUX", 0, 1);
PetscBool zero_fix_skin_flux = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-fix_skin_flux",
&zero_fix_skin_flux, PETSC_NULL);
if (zero_fix_skin_flux) {
Range faces;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, faces);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, faces, false, skin_ents);
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
if (pcomm == NULL)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Communicator not created");
CHKERR pcomm->filter_pstatus(skin_ents,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &skin_ents);
CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "FLUX",
skin_ents, 0, 1);
}
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 b : bc_map) {
MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << b.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>();
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;
};
feAxiatorRhs = boost::make_shared<DomainEle>(mField);
feAxiatorLhs = boost::make_shared<DomainEle>(mField);
auto integration_rule_axiator = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
};
feAxiatorRhs->getRuleHook = integration_rule_axiator;
feAxiatorLhs->getRuleHook = integration_rule_axiator;
feThermalRhs = boost::make_shared<DomainEle>(mField);
feThermalLhs = boost::make_shared<DomainEle>(mField);
auto integration_rule_thermal = [](int, int, int approx_order) {
return 2 * approx_order;
};
feThermalRhs->getRuleHook = integration_rule_thermal;
feThermalLhs->getRuleHook = integration_rule_thermal;
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 OpMakeHdivFromHcurl());
pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
pipeline.push_back(new OpSetInvJacL2ForFace(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()));
pipeline.push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
pipeline.push_back(new OpCalculateScalarFieldValues(
"T", commonPlasticDataPtr->getTempValPtr()));
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"T", commonPlasticDataPtr->getTempValDotPtr()));
pipeline.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
"FLUX", commonPlasticDataPtr->getTempDivFluxPtr()));
pipeline.push_back(new OpCalculateHVecVectorField<3>(
"FLUX", commonPlasticDataPtr->getTempFluxValPtr()));
};
auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
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 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;
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 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 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));
"TAU", "T", commonPlasticDataPtr));
"T", "EP", 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));
"T", 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 amplitude_cycle * sin((2 * fe_domain_rhs->ts_t - phase_shift) *
} else {
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"));
}
}
};
auto add_domain_ops_lhs_thermal = [&](auto &pipeline, auto &fe) {
auto resistance = [](const double, const double, const double) {
return (1. / heat_conductivity);
};
auto capacity = [&](const double, const double, const double) {
return -heat_capacity * fe.ts_a;
};
auto unity = []() { return 1; };
pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
pipeline.push_back(new OpHdivT("FLUX", "T", unity, true));
pipeline.push_back(new OpCapacity("T", "T", capacity));
};
auto add_domain_ops_rhs_thermal = [&](auto &pipeline) {
auto resistance = [](const double, const double, const double) {
return (1. / heat_conductivity);
};
auto capacity = [&](const double, const double, const double) {
return -heat_capacity;
};
auto unity = [](const double, const double, const double) { return 1; };
pipeline.push_back(new OpHdivFlux(
"FLUX", commonPlasticDataPtr->getTempFluxValPtr(), resistance));
pipeline.push_back(
new OpHDivTemp("FLUX", commonPlasticDataPtr->getTempValPtr(), unity));
pipeline.push_back(new OpBaseDivFlux(
"T", commonPlasticDataPtr->getTempDivFluxPtr(), unity));
// auto source = [&](const double x, const double y, const double z) {
// return 1;
// };
// pipeline.push_back(new OpHeatSource("T", source));
pipeline.push_back(new OpBaseDotT(
"T", commonPlasticDataPtr->getTempValDotPtr(), capacity));
};
auto add_boundary_ops_rhs_thermal = [&](auto &pipeline, auto &fe) {
if (SPACE_DIM == 2) {
pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
} else if (SPACE_DIM == 3) {
pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
}
const std::string block_name = "TEMPERATURE";
if (it->getName().compare(0, block_name.size(), block_name) == 0) {
Range temp_edges;
std::vector<double> attr_vec;
CHKERR it->getMeshsetIdEntitiesByDimension(
mField.get_moab(), SPACE_DIM - 1, temp_edges, true);
it->getAttributes(attr_vec);
if (attr_vec.size() != 1)
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"Should be one attribute");
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Set temerature " << attr_vec[0] << " on ents:\n"
<< temp_edges;
bcTemperatureFunctions.push_back(BcTempFun(attr_vec[0], fe));
pipeline.push_back(
new OpTemperatureBC("FLUX", bcTemperatureFunctions.back(),
boost::make_shared<Range>(temp_edges)));
}
}
};
// Mechanics
commonPlasticDataPtr->mDPtr_Deviator);
commonPlasticDataPtr->mDPtr_Deviator);
commonPlasticDataPtr->mDPtr_Deviator);
pipeline_mng->getOpBoundaryLhsPipeline());
// Axiator
commonPlasticDataPtr->mDPtr_Axiator);
// Temperature
*feThermalLhs);
// Mechanics
commonPlasticDataPtr->mDPtr_Deviator);
pipeline_mng->getOpBoundaryRhsPipeline());
// Axiator
commonPlasticDataPtr->mDPtr_Axiator);
// Temperature
*pipeline_mng->getBoundaryRhsFE());
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(
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
pipeline.push_back(
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 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(
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
postProcFe->getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"T", commonPlasticDataPtr->getTempValPtr()));
postProcFe->getOpPtrVector().push_back(new OpSymmetrizeTensor<SPACE_DIM>(
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));
};
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 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 DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feThermalLhs,
null, null);
CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feThermalRhs,
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 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();
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]
Simple *simple = m_field.getInterface<Simple>();
CHKERR simple->getOptions();
//! [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
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#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
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1, 1 > OpBaseTimesScalarField
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
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
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
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
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()
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]
double heat_conductivity
double omega_inf
double amplitude_cycle
double amplitude_shift
double phase_shift
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseDotT
Integrate Rhs base of temerature time heat capacity times heat rate (T)
int number_of_cycles_in_total_time
double fraction_of_dissipation
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpHeatSource
double omega_0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temerature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
double heat_capacity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpNormalMixVecTimesScalar< SPACE_DIM > OpTemperatureBC
long double hardening_dtemp(long double tau, double temp)
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temerature times divergenc of flux (T)
double omega_h
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* 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
*
* You should have received a copy of the GNU Lesser General Public
/** \file PlasticThermalOps.hpp
* \example PlasticThermalOps.hpp
*/
namespace PlasticThermalOps {
//! [Common data]
inline auto getTempFluxValPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &tempFluxVal);
}
inline auto getTempValDotPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &tempValDot);
}
inline auto getTempDivFluxPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &templDivFlux);
}
};
//! [Common data]
struct OpPlasticHeatProduction : public AssemblyDomainEleOp {
OpPlasticHeatProduction(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
protected:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpPlasticHeatProduction_dEP : public AssemblyDomainEleOp {
OpPlasticHeatProduction_dEP(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_dT : public AssemblyDomainEleOp {
OpCalculateContrainsLhs_dT(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;
};
const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyDomainEleOp(field_name, field_name, DomainEleOp::OPROW),
commonDataPtr(common_data_ptr) {}
OpPlasticHeatProduction::iNtegrate(EntitiesFieldData::EntData &data) {
const size_t nb_integration_pts = data.getN().size1();
const size_t nb_base_functions = data.getN().size2();
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
auto t_w = getFTensor0IntegrationWeight();
auto t_base = data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
// dot(Eps^p): D : dot(Eps^p) * fraction_of_dissipation
const double beta = (alpha * fraction_of_dissipation) *
(t_D(i, j, k, l) * t_plastic_strain_dot(k, l)) *
t_plastic_strain_dot(i, j);
size_t bb = 0;
for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
nf[bb] += beta * t_base;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_w;
++t_plastic_strain_dot;
}
}
OpPlasticHeatProduction_dEP::OpPlasticHeatProduction_dEP(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode OpPlasticHeatProduction_dEP::iNtegrate(
const size_t nb_integration_pts = row_data.getN().size1();
const size_t nb_row_base_functions = row_data.getN().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
const double ts_a = getTSa();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w * fraction_of_dissipation * ts_a * 2;
auto mat_ptr = locMat.data().begin();
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
auto t_L = symm_L_tensor();
auto t_mat =
t_D_ep_dot(L) =
((t_D(i, j, k, l) * t_plastic_strain_dot(k, l)) * t_L(i, j, L)) *
(alpha);
size_t rr = 0;
for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
t_mat(L) += (t_row_base * t_col_base) * (t_D_ep_dot(L));
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_base;
++t_w;
++t_plastic_strain_dot;
}
}
OpCalculateContrainsLhs_dT::OpCalculateContrainsLhs_dT(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode OpCalculateContrainsLhs_dT::iNtegrate(
const size_t nb_integration_pts = row_data.getN().size1();
const size_t nb_row_base_functions = row_data.getN().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
if (commonDataPtr->tempVal.size() != nb_integration_pts) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Temperature not set");
}
auto t_temp = getFTensor0FromVec(commonDataPtr->tempVal);
auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
const double t_a = getTSa();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
const double c1 =
* diff_constrain_dsigma_y(t_tau_dot, t_f, hardening(t_tau, t_temp)) *
hardening_dtemp(t_tau, t_temp);
auto mat_ptr = locMat.data().begin();
size_t rr = 0;
for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
*mat_ptr += c1 * t_row_base * t_col_base;
++mat_ptr;
++t_col_base;
}
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr)
++t_row_base;
++t_w;
++t_f;
++t_tau;
++t_temp;
++t_tau_dot;
}
}
} // namespace PlasticThermalOps
static Index< 'L', 3 > L
constexpr double alpha
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
auto diff_constrain_dsigma_y(long double tau_dot, long double f, long double hardening)
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
boost::shared_ptr< CommonData > commonDataPtr
OpCalculateContrainsLhs_dT(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)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpPlasticHeatProduction_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
OpPlasticHeatProduction(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)