v0.13.1
thermo_plastic.cpp

Thermo plasticity

/**
* \file thermo_plastic.cpp
* \example thermo_plastic.cpp
*
* Thermo plasticity
*
*/
#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;
//! [Body force]
GAUSS>::OpSource<1, SPACE_DIM>;
//! [Body force]
//! [Only used with Hooke equation (linear material model)]
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
//! [Only used for dynamics]
//! [Only used with Hencky/nonlinear material]
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
//! [Only used with Hencky/nonlinear material]
//! [Essential boundary conditions]
//! [Essential boundary conditions]
// Thermal operators
/**
* @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
*
*/
/**
* @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
* and transpose of it, i.e. (T x FLAX)
*
*/
GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
/**
* @brief Integrate Lhs base of temerature times (heat capacity) times base of
* temperature (T x T)
*
*/
/**
* @brief Integrating Rhs flux base (1/k) flux (FLUX)
*/
GAUSS>::OpBaseTimesVector<3, 3, 1>;
/**
* @brief Integrate Rhs div flux base times temperature (T)
*
*/
GAUSS>::OpMixDivTimesU<3, 1, 2>;
/**
* @brief Integrate Rhs base of temerature time heat capacity times heat rate
* (T)
*
*/
/**
* @brief Integrate Rhs base of temerature times divergenc of flux (T)
*
*/
GAUSS>::OpSource<1, 1>;
double scale = 1.;
double young_modulus = 206913;
double poisson_ratio = 0.29;
double coeff_expansion = 10e-6;
double ref_temp = 0.0;
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) {}
private:
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;
};
std::vector<BcTempFun> bcTemperatureFunctions;
};
//! [Run problem]
}
//! [Run problem]
//! [Set up problem]
// Add field
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
// Mechanical fields
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);
// Temperature
const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
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, "", "-coeff_expansion",
&coeff_expansion, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
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_h,
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) << "Coeff_expansion " << coeff_expansion;
MOFEM_LOG("EXAMPLE", Sev::inform) << "Reference_temperature " << ref_temp;
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>(
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->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();
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->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,
} else {
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
}
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
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) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
if (SPACE_DIM == 2) {
pipeline.push_back(new 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()));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
"EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
"U", commonPlasticDataPtr->mGradPtr));
pipeline.push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
pipeline.push_back(new OpCalculateScalarFieldValues(
"T", commonPlasticDataPtr->getTempValPtr()));
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"T", commonPlasticDataPtr->getTempValDotPtr()));
"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>(
"U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
"U", commonPlasticDataPtr, m_D_ptr, 1));
if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
pipeline.push_back(
};
auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
"U", "T", commonPlasticDataPtr, m_D_ptr));
"U", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
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(
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(
pipeline.push_back(new OpCalculateContrainsLhs_dEP(
"TAU", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
"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(
pipeline.push_back(
};
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(
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(
boost::make_shared<Range>(temp_edges)));
}
}
};
// Mechanics
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());
// Axiator
CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
// Temperature
CHKERR add_domain_base_ops(feThermalLhs->getOpPtrVector());
CHKERR add_domain_ops_lhs_thermal(feThermalLhs->getOpPtrVector(),
// Mechanics
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());
// Axiator
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());
// Temperature
CHKERR add_domain_base_ops(feThermalRhs->getOpPtrVector());
CHKERR add_domain_ops_rhs_thermal(feThermalRhs->getOpPtrVector());
add_boundary_ops_rhs_thermal(pipeline_mng->getOpBoundaryRhsPipeline(),
*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) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
pipeline.push_back(
"U", commonPlasticDataPtr->mGradPtr));
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
pipeline.push_back(
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(new OpSetBc("U", false, reactionMarker));
// Calculate internal forece
pipeline.push_back(
pipeline.push_back(new OpUnSetBc("U"));
}
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = integration_rule_deviator;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
}
//! [Push operators to pipeline]
//! [Solve]
auto 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);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
"U", commonPlasticDataPtr->mGradPtr));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
postProcFe->getOpPtrVector().push_back(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"T", commonPlasticDataPtr->getTempValPtr()));
postProcFe->getOpPtrVector().push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
new OpPPMap(
postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
{{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()},
{"T", commonPlasticDataPtr->getTempValPtr()}},
{{"U", u_ptr}},
{},
{{"STRAIN", commonPlasticDataPtr->mStrainPtr},
{"STRESS", commonPlasticDataPtr->mStressPtr}}
)
);
};
auto scatter_create = [&](auto D, auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
ROW, "U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
Vec v;
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
VecScatter scatter;
CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(v),
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto 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();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
std::string param_file
constexpr double third
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr double shear_modulus_G
Definition: elastic.cpp:46
constexpr double bulk_modulus_K
Definition: elastic.cpp:45
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:747
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:800
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#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
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1003
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
double A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:70
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:56
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:77
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:54
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:79
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:75
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
BcTempFun(double v, FEMethod &fe)
double operator()(const double, const double, const double)
[Example]
Definition: plastic.cpp:111
boost::shared_ptr< DomainEle > feAxiatorRhs
Definition: plastic.cpp:130
std::vector< BcTempFun > bcTemperatureFunctions
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:129
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:138
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:128
boost::shared_ptr< DomainEle > feAxiatorLhs
Definition: plastic.cpp:131
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:776
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:134
boost::shared_ptr< DomainEle > feThermalRhs
boost::shared_ptr< DomainEle > feThermalLhs
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:176
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Definition: plastic.cpp:140
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:113
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:368
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:144
MoFEM::Interface & mField
Definition: plastic.cpp:118
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:135
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:137
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:156
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:300
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:126
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:133
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:127
Simple interface for fast problem set-up.
Definition: BcManager.hpp:18
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:120
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
transform Hcurl base fluxes from reference element to physical triangle
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
Definition: PostProc.hpp:166
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PostProc.hpp:120
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
VolumeElementForcesAndSourcesCore VolEle
double young_modulus
double heat_conductivity
int main(int argc, char *argv[])
double omega_inf
double amplitude_cycle
EntitiesFieldData::EntData EntData
double Qinf
static char help[]
[Solve]
double amplitude_shift
double phase_shift
double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
#define EXECUTABLE_DIMENSION
constexpr int SPACE_DIM
double visH
double poisson_ratio
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
double scale
long double hardening(long double tau, double temp)
double H
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
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< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
long double hardening_dtau(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
int order
double ref_temp
double b_iso
double coeff_expansion
double sigmaY
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpNormalMixVecTimesScalar< SPACE_DIM > OpTemperatureBC
long double hardening_dtemp(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
double cn
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temerature times divergenc of flux (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
constexpr EntityType boundary_ent
double omega_h