SetUpSchurImpl Struct Reference

Inheritance diagram for SetUpSchurImpl:

Collaboration diagram for SetUpSchurImpl:

Classes
struct	P_MultiGridData

Public Member Functions
	SetUpSchurImpl (MoFEM::Interface &m_field)
virtual	~SetUpSchurImpl ()=default
MoFEMErrorCode	setUp (SmartPetscObj< KSP >)
	SetUpSchurImpl (MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
virtual	~SetUpSchurImpl ()
MoFEMErrorCode	setUp (TS solver)
MoFEMErrorCode	preProc ()
MoFEMErrorCode	postProc ()
	SetUpSchurImpl (MoFEM::Interface &m_field)
virtual	~SetUpSchurImpl ()
MoFEMErrorCode	setUp (SmartPetscObj< TS > solver)
	SetUpSchurImpl (MoFEM::Interface &m_field)
virtual	~SetUpSchurImpl ()
MoFEMErrorCode	setUp (SmartPetscObj< TS > solver)
	SetUpSchurImpl (MoFEM::Interface &m_field)
virtual	~SetUpSchurImpl ()=default
MoFEMErrorCode	setUp (SmartPetscObj< KSP > solver)
	SetUpSchurImpl (MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
virtual	~SetUpSchurImpl ()
MoFEMErrorCode	setUp (TS ts)
MoFEMErrorCode	preProc ()
MoFEMErrorCode	postProc ()

Private Member Functions
MoFEMErrorCode	createSubDM ()
MoFEMErrorCode	setOperator ()
MoFEMErrorCode	setPC (PC pc)
MoFEMErrorCode	setDiagonalPC (PC pc)
MoFEMErrorCode	setOperator ()
MoFEMErrorCode	setPC (PC pc)
SmartPetscObj< DM >	createSubDM ()
MoFEMErrorCode	setEntities ()
MoFEMErrorCode	createSubDM ()
MoFEMErrorCode	setOperator ()
MoFEMErrorCode	setPC (PC pc)
MoFEMErrorCode	setDiagonalPC (PC pc)
auto	getSchurFields ()
auto	getA00Fields ()

Private Attributes
MoFEM::Interface &	mField
SmartPetscObj< Mat >	S
SmartPetscObj< DM >	subDM
	field split sub dm More...
SmartPetscObj< IS >	fieldSplitIS
	IS for split Schur block. More...
SmartPetscObj< AO >	aoSchur
SmartPetscObj< DM >	schurDM
SmartPetscObj< DM >	blockDM
Range	volEnts
Range	subEnts
EshelbianCore *	epCorePtr
SmartPetscObj< AO >	aoS
boost::shared_ptr< std::vector< boost::weak_ptr< NumeredDofEntity > > >	piolaZeroDofsVec
boost::shared_ptr< std::vector< unsigned char > >	piolaZeroDofsMarker
boost::shared_ptr< P_MultiGridData >	pMGPtr

Additional Inherited Members
Static Public Member Functions inherited from SetUpSchur
static boost::shared_ptr< SetUpSchur >	createSetUpSchur (MoFEM::Interface &m_field)
static boost::shared_ptr< SetUpSchur >	createSetUpSchur (MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
	Create data structure for handling Schur complement. More...
static boost::shared_ptr< SetUpSchur >	createSetUpSchur (MoFEM::Interface &m_field)
static boost::shared_ptr< SetUpSchur >	createSetUpSchur (MoFEM::Interface &m_field)
static boost::shared_ptr< SetUpSchur >	createSetUpSchur (MoFEM::Interface &m_field)
Protected Member Functions inherited from SetUpSchur
	SetUpSchur ()=default
virtual	~SetUpSchur ()=default
	SetUpSchur ()=default
	SetUpSchur ()=default
	SetUpSchur ()=default
	SetUpSchur ()=default

Detailed Description

Examples

plastic.cpp, and test_broken_space.cpp.

Definition at line 519 of file test_broken_space.cpp.

Constructor & Destructor Documentation

SetUpSchurImpl() [1/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field )

inline

Examples

plastic.cpp, and test_broken_space.cpp.

Definition at line 521 of file test_broken_space.cpp.

: SetUpSchur(), mField(m_field) {}

~SetUpSchurImpl() [1/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

virtualdefault

Examples

plastic.cpp, and test_broken_space.cpp.

SetUpSchurImpl() [2/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field, SmartPetscObj< DM >  sub_dm, SmartPetscObj< IS >  field_split_is, SmartPetscObj< AO >  ao_up  )

inline

Definition at line 1639 of file plastic.cpp.

      : SetUpSchur(), mField(m_field), subDM(sub_dm),
        fieldSplitIS(field_split_is), aoSchur(ao_up) {
    if (S) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Is expected that schur matrix is not "
                        "allocated. This is "
                        "possible only is PC is set up twice");
    }
  }

~SetUpSchurImpl() [2/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

inlinevirtual

Definition at line 1650 of file plastic.cpp.

{ S.reset(); }

SetUpSchurImpl() [3/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field )

inline

Definition at line 988 of file contact.cpp.

: SetUpSchur(), mField(m_field) {}

~SetUpSchurImpl() [3/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

inlinevirtual

Definition at line 990 of file contact.cpp.

{}

SetUpSchurImpl() [4/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field )

inline

Definition at line 931 of file incompressible_elasticity.cpp.

: SetUpSchur(), mField(m_field) {}

~SetUpSchurImpl() [4/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

inlinevirtual

Definition at line 932 of file incompressible_elasticity.cpp.

{ S.reset(); }

SetUpSchurImpl() [5/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field )

inline

Definition at line 924 of file elastic.cpp.

                                          : SetUpSchur(), mField(m_field) {
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is PC is set up twice");
    }
  }

~SetUpSchurImpl() [5/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

virtualdefault

SetUpSchurImpl() [6/6]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface &  m_field, EshelbianCore *  ep_core_ptr  )

inline

Definition at line 11 of file SetUpSchurImpl.cpp.

      : SetUpSchur(), mField(m_field), epCorePtr(ep_core_ptr) {}

~SetUpSchurImpl() [6/6]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )

inlinevirtual

Definition at line 13 of file SetUpSchurImpl.cpp.

{}

Member Function Documentation

createSubDM() [1/3]

MoFEMErrorCode SetUpSchurImpl::createSubDM ( )

private

createSubDM() [2/3]

SmartPetscObj<DM> SetUpSchurImpl::createSubDM ( )

private

createSubDM() [3/3]

MoFEMErrorCode SetUpSchurImpl::createSubDM ( )

private

Definition at line 1082 of file contact.cpp.

                                           {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  auto create_dm = [&](const char *name, const char *field_name, auto dm_type) {
    auto dm = createDM(mField.get_comm(), dm_type);
    auto create_dm_imp = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
      CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
      CHKERR DMMoFEMAddSubFieldRow(dm, field_name);
      CHKERR DMMoFEMAddSubFieldCol(dm, field_name);
      CHKERR DMSetUp(dm);
      MoFEMFunctionReturn(0);
    };
    CHK_THROW_MESSAGE(
        create_dm_imp(),
        "Error in creating schurDM. It is possible that schurDM is "
        "already created");
    return dm;
  };
 
  // Note: here we can make block with bubbles of "U" and "SIGMA" fields. See
  // vec-0 where bubbles are added.
 
  schurDM = create_dm("SCHUR", "U", "DMMOFEM_MG");
  blockDM = create_dm("BLOCK", "SIGMA", "DMMOFEM");
 
  if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
 
    auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
      auto block_mat_data = createBlockMatStructure(
          simple->getDM(),
 
          {{
 
              simple->getDomainFEName(),
 
              {
 
                  {"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
 
              }}}
 
      );
 
      return createSchurNestedMatrixStruture(
 
          {schur_dm, block_dm}, block_mat_data,
 
          {"SIGMA"}, {nullptr}, true
 
      );
    };
 
    auto nested_mat_data = get_nested_mat_data(schurDM, blockDM);
    CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
 
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only BLOCK_SCHUR is implemented");
  }
 
  MoFEMFunctionReturn(0);
}

getA00Fields()

auto SetUpSchurImpl::getA00Fields ( )

inlineprivate

Definition at line 60 of file SetUpSchurImpl.cpp.

                      {
    std::vector<std::string> a00_field_list{
 
        epCorePtr->stretchTensor,
 
        epCorePtr->bubbleField,
 
        epCorePtr->piolaStress,
 
        epCorePtr->rotAxis,
 
        epCorePtr->spatialL2Disp
 
    };
    std::vector<boost::shared_ptr<Range>> range_list_ptr(a00_field_list.size(),
                                                         nullptr);
    return std::make_pair(a00_field_list, range_list_ptr);
  }

getSchurFields()

auto SetUpSchurImpl::getSchurFields ( )

inlineprivate

Definition at line 53 of file SetUpSchurImpl.cpp.

                        {
    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);
  };

postProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::postProc

(

)

postProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::postProc

(

)

Examples

plastic.cpp.

preProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::preProc

(

)

preProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::preProc

(

)

Examples

plastic.cpp.

setDiagonalPC() [1/2]

MoFEMErrorCode SetUpSchurImpl::setDiagonalPC ( PC  pc )

private

setDiagonalPC() [2/2]

MoFEMErrorCode SetUpSchurImpl::setDiagonalPC ( PC  pc )

private

Definition at line 1233 of file contact.cpp.

                                                  {
  MoFEMFunctionBegin;
  KSP *subksp;
  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &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) {
      CHKERR PCMGSetUpViaApproxOrders(
          pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
      CHKERR PCSetFromOptions(pc);
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
 
  CHKERR PetscFree(subksp);
  MoFEMFunctionReturn(0);
}

setEntities()

MoFEMErrorCode SetUpSchurImpl::setEntities ( )

private

Definition at line 997 of file elastic.cpp.

                                           {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
                                                     SPACE_DIM, volEnts);
  CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
                                                  subEnts);
  subEnts = subtract(subEnts, volEnts);
  MoFEMFunctionReturn(0);
};

setOperator() [1/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )

private

setOperator() [2/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )

private

setOperator() [3/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )

private

Definition at line 1149 of file contact.cpp.

                                           {
  MoFEMFunctionBegin;
 
  double eps_stab = 1e-4;
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
                               PETSC_NULL);
 
  using B = FormsIntegrators<BoundaryEleOp>::Assembly<
      BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
 
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
 
  auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
 
  // Boundary
  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
  pip->getOpBoundaryLhsPipeline().push_back(
      new OpMassStab("SIGMA", "SIGMA",
                     [eps_stab](double, double, double) { return eps_stab; }));
  pip->getOpBoundaryLhsPipeline().push_back(
 
      createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
 
  );
 
  // Domain
  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
  pip->getOpDomainLhsPipeline().push_back(
 
      createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false,
                               false)
 
  );
 
  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
 
  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
    MoFEMFunctionBegin;
    CHKERR MatZeroEntries(S);
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
    MoFEMFunctionReturn(0);
  };
 
  post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
                                              post_proc_schur_lhs_ptr]() {
    MoFEMFunctionBegin;
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
    auto print_mat_norm = [this](auto a, std::string prefix) {
      MoFEMFunctionBegin;
      double nrm;
      CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
      MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
      MoFEMFunctionReturn(0);
    };
    CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
    CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
        mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
#ifndef NDEBUG
    CHKERR print_mat_norm(S, "S");
#endif // NDEBUG
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
    MoFEMFunctionReturn(0);
  };
 
  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
 
  MoFEMFunctionReturn(0);
}

