46 std::pair<boost::shared_ptr<PostProcEle>,
47 boost::shared_ptr<SkinPostProcEle>>
52 std::array<double, SPACE_DIM> pass_field_eval_coords,
53 boost::shared_ptr<SetPtsData> pass_field_eval_data,
54 boost::shared_ptr<MatrixDouble> vec_field_ptr)
55 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
56 uZScatter(uz_scatter), fieldEvalCoords(pass_field_eval_coords),
57 fieldEvalData(pass_field_eval_data), vecFieldPtr(vec_field_ptr) {
58 postProcFe = pair_post_proc_fe.first;
59 skinPostProcFe = pair_post_proc_fe.second;
76 fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
77 simple->getDomainFEName(), fieldEvalData,
83 fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
84 simple->getDomainFEName(), fieldEvalData,
89 if (vecFieldPtr->size1()) {
90 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vecFieldPtr);
93 <<
"U_X: " << t_disp(0) <<
" U_Y: " << t_disp(1);
96 <<
"U_X: " << t_disp(0) <<
" U_Y: " << t_disp(1)
97 <<
" U_Z: " << t_disp(2);
102 auto make_vtk = [&]() {
107 CHKERR postProcFe->writeFile(
"out_incomp_elasticity" +
108 boost::lexical_cast<std::string>(ts_step) +
111 if (skinPostProcFe) {
114 CHKERR skinPostProcFe->writeFile(
115 "out_skin_incomp_elasticity_" +
116 boost::lexical_cast<std::string>(ts_step) +
".h5m");
121 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
123 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
124 INSERT_VALUES, SCATTER_FORWARD);
125 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
126 INSERT_VALUES, SCATTER_FORWARD);
128 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
129 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
131 "%s time %3.4e min %3.4e max %3.4e", msg.c_str(), ts_t, min,
137 CHKERR print_max_min(uXScatter,
"Ux");
138 CHKERR print_max_min(uYScatter,
"Uy");
140 CHKERR print_max_min(uZScatter,
"Uz");
166 return MatSetValues<AssemblyTypeSelector<AT>>(
167 op_ptr->
getKSPA(), row_data, col_data,
m, ADD_VALUES);
189 return MatSetValues<AssemblyTypeSelector<AT>>(
190 op_ptr->
getKSPB(), row_data, col_data,
m, ADD_VALUES);
198 inline static double mu;
227 boost::shared_ptr<MatrixDouble> strain_ptr,
228 boost::shared_ptr<VectorDouble> pressure_ptr)
230 stressPtr(stress_ptr), strainPtr(strain_ptr),
231 pressurePtr(pressure_ptr) {}
244 const size_t nb_gauss_pts = getGaussPts().size2();
246 stressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
247 auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(stressPtr));
248 auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(strainPtr));
251 const double l_mu = mU;
253 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
255 t_stress(
i,
j) = t_pressure *
t_kd(
i,
j) + 2. * l_mu * t_strain(
i,
j);
276 CHKERR createCommonData();
300 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
301 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
302 PetscInt choice_base_value = AINSWORTH;
304 LASBASETOPT, &choice_base_value, PETSC_NULL);
307 switch (choice_base_value) {
310 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::inform)
311 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
315 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::inform)
316 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
351 auto project_ho_geometry = [&]() {
353 return mField.loop_dofs(
"GEOMETRY", ent_method);
355 CHKERR project_ho_geometry();
364 auto get_options = [&]() {
371 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::inform)
373 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::inform)
377 const double lambda_denom =
397 auto time_scale = boost::make_shared<TimeScale>();
413 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE},
"GEOMETRY");
416 pipeline_mng->getOpBoundaryRhsPipeline(), mField,
"U", {time_scale},
417 "FORCE", Sev::inform);
421 pipeline_mng->getOpDomainRhsPipeline(), mField,
"U", {time_scale},
422 "BODY_FORCE", Sev::inform);
425 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
427 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
429 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
432 simple->getProblemName(),
"U");
447 auto add_domain_base_ops = [&](
auto &pip) {
454 auto add_domain_ops_lhs = [&](
auto &pip) {
472 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
482 auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
487 [](
const double,
const double,
const double) constexpr {
return -1.; },
489 pip.push_back(
new OpGradSymTensorGrad(
"U",
"U", mat_D_ptr));
495 pip.push_back(
new OpMassPressure(
"P",
"P", get_lambda_reciprocal));
501 pip.push_back(
new OpMassPressureStab(
502 "P",
"P", [eps_stab](
double,
double,
double) {
return eps_stab; }));
507 auto add_domain_ops_rhs = [&](
auto &pip) {
518 AT>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
522 auto pressure_ptr = boost::make_shared<VectorDouble>();
525 auto div_u_ptr = boost::make_shared<VectorDouble>();
529 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
533 auto strain_ptr = boost::make_shared<MatrixDouble>();
543 pip.push_back(
new OpDomainGradTimesTensor(
"U", strain_ptr, get_four_mu));
545 pip.push_back(
new OpDivDeltaUTimesP(
"U", pressure_ptr, minus_one));
547 pip.push_back(
new OpBaseTimesScalarValues(
"P", div_u_ptr, minus_one));
553 pip.push_back(
new OpBaseTimesScalarValues(
"P", pressure_ptr,
554 get_lambda_reciprocal));
560 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
561 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
562 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
563 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
565 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
566 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
574 static boost::shared_ptr<SetUpSchur>
589 auto set_section_monitor = [&](
auto solver) {
592 CHKERR TSGetSNES(solver, &snes);
593 PetscViewerAndFormat *vf;
594 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
595 PETSC_VIEWER_DEFAULT, &vf);
598 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
603 auto scatter_create = [&](
auto D,
auto coeff) {
605 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
606 ROW,
"U", coeff, coeff, is);
608 CHKERR ISGetLocalSize(is, &loc_size);
610 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
612 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
617 auto create_post_process_elements = [&]() {
618 auto pp_fe = boost::make_shared<PostProcEle>(mField);
620 auto push_vol_ops = [
this](
auto &pip) {
628 auto push_vol_post_proc_ops = [
this](
auto &pp_fe,
auto &&p) {
631 auto &pip = pp_fe->getOpPtrVector();
635 auto x_ptr = boost::make_shared<MatrixDouble>();
638 auto u_ptr = boost::make_shared<MatrixDouble>();
641 auto pressure_ptr = boost::make_shared<VectorDouble>();
644 auto div_u_ptr = boost::make_shared<VectorDouble>();
648 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
652 auto strain_ptr = boost::make_shared<MatrixDouble>();
655 auto stress_ptr = boost::make_shared<MatrixDouble>();
657 mu, stress_ptr, strain_ptr, pressure_ptr));
663 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
665 {{
"P", pressure_ptr}},
667 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
671 {{
"STRAIN", strain_ptr}, {
"STRESS", stress_ptr}}
680 auto vol_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops]() {
681 PetscBool post_proc_vol = PETSC_FALSE;
683 &post_proc_vol, PETSC_NULL);
684 if (post_proc_vol == PETSC_FALSE)
685 return boost::shared_ptr<PostProcEle>();
686 auto pp_fe = boost::make_shared<PostProcEle>(mField);
688 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
689 "push_vol_post_proc_ops");
693 auto skin_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops]() {
694 PetscBool post_proc_skin = PETSC_TRUE;
696 &post_proc_skin, PETSC_NULL);
697 if (post_proc_skin == PETSC_FALSE)
698 return boost::shared_ptr<SkinPostProcEle>();
701 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
704 pp_fe->getOpPtrVector().push_back(op_side);
706 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
707 "push_vol_post_proc_ops");
711 return std::make_pair(vol_post_proc(), skin_post_proc());
714 boost::shared_ptr<SetPtsData> field_eval_data;
715 boost::shared_ptr<MatrixDouble> vector_field_ptr;
717 std::array<double, SPACE_DIM> field_eval_coords;
720 field_eval_coords.data(), &coords_dim,
724 vector_field_ptr = boost::make_shared<MatrixDouble>();
729 field_eval_data,
simple->getDomainFEName());
732 field_eval_data,
simple->getDomainFEName());
735 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
736 auto no_rule = [](
int,
int,
int) {
return -1; };
737 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
738 field_eval_fe_ptr->getRuleHook = no_rule;
739 field_eval_fe_ptr->getOpPtrVector().push_back(
743 auto set_time_monitor = [&](
auto dm,
auto solver) {
745 boost::shared_ptr<MonitorIncompressible> monitor_ptr(
747 create_post_process_elements(), uXScatter,
748 uYScatter, uZScatter, field_eval_coords,
749 field_eval_data, vector_field_ptr));
750 boost::shared_ptr<ForcesAndSourcesCore>
null;
752 monitor_ptr,
null,
null);
756 auto set_essential_bc = [&]() {
760 auto pre_proc_ptr = boost::make_shared<FEMethod>();
761 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
764 auto time_scale = boost::make_shared<TimeScale>();
767 mField, pre_proc_ptr, {time_scale},
false);
768 post_proc_rhs_ptr->postProcessHook =
771 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
772 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
773 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
777 auto set_schur_pc = [&](
auto solver) {
779 CHKERR TSGetSNES(solver, &snes);
781 CHKERR SNESGetKSP(snes, &ksp);
783 CHKERR KSPGetPC(ksp, &pc);
784 PetscBool is_pcfs = PETSC_FALSE;
785 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
786 boost::shared_ptr<SetUpSchur> schur_ptr;
795 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
796 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
797 pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
799 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Zero matrices";
800 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
801 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
804 post_proc_schur_lhs_ptr->postProcessHook = [
this,
805 post_proc_schur_lhs_ptr]() {
807 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Assemble Begin";
808 *(post_proc_schur_lhs_ptr->matAssembleSwitch) =
false;
809 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
810 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
812 mField, post_proc_schur_lhs_ptr, 1.,
815 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
816 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
817 CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
818 SAME_NONZERO_PATTERN);
819 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Assemble End";
822 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
823 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
825 if (is_pcfs == PETSC_TRUE) {
828 CHK_MOAB_THROW(schur_ptr->setUp(solver),
"setup schur preconditioner");
830 auto set_pcfieldsplit_preconditioned_ts = [&](
auto solver) {
832 auto name_prb =
simple->getProblemName();
835 name_prb,
ROW,
"P", 0, 1, is_p);
836 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_p);
840 "set pcfieldsplit preconditioned");
842 return boost::make_tuple(schur_ptr,
A, B);
845 return boost::make_tuple(schur_ptr,
A, B);
848 auto dm =
simple->getDM();
851 uXScatter = scatter_create(
D, 0);
852 uYScatter = scatter_create(
D, 1);
854 uZScatter = scatter_create(
D, 2);
858 CHKERR set_essential_bc();
860 auto solver = pip_mng->createTSIM();
861 CHKERR TSSetFromOptions(solver);
863 CHKERR set_section_monitor(solver);
864 CHKERR set_time_monitor(dm, solver);
865 auto [schur_pc_ptr,
A, B] = set_schur_pc(solver);
867 CHKERR TSSetSolution(solver,
D);
869 CHKERR TSSolve(solver, NULL);
879 static char help[] =
"...\n\n";
881 int main(
int argc,
char *argv[]) {
884 const char param_file[] =
"param_file.petsc";
888 auto core_log = logging::core::get();
890 LogManager::createSink(LogManager::getStrmWorld(),
"INCOMP_ELASTICITY"));
891 LogManager::setLog(
"INCOMP_ELASTICITY");
897 DMType dm_name =
"DMMOFEM";
952 auto dm =
simple->getDM();
955 CHKERR TSGetSNES(solver, &snes);
957 CHKERR SNESGetKSP(snes, &ksp);
959 CHKERR KSPGetPC(ksp, &pc);
961 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::inform) <<
"Setup Schur pc";
965 "Is expected that schur matrix is not allocated. This is "
966 "possible only is PC is set up twice");
969 auto create_sub_dm = [&]() {
971 auto set_up = [&]() {
984 subDM = create_sub_dm();
992 pip->getOpDomainLhsPipeline().push_back(
995 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
996 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
998 pre_proc_schur_lhs_ptr->preProcessHook = [
this]() {
1004 post_proc_schur_lhs_ptr->postProcessHook = [
this, ao_up,
1005 post_proc_schur_lhs_ptr]() {
1008 CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
1009 CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
1011 mField, post_proc_schur_lhs_ptr, 1,
S, ao_up)();
1014 auto print_mat_norm = [
this](
auto a, std::string prefix) {
1017 CHKERR MatNorm(
a, NORM_FROBENIUS, &nrm);
1018 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose)
1019 << prefix <<
" norm = " << nrm;
1022 CHKERR print_mat_norm(
S,
"S");
1028 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1029 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1033 simple->getProblemName(),
ROW,
"P", 0, 1, is_p);
1034 CHKERR PCFieldSplitSetIS(pc, NULL, is_p);
1035 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
1040 boost::shared_ptr<SetUpSchur>
1042 return boost::shared_ptr<SetUpSchur>(
new SetUpSchurImpl(m_field));