SetUpSchurImpl.cpp

Go to the documentation of this file.

/** @file SetUpSchurImpl
 * @brief
 * @date 2023-05-13
 *
 * @license{This project is released under the MIT License.}
 *
 */
 
struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
      : SetUpSchur(), mField(m_field), epCorePtr(ep_core_ptr) {}
  virtual ~SetUpSchurImpl() {}
 
  MoFEMErrorCode setUp(TS ts);
  MoFEMErrorCode preProc();
  MoFEMErrorCode postProc();
 
private:
  MoFEM::Interface &mField;
  EshelbianCore *epCorePtr;
 
  SmartPetscObj<Mat> S;
  SmartPetscObj<AO> aoS;
 
  boost::shared_ptr<std::vector<boost::weak_ptr<NumeredDofEntity>>>
      piolaZeroDofsVec; //< Dofs on crack surface
  boost::shared_ptr<std::vector<unsigned char>>
      piolaZeroDofsMarker; //< marker for crack dofs on surface
 
  struct P_MultiGridData {
 
    P_MultiGridData(SmartPetscObj<DM> dm, SmartPetscObj<PC> pc,
                    SmartPetscObj<Mat> S)
        : dm(dm), pc(pc), S(S) {}
 
    MoFEMErrorCode setUP() {
      MoFEMFunctionBegin;
      CHKERR PCMGSetUpViaApproxOrders(
          pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
      CHKERR PCSetFromOptions(pc);
      MoFEMFunctionReturn(0);
    };
 
  private:
    SmartPetscObj<DM> dm;
    SmartPetscObj<PC> pc;
    SmartPetscObj<Mat> S;
  };
 
  boost::shared_ptr<P_MultiGridData> pMGPtr;
 
  auto getSchurFields() {
    std::vector<std::string> schur_field_list{epCorePtr->hybridSpatialDisp,
                                              epCorePtr->contactDisp};
    std::vector<boost::shared_ptr<Range>> dm_range_list{nullptr, nullptr};
    return std::make_pair(schur_field_list, dm_range_list);
  };
 
  auto getA00Fields() {
    auto a00_field_list =
        epCorePtr->physicalEquations->getMaterialFields(*epCorePtr);
    for (const auto &field : {epCorePtr->bubbleField, epCorePtr->piolaStress,
                              epCorePtr->rotAxis, epCorePtr->spatialL2Disp})
      if (mField.check_field(field))
        a00_field_list.push_back(field);
    std::vector<boost::shared_ptr<Range>> range_list_ptr(a00_field_list.size(),
                                                         nullptr);
    return std::make_pair(a00_field_list, range_list_ptr);
  }
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(TS ts) {
  MoFEMFunctionBegin;
 
  auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
    MoFEMFunctionBegin;
 
    dm_sub = createDM(mField.get_comm(), "DMMOFEM_MG");
    CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
    CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->skeletonElement);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
 
    int r_idx = 0;
    auto [schur_field_list, schur_range_list] = getSchurFields();
    for (auto f : schur_field_list) {
      MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
      CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, schur_range_list[r_idx]);
      CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, schur_range_list[r_idx]);
      ++r_idx;
    }
    CHKERR DMSetUp(dm_sub);
    MoFEMFunctionReturn(0);
  };
 
  auto create_a00_dm = [&](SmartPetscObj<DM> &dm_sub) {
    MoFEMFunctionBegin;
    dm_sub = createDM(mField.get_comm(), "DMMOFEM");
    CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_A00");
    CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->skeletonElement);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
 
    int r_idx = 0;
    auto [a00_field_list, a00_range_list] = getA00Fields();
    for (auto f : a00_field_list) {
      MOFEM_LOG("EP", Sev::inform) << "Add a00 field: " << f;
      CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, a00_range_list[r_idx]);
      CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, a00_range_list[r_idx]);
      ++r_idx;
    }
    CHKERR DMSetUp(dm_sub);
    MoFEMFunctionReturn(0);
  };
 
  auto get_snes = [&](TS ts) {
    SNES snes;
    CHKERR TSGetSNES(ts, &snes);
    return snes;
  };
 
  auto get_ksp = [&](SNES snes) {
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    CHKERR KSPSetFromOptions(ksp);
    return ksp;
  };
 
  auto get_pc = [&](KSP ksp) {
    PC pc;
    CHKERR KSPGetPC(ksp, &pc);
    return pc;
  };
 
  auto ksp = get_ksp(get_snes(ts));
  auto pc = get_pc(ksp);
 
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (epCorePtr->alphaOmegaPrec > 0. && !is_pcfs)
    SETERRQ(mField.get_comm(), MOFEM_INVALID_DATA,
            "-alpha_omega_prec requires the elastic PCFIELDSPLIT Schur solver");
 
  if (is_pcfs) {
    const bool use_rotation_prec = epCorePtr->alphaOmegaPrec > 0.;
    if (use_rotation_prec) {
      PCCompositeType split_type;
      CHKERR PCFieldSplitGetType(pc, &split_type);
      if (split_type != PC_COMPOSITE_SCHUR)
        SETERRQ(mField.get_comm(), MOFEM_INVALID_DATA,
                "-alpha_omega_prec requires -elastic_pc_fieldsplit_type schur");
    }
 
    MOFEM_LOG("EP", Sev::inform) << "SetUpSchurImpl::setUp: PCFIELDSPLIT";
 
    SmartPetscObj<DM> schur_dm, a00_dm;
    CHKERR create_schur_dm(schur_dm);
    CHKERR create_a00_dm(a00_dm);
 
    auto dm_elastic = epCorePtr->dmElastic;
    auto vol_elem_name = epCorePtr->elementVolumeName;
    auto skel_elem_name = epCorePtr->skeletonElement;
    auto contact_elem_name = epCorePtr->contactElement;
    auto natural_bc_element_name = epCorePtr->naturalBcElement;
 
    std::vector<std::pair<std::string, std::string>> mat_block_list = {
 
        {epCorePtr->piolaStress, epCorePtr->piolaStress},
        {epCorePtr->piolaStress, epCorePtr->spatialL2Disp},
        {epCorePtr->spatialL2Disp, epCorePtr->piolaStress},
        {epCorePtr->piolaStress, epCorePtr->rotAxis},
        {epCorePtr->rotAxis, epCorePtr->piolaStress},
        {epCorePtr->bubbleField, epCorePtr->rotAxis},
        {epCorePtr->rotAxis, epCorePtr->bubbleField},
        {epCorePtr->spatialL2Disp, epCorePtr->spatialL2Disp}
 
    };
 
    mat_block_list.push_back({epCorePtr->rotAxis, epCorePtr->rotAxis});
    auto add_pair = [&](const std::string &row, const std::string &col) {
      if (mField.check_field(row) && mField.check_field(col)) {
        mat_block_list.emplace_back(row, col);
        if (row != col)
          mat_block_list.emplace_back(col, row);
      }
    };
    for (const auto &[row, col] :
         epCorePtr->physicalEquations->getMaterialCouplings(*epCorePtr))
      add_pair(row, col);
 
    auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
      auto block_mat_data = createBlockMatStructure(
          dm_elastic,
 
          {
 
              {vol_elem_name, mat_block_list},
 
              {skel_elem_name,
 
               {
 
                   {epCorePtr->piolaStress, epCorePtr->hybridSpatialDisp},
                   {epCorePtr->hybridSpatialDisp, epCorePtr->piolaStress},
                   {epCorePtr->hybridSpatialDisp, epCorePtr->hybridSpatialDisp},
                   {epCorePtr->hybridSpatialDisp, epCorePtr->spatialL2Disp},
                   {epCorePtr->spatialL2Disp, epCorePtr->hybridSpatialDisp}
 
               }},
 
              {contact_elem_name,
 
               {
 
                   {epCorePtr->contactDisp, epCorePtr->contactDisp},
                   {epCorePtr->piolaStress, epCorePtr->contactDisp},
                   {epCorePtr->contactDisp, epCorePtr->piolaStress},
                   {epCorePtr->contactDisp, epCorePtr->spatialL2Disp},
                   {epCorePtr->spatialL2Disp, epCorePtr->contactDisp}
 
               }},
 
              {natural_bc_element_name,
 
               {
 
                   {epCorePtr->hybridSpatialDisp, epCorePtr->hybridSpatialDisp},
                   {epCorePtr->piolaStress, epCorePtr->hybridSpatialDisp},
                   {epCorePtr->hybridSpatialDisp, epCorePtr->piolaStress}
 
               }}
 
          }
 
      );
 
      auto [a00_field_list, a00_range_list] = getA00Fields();
 
      return createSchurNestedMatrixStruture(
 
          {schur_dm, a00_dm}, block_mat_data,
 
          a00_field_list,
 
          a00_range_list,
 
          use_rotation_prec
 
      );
    };
 
    auto nested_mat_data = get_nested_mat_data(schur_dm, a00_dm);
    CHKERR DMMoFEMSetNestSchurData(epCorePtr->dmElastic, nested_mat_data);
    const auto block_mat_data = std::get<2>(*nested_mat_data);
    if (use_rotation_prec)
      epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
          new OpSpatialRotationPrec(epCorePtr->rotAxis, epCorePtr->dataAtPts,
                                    epCorePtr->alphaOmegaPrec));
    CHKERR DMSetMatType(epCorePtr->dmElastic, MATSHELL);
 
    auto m = createDMBlockMat(epCorePtr->dmElastic);
    auto p = createDMNestSchurMat(epCorePtr->dmElastic);
 
    if (std::abs(epCorePtr->alphaRho) >
        std::numeric_limits<double>::epsilon()) {
      auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
                              PetscReal a, PetscReal aa, Mat A, Mat B,
                              void *ctx) {
        return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
      };
      auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
      CHKERR TSSetI2Jacobian(ts, m, p, swap_assemble, ts_ctx_ptr.get());
    } else {
      auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
                              Mat A, Mat B, void *ctx) {
        return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
      };
      auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
      CHKERR TSSetIJacobian(ts, m, p, swap_assemble, ts_ctx_ptr.get());
    }
    CHKERR KSPSetOperators(ksp, m, p);
 
    auto set_assembly = [&]() {
      MoFEMFunctionBegin;
 
      aoS = getDMSubData(schur_dm)->getSmartRowMap();
      S = createDMHybridisedL2Matrix(schur_dm);
      auto set_block_size = [&]<FieldApproximationBase Base>() {
        MoFEMFunctionBegin;
        using Orders = EshelbianCore::FieldOrders<Base>;
        CHKERR MatSetBlockSize(
            S, SPACE_DIM * NBFACETRI_L2(Orders::hybrid(epCorePtr->spaceOrder)));
        MoFEMFunctionReturn(0);
      };
      CHKERR epCorePtr->withFieldOrders(set_block_size);
      epCorePtr->S = S;
      epCorePtr->aoS = aoS;
 
      auto set_assemble = [&]() {
        MoFEMFunctionBegin;
        auto schur_asmb_pre_proc_lhs = boost::make_shared<FEMethod>();
        auto schur_asmb_pre_proc_rhs = boost::make_shared<FEMethod>();
 
        schur_asmb_pre_proc_lhs->preProcessHook = [this]() {
          MoFEMFunctionBegin;
          CHKERR MatZeroEntries(S);
          MoFEMFunctionReturn(0);
        };
 
        auto *schur_asmb_pre_proc_rhs_ptr = schur_asmb_pre_proc_rhs.get();
        schur_asmb_pre_proc_rhs->preProcessHook =
            [this, schur_asmb_pre_proc_rhs_ptr]() {
              MoFEMFunctionBegin;
              auto prb_ptr = schur_asmb_pre_proc_rhs_ptr->problemPtr;
              auto dofs_prb = prb_ptr->getNumeredRowDofsPtr();
 
              auto crack_faces = epCorePtr->crackFaces;
              piolaZeroDofsVec->clear();
              CHKERR mField.getInterface<ProblemsManager>()
                  ->getSideDofsOnBrokenSpaceEntities(
                      *piolaZeroDofsVec, prb_ptr->getName(), ROW,
                      epCorePtr->piolaStress, *crack_faces, SPACE_DIM, 0,
                      SPACE_DIM);
 
              piolaZeroDofsMarker->resize(dofs_prb->size(), 0);
              piolaZeroDofsMarker->clear();
              for (auto &dof : *piolaZeroDofsVec) {
                if (auto dof_ptr = dof.lock()) {
                  auto idx = dof_ptr->getPetscLocalDofIdx();
                  (*piolaZeroDofsMarker)[idx] = 1;
                }
              }
              // FIXME: This needs investigation, true is required for fracture
              constexpr bool hard_coded_set_bc_debug = true;
              if constexpr (hard_coded_set_bc_debug) {
 
                auto problem_name =
                    schur_asmb_pre_proc_rhs_ptr->problemPtr->getName();
                auto crack_faces = epCorePtr->crackFaces;
 
                SmartPetscObj<IS> crack_hybrid_is;
                CHKERR epCorePtr->mField.getInterface<ISManager>()
                    ->isCreateProblemFieldAndRankLocal(
                        problem_name, ROW, epCorePtr->hybridSpatialDisp, 0,
                        SPACE_DIM, crack_hybrid_is, &*crack_faces);
 
                SmartPetscObj<IS> crack_piola_is;
                CHKERR epCorePtr->mField.getInterface<ISManager>()
                    ->isCreateProblemBrokenFieldAndRankLocal(*piolaZeroDofsVec,
                                                             crack_piola_is);
 
                const double *a_x;
                CHKERR VecGetArrayRead(schur_asmb_pre_proc_rhs_ptr->x, &a_x);
                auto zero_by_is = [&](auto is) {
                  MoFEMFunctionBegin;
                  const PetscInt *is_array;
                  PetscInt is_size;
                  CHKERR ISGetLocalSize(is, &is_size);
                  CHKERR ISGetIndices(is, &is_array);
                  for (int i = 0; i != is_size; ++i) {
                    // FIXME: That is irregular, need investigate if needed
                    const_cast<double *>(a_x)[is_array[i]] = 0;
                  }
                  CHKERR ISRestoreIndices(is, &is_array);
                  MoFEMFunctionReturn(0);
                };
 
                CHKERR zero_by_is(crack_hybrid_is);
                CHKERR zero_by_is(crack_piola_is);
 
                CHKERR VecRestoreArrayRead(schur_asmb_pre_proc_rhs_ptr->x,
                                           &a_x);
 
                CHKERR epCorePtr->mField.getInterface<VecManager>()
                    ->setLocalGhostVector(problem_name, COL,
                                          schur_asmb_pre_proc_rhs_ptr->x,
                                          INSERT_VALUES, SCATTER_REVERSE);
              }
 
              MoFEMFunctionReturn(0);
            };
 
        auto schur_asmb_post_proc_lhs = boost::make_shared<FEMethod>();
        auto schur_asmb_post_proc_rhs = boost::make_shared<FEMethod>();
 
        auto schur_asmb_post_proc_rhs_weak_ptr =
            boost::weak_ptr<FEMethod>(schur_asmb_post_proc_rhs);
        schur_asmb_post_proc_rhs->postProcessHook =
            [this, schur_asmb_post_proc_rhs_weak_ptr]() {
              MoFEMFunctionBegin;
 
              if (auto schur_asmb_post_proc_rhs_ptr =
                      schur_asmb_post_proc_rhs_weak_ptr.lock()) {
 
                CHKERR VecGhostUpdateBegin(schur_asmb_post_proc_rhs_ptr->f,
                                           ADD_VALUES, SCATTER_REVERSE);
                CHKERR VecGhostUpdateEnd(schur_asmb_post_proc_rhs_ptr->f,
                                         ADD_VALUES, SCATTER_REVERSE);
                CHKERR VecAssemblyBegin(schur_asmb_post_proc_rhs_ptr->f);
                CHKERR VecAssemblyEnd(schur_asmb_post_proc_rhs_ptr->f);
                *(schur_asmb_post_proc_rhs_ptr->vecAssembleSwitch) = false;
 
                {
 
                  auto problem_name =
                      schur_asmb_post_proc_rhs_ptr->problemPtr->getName();
 
                  auto crack_faces = epCorePtr->crackFaces;
 
                  SmartPetscObj<IS> crack_hybrid_is;
                  CHKERR epCorePtr->mField.getInterface<ISManager>()
                      ->isCreateProblemFieldAndRankLocal(
                          problem_name, ROW, epCorePtr->hybridSpatialDisp, 0,
                          SPACE_DIM, crack_hybrid_is, &*crack_faces);
 
                  SmartPetscObj<IS> crack_piola_is;
                  CHKERR epCorePtr->mField.getInterface<ISManager>()
                      ->isCreateProblemBrokenFieldAndRankLocal(
                          *piolaZeroDofsVec, crack_piola_is);
 
                  double *a_f;
                  CHKERR VecGetArray(schur_asmb_post_proc_rhs_ptr->f, &a_f);
                  const double *a_x;
                  CHKERR VecGetArrayRead(schur_asmb_post_proc_rhs_ptr->x, &a_x);
                  auto zero_by_is = [&](auto is) {
                    MoFEMFunctionBegin;
                    const PetscInt *is_array;
                    PetscInt is_size;
                    CHKERR ISGetLocalSize(is, &is_size);
                    CHKERR ISGetIndices(is, &is_array);
                    for (int i = 0; i != is_size; ++i) {
                      a_f[is_array[i]] = -a_x[is_array[i]];
                    }
                    CHKERR ISRestoreIndices(is, &is_array);
                    MoFEMFunctionReturn(0);
                  };
 
                  CHKERR zero_by_is(crack_hybrid_is);
                  CHKERR zero_by_is(crack_piola_is);
 
                  CHKERR VecRestoreArray(schur_asmb_post_proc_rhs_ptr->f, &a_f);
                  CHKERR VecRestoreArrayRead(schur_asmb_post_proc_rhs_ptr->x,
                                             &a_x);
                }
 
              } else {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
                        "schur_asmb_post_proc_rhs is expired");
              }
 
              MoFEMFunctionReturn(0);
            };
 
        auto schur_asmb_post_proc_lhs_weak_ptr =
            boost::weak_ptr<FEMethod>(schur_asmb_post_proc_lhs);
        schur_asmb_post_proc_lhs->postProcessHook =
            [this, schur_asmb_post_proc_lhs_weak_ptr, block_mat_data,
             use_rotation_prec]() {
              MoFEMFunctionBegin;
 
              if (auto schur_asmb_post_proc_lhs_ptr =
                      schur_asmb_post_proc_lhs_weak_ptr.lock()) {
 
                if (pMGPtr) {
                  CHKERR pMGPtr->setUP();
                  pMGPtr.reset();
                }
 
                auto crack_faces = epCorePtr->crackFaces;
 
                // Assemble matrix
                CHKERR MatAssemblyBegin(schur_asmb_post_proc_lhs_ptr->B,
                                        MAT_FINAL_ASSEMBLY);
                CHKERR MatAssemblyEnd(schur_asmb_post_proc_lhs_ptr->B,
                                      MAT_FINAL_ASSEMBLY);
                *(schur_asmb_post_proc_lhs_ptr->matAssembleSwitch) = false;
                {
                  SmartPetscObj<IS> crack_hybrid_is;
                  CHKERR epCorePtr->mField.getInterface<ISManager>()
                      ->isCreateProblemFieldAndRank(
                          "ELASTIC_PROBLEM", ROW, epCorePtr->hybridSpatialDisp,
                          0, SPACE_DIM, crack_hybrid_is, &*crack_faces);
                  CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs_ptr->B,
                                              crack_hybrid_is, 1, PETSC_NULLPTR,
                                              PETSC_NULLPTR);
                }
                {
                  SmartPetscObj<IS> crack_piola_is;
                  CHKERR epCorePtr->mField.getInterface<ISManager>()
                      ->isCreateProblemBrokenFieldAndRank(*piolaZeroDofsVec,
                                                          crack_piola_is);
                  CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs_ptr->B,
                                              crack_piola_is, 1, PETSC_NULLPTR,
                                              PETSC_NULLPTR);
                }
 
                auto [a00_field_list, a00_range_list] = getA00Fields();
                // Use the same shifted A00 for its solve and the Schur PC.
                // Restore the exact parent operator even if assembly fails.
                if (use_rotation_prec)
                  CHKERR schurSwitchPreconditioner(block_mat_data);
                const auto schur_error = assembleBlockMatSchur(
                    epCorePtr->mField, schur_asmb_post_proc_lhs_ptr->B, S,
                    a00_field_list, a00_range_list, aoS);
                if (use_rotation_prec)
                  CHKERR schurSwitchPreconditioner(block_mat_data);
                CHKERR schur_error;
                epCorePtr->a00FieldList = a00_field_list;
                epCorePtr->a00RangeList = a00_range_list;
 
                // Apply essential constrains to Schur complement
                CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
                CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 
                SmartPetscObj<IS> crack_hybrid_is;
                CHKERR epCorePtr->mField.getInterface<ISManager>()
                    ->isCreateProblemFieldAndRank(
                        "SUB_SCHUR", ROW, epCorePtr->hybridSpatialDisp, 0,
                        SPACE_DIM, crack_hybrid_is, &*crack_faces);
                epCorePtr->crackHybridIs = crack_hybrid_is;
                CHKERR MatZeroRowsColumnsIS(S, crack_hybrid_is, 1,
                                            PETSC_NULLPTR, PETSC_NULLPTR);
 
              } else {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
                        "schur_asmb_post_proc_lhs is expired");
              }
 
              MoFEMFunctionReturn(0);
            };
 
        auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
        ts_ctx_ptr->getPreProcessIFunction().push_front(
            schur_asmb_pre_proc_rhs);
        ts_ctx_ptr->getPostProcessIFunction().push_back(
            schur_asmb_post_proc_rhs);
        ts_ctx_ptr->getPreProcessIJacobian().push_front(
            schur_asmb_pre_proc_lhs);
        ts_ctx_ptr->getPostProcessIJacobian().push_back(
            schur_asmb_post_proc_lhs);
        MoFEMFunctionReturn(0);
      };
 
      piolaZeroDofsVec =
          boost::make_shared<std::vector<boost::weak_ptr<NumeredDofEntity>>>();
      piolaZeroDofsMarker = boost::make_shared<std::vector<unsigned char>>();
      CHKERR set_assemble();
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&]() {
      MoFEMFunctionBegin;
      auto a00_is = getDMSubData(a00_dm)->getSmartRowIs();
      auto schur_is = getDMSubData(schur_dm)->getSmartRowIs();
      CHKERR PCFieldSplitSetIS(pc, NULL, a00_is);
      CHKERR PCFieldSplitSetIS(pc, NULL, schur_is);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    auto set_diagonal_pc = [&]() {
      MoFEMFunctionBegin;
      KSP *subksp;
      CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
      auto get_pc = [](auto ksp) {
        PC pc_raw;
        CHKERR KSPGetPC(ksp, &pc_raw);
        return SmartPetscObj<PC>(pc_raw, true); // bump reference
      };
      CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
 
      auto set_pc_p_mg = [&](auto dm, auto pc, auto S) {
        MoFEMFunctionBegin;
        CHKERR PCSetDM(pc, dm);
        PetscBool same = PETSC_FALSE;
        PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
        if (same) {
          auto smart_pc = SmartPetscObj<PC>(pc, true);
          pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc, S);
        }
        PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
        if (same) {
          MOFEM_LOG("EP", Sev::inform)
              << "SetUpSchurImpl::setUp: fieldsplit 1 PCKSP";
          CHKERR PCSetFromOptions(pc);
          KSP ksp;
          CHKERR PCKSPGetKSP(pc, &ksp);
          CHKERR KSPSetFromOptions(ksp);
          PC ksp_pc;
          CHKERR KSPGetPC(ksp, &ksp_pc);
          CHKERR PCSetFromOptions(ksp_pc);
          PetscObjectTypeCompare((PetscObject)ksp_pc, PCMG, &same);
          if (same) {
            auto smart_pc = SmartPetscObj<PC>(ksp_pc, true);
            pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc, S);
          }
        }
        MoFEMFunctionReturn(0);
      };
 
      CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]), S);
 
      CHKERR PetscFree(subksp);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_assembly();
    CHKERR set_pc();
    CHKERR TSSetUp(ts);
    CHKERR KSPSetUp(ksp);
    CHKERR set_diagonal_pc();
 
  } else {
    MOFEM_LOG("EP", Sev::inform) << "SetUpSchurImpl::setUp: PCLU or other";
 
    auto m = createDMMatrix(epCorePtr->dmElastic);
    // Plastic coupling may extend the elastic matrix graph on first assembly.
    if (epCorePtr->plasticVolume)
      CHKERR MatSetOption(m, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
    CHKERR TSSetIJacobian(ts, m, m, TsSetIJacobian, ts_ctx_ptr.get());
 
    epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
        createOpSchurAssembleBegin());
    epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
        createOpSchurAssembleEnd({}, {}));
    epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
        createOpSchurAssembleBegin());
    epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
        createOpSchurAssembleEnd({}, {}));
  }
 
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<EshelbianCore::SetUpSchur>
EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
                                            EshelbianCore *ep_core_ptr) {
  return boost::shared_ptr<SetUpSchur>(
      new SetUpSchurImpl(m_field, ep_core_ptr));
}