setPC() [1/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc )

private

setPC() [2/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc )

private

setPC() [3/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc )

private

Definition at line 1225 of file contact.cpp.

                                          {
  MoFEMFunctionBegin;
  auto block_is = getDMSubData(blockDM)->getSmartRowIs();
  CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
  MoFEMFunctionReturn(0);
}

setUp() [1/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< KSP >  solver )

virtual

Implements SetUpSchur.

setUp() [2/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< KSP >  ksp )

virtual

Implements SetUpSchur.

Examples

plastic.cpp, and test_broken_space.cpp.

Definition at line 532 of file test_broken_space.cpp.

                                                           {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
 
  CHKERR KSPSetFromOptions(ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
 
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
 
    MOFEM_LOG("AT", Sev::inform) << "Setup Schur pc";
 
    auto create_sub_dm = [&]() {
      auto simple = mField.getInterface<Simple>();
 
      auto create_dm = [&](
 
                           std::string problem_name,
                           std::vector<std::string> fe_names,
                           std::vector<std::string> fields,
 
                           auto dm_type
 
                       ) {
        auto dm = createDM(mField.get_comm(), dm_type);
        auto create_dm_imp = [&]() {
          MoFEMFunctionBegin;
          CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), problem_name.c_str());
          CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
          for (auto fe : fe_names) {
            CHKERR DMMoFEMAddElement(dm, fe);
          }
          CHKERR DMMoFEMAddElement(dm, simple->getSkeletonFEName());
          for (auto field : fields) {
            CHKERR DMMoFEMAddSubFieldRow(dm, field);
            CHKERR DMMoFEMAddSubFieldCol(dm, field);
          }
          CHKERR DMSetUp(dm);
          MoFEMFunctionReturn(0);
        };
        CHK_THROW_MESSAGE(
            create_dm_imp(),
            "Error in creating schurDM. It is possible that schurDM is "
            "already created");
        return dm;
      };
 
      auto schur_dm = create_dm(
 
          "SCHUR",
 
          {simple->getDomainFEName(), simple->getSkeletonFEName()},
 
          {"HYBRID"},
 
          "DMMOFEM_MG");
 
      auto block_dm = create_dm(
 
          "BLOCK",
 
          {simple->getDomainFEName(), simple->getSkeletonFEName()},
 
          {"BROKEN", "U"},
 
          "DMMOFEM");
 
      return std::make_tuple(schur_dm, block_dm);
    };
 
    auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
      auto block_mat_data = createBlockMatStructure(
          simple->getDM(),
 
          {
 
              {
 
                  simple->getDomainFEName(),
 
                  {
 
                      {"BROKEN", "BROKEN"},
                      {"U", "U"},
                      {"BROKEN", "U"},
                      {"U", "BROKEN"}
 
                  }
 
              },
 
              {
 
                  simple->getSkeletonFEName(),
 
                  {
 
                      {"BROKEN", "HYBRID"}, {"HYBRID", "BROKEN"}
 
                  }
 
              }
 
          }
 
      );
 
      return createSchurNestedMatrixStruture(
 
          {schur_dm, block_dm}, block_mat_data,
 
          {"BROKEN", "U"}, {nullptr, nullptr}, true
 
      );
    };
 
    auto set_ops = [&](auto schur_dm) {
      MoFEMFunctionBegin;
      auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
      auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
 
      boost::shared_ptr<BlockStructure> block_data;
      CHKERR DMMoFEMGetBlocMatData(simple->getDM(), block_data);
 
      if (AT == BLOCK_SCHUR) {
        pip_mng->getOpDomainLhsPipeline().push_front(
            createOpSchurAssembleBegin());
        pip_mng->getOpDomainLhsPipeline().push_back(
 
            createOpSchurAssembleEnd({"BROKEN", "U"}, {nullptr, nullptr}, ao_up,
                                     S, true, true)
 
        );
      }
 
      auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
      auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
 
      pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
        MoFEMFunctionBegin;
        CHKERR MatZeroEntries(S);
        MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble Begin";
        MoFEMFunctionReturn(0);
      };
 
      post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
                                                  post_proc_schur_lhs_ptr]() {
        MoFEMFunctionBegin;
        MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble End";
 
        if (AT == BLOCK_SCHUR) {
          CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
          CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
          if (1) {
            auto S_from_block = matDuplicate(S, MAT_SHARE_NONZERO_PATTERN);
            // Create matrix from block mat
            CHKERR assembleBlockMatSchur(mField, post_proc_schur_lhs_ptr->B,
                                 S_from_block, {"BROKEN", "U"},
                                 {nullptr, nullptr}, ao_up);
            CHKERR MatAssemblyBegin(S_from_block, MAT_FINAL_ASSEMBLY);
            CHKERR MatAssemblyEnd(S_from_block, MAT_FINAL_ASSEMBLY);
            CHKERR MatAYPX(S_from_block, -1, S, DIFFERENT_NONZERO_PATTERN);
            double norm;
            CHKERR MatNorm(S_from_block, NORM_FROBENIUS, &norm);
            MOFEM_LOG("AT", Sev::inform) << "Norm of difference: " << norm;
            if (norm > 1e-6)
              SETERRQ(
                  PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Norm of difference between Schur and block matrix is larger "
                  "than accepted");
          }
        } else {
          CHKERR assembleBlockMatSchur(mField, post_proc_schur_lhs_ptr->B, S,
                               {"BROKEN", "U"}, {nullptr, nullptr}, ao_up);
          CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
          CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
        }
 
        MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble Finish";
        MoFEMFunctionReturn(0);
      };
 
      auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
      ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
      ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&](auto pc, auto block_dm) {
      MoFEMFunctionBegin;
      auto block_is = getDMSubData(block_dm)->getSmartRowIs();
      CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    auto set_diagonal_pc = [&](auto pc, auto schur_dm) {
      MoFEMFunctionBegin;
 
      if (AT == BLOCK_SCHUR || AT == BLOCK_MAT) {
        auto A = createDMBlockMat(simple->getDM());
        auto P = createDMNestSchurMat(simple->getDM());
        CHKERR PCSetOperators(pc, A, P);
      }
 
      KSP *subksp;
      CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
      auto get_pc = [](auto ksp) {
        PC pc_raw;
        CHKERR KSPGetPC(ksp, &pc_raw);
        return pc_raw;
      };
      CHKERR setSchurA00MatSolvePC(SmartPetscObj<PC>(get_pc(subksp[0]), true));
 
      auto set_pc_p_mg = [&](auto dm, auto pc) {
        MoFEMFunctionBegin;
 
        CHKERR PCSetDM(pc, dm);
        PetscBool same = PETSC_FALSE;
        PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
        if (same) {
          // By default do not use shell mg mat. Implementation of SOR is slow.
          CHKERR PCMGSetUpViaApproxOrders(
              pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, false));
          CHKERR PCSetFromOptions(pc);
        }
        MoFEMFunctionReturn(0);
      };
 
      CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]));
 
      CHKERR PetscFree(subksp);
      MoFEMFunctionReturn(0);
    };
 
    auto [schur_dm, block_dm] = create_sub_dm();
    if (AT == BLOCK_SCHUR || AT == BLOCK_MAT) {
      auto nested_mat_data = get_nested_mat_data(schur_dm, block_dm);
      CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
    }
    S = createDMHybridisedL2Matrix(schur_dm);
    CHKERR MatSetDM(S, PETSC_NULL);
 
    int bs = (SPACE_DIM == 2) ? NBEDGE_L2(approx_order - 1)
                              : NBFACETRI_L2(approx_order - 1);
    CHKERR MatSetBlockSize(S, bs);
 
    CHKERR set_ops(schur_dm);
    CHKERR set_pc(pc, block_dm);
    DM solver_dm;
    CHKERR KSPGetDM(ksp, &solver_dm);
    if (AT == BLOCK_SCHUR || AT == BLOCK_MAT)
      CHKERR DMSetMatType(solver_dm, MATSHELL);
 
    CHKERR KSPSetUp(ksp);
    if (AT == BLOCK_SCHUR || AT == BLOCK_MAT)
      CHKERR set_diagonal_pc(pc, schur_dm);
 
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "PC is not set to PCFIELDSPLIT");
  }
  MoFEMFunctionReturn(0);
}

