24 SmartPetscObj<AO>
aoS;
26 boost::shared_ptr<std::vector<boost::weak_ptr<NumeredDofEntity>>>
28 boost::shared_ptr<std::vector<unsigned char>>
39 CHKERR PCMGSetUpViaApproxOrders(
40 pc, createPCMGSetUpViaApproxOrdersCtx(
dm,
S,
true));
51 boost::shared_ptr<P_MultiGridData>
pMGPtr;
56 std::vector<boost::shared_ptr<Range>> dm_range_list{
nullptr,
nullptr};
57 return std::make_pair(schur_field_list, dm_range_list);
61 std::vector<std::string> a00_field_list{
74 std::vector<boost::shared_ptr<Range>> range_list_ptr(a00_field_list.size(),
76 return std::make_pair(a00_field_list, range_list_ptr);
83 auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
95 for (
auto f : schur_field_list) {
96 MOFEM_LOG(
"EP", Sev::inform) <<
"Add schur field: " <<
f;
105 auto create_a00_dm = [&](SmartPetscObj<DM> &dm_sub) {
117 for (
auto f : a00_field_list) {
118 MOFEM_LOG(
"EP", Sev::inform) <<
"Add a00 field: " <<
f;
127 auto get_snes = [&](TS ts) {
129 CHKERR TSGetSNES(ts, &snes);
133 auto get_ksp = [&](SNES snes) {
135 CHKERR SNESGetKSP(snes, &ksp);
136 CHKERR KSPSetFromOptions(ksp);
140 auto get_pc = [&](KSP ksp) {
142 CHKERR KSPGetPC(ksp, &pc);
146 auto ksp = get_ksp(get_snes(ts));
147 auto pc = get_pc(ksp);
149 PetscBool is_pcfs = PETSC_FALSE;
150 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
153 MOFEM_LOG(
"EP", Sev::inform) <<
"SetUpSchurImpl::setUp: PCFIELDSPLIT";
155 SmartPetscObj<DM> schur_dm, a00_dm;
156 CHKERR create_schur_dm(schur_dm);
157 CHKERR create_a00_dm(a00_dm);
165 std::vector<std::pair<std::string, std::string>> mat_block_list = {
185 mat_block_list.push_back(
187 mat_block_list.push_back(
189 mat_block_list.push_back(
200 auto get_nested_mat_data = [&](
auto schur_dm,
auto block_dm) {
206 {vol_elem_name, mat_block_list},
227 {natural_bc_element_name,
245 {schur_dm, a00_dm}, block_mat_data,
256 auto nested_mat_data = get_nested_mat_data(schur_dm, a00_dm);
264 std::numeric_limits<double>::epsilon()) {
265 auto swap_assemble = [](TS ts, PetscReal
t, Vec u, Vec u_t, Vec utt,
266 PetscReal
a, PetscReal aa, Mat
A, Mat
B,
271 CHKERR TSSetI2Jacobian(ts,
m, p, swap_assemble, ts_ctx_ptr.get());
273 auto swap_assemble = [](TS ts, PetscReal
t, Vec u, Vec u_t, PetscReal
a,
274 Mat
A, Mat
B,
void *ctx) {
278 CHKERR TSSetIJacobian(ts,
m, p, swap_assemble, ts_ctx_ptr.get());
280 CHKERR KSPSetOperators(ksp,
m, p);
282 auto set_assembly = [&]() {
292 auto set_assemble = [&]() {
294 auto schur_asmb_pre_proc_lhs = boost::make_shared<FEMethod>();
295 auto schur_asmb_pre_proc_rhs = boost::make_shared<FEMethod>();
297 schur_asmb_pre_proc_lhs->preProcessHook = [
this]() {
303 schur_asmb_pre_proc_rhs->preProcessHook = [
this,
304 schur_asmb_pre_proc_rhs]() {
306 auto prb_ptr = schur_asmb_pre_proc_rhs->problemPtr;
307 auto dofs_prb = prb_ptr->getNumeredRowDofsPtr();
312 ->getSideDofsOnBrokenSpaceEntities(
320 if (
auto dof_ptr = dof.lock()) {
321 auto idx = dof_ptr->getPetscLocalDofIdx();
322 (*piolaZeroDofsMarker)[idx] = 1;
326 constexpr bool hard_coded_set_bc_debug =
false;
327 if constexpr (hard_coded_set_bc_debug) {
329 auto problem_name = schur_asmb_pre_proc_rhs->problemPtr->getName();
332 SmartPetscObj<IS> crack_hybrid_is;
334 ->isCreateProblemFieldAndRankLocal(
336 SPACE_DIM, crack_hybrid_is, &*crack_faces);
338 SmartPetscObj<IS> crack_piola_is;
344 CHKERR VecGetArrayRead(schur_asmb_pre_proc_rhs->x, &a_x);
345 auto zero_by_is = [&](
auto is) {
347 const PetscInt *is_array;
349 CHKERR ISGetLocalSize(is, &is_size);
350 CHKERR ISGetIndices(is, &is_array);
351 for (
int i = 0;
i != is_size; ++
i) {
353 const_cast<double *
>(a_x)[is_array[
i]] = 0;
355 CHKERR ISRestoreIndices(is, &is_array);
359 CHKERR zero_by_is(crack_hybrid_is);
360 CHKERR zero_by_is(crack_piola_is);
362 CHKERR VecRestoreArrayRead(schur_asmb_pre_proc_rhs->x, &a_x);
365 ->setLocalGhostVector(problem_name,
COL,
366 schur_asmb_pre_proc_rhs->x, INSERT_VALUES,
373 auto schur_asmb_post_proc_lhs = boost::make_shared<FEMethod>();
374 auto schur_asmb_post_proc_rhs = boost::make_shared<FEMethod>();
376 schur_asmb_post_proc_rhs->postProcessHook =
377 [
this, schur_asmb_post_proc_rhs]() {
380 CHKERR VecGhostUpdateBegin(schur_asmb_post_proc_rhs->f,
381 ADD_VALUES, SCATTER_REVERSE);
382 CHKERR VecGhostUpdateEnd(schur_asmb_post_proc_rhs->f, ADD_VALUES,
384 CHKERR VecAssemblyBegin(schur_asmb_post_proc_rhs->f);
385 CHKERR VecAssemblyEnd(schur_asmb_post_proc_rhs->f);
386 *(schur_asmb_post_proc_rhs->vecAssembleSwitch) =
false;
391 schur_asmb_post_proc_rhs->problemPtr->getName();
395 SmartPetscObj<IS> crack_hybrid_is;
397 ->isCreateProblemFieldAndRankLocal(
399 SPACE_DIM, crack_hybrid_is, &*crack_faces);
401 SmartPetscObj<IS> crack_piola_is;
407 CHKERR VecGetArray(schur_asmb_post_proc_rhs->f, &a_f);
409 CHKERR VecGetArrayRead(schur_asmb_post_proc_rhs->x, &a_x);
410 auto zero_by_is = [&](
auto is) {
412 const PetscInt *is_array;
414 CHKERR ISGetLocalSize(is, &is_size);
415 CHKERR ISGetIndices(is, &is_array);
416 for (
int i = 0;
i != is_size; ++
i) {
417 a_f[is_array[
i]] = -a_x[is_array[
i]];
419 CHKERR ISRestoreIndices(is, &is_array);
423 CHKERR zero_by_is(crack_hybrid_is);
424 CHKERR zero_by_is(crack_piola_is);
426 CHKERR VecRestoreArray(schur_asmb_post_proc_rhs->f, &a_f);
427 CHKERR VecRestoreArrayRead(schur_asmb_post_proc_rhs->x, &a_x);
433 schur_asmb_post_proc_lhs->postProcessHook =
434 [
this, schur_asmb_post_proc_lhs]() {
445 CHKERR MatAssemblyBegin(schur_asmb_post_proc_lhs->B,
447 CHKERR MatAssemblyEnd(schur_asmb_post_proc_lhs->B,
449 *(schur_asmb_post_proc_lhs->matAssembleSwitch) =
false;
451 SmartPetscObj<IS> crack_hybrid_is;
453 ->isCreateProblemFieldAndRank(
455 SPACE_DIM, crack_hybrid_is, &*crack_faces);
456 CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
457 crack_hybrid_is, 1, PETSC_NULLPTR,
461 SmartPetscObj<IS> crack_piola_is;
465 CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
466 crack_piola_is, 1, PETSC_NULLPTR,
472 schur_asmb_post_proc_lhs->B,
S,
473 a00_field_list, a00_range_list,
aoS);
478 CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
479 CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
481 SmartPetscObj<IS> crack_hybrid_is;
483 ->isCreateProblemFieldAndRank(
485 SPACE_DIM, crack_hybrid_is, &*crack_faces);
487 CHKERR MatZeroRowsColumnsIS(
S, crack_hybrid_is, 1, PETSC_NULLPTR,
494 ts_ctx_ptr->getPreProcessIFunction().push_front(
495 schur_asmb_pre_proc_rhs);
496 ts_ctx_ptr->getPostProcessIFunction().push_back(
497 schur_asmb_post_proc_rhs);
498 ts_ctx_ptr->getPreProcessIJacobian().push_front(
499 schur_asmb_pre_proc_lhs);
500 ts_ctx_ptr->getPostProcessIJacobian().push_back(
501 schur_asmb_post_proc_lhs);
506 boost::make_shared<std::vector<boost::weak_ptr<NumeredDofEntity>>>();
513 auto set_pc = [&]() {
516 auto schur_is =
getDMSubData(schur_dm)->getSmartRowIs();
517 CHKERR PCFieldSplitSetIS(pc, NULL, a00_is);
518 CHKERR PCFieldSplitSetIS(pc, NULL, schur_is);
519 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
523 auto set_diagonal_pc = [&]() {
526 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
527 auto get_pc = [](
auto ksp) {
529 CHKERR KSPGetPC(ksp, &pc_raw);
530 return SmartPetscObj<PC>(pc_raw,
true);
534 auto set_pc_p_mg = [&](
auto dm,
auto pc,
auto S) {
537 PetscBool same = PETSC_FALSE;
538 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
540 auto smart_pc = SmartPetscObj<PC>(pc,
true);
541 pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc,
S);
543 PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
546 <<
"SetUpSchurImpl::setUp: fieldsplit 1 PCKSP";
547 CHKERR PCSetFromOptions(pc);
549 CHKERR PCKSPGetKSP(pc, &ksp);
550 CHKERR KSPSetFromOptions(ksp);
552 CHKERR KSPGetPC(ksp, &ksp_pc);
553 CHKERR PCSetFromOptions(ksp_pc);
554 PetscObjectTypeCompare((PetscObject)ksp_pc, PCMG, &same);
556 auto smart_pc = SmartPetscObj<PC>(ksp_pc,
true);
557 pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc,
S);
563 CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]),
S);
576 MOFEM_LOG(
"EP", Sev::inform) <<
"SetUpSchurImpl::setUp: PCLU or other";
591boost::shared_ptr<EshelbianCore::SetUpSchur>
594 return boost::shared_ptr<SetUpSchur>(
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
#define MOFEM_LOG(channel, severity)
Log.
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
auto createDMHybridisedL2Matrix(DM dm)
Get smart hybridised L2 matrix from DM.
auto getDMSubData(DM dm)
Get sub problem data structure.
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
auto createDMNestSchurMat(DM dm)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao)
Assemble Schur matrix.
OpSchurAssembleBase * createOpSchurAssembleBegin()
auto createDMBlockMat(DM dm)
constexpr double t
plate stiffness
FTensor::Index< 'm', 3 > m
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
std::vector< boost::shared_ptr< Range > > a00RangeList
const std::string skeletonElement
MoFEM::Interface & mField
const std::string spatialL2Disp
SmartPetscObj< IS > crackHybridIs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
const std::string elementVolumeName
const std::string piolaStress
std::vector< std::string > a00FieldList
const std::string bubbleField
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
static PetscBool noStretch
const std::string rotAxis
const std::string contactDisp
const std::string naturalBcElement
boost::shared_ptr< Range > crackFaces
const std::string hybridSpatialDisp
SmartPetscObj< DM > dmElastic
Elastic problem.
const std::string stretchTensor
const std::string contactElement
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
P_MultiGridData(SmartPetscObj< DM > dm, SmartPetscObj< PC > pc, SmartPetscObj< Mat > S)
EshelbianCore * epCorePtr
boost::shared_ptr< std::vector< unsigned char > > piolaZeroDofsMarker
boost::shared_ptr< std::vector< boost::weak_ptr< NumeredDofEntity > > > piolaZeroDofsVec
boost::shared_ptr< P_MultiGridData > pMGPtr
MoFEMErrorCode setUp(TS ts)
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
virtual ~SetUpSchurImpl()
MoFEMErrorCode postProc()
MoFEM::Interface & mField
SetUpSchurImpl(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
[Push operators to pipeline]