13 #ifndef EXECUTABLE_DIMENSION
14 #define EXECUTABLE_DIMENSION 3
17 #ifndef SCHUR_ASSEMBLE
18 #define SCHUR_ASSEMBLE 0
25 #include <boost/python.hpp>
26 #include <boost/python/def.hpp>
27 #include <boost/python/numpy.hpp>
28 namespace bp = boost::python;
29 namespace np = boost::python::numpy;
32 using namespace MoFEM;
76 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
77 op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
86 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
87 op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
109 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
110 op_ptr->
getKSPB(), row_data, col_data,
m, ADD_VALUES);
120 IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
153 #ifdef WITH_MODULE_MFRONT_INTERFACE
154 #include <MFrontMoFEMInterface.hpp>
193 boost::shared_ptr<SDFPython> sdfPythonPtr;
208 CHKERR createCommonData();
221 PetscBool use_mfront = PETSC_FALSE;
241 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
242 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
243 PetscInt choice_base_value = AINSWORTH;
245 LASBASETOPT, &choice_base_value, PETSC_NULL);
248 switch (choice_base_value) {
252 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
257 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
277 auto get_skin = [&]() {
279 CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
280 Skinner skin(&mField.get_moab());
282 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
286 auto filter_blocks = [&](
auto skin) {
287 bool is_contact_block =
false;
292 (boost::format(
"%s(.*)") %
"CONTACT").str()
301 <<
"Find contact block set: " <<
m->getName();
302 auto meshset =
m->getMeshset();
303 Range contact_meshset_range;
304 CHKERR mField.get_moab().get_entities_by_dimension(
305 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
308 contact_meshset_range);
309 contact_range.merge(contact_meshset_range);
311 if (is_contact_block) {
313 <<
"Nb entities in contact surface: " << contact_range.size();
315 skin = intersect(skin, contact_range);
320 auto filter_true_skin = [&](
auto skin) {
322 ParallelComm *pcomm =
324 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
325 PSTATUS_NOT, -1, &boundary_ents);
326 return boundary_ents;
329 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
332 CHKERR simple->setFieldOrder(
"SIGMA", sigma_order, &boundary_ents);
337 CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1,
false, ho_ents,
338 moab::Interface::UNION);
340 ho_ents = boundary_ents;
350 "Use executable contact_2d with axisymmetric model");
354 "Axisymmetric model is only available with MFront (set "
357 MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using axisymmetric model";
362 MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using plane strain model";
369 #ifndef WITH_MODULE_MFRONT_INTERFACE
372 "MFrontInterface module was not found while use_mfront was set to 1");
376 "MFrontInterface module is not compatible with Schur assembly");
380 boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
385 boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
388 mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
394 CHKERR mfrontInterface->getCommandLineParameters();
395 CHKERR mfrontInterface->addElementFields();
396 CHKERR mfrontInterface->createElements();
403 CHKERR mField.build_finite_elements(
"MFRONT_EL");
404 CHKERR mfrontInterface->addElementsToDM(
simple->getDM());
409 auto dm =
simple->getDM();
413 mfrontInterface->setMonitorPtr(monitorPtr);
416 auto project_ho_geometry = [&]() {
418 return mField.loop_dofs(
"GEOMETRY", ent_method);
421 PetscBool project_geometry = PETSC_TRUE;
423 &project_geometry, PETSC_NULL);
424 if (project_geometry){
425 CHKERR project_ho_geometry();
435 auto get_options = [&]() {
452 if (!mfrontInterface) {
456 MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using MFront for material model";
468 PetscBool use_scale = PETSC_FALSE;
488 char sdf_file_name[255];
490 sdf_file_name, 255, PETSC_NULL);
492 sdfPythonPtr = boost::make_shared<SDFPython>();
493 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
494 sdfPythonWeakPtr = sdfPythonPtr;
504 auto bc_mng = mField.getInterface<
BcManager>();
507 for (
auto f : {
"U",
"SIGMA"}) {
508 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
509 "REMOVE_X",
f, 0, 0);
510 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
511 "REMOVE_Y",
f, 1, 1);
512 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
513 "REMOVE_Z",
f, 2, 2);
514 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
515 "REMOVE_ALL",
f, 0, 3);
518 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
519 "SIGMA", 0, 0,
false,
true);
520 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
521 "SIGMA", 1, 1,
false,
true);
522 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
523 "SIGMA", 2, 2,
false,
true);
524 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
525 "SIGMA", 0, 3,
false,
true);
526 CHKERR bc_mng->removeBlockDOFsOnEntities(
527 simple->getProblemName(),
"NO_CONTACT",
"SIGMA", 0, 3,
false,
true);
532 simple->getProblemName(),
"U");
544 auto time_scale = boost::make_shared<ScaledTimeScale>();
545 auto body_force_time_scale =
546 boost::make_shared<ScaledTimeScale>(
"body_force_hist.txt");
555 auto add_domain_base_ops = [&](
auto &pip) {
562 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
563 henky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
564 henky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
566 auto add_domain_ops_lhs = [&](
auto &pip) {
576 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
578 auto get_inertia_and_mass_damping =
580 return (
rho *
scale) * fe_domain_lhs->ts_aa +
583 pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
587 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
589 auto get_inertia_and_mass_damping =
593 pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
597 if (!mfrontInterface) {
598 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
599 mField, pip,
"U",
"MAT_ELASTIC", Sev::verbose,
scale);
605 auto add_domain_ops_rhs = [&](
auto &pip) {
609 pip, mField,
"U", {body_force_time_scale}, Sev::inform);
613 AT>::LinearForm<IT>::OpBaseTimesVector<1,
SPACE_DIM, 1>;
618 auto mat_acceleration = boost::make_shared<MatrixDouble>();
620 "U", mat_acceleration));
622 new OpInertiaForce(
"U", mat_acceleration, [](
double,
double,
double) {
630 auto mat_velocity = boost::make_shared<MatrixDouble>();
634 new OpInertiaForce(
"U", mat_velocity, [](
double,
double,
double) {
639 if (!mfrontInterface) {
640 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
641 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
644 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
650 auto add_boundary_base_ops = [&](
auto &pip) {
660 auto add_boundary_ops_lhs = [&](
auto &pip) {
670 pip, mField,
"U", Sev::inform);
675 auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
680 [
this, fe_boundary_lhs](
double,
double,
double) {
689 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
693 mField, pip,
simple->getDomainFEName(),
"SIGMA",
"U",
"GEOMETRY",
699 auto add_boundary_ops_rhs = [&](
auto &pip) {
704 AT>::LinearForm<IT>::OpBaseTimesVector<1,
SPACE_DIM, 1>;
709 pip, mField,
"U", {time_scale}, Sev::inform);
712 auto u_disp = boost::make_shared<MatrixDouble>();
713 auto dot_u_disp = boost::make_shared<MatrixDouble>();
718 new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
722 new OpSpringRhs(
"U", dot_u_disp, [
this](
double,
double,
double) {
728 ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
734 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
735 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
736 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
737 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
739 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
740 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
741 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
742 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
744 if (mfrontInterface) {
749 CHKERR mfrontInterface->setOperators();
750 CHKERR mfrontInterface->setupSolverFunctionTS(t_type);
751 CHKERR mfrontInterface->setupSolverJacobianTS(t_type);
754 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
755 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
756 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
757 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
765 static boost::shared_ptr<SetUpSchur>
770 SetUpSchur() =
default;
780 auto set_section_monitor = [&](
auto solver) {
783 CHKERR TSGetSNES(solver, &snes);
784 PetscViewerAndFormat *vf;
785 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
786 PETSC_VIEWER_DEFAULT, &vf);
789 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
794 auto scatter_create = [&](
auto D,
auto coeff) {
796 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
797 ROW,
"U", coeff, coeff, is);
799 CHKERR ISGetLocalSize(is, &loc_size);
801 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
803 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
808 auto set_time_monitor = [&](
auto dm,
auto solver) {
810 monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
811 boost::shared_ptr<ForcesAndSourcesCore>
null;
813 monitorPtr,
null,
null);
817 auto set_essential_bc = [&]() {
821 auto pre_proc_ptr = boost::make_shared<FEMethod>();
822 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
823 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
826 auto time_scale = boost::make_shared<TimeScale>();
828 auto get_bc_hook_rhs = [&]() {
830 {time_scale},
false);
833 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
835 auto get_post_proc_hook_rhs = [&]() {
837 mField, post_proc_rhs_ptr, 1.);
839 auto get_post_proc_hook_lhs = [&]() {
841 mField, post_proc_lhs_ptr, 1.);
843 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
845 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
846 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
847 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
849 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
850 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
855 auto set_schur_pc = [&](
auto solver) {
856 boost::shared_ptr<SetUpSchur> schur_ptr;
864 auto dm =
simple->getDM();
869 uXScatter = scatter_create(
D, 0);
870 uYScatter = scatter_create(
D, 1);
872 uZScatter = scatter_create(
D, 2);
876 CHKERR set_essential_bc();
879 auto solver = pip_mng->createTSIM();
880 CHKERR TSSetFromOptions(solver);
883 auto schur_pc_ptr = set_schur_pc(solver);
885 CHKERR set_section_monitor(solver);
886 CHKERR set_time_monitor(dm, solver);
887 CHKERR TSSetSolution(solver,
D);
889 CHKERR TSSolve(solver, NULL);
891 auto solver = pip_mng->createTSIM2();
892 CHKERR TSSetFromOptions(solver);
894 auto dm =
simple->getDM();
897 auto schur_pc_ptr = set_schur_pc(solver);
899 CHKERR set_section_monitor(solver);
900 CHKERR set_time_monitor(dm, solver);
901 CHKERR TS2SetSolution(solver,
D, DD);
903 CHKERR TSSolve(solver, NULL);
913 if (
atom_test == 1 && !mField.get_comm_rank()) {
916 double hertz_tract = 158.73;
918 if (fabs(t_ptr[1] - hertz_tract) / hertz_tract >
tol) {
920 "atom test %d diverged! %3.4e != %3.4e",
atom_test, t_ptr[1],
932 static char help[] =
"...\n\n";
934 int main(
int argc,
char *argv[]) {
942 const char param_file[] =
"param_file.petsc";
946 auto core_log = logging::core::get();
955 DMType dm_name =
"DMMOFEM";
985 if (Py_FinalizeEx() < 0) {
1025 auto dm =
simple->getDM();
1028 CHKERR TSGetSNES(solver, &snes);
1030 CHKERR SNESGetKSP(snes, &ksp);
1031 CHKERR KSPSetFromOptions(ksp);
1034 CHKERR KSPGetPC(ksp, &pc);
1036 PetscBool is_pcfs = PETSC_FALSE;
1037 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1040 MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Setup Schur pc";
1045 "It is expected that Schur matrix is not allocated. This is "
1046 "possible only if PC is set up twice");
1051 subDM = createSubDM();
1062 MOFEM_LOG(
"CONTACT", Sev::inform) <<
"No Schur pc";
1065 pip->getOpBoundaryLhsPipeline().push_back(
1068 pip->getOpDomainLhsPipeline().push_back(
1071 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1072 post_proc_schur_lhs_ptr->postProcessHook = [
this,
1073 post_proc_schur_lhs_ptr]() {
1076 mField, post_proc_schur_lhs_ptr, 1.)();
1080 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1087 auto sub_dm =
createDM(mField.get_comm(),
"DMMOFEM");
1088 auto set_up = [&]() {
1104 double eps_stab = 1e-4;
1110 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1118 pip->getOpBoundaryLhsPipeline().push_back(
1119 new OpMassStab(
"SIGMA",
"SIGMA",
1120 [eps_stab](
double,
double,
double) {
return eps_stab; }));
1122 {
"SIGMA"}, {
nullptr}, {ao_up}, {S}, {
false},
false));
1127 {
"SIGMA"}, {
nullptr}, {ao_up}, {S}, {
false},
false));
1129 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1130 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1132 pre_proc_schur_lhs_ptr->preProcessHook = [
this]() {
1136 CHKERR MatZeroEntries(S);
1137 MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble Begin";
1141 post_proc_schur_lhs_ptr->postProcessHook = [
this, ao_up,
1142 post_proc_schur_lhs_ptr]() {
1144 MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble End";
1146 *post_proc_schur_lhs_ptr->matAssembleSwitch =
false;
1148 auto print_mat_norm = [
this](
auto a, std::string prefix) {
1151 CHKERR MatNorm(
a, NORM_FROBENIUS, &nrm);
1152 MOFEM_LOG(
"CONTACT", Sev::noisy) << prefix <<
" norm = " << nrm;
1156 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
1157 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
1159 mField, post_proc_schur_lhs_ptr, 1,
A)();
1161 CHKERR MatAssemblyBegin(
P, MAT_FINAL_ASSEMBLY);
1162 CHKERR MatAssemblyEnd(
P, MAT_FINAL_ASSEMBLY);
1163 CHKERR MatAXPY(
P, 1,
A, SAME_NONZERO_PATTERN);
1165 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1166 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1169 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1172 CHKERR print_mat_norm(
A,
"A");
1173 CHKERR print_mat_norm(
P,
"P");
1174 CHKERR print_mat_norm(S,
"S");
1177 MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble Finish";
1183 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1184 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1193 mField.getInterface<
ISManager>()->isCreateProblemFieldAndRank(
1195 CHKERR PCFieldSplitSetIS(pc, NULL, is);
1196 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1200 boost::shared_ptr<SetUpSchur>
1202 return boost::shared_ptr<SetUpSchur>(
new SetUpSchurImpl(m_field));