setUp() [3/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver )

virtual

Implements SetUpSchur.

setUp() [4/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver )

virtual

Implements SetUpSchur.

Definition at line 1007 of file contact.cpp.

                                                             {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
 
  SNES snes;
  CHKERR TSGetSNES(solver, &snes);
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  CHKERR KSPSetFromOptions(ksp);
 
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
 
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
 
    MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
 
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "It is expected that Schur matrix is not allocated. This is "
          "possible only if PC is set up twice");
    }
 
    CHKERR createSubDM();
 
    // Add data to DM storage
    S = createDMMatrix(schurDM);
    CHKERR MatSetBlockSize(S, SPACE_DIM);
    // CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
 
    // Set DM to use shell block matrix
    DM solver_dm;
    CHKERR TSGetDM(solver, &solver_dm);
    CHKERR DMSetMatType(solver_dm, MATSHELL);
 
    auto ts_ctx_ptr = getDMTsCtx(solver_dm);
    auto A = createDMBlockMat(simple->getDM());
    auto P = createDMNestSchurMat(simple->getDM());
 
    if (is_quasi_static == PETSC_TRUE) {
      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);
      };
      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) {
        return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
      };
      CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
    }
    CHKERR KSPSetOperators(ksp, A, P);
 
    CHKERR setOperator();
    CHKERR setPC(pc);
    CHKERR TSSetUp(solver);
    CHKERR KSPSetUp(ksp);
    CHKERR setDiagonalPC(pc);
 
  } else {
    MOFEM_LOG("CONTACT", Sev::inform) << "No Schur PC";
    pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
    pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
  }
  MoFEMFunctionReturn(0);
}

