12 #ifndef EXECUTABLE_DIMENSION
13 #define EXECUTABLE_DIMENSION 3
22 using namespace MoFEM;
75 return MatSetValues<AssemblyTypeSelector<AT>>(
76 op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
85 return MatSetValues<AssemblyTypeSelector<AT>>(
86 op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
108 return MatSetValues<AssemblyTypeSelector<AT>>(
109 op_ptr->
getKSPB(), row_data, col_data,
m, ADD_VALUES);
112 #endif // ADD_CONTACT
116 std::max(
static_cast<double>(std::numeric_limits<float>::min_exponent10),
130 auto r = [&](
auto tau) {
133 constexpr
double eps = 1e-12;
134 return std::max(
r(tau),
eps *
r(0));
140 template <
typename T,
int DIM>
147 if (
C1_k < std::numeric_limits<double>::epsilon()) {
151 t_alpha(
i,
j) =
C1_k * t_plastic_strain(
i,
j);
200 #include <boost/python.hpp>
201 #include <boost/python/def.hpp>
202 #include <boost/python/numpy.hpp>
203 namespace bp = boost::python;
204 namespace np = boost::python::numpy;
214 #endif // ADD_CONTACT
243 boost::shared_ptr<DomainEle> reactionFe;
258 boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
260 #endif // ADD_CONTACT
266 CHKERR createCommonData();
281 CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, domain_ents,
283 auto get_ents_by_dim = [&](
const auto dim) {
289 CHKERR mField.get_moab().get_connectivity(domain_ents, ents,
true);
291 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents,
true);
296 auto get_base = [&]() {
297 auto domain_ents = get_ents_by_dim(
SPACE_DIM);
298 if (domain_ents.empty())
316 const auto base = get_base();
327 PetscBool order_edge = PETSC_FALSE;
330 PetscBool order_face = PETSC_FALSE;
333 PetscBool order_volume = PETSC_FALSE;
337 if (order_edge || order_face || order_volume) {
339 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order edge " << order_edge
342 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order face " << order_face
345 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Order volume " << order_volume
349 auto ents = get_ents_by_dim(0);
351 ents.merge(get_ents_by_dim(1));
353 ents.merge(get_ents_by_dim(2));
355 ents.merge(get_ents_by_dim(3));
371 auto get_skin = [&]() {
373 CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
374 Skinner skin(&mField.get_moab());
376 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
380 auto filter_blocks = [&](
auto skin) {
381 bool is_contact_block =
false;
386 (boost::format(
"%s(.*)") %
"CONTACT").str()
395 <<
"Find contact block set: " <<
m->getName();
396 auto meshset =
m->getMeshset();
397 Range contact_meshset_range;
398 CHKERR mField.get_moab().get_entities_by_dimension(
399 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
402 contact_meshset_range);
403 contact_range.merge(contact_meshset_range);
405 if (is_contact_block) {
407 <<
"Nb entities in contact surface: " << contact_range.size();
409 skin = intersect(skin, contact_range);
414 auto filter_true_skin = [&](
auto skin) {
416 ParallelComm *pcomm =
418 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
419 PSTATUS_NOT, -1, &boundary_ents);
420 return boundary_ents;
423 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
431 auto project_ho_geometry = [&]() {
433 return mField.loop_dofs(
"GEOMETRY", ent_method);
435 PetscBool project_geometry = PETSC_TRUE;
437 &project_geometry, PETSC_NULL);
438 if (project_geometry){
439 CHKERR project_ho_geometry();
450 auto get_command_line_parameters = [&]() {
475 PetscBool tau_order_is_set;
478 PetscBool ep_order_is_set;
491 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Hardening " <<
H;
492 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Viscous hardening " <<
visH;
493 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation yield stress " <<
Qinf;
494 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation exponent " <<
b_iso;
495 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Kinematic hardening " <<
C1_k;
500 if (tau_order_is_set == PETSC_FALSE)
502 if (ep_order_is_set == PETSC_FALSE)
505 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Approximation order " <<
order;
507 <<
"Ep approximation order " <<
ep_order;
509 <<
"Tau approximation order " <<
tau_order;
511 <<
"Geometry approximation order " <<
geom_order;
513 MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Density " <<
rho;
516 PetscBool is_scale = PETSC_TRUE;
530 #endif // ADD_CONTACT
540 CHKERR get_command_line_parameters();
544 sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
545 CHKERR sdfPythonPtr->sdfInit(
"sdf.py");
546 ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
548 #endif // ADD_CONTACT
562 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
564 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
566 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
568 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
569 "REMOVE_ALL",
"U", 0, 3);
572 for (
auto b : {
"FIX_X",
"REMOVE_X"})
573 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), b,
574 "SIGMA", 0, 0,
false,
true);
575 for (
auto b : {
"FIX_Y",
"REMOVE_Y"})
576 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), b,
577 "SIGMA", 1, 1,
false,
true);
578 for (
auto b : {
"FIX_Z",
"REMOVE_Z"})
579 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), b,
580 "SIGMA", 2, 2,
false,
true);
581 for (
auto b : {
"FIX_ALL",
"REMOVE_ALL"})
582 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), b,
583 "SIGMA", 0, 3,
false,
true);
584 CHKERR bc_mng->removeBlockDOFsOnEntities(
585 simple->getProblemName(),
"NO_CONTACT",
"SIGMA", 0, 3,
false,
true);
589 simple->getProblemName(),
"U");
591 auto &bc_map = bc_mng->getBcMapByBlockName();
592 for (
auto bc : bc_map)
593 MOFEM_LOG(
"PLASTICITY", Sev::verbose) <<
"Marker " << bc.first;
606 auto integration_rule_bc = [](
int,
int,
int ao) {
return 2 * ao; };
610 auto add_boundary_ops_lhs_mechanical = [&](
auto &pip) {
619 pip, mField,
"U", Sev::inform);
623 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
626 ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
627 mField, pip,
simple->getDomainFEName(),
"SIGMA",
"U",
"GEOMETRY",
629 #endif // ADD_CONTACT
634 auto add_boundary_ops_rhs_mechanical = [&](
auto &pip) {
643 pip, mField,
"U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
646 CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
648 #endif // ADD_CONTACT
653 auto add_domain_ops_lhs = [
this](
auto &pip) {
665 auto get_inertia_and_mass_damping = [
this](
const double,
const double,
669 return (
rho /
scale) * fe_domain_lhs->ts_aa +
672 pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
675 CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
676 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
681 auto add_domain_ops_rhs = [
this](
auto &pip) {
689 {boost::make_shared<ScaledTimeScale>(
"body_force_hist.txt")},
697 AT>::LinearForm<IT>::OpBaseTimesVector<1,
SPACE_DIM, 1>;
700 auto mat_acceleration = boost::make_shared<MatrixDouble>();
702 "U", mat_acceleration));
704 new OpInertiaForce(
"U", mat_acceleration, [](
double,
double,
double) {
708 auto mat_velocity = boost::make_shared<MatrixDouble>();
712 new OpInertiaForce(
"U", mat_velocity, [](
double,
double,
double) {
718 CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
719 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
722 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
724 #endif // ADD_CONTACT
729 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
730 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
733 CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
734 CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
736 CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
737 CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
739 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
740 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
742 auto create_reaction_pipeline = [&](
auto &pip) {
746 CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
747 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
751 reactionFe = boost::make_shared<DomainEle>(mField);
752 reactionFe->getRuleHook = vol_rule;
753 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
754 reactionFe->postProcessHook =
773 static boost::shared_ptr<SetUpSchur> createSetUpSchur(
794 auto set_section_monitor = [&](
auto solver) {
797 CHKERR TSGetSNES(solver, &snes);
798 CHKERR SNESMonitorSet(snes,
801 (
void *)(snes_ctx_ptr.get()),
nullptr);
805 auto create_post_process_elements = [&]() {
806 auto pp_fe = boost::make_shared<PostProcEle>(mField);
807 auto &pip = pp_fe->getOpPtrVector();
809 auto push_vol_ops = [
this](
auto &pip) {
813 auto [common_plastic_ptr, common_henky_ptr] =
814 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
815 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
817 if (common_henky_ptr) {
818 if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
822 return std::make_pair(common_plastic_ptr, common_henky_ptr);
825 auto push_vol_post_proc_ops = [
this](
auto &pp_fe,
auto &&p) {
828 auto &pip = pp_fe->getOpPtrVector();
830 auto [common_plastic_ptr, common_henky_ptr] = p;
834 auto x_ptr = boost::make_shared<MatrixDouble>();
837 auto u_ptr = boost::make_shared<MatrixDouble>();
846 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
849 common_plastic_ptr->getPlasticSurfacePtr()},
850 {
"PLASTIC_MULTIPLIER",
851 common_plastic_ptr->getPlasticTauPtr()}},
853 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
855 {{
"GRAD", common_henky_ptr->matGradPtr},
856 {
"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
858 {{
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
859 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
871 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
874 common_plastic_ptr->getPlasticSurfacePtr()},
875 {
"PLASTIC_MULTIPLIER",
876 common_plastic_ptr->getPlasticTauPtr()}},
878 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
882 {{
"STRAIN", common_plastic_ptr->mStrainPtr},
883 {
"STRESS", common_plastic_ptr->mStressPtr},
884 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
885 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
895 auto vol_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops]() {
896 PetscBool post_proc_vol = PETSC_FALSE;
898 &post_proc_vol, PETSC_NULL);
899 if (post_proc_vol == PETSC_FALSE)
900 return boost::shared_ptr<PostProcEle>();
901 auto pp_fe = boost::make_shared<PostProcEle>(mField);
903 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
904 "push_vol_post_proc_ops");
908 auto skin_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops]() {
909 PetscBool post_proc_skin = PETSC_TRUE;
911 &post_proc_skin, PETSC_NULL);
912 if (post_proc_skin == PETSC_FALSE)
913 return boost::shared_ptr<SkinPostProcEle>();
916 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
919 pp_fe->getOpPtrVector().push_back(op_side);
921 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
922 "push_vol_post_proc_ops");
926 return std::make_pair(vol_post_proc(), skin_post_proc());
929 auto scatter_create = [&](
auto D,
auto coeff) {
931 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
932 ROW,
"U", coeff, coeff, is);
934 CHKERR ISGetLocalSize(is, &loc_size);
936 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
938 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
943 auto set_time_monitor = [&](
auto dm,
auto solver) {
945 boost::shared_ptr<Monitor> monitor_ptr(
946 new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
947 uYScatter, uZScatter));
948 boost::shared_ptr<ForcesAndSourcesCore>
null;
950 monitor_ptr,
null,
null);
954 auto set_schur_pc = [&](
auto solver,
955 boost::shared_ptr<SetUpSchur> &schur_ptr) {
959 auto name_prb =
simple->getProblemName();
965 dm_sub =
createDM(mField.get_comm(),
"DMMOFEM");
970 for (
auto f : {
"U"}) {
989 CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
993 auto add_sigma_to_is = [&](
auto is_union) {
995 auto add_sigma_to_is_impl = [&]() {
1001 IS is_union_raw_sigma;
1002 CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);
1007 return is_union_sigma;
1009 is_union = add_sigma_to_is(is_union);
1010 #endif // ADD_CONTACT
1021 CHKERR schur_ptr->setUp(solver);
1027 auto dm =
simple->getDM();
1030 uXScatter = scatter_create(
D, 0);
1031 uYScatter = scatter_create(
D, 1);
1033 uZScatter = scatter_create(
D, 2);
1035 auto create_solver = [pip_mng]() {
1037 return pip_mng->createTSIM();
1039 return pip_mng->createTSIM2();
1042 auto solver = create_solver();
1044 auto active_pre_lhs = []() {
1051 auto active_post_lhs = [&]() {
1053 auto get_iter = [&]() {
1058 "Can not get iter");
1062 auto iter = get_iter();
1065 std::array<int, 5> activity_data;
1066 std::fill(activity_data.begin(), activity_data.end(), 0);
1068 activity_data.data(), activity_data.size(), MPI_INT,
1069 MPI_SUM, mField.get_comm());
1071 int &active_points = activity_data[0];
1072 int &avtive_full_elems = activity_data[1];
1073 int &avtive_elems = activity_data[2];
1074 int &nb_points = activity_data[3];
1075 int &nb_elements = activity_data[4];
1079 double proc_nb_points =
1080 100 *
static_cast<double>(active_points) / nb_points;
1081 double proc_nb_active =
1082 100 *
static_cast<double>(avtive_elems) / nb_elements;
1083 double proc_nb_full_active = 100;
1085 proc_nb_full_active =
1086 100 *
static_cast<double>(avtive_full_elems) / avtive_elems;
1089 "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "
1091 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1092 iter, nb_points, active_points, proc_nb_points,
1093 avtive_elems, proc_nb_active, avtive_full_elems,
1094 proc_nb_full_active, iter);
1101 auto add_active_dofs_elem = [&](
auto dm) {
1103 auto fe_pre_proc = boost::make_shared<FEMethod>();
1104 fe_pre_proc->preProcessHook = active_pre_lhs;
1105 auto fe_post_proc = boost::make_shared<FEMethod>();
1106 fe_post_proc->postProcessHook = active_post_lhs;
1108 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1109 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1113 auto set_essential_bc = [&](
auto dm,
auto solver) {
1117 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1118 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1119 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1122 auto disp_time_scale = boost::make_shared<TimeScale>();
1124 auto get_bc_hook_rhs = [
this, pre_proc_ptr, disp_time_scale]() {
1126 {disp_time_scale},
false);
1129 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1131 auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
1134 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1136 mField, post_proc_rhs_ptr, 1.)();
1139 auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
1141 mField, post_proc_lhs_ptr, 1.);
1143 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1146 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1147 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1148 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1151 CHKERR TSGetSNES(solver, &snes);
1153 CHKERR SNESGetKSP(snes, &ksp);
1155 CHKERR KSPGetPC(ksp, &pc);
1156 PetscBool is_pcfs = PETSC_FALSE;
1157 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1159 if (is_pcfs == PETSC_FALSE) {
1160 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1161 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1167 CHKERR TSSetSolution(solver,
D);
1169 CHKERR TS2SetSolution(solver,
D, DD);
1172 CHKERR set_section_monitor(solver);
1173 CHKERR set_time_monitor(dm, solver);
1174 CHKERR TSSetFromOptions(solver);
1177 CHKERR add_active_dofs_elem(dm);
1178 boost::shared_ptr<SetUpSchur> schur_ptr;
1179 CHKERR set_schur_pc(solver, schur_ptr);
1180 CHKERR set_essential_bc(dm, solver);
1185 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
1186 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1188 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1189 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1190 CHKERR TSSolve(solver, NULL);
1191 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1206 #endif // ADD_CONTACT
1209 const char param_file[] =
"param_file.petsc";
1213 auto core_log = logging::core::get();
1226 #endif // ADD_CONTACT
1231 DMType dm_name =
"DMMOFEM";
1262 if (Py_FinalizeEx() < 0) {
1266 #endif // ADD_CONTACT
1275 :
SetUpSchur(), mField(m_field), subDM(sub_dm),
1276 fieldSplitIS(field_split_is), aoUp(ao_up) {
1280 "Is expected that schur matrix is not allocated. This is "
1281 "possible only is PC is set up twice");
1288 #endif // ADD_CONTACT
1300 #endif // ADD_CONTACT
1315 CHKERR TSGetSNES(solver, &snes);
1317 CHKERR SNESGetKSP(snes, &ksp);
1318 CHKERR KSPSetFromOptions(ksp);
1321 CHKERR KSPSetFromOptions(ksp);
1322 CHKERR KSPGetPC(ksp, &pc);
1323 PetscBool is_pcfs = PETSC_FALSE;
1324 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1329 "Is expected that schur matrix is not allocated. This is "
1330 "possible only is PC is set up twice");
1338 #endif // ADD_CONTACT
1342 auto set_ops = [&]() {
1348 pip_mng->getOpBoundaryLhsPipeline().push_front(
1350 pip_mng->getOpBoundaryLhsPipeline().push_back(
1358 pip_mng->getOpDomainLhsPipeline().push_front(
1360 pip_mng->getOpDomainLhsPipeline().push_back(
1369 double eps_stab = 1e-4;
1375 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1378 pip_mng->getOpBoundaryLhsPipeline().push_front(
1380 pip_mng->getOpBoundaryLhsPipeline().push_back(
1381 new OpMassStab(
"SIGMA",
"SIGMA", [eps_stab](
double,
double,
double) {
1386 {
"SIGMA",
"EP",
"TAU"}, {
nullptr,
nullptr,
nullptr},
1389 {
false,
false,
false}, false
1393 pip_mng->getOpDomainLhsPipeline().push_front(
1397 {
"SIGMA",
"EP",
"TAU"}, {
nullptr,
nullptr,
nullptr},
1400 {
false,
false,
false}, false
1403 #endif // ADD_CONTACT
1407 auto set_assemble_elems = [&]() {
1409 auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1410 schur_asmb_pre_proc->preProcessHook = [
this]() {
1415 #endif // ADD_CONTACT
1416 CHKERR MatZeroEntries(S);
1417 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble Begin";
1420 auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1422 schur_asmb_post_proc->postProcessHook = [
this, schur_asmb_post_proc]() {
1424 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble End";
1428 mField, schur_asmb_post_proc, 1)();
1429 #else // ADD_CONTACT
1430 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
1431 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
1434 mField, schur_asmb_post_proc, 1,
A)();
1435 CHKERR MatAssemblyBegin(
P, MAT_FINAL_ASSEMBLY);
1436 CHKERR MatAssemblyEnd(
P, MAT_FINAL_ASSEMBLY);
1437 CHKERR MatAXPY(
P, 1,
A, SAME_NONZERO_PATTERN);
1438 #endif // ADD_CONTACT
1441 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1442 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1444 mField, schur_asmb_post_proc, 1, S, aoUp)();
1449 ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1450 ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1454 auto set_pc = [&]() {
1456 CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1457 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1463 CHKERR set_assemble_elems();
1466 pip_mng->getOpBoundaryLhsPipeline().push_front(
1468 pip_mng->getOpBoundaryLhsPipeline().push_back(
1471 pip_mng->getOpDomainLhsPipeline().push_back(
1477 fieldSplitIS.reset();
1482 boost::shared_ptr<SetUpSchur>
1486 return boost::shared_ptr<SetUpSchur>(