v0.14.0
SetUpSchurImpl.cpp
Go to the documentation of this file.
1 /** @file SetUpSchurImpl
2  * @brief
3  * @date 2023-05-13
4  *
5  * @license{This project is released under the MIT License.}
6  *
7  */
8 
9 struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
10 
11  SetUpSchurImpl(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
12  : SetUpSchur(), mField(m_field), epCorePtr(ep_core_ptr) {}
13  virtual ~SetUpSchurImpl() {}
14 
15  MoFEMErrorCode setUp(TS ts);
16  MoFEMErrorCode preProc();
17  MoFEMErrorCode postProc();
18 
19 private:
20  MoFEM::Interface &mField;
21  EshelbianCore *epCorePtr;
22 
23  SmartPetscObj<Mat> S;
24  SmartPetscObj<AO> aoS;
25 
26  boost::shared_ptr<std::vector<boost::weak_ptr<NumeredDofEntity>>>
27  piolaZeroDofsVec; //< Dofs on crack surface
28  boost::shared_ptr<std::vector<unsigned char>>
29  piolaZeroDofsMarker; //< marker for crack dofs on surface
30 
31  struct P_MultiGridData {
32 
33  P_MultiGridData(SmartPetscObj<DM> dm, SmartPetscObj<PC> pc,
34  SmartPetscObj<Mat> S)
35  : dm(dm), pc(pc), S(S) {}
36 
37  MoFEMErrorCode setUP() {
38  MoFEMFunctionBegin;
39  CHKERR PCMGSetUpViaApproxOrders(
40  pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
41  CHKERR PCSetFromOptions(pc);
42  MoFEMFunctionReturn(0);
43  };
44 
45  private:
46  SmartPetscObj<DM> dm;
47  SmartPetscObj<PC> pc;
48  SmartPetscObj<Mat> S;
49  };
50 
51  boost::shared_ptr<P_MultiGridData> pMGPtr;
52 
53  auto getSchurFields() {
54  std::vector<std::string> schur_field_list{epCorePtr->hybridSpatialDisp,
55  epCorePtr->contactDisp};
56  std::vector<boost::shared_ptr<Range>> dm_range_list{nullptr, nullptr};
57  return std::make_pair(schur_field_list, dm_range_list);
58  };
59 
60  auto getA00Fields() {
61  std::vector<std::string> a00_field_list{
62 
63  epCorePtr->stretchTensor,
64 
65  epCorePtr->bubbleField,
66 
67  epCorePtr->piolaStress,
68 
69  epCorePtr->rotAxis,
70 
71  epCorePtr->spatialL2Disp
72 
73  };
74  std::vector<boost::shared_ptr<Range>> range_list_ptr(a00_field_list.size(),
75  nullptr);
76  return std::make_pair(a00_field_list, range_list_ptr);
77  }
78 };
79 
80 MoFEMErrorCode SetUpSchurImpl::setUp(TS ts) {
81  MoFEMFunctionBegin;
82 
83  auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
84  MoFEMFunctionBegin;
85 
86  dm_sub = createDM(mField.get_comm(), "DMMOFEM_MG");
87  CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
88  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
89  CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
90  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->skeletonElement);
91  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
92 
93  int r_idx = 0;
94  auto [schur_field_list, schur_range_list] = getSchurFields();
95  for (auto f : schur_field_list) {
96  MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
97  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, schur_range_list[r_idx]);
98  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, schur_range_list[r_idx]);
99  ++r_idx;
100  }
101  CHKERR DMSetUp(dm_sub);
102  MoFEMFunctionReturn(0);
103  };
104 
105  auto create_a00_dm = [&](SmartPetscObj<DM> &dm_sub) {
106  MoFEMFunctionBegin;
107  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
108  CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_A00");
109  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
110  CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
111  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
112  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->skeletonElement);
113  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
114 
115  int r_idx = 0;
116  auto [a00_field_list, a00_range_list] = getA00Fields();
117  for (auto f : a00_field_list) {
118  MOFEM_LOG("EP", Sev::inform) << "Add a00 field: " << f;
119  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, a00_range_list[r_idx]);
120  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, a00_range_list[r_idx]);
121  ++r_idx;
122  }
123  CHKERR DMSetUp(dm_sub);
124  MoFEMFunctionReturn(0);
125  };
126 
127  auto get_snes = [&](TS ts) {
128  SNES snes;
129  CHKERR TSGetSNES(ts, &snes);
130  return snes;
131  };
132 
133  auto get_ksp = [&](SNES snes) {
134  KSP ksp;
135  CHKERR SNESGetKSP(snes, &ksp);
136  CHKERR KSPSetFromOptions(ksp);
137  return ksp;
138  };
139 
140  auto get_pc = [&](KSP ksp) {
141  PC pc;
142  CHKERR KSPGetPC(ksp, &pc);
143  return pc;
144  };
145 
146  auto ksp = get_ksp(get_snes(ts));
147  auto pc = get_pc(ksp);
148 
149  PetscBool is_pcfs = PETSC_FALSE;
150  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
151  if (is_pcfs) {
152 
153  MOFEM_LOG("EP", Sev::inform) << "SetUpSchurImpl::setUp: PCFIELDSPLIT";
154 
155  SmartPetscObj<DM> schur_dm, a00_dm;
156  CHKERR create_schur_dm(schur_dm);
157  CHKERR create_a00_dm(a00_dm);
158 
159  auto dm_elastic = epCorePtr->dmElastic;
160  auto vol_elem_name = epCorePtr->elementVolumeName;
161  auto skel_elem_name = epCorePtr->skeletonElement;
162  auto contact_elem_name = epCorePtr->contactElement;
163 
164  std::vector<std::pair<std::string, std::string>> mat_block_list = {
165 
166  {epCorePtr->piolaStress, epCorePtr->piolaStress},
167  {epCorePtr->stretchTensor, epCorePtr->stretchTensor},
168  {epCorePtr->piolaStress, epCorePtr->stretchTensor},
169  {epCorePtr->stretchTensor, epCorePtr->piolaStress},
170  {epCorePtr->bubbleField, epCorePtr->stretchTensor},
171  {epCorePtr->stretchTensor, epCorePtr->bubbleField},
172 
173  {epCorePtr->piolaStress, epCorePtr->spatialL2Disp},
174  {epCorePtr->spatialL2Disp, epCorePtr->piolaStress},
175  {epCorePtr->piolaStress, epCorePtr->rotAxis},
176  {epCorePtr->rotAxis, epCorePtr->piolaStress},
177  {epCorePtr->bubbleField, epCorePtr->rotAxis},
178  {epCorePtr->rotAxis, epCorePtr->bubbleField},
179  {epCorePtr->spatialL2Disp, epCorePtr->spatialL2Disp}
180 
181  };
182 
183  if(epCorePtr->noStretch) {
184  mat_block_list.push_back(
185  {epCorePtr->bubbleField, epCorePtr->bubbleField});
186  mat_block_list.push_back(
187  {epCorePtr->bubbleField, epCorePtr->piolaStress});
188  mat_block_list.push_back(
189  {epCorePtr->piolaStress, epCorePtr->bubbleField});
190  }
191 
192  // if (epCorePtr->gradApproximator > EshelbianPlasticity::SMALL_ROT ||
193  // epCorePtr->rotSelector > EshelbianPlasticity::SMALL_ROT) {
194  mat_block_list.push_back({epCorePtr->rotAxis, epCorePtr->rotAxis});
195  mat_block_list.push_back({epCorePtr->stretchTensor, epCorePtr->rotAxis});
196  mat_block_list.push_back({epCorePtr->rotAxis, epCorePtr->stretchTensor});
197  // }
198 
199  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
200  auto block_mat_data = createBlockMatStructure(
201  dm_elastic,
202 
203  {
204 
205  {vol_elem_name, mat_block_list},
206 
207  {skel_elem_name,
208 
209  {
210 
211  {epCorePtr->piolaStress, epCorePtr->hybridSpatialDisp},
212  {epCorePtr->hybridSpatialDisp, epCorePtr->piolaStress}
213 
214  }},
215 
216  {contact_elem_name,
217 
218  {
219 
220  {epCorePtr->contactDisp, epCorePtr->contactDisp},
221  {epCorePtr->piolaStress, epCorePtr->contactDisp},
222  {epCorePtr->contactDisp, epCorePtr->piolaStress}
223 
224  }}
225 
226  }
227 
228  );
229 
230  auto [a00_field_list, a00_range_list] = getA00Fields();
231 
232  return createSchurNestedMatrixStruture(
233 
234  {schur_dm, a00_dm}, block_mat_data,
235 
236  a00_field_list,
237 
238  a00_range_list,
239 
240  false
241 
242  );
243  };
244 
245  auto nested_mat_data = get_nested_mat_data(schur_dm, a00_dm);
246  CHKERR DMMoFEMSetNestSchurData(epCorePtr->dmElastic, nested_mat_data);
247  CHKERR DMSetMatType(epCorePtr->dmElastic, MATSHELL);
248 
249  auto m = createDMBlockMat(epCorePtr->dmElastic);
250  auto p = createDMNestSchurMat(epCorePtr->dmElastic);
251 
252  if (std::abs(epCorePtr->alphaRho) >
253  std::numeric_limits<double>::epsilon()) {
254  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
255  PetscReal a, PetscReal aa, Mat A, Mat B,
256  void *ctx) {
257  return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
258  };
259  auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
260  CHKERR TSSetI2Jacobian(ts, m, p, swap_assemble, ts_ctx_ptr.get());
261  } else {
262  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
263  Mat A, Mat B, void *ctx) {
264  return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
265  };
266  auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
267  CHKERR TSSetIJacobian(ts, m, p, swap_assemble, ts_ctx_ptr.get());
268  }
269  CHKERR KSPSetOperators(ksp, m, p);
270 
271  auto set_assembly = [&]() {
272  MoFEMFunctionBegin;
273 
274  aoS = getDMSubData(schur_dm)->getSmartRowMap();
275  S = createDMHybridisedL2Matrix(schur_dm);
276  CHKERR MatSetBlockSize(S, SPACE_DIM *
277  NBFACETRI_L2(epCorePtr->spaceOrder - 1));
278  epCorePtr->S = S;
279  epCorePtr->aoS = aoS;
280 
281  auto set_assemble = [&]() {
282  MoFEMFunctionBegin;
283  auto schur_asmb_pre_proc_lhs = boost::make_shared<FEMethod>();
284  auto schur_asmb_pre_proc_rhs = boost::make_shared<FEMethod>();
285 
286  schur_asmb_pre_proc_lhs->preProcessHook = [this]() {
287  MoFEMFunctionBegin;
288  CHKERR MatZeroEntries(S);
289  MoFEMFunctionReturn(0);
290  };
291 
292  schur_asmb_pre_proc_rhs->preProcessHook =
293  [this, schur_asmb_pre_proc_rhs]() {
294  MoFEMFunctionBegin;
295  auto prb_ptr = schur_asmb_pre_proc_rhs->problemPtr;
296  auto dofs_prb = prb_ptr->getNumeredRowDofsPtr();
297 
298  auto crack_faces = epCorePtr->crackFaces;
299  piolaZeroDofsVec->clear();
300  CHKERR mField.getInterface<ProblemsManager>()
301  ->getSideDofsOnBrokenSpaceEntities(
302  *piolaZeroDofsVec, prb_ptr->getName(), ROW,
303  epCorePtr->piolaStress, *crack_faces, SPACE_DIM, 0,
304  SPACE_DIM);
305 
306  piolaZeroDofsMarker->resize(dofs_prb->size(), 0);
307  piolaZeroDofsMarker->clear();
308  for (auto &dof : *piolaZeroDofsVec) {
309  if (auto dof_ptr = dof.lock()) {
310  auto idx = dof_ptr->getPetscLocalDofIdx();
311  (*piolaZeroDofsMarker)[idx] = 1;
312  }
313  }
314 
315  {
316 
317  auto problem_name =
318  schur_asmb_pre_proc_rhs->problemPtr->getName();
319 
320  auto crack_faces = epCorePtr->crackFaces;
321 
322  SmartPetscObj<IS> crack_hybrid_is;
323  CHKERR epCorePtr->mField.getInterface<ISManager>()
324  ->isCreateProblemFieldAndRankLocal(
325  problem_name, ROW, epCorePtr->hybridSpatialDisp, 0,
326  SPACE_DIM, crack_hybrid_is, &*crack_faces);
327 
328  SmartPetscObj<IS> crack_piola_is;
329  CHKERR epCorePtr->mField.getInterface<ISManager>()
330  ->isCreateProblemBrokenFieldAndRankLocal(*piolaZeroDofsVec,
331  crack_piola_is);
332 
333  double *a_x;
334  CHKERR VecGetArray(schur_asmb_pre_proc_rhs->x, &a_x);
335  auto zero_by_is = [&](auto is) {
336  MoFEMFunctionBegin;
337  const PetscInt *is_array;
338  PetscInt is_size;
339  CHKERR ISGetLocalSize(is, &is_size);
340  CHKERR ISGetIndices(is, &is_array);
341  for (int i = 0; i != is_size; ++i) {
342  a_x[is_array[i]] = 0;
343  }
344  CHKERR ISRestoreIndices(is, &is_array);
345  MoFEMFunctionReturn(0);
346  };
347 
348  CHKERR zero_by_is(crack_hybrid_is);
349  CHKERR zero_by_is(crack_piola_is);
350 
351  CHKERR VecRestoreArray(schur_asmb_pre_proc_rhs->x, &a_x);
352 
353  CHKERR epCorePtr->mField.getInterface<VecManager>()
354  ->setLocalGhostVector(problem_name, COL,
355  schur_asmb_pre_proc_rhs->x,
356  INSERT_VALUES, SCATTER_REVERSE);
357  }
358 
359  MoFEMFunctionReturn(0);
360  };
361 
362  auto schur_asmb_post_proc_lhs = boost::make_shared<FEMethod>();
363  auto schur_asmb_post_proc_rhs = boost::make_shared<FEMethod>();
364 
365  schur_asmb_post_proc_rhs->postProcessHook =
366  [this, schur_asmb_post_proc_rhs]() {
367  MoFEMFunctionBegin;
368 
369  CHKERR VecGhostUpdateBegin(schur_asmb_post_proc_rhs->f,
370  ADD_VALUES, SCATTER_REVERSE);
371  CHKERR VecGhostUpdateEnd(schur_asmb_post_proc_rhs->f, ADD_VALUES,
372  SCATTER_REVERSE);
373  CHKERR VecAssemblyBegin(schur_asmb_post_proc_rhs->f);
374  CHKERR VecAssemblyEnd(schur_asmb_post_proc_rhs->f);
375  *(schur_asmb_post_proc_rhs->vecAssembleSwitch) = false;
376 
377  {
378 
379  auto problem_name =
380  schur_asmb_post_proc_rhs->problemPtr->getName();
381 
382  auto crack_faces = epCorePtr->crackFaces;
383 
384  SmartPetscObj<IS> crack_hybrid_is;
385  CHKERR epCorePtr->mField.getInterface<ISManager>()
386  ->isCreateProblemFieldAndRankLocal(
387  problem_name, ROW, epCorePtr->hybridSpatialDisp, 0,
388  SPACE_DIM, crack_hybrid_is, &*crack_faces);
389 
390  SmartPetscObj<IS> crack_piola_is;
391  CHKERR epCorePtr->mField.getInterface<ISManager>()
392  ->isCreateProblemBrokenFieldAndRankLocal(*piolaZeroDofsVec,
393  crack_piola_is);
394 
395  double *a_f;
396  CHKERR VecGetArray(schur_asmb_post_proc_rhs->f, &a_f);
397  auto zero_by_is = [&](auto is) {
398  MoFEMFunctionBegin;
399  const PetscInt *is_array;
400  PetscInt is_size;
401  CHKERR ISGetLocalSize(is, &is_size);
402  CHKERR ISGetIndices(is, &is_array);
403  for (int i = 0; i != is_size; ++i) {
404  a_f[is_array[i]] = 0;
405  }
406  CHKERR ISRestoreIndices(is, &is_array);
407  MoFEMFunctionReturn(0);
408  };
409 
410  CHKERR zero_by_is(crack_hybrid_is);
411  CHKERR zero_by_is(crack_piola_is);
412 
413  CHKERR VecRestoreArray(schur_asmb_post_proc_rhs->f, &a_f);
414 
415  }
416 
417  MoFEMFunctionReturn(0);
418  };
419 
420  schur_asmb_post_proc_lhs->postProcessHook =
421  [this, schur_asmb_post_proc_lhs]() {
422  MoFEMFunctionBegin;
423 
424  if (pMGPtr) {
425  CHKERR pMGPtr->setUP();
426  pMGPtr.reset();
427  }
428 
429  auto crack_faces = epCorePtr->crackFaces;
430 
431  // Assemble matrix
432  CHKERR MatAssemblyBegin(schur_asmb_post_proc_lhs->B,
433  MAT_FINAL_ASSEMBLY);
434  CHKERR MatAssemblyEnd(schur_asmb_post_proc_lhs->B,
435  MAT_FINAL_ASSEMBLY);
436  *(schur_asmb_post_proc_lhs->matAssembleSwitch) = false;
437  {
438  SmartPetscObj<IS> crack_hybrid_is;
439  CHKERR epCorePtr->mField.getInterface<ISManager>()
440  ->isCreateProblemFieldAndRank(
441  "ELASTIC_PROBLEM", ROW, epCorePtr->hybridSpatialDisp, 0,
442  SPACE_DIM, crack_hybrid_is, &*crack_faces);
443  CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
444  crack_hybrid_is, 1, PETSC_NULL,
445  PETSC_NULL);
446  }
447  {
448  SmartPetscObj<IS> crack_piola_is;
449  CHKERR epCorePtr->mField.getInterface<ISManager>()
450  ->isCreateProblemBrokenFieldAndRank(*piolaZeroDofsVec,
451  crack_piola_is);
452  CHKERR MatZeroRowsColumnsIS(schur_asmb_post_proc_lhs->B,
453  crack_piola_is, 1, PETSC_NULL,
454  PETSC_NULL);
455  }
456 
457  auto [a00_field_list, a00_range_list] = getA00Fields();
458  CHKERR assembleBlockMatSchur(epCorePtr->mField,
459  schur_asmb_post_proc_lhs->B, S,
460  a00_field_list, a00_range_list, aoS);
461  epCorePtr->a00FieldList = a00_field_list;
462  epCorePtr->a00RangeList = a00_range_list;
463 
464  // Apply essential constrains to Schur complement
465  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
466  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
467 
468  SmartPetscObj<IS> crack_hybrid_is;
469  CHKERR epCorePtr->mField.getInterface<ISManager>()
470  ->isCreateProblemFieldAndRank(
471  "SUB_SCHUR", ROW, epCorePtr->hybridSpatialDisp, 0,
472  SPACE_DIM, crack_hybrid_is, &*crack_faces);
473  epCorePtr->crackHybridIs = crack_hybrid_is;
474  CHKERR MatZeroRowsColumnsIS(S, crack_hybrid_is, 1, PETSC_NULL,
475  PETSC_NULL);
476 
477  MoFEMFunctionReturn(0);
478  };
479 
480  auto ts_ctx_ptr = getDMTsCtx(epCorePtr->dmElastic);
481  ts_ctx_ptr->getPreProcessIFunction().push_front(
482  schur_asmb_pre_proc_rhs);
483  ts_ctx_ptr->getPostProcessIFunction().push_back(
484  schur_asmb_post_proc_rhs);
485  ts_ctx_ptr->getPreProcessIJacobian().push_front(
486  schur_asmb_pre_proc_lhs);
487  ts_ctx_ptr->getPostProcessIJacobian().push_back(
488  schur_asmb_post_proc_lhs);
489  MoFEMFunctionReturn(0);
490  };
491 
492  piolaZeroDofsVec =
493  boost::make_shared<std::vector<boost::weak_ptr<NumeredDofEntity>>>();
494  piolaZeroDofsMarker = boost::make_shared<std::vector<unsigned char>>();
495  CHKERR set_assemble();
496 
497  MoFEMFunctionReturn(0);
498  };
499 
500  auto set_pc = [&]() {
501  MoFEMFunctionBegin;
502  auto a00_is = getDMSubData(a00_dm)->getSmartRowIs();
503  auto schur_is = getDMSubData(schur_dm)->getSmartRowIs();
504  CHKERR PCFieldSplitSetIS(pc, NULL, a00_is);
505  CHKERR PCFieldSplitSetIS(pc, NULL, schur_is);
506  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
507  MoFEMFunctionReturn(0);
508  };
509 
510  auto set_diagonal_pc = [&]() {
511  MoFEMFunctionBegin;
512  KSP *subksp;
513  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
514  auto get_pc = [](auto ksp) {
515  PC pc_raw;
516  CHKERR KSPGetPC(ksp, &pc_raw);
517  return SmartPetscObj<PC>(pc_raw, true); // bump reference
518  };
519  CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
520 
521  auto set_pc_p_mg = [&](auto dm, auto pc, auto S) {
522  MoFEMFunctionBegin;
523  CHKERR PCSetDM(pc, dm);
524  PetscBool same = PETSC_FALSE;
525  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
526  if (same) {
527  auto smart_pc = SmartPetscObj<PC>(pc, true);
528  pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc, S);
529  }
530  PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
531  if (same) {
532  MOFEM_LOG("EP", Sev::inform)
533  << "SetUpSchurImpl::setUp: fieldsplit 1 PCKSP";
534  CHKERR PCSetFromOptions(pc);
535  KSP ksp;
536  CHKERR PCKSPGetKSP(pc, &ksp);
537  PC ksp_pc;
538  CHKERR KSPGetPC(ksp, &ksp_pc);
539  CHKERR PCSetFromOptions(ksp_pc);
540  PetscObjectTypeCompare((PetscObject)ksp_pc, PCMG, &same);
541  if (same) {
542  auto smart_pc = SmartPetscObj<PC>(ksp_pc, true);
543  pMGPtr = boost::make_shared<P_MultiGridData>(dm, smart_pc, S);
544  }
545  }
546  MoFEMFunctionReturn(0);
547  };
548 
549  CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]), S);
550 
551  CHKERR PetscFree(subksp);
552  MoFEMFunctionReturn(0);
553  };
554 
555  CHKERR set_assembly();
556  CHKERR set_pc();
557  CHKERR TSSetUp(ts);
558  CHKERR KSPSetUp(ksp);
559  CHKERR set_diagonal_pc();
560 
561  } else {
562  MOFEM_LOG("EP", Sev::inform) << "SetUpSchurImpl::setUp: PCLU or other";
563 
564  epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
565  createOpSchurAssembleBegin());
566  epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
567  createOpSchurAssembleEnd({}, {}));
568  epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
569  createOpSchurAssembleBegin());
570  epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
571  createOpSchurAssembleEnd({}, {}));
572  }
573 
574  MoFEMFunctionReturn(0);
575 }
576 
577 boost::shared_ptr<EshelbianCore::SetUpSchur>
578 EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
579  EshelbianCore *ep_core_ptr) {
580  return boost::shared_ptr<SetUpSchur>(
581  new SetUpSchurImpl(m_field, ep_core_ptr));
582 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SetUpSchurImpl::P_MultiGridData::P_MultiGridData
P_MultiGridData(SmartPetscObj< DM > dm, SmartPetscObj< PC > pc, SmartPetscObj< Mat > S)
Definition: SetUpSchurImpl.cpp:33
SetUpSchurImpl
Definition: test_broken_space.cpp:519
SetUpSchurImpl::aoS
SmartPetscObj< AO > aoS
Definition: SetUpSchurImpl.cpp:24
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
EshelbianCore::a00FieldList
std::vector< std::string > a00FieldList
Definition: EshelbianCore.hpp:315
EshelbianCore::S
Mat S
Definition: EshelbianCore.hpp:312
EshelbianCore::spatialL2Disp
const std::string spatialL2Disp
Definition: EshelbianCore.hpp:102
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:532
EshelbianCore::SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
Definition: SetUpSchurImpl.cpp:578
EshelbianCore::elasticFeLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
Definition: EshelbianCore.hpp:90
EshelbianCore::noStretch
static PetscBool noStretch
Definition: EshelbianCore.hpp:18
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
SetUpSchurImpl::epCorePtr
EshelbianCore * epCorePtr
Definition: SetUpSchurImpl.cpp:21
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:630
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: SetUpSchurImpl.cpp:13
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::TsSetI2Jacobian
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.
Definition: TsCtx.cpp:519
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
EshelbianCore::spaceOrder
int spaceOrder
Definition: EshelbianCore.hpp:124
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
EshelbianCore::crackHybridIs
SmartPetscObj< IS > crackHybridIs
Definition: EshelbianCore.hpp:314
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition: DMMoFEM.hpp:1076
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EshelbianCore::contactDisp
const std::string contactDisp
Definition: EshelbianCore.hpp:108
ROW
@ ROW
Definition: definitions.h:136
EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianCore.hpp:96
SetUpSchurImpl::P_MultiGridData::pc
SmartPetscObj< PC > pc
Definition: SetUpSchurImpl.cpp:47
MoFEM::createOpSchurAssembleEnd
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.
Definition: Schur.cpp:2585
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:529
EshelbianCore::alphaRho
double alphaRho
Definition: EshelbianCore.hpp:130
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
Definition: SetUpSchurImpl.cpp:11
SetUpSchurImpl::pMGPtr
boost::shared_ptr< P_MultiGridData > pMGPtr
Definition: SetUpSchurImpl.cpp:51
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
EshelbianCore::contactElement
const std::string contactElement
Definition: EshelbianCore.hpp:117
SetUpSchurImpl::piolaZeroDofsMarker
boost::shared_ptr< std::vector< unsigned char > > piolaZeroDofsMarker
Definition: SetUpSchurImpl.cpp:29
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
SetUpSchurImpl::getSchurFields
auto getSchurFields()
Definition: SetUpSchurImpl.cpp:53
SetUpSchurImpl::P_MultiGridData
Definition: SetUpSchurImpl.cpp:31
EshelbianCore::rotAxis
const std::string rotAxis
Definition: EshelbianCore.hpp:110
MoFEM::createDMHybridisedL2Matrix
auto createDMHybridisedL2Matrix(DM dm)
Get smart hybridised L2 matrix from DM.
Definition: DMMoFEM.hpp:1069
SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:528
EshelbianCore::elementVolumeName
const std::string elementVolumeName
Definition: EshelbianCore.hpp:113
COL
@ COL
Definition: definitions.h:136
EshelbianCore::aoS
AO aoS
Definition: EshelbianCore.hpp:313
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2580
EshelbianCore::mField
MoFEM::Interface & mField
Definition: EshelbianCore.hpp:84
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
MoFEM::TsSetIJacobian
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.
Definition: TsCtx.cpp:169
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EshelbianCore::hybridSpatialDisp
const std::string hybridSpatialDisp
Definition: EshelbianCore.hpp:106
EshelbianCore::stretchTensor
const std::string stretchTensor
Definition: EshelbianCore.hpp:109
MoFEM::assembleBlockMatSchur
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.
Definition: Schur.cpp:1817
EshelbianCore::a00RangeList
std::vector< boost::shared_ptr< Range > > a00RangeList
Definition: EshelbianCore.hpp:317
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2627
EshelbianCore::skeletonElement
const std::string skeletonElement
Definition: EshelbianCore.hpp:116
SetUpSchurImpl::piolaZeroDofsVec
boost::shared_ptr< std::vector< boost::weak_ptr< NumeredDofEntity > > > piolaZeroDofsVec
Definition: SetUpSchurImpl.cpp:27
EshelbianCore::SetUpSchur
Definition: EshelbianCore.hpp:264
SetUpSchurImpl::P_MultiGridData::S
SmartPetscObj< Mat > S
Definition: SetUpSchurImpl.cpp:48
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SetUpSchurImpl::P_MultiGridData::setUP
MoFEMErrorCode setUP()
Definition: SetUpSchurImpl.cpp:37
EshelbianCore::bubbleField
const std::string bubbleField
Definition: EshelbianCore.hpp:111
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
NBFACETRI_L2
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
Definition: h1_hdiv_hcurl_l2.h:42
MoFEM::createSchurNestedMatrixStruture
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.
Definition: Schur.cpp:2343
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
EshelbianCore::elasticBcLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
Definition: EshelbianCore.hpp:91
SetUpSchurImpl::P_MultiGridData::dm
SmartPetscObj< DM > dm
Definition: SetUpSchurImpl.cpp:43
EshelbianCore::piolaStress
const std::string piolaStress
Definition: EshelbianCore.hpp:100
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:44
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1082
MoFEM::PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:634
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
EshelbianCore::crackFaces
boost::shared_ptr< Range > crackFaces
Definition: EshelbianCore.hpp:284
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
SetUpSchurImpl::getA00Fields
auto getA00Fields()
Definition: SetUpSchurImpl.cpp:60
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
EshelbianCore
Definition: EshelbianCore.hpp:12
Generated by Doxygen 1.8.17 and hosted at University of Glasgow