setUp() [5/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( TS  solver )

virtual

Implements SetUpSchur.

Definition at line 1665 of file plastic.cpp.

                                              {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
 
  SNES snes;
  CHKERR TSGetSNES(solver, &snes);
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  CHKERR KSPSetFromOptions(ksp);
 
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
    if (S) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Is expected that schur matrix is not "
                        "allocated. This is "
                        "possible only is PC is set up twice");
    }
 
    S = createDMMatrix(subDM);
    CHKERR MatSetBlockSize(S, SPACE_DIM);
 
    // Set DM to use shell block matrix
    DM solver_dm;
    CHKERR TSGetDM(solver, &solver_dm);
    CHKERR DMSetMatType(solver_dm, MATSHELL);
 
    auto ts_ctx_ptr = getDMTsCtx(solver_dm);
    auto A = createDMBlockMat(simple->getDM());
    auto P = createDMNestSchurMat(simple->getDM());
 
    if (is_quasi_static == PETSC_TRUE) {
      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);
      };
      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) {
        return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
      };
      CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
    }
    CHKERR KSPSetOperators(ksp, A, P);
 
    auto set_ops = [&]() {
      MoFEMFunctionBegin;
      auto pip_mng = mField.getInterface<PipelineManager>();
 
#ifndef ADD_CONTACT
      // Boundary
      pip_mng->getOpBoundaryLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
 
          ));
      // Domain
      pip_mng->getOpDomainLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
 
          ));
