#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
};
};
(SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
: AssemblyType::PETSC;
IntegrationType::GAUSS;
#ifdef ADD_CONTACT
template <>
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
};
template <>
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
};
using BoundaryEleOp::BoundaryEleOp;
};
template <>
return MatSetValues<AssemblyTypeSelector<AT>>(
op_ptr->
getKSPB(), row_data, col_data,
m, ADD_VALUES);
};
#endif
return std::exp(
std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
}
}
};
constexpr double eps = 1e-12;
return std::max(
r(tau),
eps *
r(0));
}
template <typename T, int DIM>
inline auto
if (
C1_k < std::numeric_limits<double>::epsilon()) {
return t_alpha;
}
t_alpha(
i,
j) =
C1_k * t_plastic_strain(
i,
j);
return t_alpha;
}
template <int DIM>
return t_diff;
}
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
#endif
};
#include <ContactOps.hpp>
#endif
private:
};
};
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif
};
}
true);
auto get_ents_by_dim = [&](
const auto dim) {
return domain_ents;
} else {
else
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(
SPACE_DIM);
if (domain_ents.empty())
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
}
};
const auto base = get_base();
PetscBool order_edge = PETSC_FALSE;
PETSC_NULL);
PetscBool order_face = PETSC_FALSE;
PETSC_NULL);
PetscBool order_volume = PETSC_FALSE;
PETSC_NULL);
if (order_edge || order_face || order_volume) {
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order edge " << order_edge
? "true"
: "false";
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order face " << order_face
? "true"
: "false";
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order volume " << order_volume
? "true"
: "false";
auto ents = get_ents_by_dim(0);
if (order_edge)
ents.merge(get_ents_by_dim(1));
if (order_face)
ents.merge(get_ents_by_dim(2));
if (order_volume)
ents.merge(get_ents_by_dim(3));
} else {
}
#ifdef ADD_CONTACT
auto get_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_blocks = [&](auto skin) {
bool is_contact_block = false;
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true;
<<
"Find contact block set: " <<
m->getName();
auto meshset =
m->getMeshset();
Range contact_meshset_range;
meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
contact_meshset_range);
contact_range.merge(contact_meshset_range);
}
if (is_contact_block) {
<< "Nb entities in contact surface: " << contact_range.size();
skin = intersect(skin, contact_range);
}
return skin;
};
auto filter_true_skin = [&](auto skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
#endif
auto project_ho_geometry = [&]() {
};
}
auto get_command_line_parameters = [&]() {
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PetscBool tau_order_is_set;
&tau_order_is_set);
PetscBool ep_order_is_set;
&ep_order_is_set);
PETSC_NULL);
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Hardening " <<
H;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Viscous hardening " <<
visH;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation yield stress " <<
Qinf;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation exponent " <<
b_iso;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Kinematic hardening " <<
C1_k;
if (tau_order_is_set == PETSC_FALSE)
if (ep_order_is_set == PETSC_FALSE)
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Approximation order " <<
order;
<<
"Ep approximation order " <<
ep_order;
PetscBool is_scale = PETSC_TRUE;
PETSC_NULL);
if (is_scale) {
}
#ifdef ADD_CONTACT
#endif
};
CHKERR get_command_line_parameters();
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdfPythonPtr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
#ifdef ADD_CONTACT
for (auto b : {"FIX_X", "REMOVE_X"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
b,
"SIGMA", 0, 0, false, true);
for (auto b : {"FIX_Y", "REMOVE_Y"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
b,
"SIGMA", 1, 1, false, true);
for (auto b : {"FIX_Z", "REMOVE_Z"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
b,
"SIGMA", 2, 2, false, true);
for (auto b : {"FIX_ALL", "REMOVE_ALL"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
b,
"SIGMA", 0, 3, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT",
"SIGMA", 0, 3,
false,
true);
#endif
simple->getProblemName(),
"U");
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map)
MOFEM_LOG(
"PLASTICITY", Sev::verbose) <<
"Marker " << bc.first;
}
auto integration_rule_bc = [](
int,
int,
int ao) {
return 2 * ao; };
auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
"GEOMETRY");
pip,
mField,
"U", Sev::inform);
#ifdef ADD_CONTACT
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U");
ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
mField, pip,
simple->getDomainFEName(),
"SIGMA",
"U",
"GEOMETRY",
vol_rule);
#endif
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
"GEOMETRY");
pip,
mField,
"U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
pip, "SIGMA", "U");
#endif
};
auto add_domain_ops_lhs = [this](auto &pip) {
"GEOMETRY");
auto get_inertia_and_mass_damping = [
this](
const double,
const double,
return (
rho /
scale) * fe_domain_lhs->ts_aa +
};
pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
}
CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
};
auto add_domain_ops_rhs = [this](auto &pip) {
"GEOMETRY");
{boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
Sev::inform);
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(
new OpInertiaForce(
"U", mat_acceleration, [](
double,
double,
double) {
}));
auto mat_velocity = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
}));
}
}
CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
pip, "SIGMA", "U");
#endif
};
CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto create_reaction_pipeline = [&](auto &pip) {
"GEOMETRY");
CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
};
}
);
protected:
};
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(
mField);
auto &pip = pp_fe->getOpPtrVector();
auto push_vol_ops = [this](auto &pip) {
"GEOMETRY");
auto [common_plastic_ptr, common_henky_ptr] =
PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
if (common_henky_ptr) {
if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
}
return std::make_pair(common_plastic_ptr, common_henky_ptr);
};
auto push_vol_post_proc_ops = [
this](
auto &pp_fe,
auto &&
p) {
auto &pip = pp_fe->getOpPtrVector();
auto [common_plastic_ptr, common_henky_ptr] =
p;
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()}},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{{"GRAD", common_henky_ptr->matGradPtr},
{"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
{{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
} else {
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()}},
{{"U", u_ptr}, {"GEOMETRY", x_ptr}},
{},
{{"STRAIN", common_plastic_ptr->mStrainPtr},
{"STRESS", common_plastic_ptr->mStressPtr},
{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
}
};
auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
PetscBool post_proc_vol = PETSC_FALSE;
&post_proc_vol, PETSC_NULL);
if (post_proc_vol == PETSC_FALSE)
return boost::shared_ptr<PostProcEle>();
auto pp_fe = boost::make_shared<PostProcEle>(mField);
push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
PetscBool post_proc_skin = PETSC_TRUE;
&post_proc_skin, PETSC_NULL);
if (post_proc_skin == PETSC_FALSE)
return boost::shared_ptr<SkinPostProcEle>();
auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
pp_fe->getOpPtrVector().push_back(op_side);
pp_fe, push_vol_ops(op_side->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
return std::make_pair(vol_post_proc(), skin_post_proc());
};
auto scatter_create = [&](
auto D,
auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
auto set_schur_pc = [&](auto solver,
boost::shared_ptr<SetUpSchur> &schur_ptr) {
auto name_prb =
simple->getProblemName();
dm_sub =
createDM(mField.get_comm(),
"DMMOFEM");
for (auto f : {"U"}) {
}
};
if constexpr (
AT == AssemblyType::SCHUR) {
IS is_union_raw;
CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
#ifdef ADD_CONTACT
auto add_sigma_to_is = [&](auto is_union) {
auto add_sigma_to_is_impl = [&]() {
is_sigma);
IS is_union_raw_sigma;
CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);
};
return is_union_sigma;
};
is_union = add_sigma_to_is(is_union);
#endif
schur_ptr =
CHKERR schur_ptr->setUp(solver);
}
};
uXScatter = scatter_create(
D, 0);
uYScatter = scatter_create(
D, 1);
uZScatter = scatter_create(
D, 2);
auto create_solver = [pip_mng]() {
return pip_mng->createTSIM();
else
return pip_mng->createTSIM2();
};
auto solver = create_solver();
auto active_pre_lhs = []() {
};
auto active_post_lhs = [&]() {
auto get_iter = [&]() {
SNES snes;
int iter;
"Can not get iter");
return iter;
};
auto iter = get_iter();
if (iter >= 0) {
std::array<int, 5> activity_data;
std::fill(activity_data.begin(), activity_data.end(), 0);
activity_data.data(), activity_data.size(), MPI_INT,
MPI_SUM, mField.get_comm());
int &active_points = activity_data[0];
int &avtive_full_elems = activity_data[1];
int &avtive_elems = activity_data[2];
int &nb_points = activity_data[3];
int &nb_elements = activity_data[4];
if (nb_points) {
double proc_nb_points =
100 * static_cast<double>(active_points) / nb_points;
double proc_nb_active =
100 * static_cast<double>(avtive_elems) / nb_elements;
double proc_nb_full_active = 100;
if (avtive_elems)
proc_nb_full_active =
100 * static_cast<double>(avtive_full_elems) / avtive_elems;
"Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "
"elements %d "
"(%3.3f\%) nb full active elems %d (%3.3f\%)",
iter, nb_points, active_points, proc_nb_points,
avtive_elems, proc_nb_active, avtive_full_elems,
proc_nb_full_active, iter);
}
}
};
auto add_active_dofs_elem = [&](auto dm) {
auto fe_pre_proc = boost::make_shared<FEMethod>();
fe_pre_proc->preProcessHook = active_pre_lhs;
auto fe_post_proc = boost::make_shared<FEMethod>();
fe_post_proc->postProcessHook = active_post_lhs;
ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
};
auto set_essential_bc = [&](auto dm, auto solver) {
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
auto disp_time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
{disp_time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_FALSE) {
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
}
};
} else {
CHKERR TS2SetSolution(solver,
D, DD);
}
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR add_active_dofs_elem(dm);
boost::shared_ptr<SetUpSchur> schur_ptr;
CHKERR set_schur_pc(solver, schur_ptr);
CHKERR set_essential_bc(dm, solver);
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
Py_Initialize();
#endif
#endif
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
LogManager::setLog("PLASTICITY");
#ifdef ADD_CONTACT
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
#endif
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
#ifdef ADD_CONTACT
#ifdef PYTHON_SFD
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif
return 0;
}
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
}
#ifdef ADD_CONTACT
#endif
}
private:
#ifdef ADD_CONTACT
#endif
};
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
CHKERR KSPSetFromOptions(ksp);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
#ifdef ADD_CONTACT
CHKERR TSSetIJacobian(solver,
A,
P, TsSetIJacobian, ts_ctx_ptr.get());
#endif
auto set_ops = [&]() {
#ifndef ADD_CONTACT
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
));
pip_mng->getOpDomainLhsPipeline().push_back(
));
#else
double eps_stab = 1e-4;
PETSC_NULL);
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
return eps_stab;
}));
pip_mng->getOpBoundaryLhsPipeline().push_back(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
{false, false, false}
));
pip_mng->getOpDomainLhsPipeline().push_back(
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
{false, false, false}
));
#endif
};
auto set_assemble_elems = [&]() {
auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
schur_asmb_pre_proc->preProcessHook = [this]() {
#ifdef ADD_CONTACT
#endif
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble Begin";
};
auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble End";
#ifndef ADD_CONTACT
mField, schur_asmb_post_proc, 1)();
#else
CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
mField, schur_asmb_post_proc, 1,
A)();
CHKERR MatAssemblyBegin(
P, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
P, MAT_FINAL_ASSEMBLY);
CHKERR MatAXPY(
P, 1,
A, SAME_NONZERO_PATTERN);
#endif
CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
};
ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
};
auto set_pc = [&]() {
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
};
} else {
pip_mng->getOpBoundaryLhsPipeline().push_back(
pip_mng->getOpDomainLhsPipeline().push_back(
}
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(
}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
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
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
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.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getDMSubData(DM dm)
Get sub problem data structure.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
constexpr AssemblyType AT
double C1_k
Kinematic hardening.
double Qinf
Saturation yield stress.
constexpr IntegrationType IT
#define EXECUTABLE_DIMENSION
PetscBool is_quasi_static
double visH
Viscous hardening.
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
PetscBool set_timer
Set timer.
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
double zeta
Viscous hardening.
int tau_order
Order of tau files.
double iso_hardening_exp(double tau, double b_iso)
double b_iso
Saturation exponent.
PetscBool is_large_strains
Large strains.
int geom_order
Order if fixed.
double sigmaY
Yield stress.
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
auto kinematic_hardening_dplastic_strain(double C1_k)
ElementsAndOps< SPACE_DIM >::SideEle SideEle
int ep_order
Order of ep files.
constexpr FieldSpace CONTACT_SPACE
FTensor::Index< 'm', 3 > m
Element used to specialise assembly.
double getScale(const double time)
Get scaling at a given time.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
FieldApproximationBase base
MoFEMErrorCode createCommonData()
[Set up problem]
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode bC()
[Create common data]
boost::shared_ptr< DomainEle > reactionFe
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to enforce essential constrains.
Class (Function) to calculate residual side diagonal.
default operator for TRI element
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Clear Schur complement internal data.
Assemble Schur complement.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
[Push operators to pipeline]
static std::array< int, 5 > activityData
[Push operators to pipeline]
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
virtual MoFEMErrorCode setUp(TS solver)=0
SmartPetscObj< DM > subDM
field split sub dm
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
virtual ~SetUpSchurImpl()
MoFEMErrorCode postProc()
MoFEM::Interface & mField
struct CommonData :
public boost::enable_shared_from_this<CommonData> {
};
return boost::shared_ptr<BlockParams>(shared_from_this(), &
blockParams);
};
boost::shared_ptr<MatrixDouble>
mDPtr;
boost::shared_ptr<MatrixDouble>
mGradPtr;
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticSurface);
}
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticTau);
}
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticTauDot);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
plasticStrain);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
plasticFlow);
}
};
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticSurfaceImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculatePlasticityImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpPlasticStressImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowRhsImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dEPImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dTAUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_LogStrain_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dEPImpl;
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dTAUImpl;
template <typename DomainEleOp> struct PlasticityIntegrators {
template <int DIM, IntegrationType I>
OpCalculatePlasticSurfaceImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
template <int DIM, IntegrationType I>
template <AssemblyType A> struct Assembly {
typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowRhsImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsRhsImpl<I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl<
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculatePlasticFlowLhs_dTAUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
template <int DIM, IntegrationType I>
OpCalculateConstraintsLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
template <IntegrationType I>
OpCalculateConstraintsLhs_dTAUImpl<I, AssemblyDomainEleOp>;
};
};
};
using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
template <int DIM>
boost::shared_ptr<MatrixDouble> mat_D_Ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
double scale_value, Sev sev) {
OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
"Can not get data from block");
}
EntitiesFieldData::EntData &data) {
auto getK = [](
auto &
p) {
};
auto getG = [](
auto &
p) {
};
auto scale_fun = [
this](
auto &
p) {
p[CommonData::YOUNG_MODULUS] *= scaleVal;
p[CommonData::SIGMA_Y] *= scaleVal;
p[CommonData::H] *= scaleVal;
p[CommonData::VIS_H] *= scaleVal;
p[CommonData::QINF] *= scaleVal;
p[CommonData::C1_k] *= scaleVal;
};
for (auto &b : blockData) {
if (
b.blockEnts.find(getFEEntityHandle()) !=
b.blockEnts.end()) {
*matParamsPtr =
b.bParams;
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
}
scale_fun(*matParamsPtr);
CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
const double scaleVal;
std::array<double, CommonData::LAST_PARAM> bParams;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != CommonData::LAST_PARAM) {
"Wrong number of block attribute");
}
auto get_block_ents = [&]() {
true);
return ents;
};
CommonData::BlockParams block_params;
for (
auto i = 0;
i != CommonData::LAST_PARAM; ++
i) {
block_params[
i] = block_data[
i];
}
<< "E = " << block_params[CommonData::YOUNG_MODULUS]
<< " nu = " << block_params[CommonData::POISSON_RATIO];
<< std::endl
<< "sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
<< "h = " << block_params[CommonData::H] << std::endl
<< "vis_h = " << block_params[CommonData::VIS_H] << std::endl
<< "qinf = " << block_params[CommonData::QINF] << std::endl
<< "biso = " << block_params[CommonData::BISO] << std::endl
<< "C1_k = " << block_params[CommonData::C1_k] << std::endl;
blockData.push_back({block_params, get_block_ents()});
}
}
auto set_material_stiffness = [&]() {
: 1;
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
};
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
set_material_stiffness();
}
};
pip.push_back(new OpMatBlocks(
mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
));
}
template <int DIM, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string u, std::string ep, std::string tau,
double scale, Sev sev) {
using P = PlasticityIntegrators<DomainEleOp>;
auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
auto make_d_mat = []() {
};
common_plastic_ptr->mDPtr = make_d_mat();
common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
auto m_D_ptr = common_plastic_ptr->mDPtr;
common_plastic_ptr->getParamsPtr(),
"add mat block ops");
pip.push_back(new OpCalculateScalarFieldValues(
tau, common_plastic_ptr->getPlasticTauPtr()));
pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
ep, common_plastic_ptr->getPlasticStrainPtr()));
u, common_plastic_ptr->mGradPtr));
common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
common_henky_ptr->matLogCPlastic =
common_plastic_ptr->getPlasticStrainPtr();
common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
u, common_henky_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
u, common_henky_ptr));
pip.push_back(new
typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
u, common_henky_ptr));
} else {
pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
u, common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
u, common_plastic_ptr));
return std::make_tuple(common_plastic_ptr, common_henky_ptr);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string u, std::string ep, std::string tau) {
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau,
scale, Sev::inform);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
if (common_henky_ptr) {
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
DIM,
I>(ep, common_plastic_ptr, m_D_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string u, std::string ep, std::string tau) {
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
using OpKCauchy =
typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
using P = PlasticityIntegrators<DomainEleOp>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau,
scale, Sev::verbose);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr));
if (common_henky_ptr) {
pip.push_back(new typename H::template OpHenckyTangent<DIM, I>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(
new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
u, ep, common_plastic_ptr, m_D_ptr));
}
if (common_henky_ptr) {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
} else {
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
DIM,
I>(tau, u, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
DIM,
I>(ep, u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
DIM,
I>(ep, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
DIM,
I>(ep, tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
DIM,
I>(tau, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
I>(tau, tau, common_plastic_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string block_name,
Pip &pip,
std::string u, std::string ep,
std::string tau) {
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
auto [common_plastic_ptr, common_henky_ptr] =
createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
ep, tau, 1., Sev::inform);
if (common_henky_ptr) {
u, common_henky_ptr->getMatFirstPiolaStress()));
} else {
}
}
}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
FTensor::Index< 'M', 3 > M
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
FTensor::Index< 'N', 3 > N
FTensor::Index< 'I', 3 > I
[Common data]
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'J', 3 > J
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
boost::shared_ptr< MatrixDouble > mStrainPtr
MatrixDouble resFlowDstrainDot
auto getPlasticStrainPtr()
auto getPlasticStrainDotPtr()
boost::shared_ptr< MatrixDouble > mGradPtr
MatrixDouble resCdPlasticStrain
VectorDouble plasticTauDot
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
VectorDouble plasticSurface
[Common data set externally]
auto getPlasticTauDotPtr()
MatrixDouble plasticStrain
boost::shared_ptr< MatrixDouble > mStressPtr
MatrixDouble plasticStrainDot
MatrixDouble resFlowDstrain
auto getPlasticSurfacePtr()
std::array< double, LAST_PARAM > BlockParams
OpCalculatePlasticFlowRhsImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowRhs
OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_LogStrain_dU
OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_LogStrain_dU
OpCalculatePlasticFlowLhs_dTAUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dTAU
OpCalculateConstraintsLhs_dTAUImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dTAU
OpCalculatePlasticInternalForceLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_dEP
OpCalculatePlasticFlowLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dU
OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticInternalForceLhs_LogStrain_dEP
OpCalculatePlasticFlowLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculatePlasticFlowLhs_dEP
OpCalculateConstraintsRhsImpl< I, AssemblyDomainEleOp > OpCalculateConstraintsRhs
OpCalculateConstraintsLhs_dUImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dU
OpCalculateConstraintsLhs_dEPImpl< DIM, I, AssemblyDomainEleOp > OpCalculateConstraintsLhs_dEP
OpPlasticStressImpl< DIM, I, DomainEleOp > OpPlasticStress
OpCalculatePlasticSurfaceImpl< DIM, I, DomainEleOp > OpCalculatePlasticSurface
OpCalculatePlasticityImpl< DIM, I, DomainEleOp > OpCalculatePlasticity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]