#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#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_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
};
#endif // ADD_CONTACT
private:
boost::shared_ptr<DomainEle> reactionFe;
double getScale(const double time) {
};
};
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif // ADD_CONTACT
};
}
CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, domain_ents,
true);
auto get_ents_by_dim = [&](const auto dim) {
return domain_ents;
} else {
if (dim == 0)
CHKERR mField.get_moab().get_connectivity(domain_ents, ents,
true);
else
CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents,
true);
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(
SPACE_DIM);
if (domain_ents.empty())
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 mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
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;
CHKERR mField.get_moab().get_entities_by_dimension(
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 = [&]() {
return mField.loop_dofs("GEOMETRY", ent_method);
};
PetscBool project_geometry = PETSC_TRUE;
&project_geometry, PETSC_NULL);
if (project_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 // ADD_CONTACT
};
CHKERR get_command_line_parameters();
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdfPythonPtr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif // ADD_CONTACT
}
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 // ADD_CONTACT
};
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 // ADD_CONTACT
};
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");
pip, mField, "U",
{boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
Sev::inform);
AT>::LinearForm<IT>::OpBaseTimesVector<1,
SPACE_DIM, 1>;
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 // ADD_CONTACT
};
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");
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = vol_rule;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
reactionFe->postProcessHook =
}
static boost::shared_ptr<SetUpSchur> createSetUpSchur(
);
protected:
};
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
(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");
}
};
dm_sub =
createDM(mField.get_comm(),
"DMMOFEM");
#ifdef ADD_CONTACT
for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
}
#else
for (
auto f : {
"EP",
"TAU"}) {
}
#endif
};
#ifdef ADD_CONTACT
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
{
{simple->getDomainFEName(),
{{"U", "U"},
{"SIGMA", "SIGMA"},
{"U", "SIGMA"},
{"SIGMA", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "EP"},
{"TAU", "U"}
}},
{simple->getBoundaryFEName(),
{{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
}}
}
);
{dm_schur, dm_block}, block_mat_data,
{"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
);
};
#else
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
auto block_mat_data =
{{simple->getDomainFEName(),
{{"U", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "U"},
{"TAU", "EP"}
}}}
);
{dm_schur, dm_block}, block_mat_data,
{"EP", "TAU"}, {nullptr, nullptr}, false
);
};
#endif
auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
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);
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_SDF
Py_Initialize();
np::initialize();
#endif
#endif // ADD_CONTACT
const char param_file[] = "param_file.petsc";
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 // ADD_CONTACT
try {
DMType dm_name = "DMMOFEM";
}
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif // ADD_CONTACT
return 0;
}
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
}
private:
};
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
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");
}
DM solver_dm;
CHKERR TSGetDM(solver, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t, PetscReal
a,
Mat
A, Mat B,
void *ctx) {
};
CHKERR TSSetIJacobian(solver,
A,
P, swap_assemble, ts_ctx_ptr.get());
} else {
auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t,
Vec utt,
PetscReal
a, PetscReal aa, Mat
A, Mat B,
void *ctx) {
};
CHKERR TSSetI2Jacobian(solver,
A,
P, swap_assemble, ts_ctx_ptr.get());
}
auto set_ops = [&]() {
#ifndef ADD_CONTACT
pip_mng->getOpBoundaryLhsPipeline().push_front(
{
"EP",
"TAU"}, {
nullptr,
nullptr},
aoSchur,
S,
false, false
));
pip_mng->getOpDomainLhsPipeline().push_front(
{
"EP",
"TAU"}, {
nullptr,
nullptr},
aoSchur,
S,
false, false
));
#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;
}));
{
"SIGMA",
"EP",
"TAU"}, {
nullptr,
nullptr,
nullptr},
aoSchur,
S,
false, false
));
pip_mng->getOpDomainLhsPipeline().push_front(
{
"SIGMA",
"EP",
"TAU"}, {
nullptr,
nullptr,
nullptr},
aoSchur,
S,
false, false
));
#endif // ADD_CONTACT
};
auto set_assemble_elems = [&]() {
auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
schur_asmb_pre_proc->preProcessHook = [this]() {
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";
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);
};
auto set_diagonal_pc = [&]() {
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
};
};
} else {
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
pip_mng->getOpDomainLhsPipeline().push_back(
}
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(
}
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 {
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");
}
auto getK = [](auto &p) {
};
auto getG = [](auto &p) {
};
auto scale_fun = [this](auto &p) {
};
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,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
"Wrong number of block attribute");
}
auto get_block_ents = [&]() {
true);
return ents;
};
block_params[
i] = block_data[
i];
}
<< 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, 0>(
u, common_henky_ptr, m_D_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
u, common_henky_ptr));
} else {
pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
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<
A>::template LinearForm<I>;
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 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, 0>(
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) {
A>::template LinearForm<I>;
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 {
}
}
}