#else
 
      double eps_stab = 1e-4;
      CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
                                   PETSC_NULL);
 
      using B = FormsIntegrators<BoundaryEleOp>::Assembly<
          BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
      using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
 
      // Boundary
      pip_mng->getOpBoundaryLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpBoundaryLhsPipeline().push_back(
          new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
            return eps_stab;
          }));
      pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
          false, false
 
          ));
      // Domain
      pip_mng->getOpDomainLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
          false, false
 
          ));
#endif // ADD_CONTACT
      MoFEMFunctionReturn(0);
    };
 
    auto set_assemble_elems = [&]() {
      MoFEMFunctionBegin;
      auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
      schur_asmb_pre_proc->preProcessHook = [this]() {
        MoFEMFunctionBegin;
        CHKERR MatZeroEntries(S);
        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
        MoFEMFunctionReturn(0);
      };
      auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
 
      schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
        MoFEMFunctionBegin;
        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
 
        // Apply essential constrains to Schur complement
        CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
            mField, schur_asmb_post_proc, 1, S, aoSchur)();
 
        MoFEMFunctionReturn(0);
      };
      auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
      ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
      ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&]() {
      MoFEMFunctionBegin;
      CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    auto set_diagonal_pc = [&]() {
      MoFEMFunctionBegin;
      KSP *subksp;
      CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &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]));
      CHKERR PetscFree(subksp);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_ops();
    CHKERR set_pc();
    CHKERR set_assemble_elems();
 
    CHKERR TSSetUp(solver);
    CHKERR KSPSetUp(ksp);
    CHKERR set_diagonal_pc();
 
  } else {
    pip_mng->getOpBoundaryLhsPipeline().push_front(
        createOpSchurAssembleBegin());
    pip_mng->getOpBoundaryLhsPipeline().push_back(
        createOpSchurAssembleEnd({}, {}));
    pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip_mng->getOpDomainLhsPipeline().push_back(
        createOpSchurAssembleEnd({}, {}));
  }
 
  // fieldSplitIS.reset();
  // aoSchur.reset();
  MoFEMFunctionReturn(0);
}

