25 boost::shared_ptr<std::vector<boost::weak_ptr<NumeredDofEntity>>>
27 boost::shared_ptr<std::vector<unsigned char>>
34 auto get_schur_fields = [&]() {
35 std::vector<std::string> schur_field_list{
epCorePtr->hybridSpatialDisp,
37 std::vector<boost::shared_ptr<Range>> dm_range_list{
nullptr,
nullptr};
38 return std::make_pair(schur_field_list, dm_range_list);
41 auto get_a00_fields = [&]() {
42 std::vector<std::string> a00_field_list{
55 std::vector<boost::shared_ptr<Range>> range_list_ptr(a00_field_list.size(),
57 return std::make_pair(a00_field_list, range_list_ptr);
60 auto schur_data = get_schur_fields();
61 auto [schur_field_list, schur_range_list] = schur_data;
62 auto a00_data = get_a00_fields();
63 auto [a00_field_list, a00_range_list] = a00_data;
65 auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
76 auto [schur_field_list, schur_range_list] = schur_data;
77 for (
auto f : schur_field_list) {
78 MOFEM_LOG(
"EP", Sev::inform) <<
"Add schur field: " <<
f;
87 auto create_a00_dm = [&](SmartPetscObj<DM> &dm_sub) {
98 auto [a00_field_list, a00_range_list] = a00_data;
99 for (
auto f : a00_field_list) {
100 MOFEM_LOG(
"EP", Sev::inform) <<
"Add a00 field: " <<
f;
109 auto get_snes = [&](TS ts) {
111 CHKERR TSGetSNES(ts, &snes);
115 auto get_ksp = [&](SNES snes) {
117 CHKERR SNESGetKSP(snes, &ksp);
118 CHKERR KSPSetFromOptions(ksp);
122 auto get_pc = [&](KSP ksp) {
124 CHKERR KSPGetPC(ksp, &pc);
128 auto ksp = get_ksp(get_snes(ts));
129 auto pc = get_pc(ksp);
131 PetscBool is_pcfs = PETSC_FALSE;
132 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
135 MOFEM_LOG(
"EP", Sev::inform) <<
"SetUpSchurImpl::setUp: PCFIELDSPLIT";
137 SmartPetscObj<DM> schur_dm, a00_dm;
138 CHKERR create_schur_dm(schur_dm);
139 CHKERR create_a00_dm(a00_dm);
142 auto vol_elem_name =
epCorePtr->elementVolumeName;
143 auto skel_elem_name =
epCorePtr->skeletonElement;
144 auto contact_elem_name =
epCorePtr->contactElement;
146 std::vector<std::pair<std::string, std::string>> mat_block_list = {
166 mat_block_list.push_back(
168 mat_block_list.push_back(
170 mat_block_list.push_back(
181 auto get_nested_mat_data = [&](
auto schur_dm,
auto block_dm) {
187 {vol_elem_name, mat_block_list},
212 auto [a00_field_list, a00_range_list] = a00_data;
216 {schur_dm, a00_dm}, block_mat_data,
227 auto nested_mat_data = get_nested_mat_data(schur_dm, a00_dm);
235 std::numeric_limits<double>::epsilon()) {
236 auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t,
Vec utt,
237 PetscReal
a, PetscReal aa, Mat
A, Mat B,
242 CHKERR TSSetI2Jacobian(ts,
m, p, swap_assemble, ts_ctx_ptr.get());
244 auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t, PetscReal
a,
245 Mat
A, Mat B,
void *ctx) {
249 CHKERR TSSetIJacobian(ts,
m, p, swap_assemble, ts_ctx_ptr.get());
251 CHKERR KSPSetOperators(ksp,
m, p);
253 auto set_assembly = [&]() {
255 auto [a00_field_list, a00_range_list] = a00_data;
256 std::vector<boost::shared_ptr<Range>> ranges_list(a00_field_list.size(),
259 auto schur_ao =
getDMSubData(schur_dm)->getSmartRowMap();
264 auto set_ops = [&]() {
267 auto [a00_field_list, a00_range_list] = a00_data;
268 const bool symmetric_schur =
false;
269 const bool symmetric_block =
false;
270 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
272 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
274 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
276 S, symmetric_schur, symmetric_block));
277 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
279 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
281 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
283 symmetric_schur, symmetric_block));
287 auto set_assemble = [&]() {
289 auto schur_asmb_pre_proc_lhs = boost::make_shared<FEMethod>();
290 auto schur_asmb_pre_proc_rhs = boost::make_shared<FEMethod>();
292 schur_asmb_pre_proc_lhs->preProcessHook = [
this]() {
298 schur_asmb_pre_proc_rhs->preProcessHook =
299 [
this, schur_asmb_pre_proc_rhs]() {
301 auto prb_ptr = schur_asmb_pre_proc_rhs->problemPtr;
302 auto dofs_prb = prb_ptr->getNumeredRowDofsPtr();
304 auto crack_faces =
epCorePtr->crackFaces;
307 ->getSideDofsOnBrokenSpaceEntities(
315 if (
auto dof_ptr = dof.lock()) {
316 auto idx = dof_ptr->getPetscLocalDofIdx();
317 (*piolaZeroDofsMarker)[idx] = 1;
324 auto schur_asmb_post_proc_lhs = boost::make_shared<FEMethod>();
325 auto schur_asmb_post_proc_rhs = boost::make_shared<FEMethod>();
327 schur_asmb_post_proc_rhs->postProcessHook =
328 [
this, schur_asmb_post_proc_rhs]() {
331 CHKERR VecGhostUpdateBegin(schur_asmb_post_proc_rhs->f,
332 ADD_VALUES, SCATTER_REVERSE);
333 CHKERR VecGhostUpdateEnd(schur_asmb_post_proc_rhs->f, ADD_VALUES,
335 CHKERR VecAssemblyBegin(schur_asmb_post_proc_rhs->f);
336 CHKERR VecAssemblyEnd(schur_asmb_post_proc_rhs->f);
337 *(schur_asmb_post_proc_rhs->vecAssembleSwitch) =
false;
342 schur_asmb_post_proc_rhs->problemPtr->getName();
344 auto crack_faces =
epCorePtr->crackFaces;
346 SmartPetscObj<IS> crack_hybrid_is;
348 ->isCreateProblemFieldAndRank(
350 SPACE_DIM, crack_hybrid_is, &*crack_faces);
352 SmartPetscObj<IS> crack_piola_is;
358 CHKERR VecGetArray(schur_asmb_post_proc_rhs->f, &a_f);
359 auto zero_by_is = [&](
auto is) {
361 const PetscInt *is_array;
363 CHKERR ISGetLocalSize(is, &is_size);
364 CHKERR ISGetIndices(is, &is_array);
365 for (
int i = 0;
i != is_size; ++
i) {
366 a_f[is_array[
i]] = 0;
368 CHKERR ISRestoreIndices(is, &is_array);
372 CHKERR zero_by_is(crack_hybrid_is);
373 CHKERR zero_by_is(crack_piola_is);
375 CHKERR VecRestoreArray(schur_asmb_post_proc_rhs->f, &a_f);
382 schur_asmb_post_proc_lhs->postProcessHook =
383 [
this, schur_asmb_post_proc_lhs]() {
387 CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
388 CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
391 auto crack_faces =
epCorePtr->crackFaces;
393 SmartPetscObj<IS> crack_hybrid_is;
395 ->isCreateProblemFieldAndRank(
397 SPACE_DIM, crack_hybrid_is, &*crack_faces);
398 CHKERR MatZeroRowsColumnsIS(
S, crack_hybrid_is, 1, PETSC_NULL,
403 CHKERR MatAssemblyBegin(schur_asmb_post_proc_lhs->B,
405 CHKERR MatAssemblyEnd(schur_asmb_post_proc_lhs->B,
407 *(schur_asmb_post_proc_lhs->matAssembleSwitch) =
false;
409 SmartPetscObj<IS> crack_hybrid_is;
411 ->isCreateProblemFieldAndRank(
412 "ELASTIC_PROBLEM",
ROW,
epCorePtr->hybridSpatialDisp, 0,
413 SPACE_DIM, crack_hybrid_is, &*crack_faces);
414 CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
415 crack_hybrid_is, 1, PETSC_NULL,
419 SmartPetscObj<IS> crack_piola_is;
423 CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
424 crack_piola_is, 1, PETSC_NULL,
432 ts_ctx_ptr->getPreProcessIFunction().push_front(
433 schur_asmb_pre_proc_rhs);
434 ts_ctx_ptr->getPostProcessIFunction().push_back(
435 schur_asmb_post_proc_rhs);
436 ts_ctx_ptr->getPreProcessIJacobian().push_front(
437 schur_asmb_pre_proc_lhs);
438 ts_ctx_ptr->getPostProcessIJacobian().push_back(
439 schur_asmb_post_proc_lhs);
444 boost::make_shared<std::vector<boost::weak_ptr<NumeredDofEntity>>>();
452 auto set_pc = [&]() {
455 auto schur_is =
getDMSubData(schur_dm)->getSmartRowIs();
456 CHKERR PCFieldSplitSetIS(pc, NULL, a00_is);
457 CHKERR PCFieldSplitSetIS(pc, NULL, schur_is);
458 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
462 auto set_diagonal_pc = [&]() {
465 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
466 auto get_pc = [](
auto ksp) {
468 CHKERR KSPGetPC(ksp, &pc_raw);
469 return SmartPetscObj<PC>(pc_raw,
true);
473 auto set_pc_p_mg = [](
auto dm,
auto pc,
auto S) {
476 PetscBool same = PETSC_FALSE;
477 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
481 CHKERR PCSetFromOptions(pc);
486 CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]),
S);
499 MOFEM_LOG(
"EP", Sev::inform) <<
"SetUpSchurImpl::setUp: PCLU or other";
501 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
503 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
505 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
507 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
514 boost::shared_ptr<EshelbianCore::SetUpSchur>
517 return boost::shared_ptr<SetUpSchur>(