setUp() [6/6]

MoFEMErrorCode SetUpSchurImpl::setUp ( TS  ts )

virtual

Implements SetUpSchur.

Member Data Documentation

aoS

SmartPetscObj<AO> SetUpSchurImpl::aoS

private

Definition at line 24 of file SetUpSchurImpl.cpp.

aoSchur

SmartPetscObj<AO> SetUpSchurImpl::aoSchur

private

Examples

plastic.cpp.

Definition at line 1662 of file plastic.cpp.

blockDM

SmartPetscObj< DM > SetUpSchurImpl::blockDM

private

Definition at line 1004 of file contact.cpp.

epCorePtr

EshelbianCore* SetUpSchurImpl::epCorePtr

private

Definition at line 21 of file SetUpSchurImpl.cpp.

fieldSplitIS

SmartPetscObj<IS> SetUpSchurImpl::fieldSplitIS

private

IS for split Schur block.

Examples

plastic.cpp.

Definition at line 1661 of file plastic.cpp.

mField

MoFEM::Interface & SetUpSchurImpl::mField

private

Examples

plastic.cpp, and test_broken_space.cpp.

Definition at line 528 of file test_broken_space.cpp.

piolaZeroDofsMarker

boost::shared_ptr<std::vector<unsigned char> > SetUpSchurImpl::piolaZeroDofsMarker

private

Definition at line 29 of file SetUpSchurImpl.cpp.

piolaZeroDofsVec

boost::shared_ptr<std::vector<boost::weak_ptr<NumeredDofEntity> > > SetUpSchurImpl::piolaZeroDofsVec

private

Definition at line 27 of file SetUpSchurImpl.cpp.

pMGPtr

boost::shared_ptr<P_MultiGridData> SetUpSchurImpl::pMGPtr

private

Definition at line 51 of file SetUpSchurImpl.cpp.

S

SmartPetscObj< Mat > SetUpSchurImpl::S

private

Examples

plastic.cpp, and test_broken_space.cpp.

Definition at line 529 of file test_broken_space.cpp.

schurDM

SmartPetscObj< DM > SetUpSchurImpl::schurDM

private

Definition at line 1003 of file contact.cpp.

subDM

SmartPetscObj< DM > SetUpSchurImpl::subDM

private

field split sub dm

Examples

plastic.cpp.

Definition at line 1660 of file plastic.cpp.

subEnts

Range SetUpSchurImpl::subEnts

private

Definition at line 950 of file elastic.cpp.

volEnts

Range SetUpSchurImpl::volEnts

private

Definition at line 949 of file elastic.cpp.

---

The documentation for this struct was generated from the following files:

• atom_tests/test_broken_space.cpp
• tutorials/adv-0/plastic.cpp
• tutorials/adv-1/contact.cpp
• tutorials/mix-2/incompressible_elasticity.cpp
• tutorials/vec-0/elastic.cpp
• users_modules/eshelbian_plasticit/src/impl/SetUpSchurImpl.cpp