v0.14.0
elastic.cpp
Go to the documentation of this file.
1 /**
2  * @file elastic.cpp
3  * @brief elastic example
4  * @version 0.13.2
5  * @date 2022-09-19
6  *
7  * @copyright Copyright (c) 2022
8  *
9  */
10 
11 #include <MoFEM.hpp>
12 
13 using namespace MoFEM;
14 
15 constexpr int BASE_DIM = 1; //< Dimension of the base functions
16 
17 //! [Define dimension]
18 constexpr int SPACE_DIM =
19  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
20 //! [Define dimension]
21 constexpr AssemblyType A = (SCHUR_ASSEMBLE)
22  ? AssemblyType::BLOCK_SCHUR
23  : AssemblyType::PETSC; //< selected assembly type
24 
25 constexpr IntegrationType I =
26  IntegrationType::GAUSS; //< selected integration type
27 
28 //! [Define entities]
29 using EntData = EntitiesFieldData::EntData;
30 using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
31 using BoundaryEle =
32  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
33 using DomainEleOp = DomainEle::UserDataOperator;
34 using BoundaryEleOp = BoundaryEle::UserDataOperator;
35 //! [Define entities]
36 
37 //! [OpK]
38 using OpK = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
39  I>::OpGradSymTensorGrad<BASE_DIM, SPACE_DIM, SPACE_DIM, 0>;
40 //! [OpK]
41 //! [OpInternalForce]
42 using OpInternalForce = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
43  I>::OpGradTimesSymTensor<BASE_DIM, SPACE_DIM, SPACE_DIM>;
44 //! [OpInternalForce]
45 struct DomainBCs {};
46 struct BoundaryBCs {};
47 
48 using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<I>;
49 using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
50 using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>;
51 using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
52 using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<I>;
53 using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
54 
55 template <int DIM> struct PostProcEleByDim;
56 
57 template <> struct PostProcEleByDim<2> {
58  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
59  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
60  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
61 };
62 
63 template <> struct PostProcEleByDim<3> {
64  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
65  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
66  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
67 };
68 
69 using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
70 using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
71 using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
72 
73 #include <ElasticSpring.hpp>
74 #include <FluidLevel.hpp>
75 #include <CalculateTraction.hpp>
76 #include <NaturalDomainBC.hpp>
77 #include <NaturalBoundaryBC.hpp>
78 
79 constexpr double young_modulus = 1;
80 constexpr double poisson_ratio = 0.3;
81 constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
82 constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
83 
84 PetscBool is_plane_strain = PETSC_FALSE;
85 
86 struct Example {
87 
88  Example(MoFEM::Interface &m_field) : mField(m_field) {}
89 
90  MoFEMErrorCode runProblem();
91 
92 private:
93  MoFEM::Interface &mField;
94 
95  PetscBool doEvalField;
96  std::array<double, SPACE_DIM> fieldEvalCoords;
97  boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
98  boost::shared_ptr<MatrixDouble> vectorFieldPtr;
99 
100  MoFEMErrorCode readMesh();
101  MoFEMErrorCode setupProblem();
102  MoFEMErrorCode boundaryCondition();
103  MoFEMErrorCode assembleSystem();
104  MoFEMErrorCode solveSystem();
105  MoFEMErrorCode outputResults();
106  MoFEMErrorCode checkResults();
107 
108  MoFEMErrorCode addMatBlockOps(
109  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
110  std::string field_name, std::string block_name,
111  boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev);
112 };
113 
114 MoFEMErrorCode Example::addMatBlockOps(
115  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
116  std::string field_name, std::string block_name,
117  boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev) {
118  MoFEMFunctionBegin;
119 
120  struct OpMatBlocks : public DomainEleOp {
121  OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
122  double bulk_modulus_K, double shear_modulus_G,
123  MoFEM::Interface &m_field, Sev sev,
124  std::vector<const CubitMeshSets *> meshset_vec_ptr)
125  : DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
126  bulkModulusKDefault(bulk_modulus_K),
127  shearModulusGDefault(shear_modulus_G) {
128  std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
129  CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
130  "Can not get data from block");
131  }
132 
133  MoFEMErrorCode doWork(int side, EntityType type,
134  EntitiesFieldData::EntData &data) {
135  MoFEMFunctionBegin;
136 
137  for (auto &b : blockData) {
138 
139  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
140  CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
141  MoFEMFunctionReturnHot(0);
142  }
143  }
144 
145  CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
146  MoFEMFunctionReturn(0);
147  }
148 
149  private:
150  boost::shared_ptr<MatrixDouble> matDPtr;
151 
152  struct BlockData {
153  double bulkModulusK;
154  double shearModulusG;
155  Range blockEnts;
156  };
157 
158  double bulkModulusKDefault;
159  double shearModulusGDefault;
160  std::vector<BlockData> blockData;
161 
162  MoFEMErrorCode
163  extractBlockData(MoFEM::Interface &m_field,
164  std::vector<const CubitMeshSets *> meshset_vec_ptr,
165  Sev sev) {
166  MoFEMFunctionBegin;
167 
168  for (auto m : meshset_vec_ptr) {
169  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
170  std::vector<double> block_data;
171  CHKERR m->getAttributes(block_data);
172  if (block_data.size() < 2) {
173  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
174  "Expected that block has two attributes");
175  }
176  auto get_block_ents = [&]() {
177  Range ents;
178  CHKERR
179  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
180  return ents;
181  };
182 
183  double young_modulus = block_data[0];
184  double poisson_ratio = block_data[1];
185  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
186  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
187 
188  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
189  << "E = " << young_modulus << " nu = " << poisson_ratio;
190 
191  blockData.push_back(
192  {bulk_modulus_K, shear_modulus_G, get_block_ents()});
193  }
194  MOFEM_LOG_CHANNEL("WORLD");
195  MoFEMFunctionReturn(0);
196  }
197 
198  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
199  double bulk_modulus_K, double shear_modulus_G) {
200  MoFEMFunctionBegin;
201  //! [Calculate elasticity tensor]
202  auto set_material_stiffness = [&]() {
203  FTensor::Index<'i', SPACE_DIM> i;
204  FTensor::Index<'j', SPACE_DIM> j;
205  FTensor::Index<'k', SPACE_DIM> k;
206  FTensor::Index<'l', SPACE_DIM> l;
207  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
208  double A = 1.;
209  if (SPACE_DIM == 2 && !is_plane_strain) {
210  A = 2 * shear_modulus_G /
211  (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
212  }
213  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
214  t_D(i, j, k, l) =
215  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
216  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
217  t_kd(k, l);
218  };
219  //! [Calculate elasticity tensor]
220  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
221  mat_D_ptr->resize(size_symm * size_symm, 1);
222  set_material_stiffness();
223  MoFEMFunctionReturn(0);
224  }
225  };
226 
227  pipeline.push_back(new OpMatBlocks(
228  field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, mField, sev,
229 
230  // Get blockset using regular expression
231  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
232 
233  (boost::format("%s(.*)") % block_name).str()
234 
235  ))
236 
237  ));
238 
239  MoFEMFunctionReturn(0);
240 }
241 
242 //! [Run problem]
243 MoFEMErrorCode Example::runProblem() {
244  MoFEMFunctionBegin;
245  CHKERR readMesh();
246  CHKERR setupProblem();
247  CHKERR boundaryCondition();
248  CHKERR assembleSystem();
249  CHKERR solveSystem();
250  CHKERR outputResults();
251  CHKERR checkResults();
252  MoFEMFunctionReturn(0);
253 }
254 //! [Run problem]
255 
256 //! [Read mesh]
257 MoFEMErrorCode Example::readMesh() {
258  MoFEMFunctionBegin;
259  auto simple = mField.getInterface<Simple>();
260  CHKERR simple->getOptions();
261  CHKERR simple->loadFile();
262  MoFEMFunctionReturn(0);
263 }
264 //! [Read mesh]
265 
266 //! [Set up problem]
267 MoFEMErrorCode Example::setupProblem() {
268  MoFEMFunctionBegin;
269  Simple *simple = mField.getInterface<Simple>();
270 
271  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
272  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
273  PetscInt choice_base_value = AINSWORTH;
274  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
275  LASBASETOPT, &choice_base_value, PETSC_NULL);
276 
277  FieldApproximationBase base;
278  switch (choice_base_value) {
279  case AINSWORTH:
280  base = AINSWORTH_LEGENDRE_BASE;
281  MOFEM_LOG("WORLD", Sev::inform)
282  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
283  break;
284  case DEMKOWICZ:
285  base = DEMKOWICZ_JACOBI_BASE;
286  MOFEM_LOG("WORLD", Sev::inform)
287  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
288  break;
289  default:
290  base = LASTBASE;
291  break;
292  }
293 
294  // Add field
295  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
296  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
297  int order = 3;
298  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
299 
300  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
301 
302  CHKERR simple->setFieldOrder("U", order);
303  CHKERR simple->setFieldOrder("GEOMETRY", 2);
304  CHKERR simple->setUp();
305 
306  auto project_ho_geometry = [&]() {
307  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
308  return mField.loop_dofs("GEOMETRY", ent_method);
309  };
310  CHKERR project_ho_geometry();
311 
312  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-plane_strain", &is_plane_strain,
313  PETSC_NULL);
314 
315  int coords_dim = SPACE_DIM;
316  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
317  fieldEvalCoords.data(), &coords_dim,
318  &doEvalField);
319 
320  if (doEvalField) {
321  vectorFieldPtr = boost::make_shared<MatrixDouble>();
322  fieldEvalData =
323  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
324 
325  if constexpr (SPACE_DIM == 3) {
326  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
327  fieldEvalData, simple->getDomainFEName());
328  } else {
329  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
330  fieldEvalData, simple->getDomainFEName());
331  }
332 
333  fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
334  auto no_rule = [](int, int, int) { return -1; };
335  auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
336  field_eval_fe_ptr->getRuleHook = no_rule;
337 
338  field_eval_fe_ptr->getOpPtrVector().push_back(
339  new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
340  }
341 
342  MoFEMFunctionReturn(0);
343 }
344 //! [Set up problem]
345 
346 //! [Boundary condition]
347 MoFEMErrorCode Example::boundaryCondition() {
348  MoFEMFunctionBegin;
349  auto simple = mField.getInterface<Simple>();
350  auto bc_mng = mField.getInterface<BcManager>();
351 
352  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
353  "U", 0, 0);
354  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
355  "U", 1, 1);
356  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
357  "U", 2, 2);
358  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
359  "REMOVE_ALL", "U", 0, 3);
360  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
361  simple->getProblemName(), "U");
362 
363  // adding MPCs
364  CHKERR bc_mng->addBlockDOFsToMPCs(simple->getProblemName(), "U");
365 
366  MoFEMFunctionReturn(0);
367 }
368 //! [Boundary condition]
369 
370 //! [Push operators to pipeline]
371 MoFEMErrorCode Example::assembleSystem() {
372  MoFEMFunctionBegin;
373  auto pip = mField.getInterface<PipelineManager>();
374 
375  //! [Integration rule]
376  auto integration_rule = [](int, int, int approx_order) {
377  return 2 * approx_order + 1;
378  };
379 
380  auto integration_rule_bc = [](int, int, int approx_order) {
381  return 2 * approx_order + 1;
382  };
383 
384  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
385  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
386  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
387  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
388  //! [Integration rule]
389 
390  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
391  pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
392  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
393  pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
394  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
395  pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
396  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
397  pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
398 
399  //! [Push domain stiffness matrix to pipeline]
400  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
401 
402  // Assemble domain stiffness matrix
403  CHKERR addMatBlockOps(pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC",
404  mat_D_ptr, Sev::verbose);
405  pip->getOpDomainLhsPipeline().push_back(new OpK("U", "U", mat_D_ptr));
406  //! [Push domain stiffness matrix to pipeline]
407 
408  //! [Push Internal forces]
409  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
410  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
411  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
412  pip->getOpDomainRhsPipeline().push_back(
413  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
414  mat_grad_ptr));
415  CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
416  mat_D_ptr, Sev::inform);
417  pip->getOpDomainRhsPipeline().push_back(
418  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
419  pip->getOpDomainRhsPipeline().push_back(
420  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
421  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
422 
423  pip->getOpDomainRhsPipeline().push_back(
424  new OpInternalForce("U", mat_stress_ptr,
425  [](double, double, double) constexpr { return -1; }));
426  //! [Push Internal forces]
427 
428  //! [Push Body forces]
429  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
430  pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
431  //! [Push Body forces]
432 
433  //! [Push natural boundary conditions]
434  // Add force boundary condition
435  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
436  pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
437  // Add case for mix type of BCs
438  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
439  pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
440  //! [Push natural boundary conditions]
441  MoFEMFunctionReturn(0);
442 }
443 //! [Push operators to pipeline]
444 
445 struct SetUpSchur {
446  static boost::shared_ptr<SetUpSchur>
447  createSetUpSchur(MoFEM::Interface &m_field);
448  virtual MoFEMErrorCode setUp(SmartPetscObj<KSP> solver) = 0;
449 
450 protected:
451  SetUpSchur() = default;
452 };
453 
454 //! [Solve]
455 MoFEMErrorCode Example::solveSystem() {
456  MoFEMFunctionBegin;
457  auto simple = mField.getInterface<Simple>();
458  auto pip = mField.getInterface<PipelineManager>();
459  auto solver = pip->createKSP();
460  CHKERR KSPSetFromOptions(solver);
461 
462  auto dm = simple->getDM();
463  auto D = createDMVector(dm);
464  auto F = vectorDuplicate(D);
465 
466  auto set_essential_bc = [&]() {
467  MoFEMFunctionBegin;
468  // This is low level pushing finite elements (pipelines) to solver
469  auto ksp_ctx_ptr = getDMKspCtx(dm);
470 
471  auto pre_proc_rhs = boost::make_shared<FEMethod>();
472  auto post_proc_rhs = boost::make_shared<FEMethod>();
473  auto post_proc_lhs = boost::make_shared<FEMethod>();
474 
475  auto get_pre_proc_hook = [&]() {
476  return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
477  {});
478  };
479  pre_proc_rhs->preProcessHook = get_pre_proc_hook();
480 
481  auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
482  MoFEMFunctionBeginHot;
483 
484  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
485  post_proc_rhs, 1.)();
486  CHKERR EssentialPostProcRhs<MPCsType>(mField, post_proc_rhs, 1.)();
487  MoFEMFunctionReturnHot(0);
488  };
489 
490  auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
491  MoFEMFunctionBeginHot;
492 
493  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
494  post_proc_lhs, 1.)();
495  CHKERR EssentialPostProcLhs<MPCsType>(mField, post_proc_lhs)();
496  MoFEMFunctionReturnHot(0);
497  };
498 
499  post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
500  post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
501 
502  ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
503  ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
504  ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
505  MoFEMFunctionReturn(0);
506  };
507 
508  auto setup_and_solve = [&]() {
509  MoFEMFunctionBegin;
510  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
511  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
512  CHKERR KSPSetUp(solver);
513  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
514  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
515  CHKERR KSPSolve(solver, F, D);
516  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
517  MoFEMFunctionReturn(0);
518  };
519 
520  MOFEM_LOG_CHANNEL("TIMER");
521  MOFEM_LOG_TAG("TIMER", "timer");
522 
523  CHKERR set_essential_bc();
524 
525  if (A == AssemblyType::BLOCK_SCHUR || A == AssemblyType::SCHUR) {
526  auto schur_ptr = SetUpSchur::createSetUpSchur(mField);
527  CHKERR schur_ptr->setUp(solver);
528  CHKERR setup_and_solve();
529  } else {
530  CHKERR setup_and_solve();
531  }
532 
533  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
534  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
535  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
536 
537  if (doEvalField) {
538  if constexpr (SPACE_DIM == 3) {
539  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
540  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
541  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
542  mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
543  } else {
544  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
545  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
546  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
547  mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
548  }
549 
550  if (vectorFieldPtr->size1()) {
551  auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
552  if constexpr (SPACE_DIM == 2)
553  MOFEM_LOG("FieldEvaluator", Sev::inform)
554  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
555  else
556  MOFEM_LOG("FieldEvaluator", Sev::inform)
557  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
558  << " U_Z: " << t_disp(2);
559  }
560 
561  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
562  }
563 
564  MoFEMFunctionReturn(0);
565 }
566 //! [Solve]
567 
568 //! [Postprocess results]
569 MoFEMErrorCode Example::outputResults() {
570  MoFEMFunctionBegin;
571  auto simple = mField.getInterface<Simple>();
572  auto pip = mField.getInterface<PipelineManager>();
573  auto det_ptr = boost::make_shared<VectorDouble>();
574  auto jac_ptr = boost::make_shared<MatrixDouble>();
575  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
576  //! [Postprocess clean]
577  pip->getDomainRhsFE().reset();
578  pip->getDomainLhsFE().reset();
579  pip->getBoundaryRhsFE().reset();
580  pip->getBoundaryLhsFE().reset();
581  //! [Postprocess clean]
582 
583  //! [Postprocess initialise]
584  auto post_proc_mesh = boost::make_shared<moab::Core>();
585  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
586  mField, post_proc_mesh);
587  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
588  mField, post_proc_mesh);
589  //! [Postprocess initialise]
590 
591  auto calculate_stress_ops = [&](auto &pip) {
592  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
593  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
594  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
595  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
596  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
597  "U", mat_grad_ptr));
598  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
599  CHKERR addMatBlockOps(pip, "U", "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
600  pip.push_back(
601  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
602  pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
603  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
604  auto u_ptr = boost::make_shared<MatrixDouble>();
605  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
606  auto x_ptr = boost::make_shared<MatrixDouble>();
607  pip.push_back(
608  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
609  return boost::make_tuple(u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr);
610  };
611 
612  auto get_tag_id_on_pmesh = [&](bool post_proc_skin)
613  {
614  int def_val_int = 0;
615  Tag tag_mat;
616  CHKERR mField.get_moab().tag_get_handle(
617  "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
618  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
619  auto meshset_vec_ptr =
620  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
621  std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
622 
623  Range skin_ents;
624  std::unique_ptr<Skinner> skin_ptr;
625  if (post_proc_skin) {
626  skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
627  auto boundary_meshset = simple->getBoundaryMeshSet();
628  CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
629  skin_ents, true);
630  }
631 
632  for (auto m : meshset_vec_ptr) {
633  Range ents_3d;
634  CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
635  true);
636  int const id = m->getMeshsetId();
637  ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
638  if (post_proc_skin) {
639  Range skin_faces;
640  CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
641  ents_3d = intersect(skin_ents, skin_faces);
642  }
643  CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
644  }
645 
646  return tag_mat;
647  };
648 
649  auto post_proc_domain = [&](auto post_proc_mesh) {
650  auto post_proc_fe =
651  boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
652  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
653 
654  auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
655  calculate_stress_ops(post_proc_fe->getOpPtrVector());
656 
657  post_proc_fe->getOpPtrVector().push_back(
658 
659  new OpPPMap(
660 
661  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
662 
663  {},
664 
665  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
666 
667  {},
668 
669  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
670 
671  )
672 
673  );
674 
675  post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
676  return post_proc_fe;
677  };
678 
679  auto post_proc_boundary = [&](auto post_proc_mesh) {
680  auto post_proc_fe =
681  boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
682  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
683  post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
684  auto op_loop_side =
685  new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
686  // push ops to side element, through op_loop_side operator
687  auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
688  calculate_stress_ops(op_loop_side->getOpPtrVector());
689  post_proc_fe->getOpPtrVector().push_back(op_loop_side);
690  auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
691  post_proc_fe->getOpPtrVector().push_back(
692  new ElasticExample::OpCalculateTraction(mat_stress_ptr,
693  mat_traction_ptr));
694 
695  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
696 
697  post_proc_fe->getOpPtrVector().push_back(
698 
699  new OpPPMap(
700 
701  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
702 
703  {},
704 
705  {{"U", u_ptr}, {"GEOMETRY", x_ptr}, {"T", mat_traction_ptr}},
706 
707  {},
708 
709  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
710 
711  )
712 
713  );
714 
715  post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
716  return post_proc_fe;
717  };
718 
719  PetscBool post_proc_skin_only = PETSC_FALSE;
720  if (SPACE_DIM == 3) {
721  post_proc_skin_only = PETSC_TRUE;
722  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin_only",
723  &post_proc_skin_only, PETSC_NULL);
724  }
725  if (post_proc_skin_only == PETSC_FALSE) {
726  pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
727  } else {
728  pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
729  }
730 
731  CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
732  post_proc_begin->getFEMethod());
733  CHKERR pip->loopFiniteElements();
734  CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
735  post_proc_end->getFEMethod());
736 
737  CHKERR post_proc_end->writeFile("out_elastic.h5m");
738  MoFEMFunctionReturn(0);
739 }
740 //! [Postprocess results]
741 
742 //! [Check]
743 MoFEMErrorCode Example::checkResults() {
744  MOFEM_LOG_CHANNEL("WORLD");
745  auto simple = mField.getInterface<Simple>();
746  auto pip = mField.getInterface<PipelineManager>();
747  MoFEMFunctionBegin;
748  pip->getDomainRhsFE().reset();
749  pip->getDomainLhsFE().reset();
750  pip->getBoundaryRhsFE().reset();
751  pip->getBoundaryLhsFE().reset();
752 
753  auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
754  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
755  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule);
756 
757  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
758  pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
759  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
760  pip->getOpBoundaryRhsPipeline(), {}, "GEOMETRY");
761 
762  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
763  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
764  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
765 
766  pip->getOpDomainRhsPipeline().push_back(
767  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
768  mat_grad_ptr));
769  pip->getOpDomainRhsPipeline().push_back(
770  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
771 
772  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
773  CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
774  mat_D_ptr, Sev::verbose);
775  pip->getOpDomainRhsPipeline().push_back(
776  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
777  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
778 
779  pip->getOpDomainRhsPipeline().push_back(
780  new OpInternalForce("U", mat_stress_ptr));
781 
782  pip->getOpBoundaryRhsPipeline().push_back(
783  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
784  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
785  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
786  pip->getOpDomainRhsPipeline(), mField, "U", Sev::verbose);
787  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
788  pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::verbose);
789 
790  auto dm = simple->getDM();
791  auto res = createDMVector(dm);
792  CHKERR VecSetDM(res, PETSC_NULL);
793 
794  pip->getDomainRhsFE()->f = res;
795  pip->getBoundaryRhsFE()->f = res;
796 
797  CHKERR VecZeroEntries(res);
798 
799  CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
800  CHKERR pip->loopFiniteElements();
801  CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
802 
803  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
804  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
805  CHKERR VecAssemblyBegin(res);
806  CHKERR VecAssemblyEnd(res);
807 
808  auto zero_residual_at_constrains = [&]() {
809  MoFEMFunctionBegin;
810  auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
811  auto get_post_proc_hook_rhs = [&]() {
812  MoFEMFunctionBegin;
813  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
814  mField, fe_post_proc_ptr, res)();
815  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
816  mField, fe_post_proc_ptr, 0, res)();
817  CHKERR EssentialPostProcRhs<MPCsType>(mField, fe_post_proc_ptr, 0, res)();
818  MoFEMFunctionReturn(0);
819  };
820  fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
821  CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
822  MoFEMFunctionReturn(0);
823  };
824 
825  CHKERR zero_residual_at_constrains();
826 
827  double nrm2;
828  CHKERR VecNorm(res, NORM_2, &nrm2);
829  MOFEM_LOG_CHANNEL("WORLD");
830  MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
831 
832  PetscBool test = PETSC_FALSE;
833  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
834  if (test == PETSC_TRUE) {
835 
836  auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
837  MoFEMFunctionBegin;
838  auto post_proc_fe =
839  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
840  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
841  auto u_vec = boost::make_shared<MatrixDouble>();
842  post_proc_fe->getOpPtrVector().push_back(
843  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
844  post_proc_fe->getOpPtrVector().push_back(
845 
846  new OpPPMap(
847 
848  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
849 
850  {},
851 
852  {{"RES", u_vec}},
853 
854  {}, {})
855 
856  );
857 
858  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
859  post_proc_fe);
860  post_proc_fe->writeFile(out_name);
861  MoFEMFunctionReturn(0);
862  };
863 
864  CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
865 
866  constexpr double eps = 1e-8;
867  if (nrm2 > eps)
868  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
869  "Residual is not zero");
870  }
871 
872  MoFEMFunctionReturn(0);
873 }
874 //! [Check]
875 
876 static char help[] = "...\n\n";
877 
878 int main(int argc, char *argv[]) {
879 
880  // Initialisation of MoFEM/PETSc and MOAB data structures
881  const char param_file[] = "param_file.petsc";
882  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
883 
884  auto core_log = logging::core::get();
885  core_log->add_sink(
886  LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
887 
888  core_log->add_sink(
889  LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
890  LogManager::setLog("FieldEvaluator");
891  MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
892 
893  try {
894 
895  //! [Register MoFEM discrete manager in PETSc]
896  DMType dm_name = "DMMOFEM";
897  CHKERR DMRegister_MoFEM(dm_name);
898  DMType dm_name_mg = "DMMOFEM_MG";
899  CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
900  //! [Register MoFEM discrete manager in PETSc
901 
902  //! [Create MoAB]
903  moab::Core mb_instance; ///< mesh database
904  moab::Interface &moab = mb_instance; ///< mesh database interface
905  //! [Create MoAB]
906 
907  //! [Create MoFEM]
908  MoFEM::Core core(moab); ///< finite element database
909  MoFEM::Interface &m_field = core; ///< finite element database interface
910  //! [Create MoFEM]
911 
912  //! [Example]
913  Example ex(m_field);
914  CHKERR ex.runProblem();
915  //! [Example]
916  }
917  CATCH_ERRORS;
918 
919  CHKERR MoFEM::Core::Finalize();
920 }
921 
922 struct SetUpSchurImpl : public SetUpSchur {
923 
924  SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {
925  if (S) {
926  CHK_THROW_MESSAGE(
927  MOFEM_DATA_INCONSISTENCY,
928  "Is expected that schur matrix is not allocated. This is "
929  "possible only is PC is set up twice");
930  }
931  }
932  virtual ~SetUpSchurImpl() = default;
933 
934  MoFEMErrorCode setUp(SmartPetscObj<KSP> solver);
935 
936 private:
937  MoFEMErrorCode setEntities();
938  MoFEMErrorCode createSubDM();
939  MoFEMErrorCode setOperator();
940  MoFEMErrorCode setPC(PC pc);
941  MoFEMErrorCode setDiagonalPC(PC pc);
942 
943  SmartPetscObj<Mat> S;
944 
945  MoFEM::Interface &mField;
946 
947  SmartPetscObj<DM> schurDM;
948  SmartPetscObj<DM> blockDM;
949  Range volEnts;
950  Range subEnts;
951 };
952 
953 MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<KSP> solver) {
954  MoFEMFunctionBegin;
955  auto pip = mField.getInterface<PipelineManager>();
956  PC pc;
957  CHKERR KSPGetPC(solver, &pc);
958  PetscBool is_pcfs = PETSC_FALSE;
959  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
960  if (is_pcfs) {
961  if (S) {
962  CHK_THROW_MESSAGE(
963  MOFEM_DATA_INCONSISTENCY,
964  "Is expected that schur matrix is not allocated. This is "
965  "possible only is PC is set up twice");
966  }
967  CHKERR setEntities();
968  CHKERR createSubDM();
969 
970  // Add data to DM storage
971  S = createDMMatrix(schurDM);
972  CHKERR MatSetDM(S, PETSC_NULL);
973  CHKERR MatSetBlockSize(S, SPACE_DIM);
974  CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
975 
976  CHKERR setOperator();
977  CHKERR setPC(pc);
978 
979  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
980  // Set DM to use shell block matrix
981  DM solver_dm;
982  CHKERR KSPGetDM(solver, &solver_dm);
983  CHKERR DMSetMatType(solver_dm, MATSHELL);
984  }
985  CHKERR KSPSetUp(solver);
986  CHKERR setDiagonalPC(pc);
987 
988  } else {
989  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
990  pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
991  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
992  pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
993  }
994  MoFEMFunctionReturn(0);
995 }
996 
997 MoFEMErrorCode SetUpSchurImpl::setEntities() {
998  MoFEMFunctionBegin;
999  auto simple = mField.getInterface<Simple>();
1000  CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
1001  SPACE_DIM, volEnts);
1002  CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
1003  subEnts);
1004  subEnts = subtract(subEnts, volEnts);
1005  MoFEMFunctionReturn(0);
1006 };
1007 
1008 MoFEMErrorCode SetUpSchurImpl::createSubDM() {
1009  MoFEMFunctionBegin;
1010  auto simple = mField.getInterface<Simple>();
1011 
1012  auto create_dm = [&](const char *name, auto &ents, auto dm_type) {
1013  auto dm = createDM(mField.get_comm(), dm_type);
1014  auto create_dm_imp = [&]() {
1015  MoFEMFunctionBegin;
1016  CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
1017  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
1018  CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
1019  auto sub_ents_ptr = boost::make_shared<Range>(ents);
1020  CHKERR DMMoFEMAddSubFieldRow(dm, "U", sub_ents_ptr);
1021  CHKERR DMMoFEMAddSubFieldCol(dm, "U", sub_ents_ptr);
1022  CHKERR DMSetUp(dm);
1023  MoFEMFunctionReturn(0);
1024  };
1025  CHK_THROW_MESSAGE(create_dm_imp(),
1026  "Error in creating schurDM. It is possible that schurDM is "
1027  "already created");
1028  return dm;
1029  };
1030 
1031  schurDM = create_dm("SCHUR", subEnts, "DMMOFEM_MG");
1032  blockDM = create_dm("BLOCK", volEnts, "DMMOFEM");
1033 
1034  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
1035 
1036  auto get_nested_mat_data = [&]() -> boost::shared_ptr<NestSchurData> {
1037  auto block_mat_data =
1038  createBlockMatStructure(simple->getDM(),
1039 
1040  {{
1041 
1042  simple->getDomainFEName(),
1043 
1044  {{"U", "U"}
1045 
1046  }}}
1047 
1048  );
1049 
1050  return createSchurNestedMatrixStruture(
1051 
1052  {schurDM, blockDM}, block_mat_data,
1053 
1054  {"U"}, {boost::make_shared<Range>(volEnts)}
1055 
1056  );
1057  };
1058 
1059  auto nested_mat_data = get_nested_mat_data();
1060 
1061  CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1062  }
1063 
1064  MoFEMFunctionReturn(0);
1065 }
1066 
1067 MoFEMErrorCode SetUpSchurImpl::setOperator() {
1068  MoFEMFunctionBegin;
1069  auto simple = mField.getInterface<Simple>();
1070  auto pip = mField.getInterface<PipelineManager>();
1071 
1072  // Boundary
1073  auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
1074  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
1075  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1076  pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1077  {"U"}, {boost::make_shared<Range>(volEnts)}, ao_up, S, true, true));
1078  // Domain
1079  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1080  pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1081  {"U"}, {boost::make_shared<Range>(volEnts)}, ao_up, S, true, true));
1082 
1083  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1084  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1085 
1086  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1087  MoFEMFunctionBegin;
1088  if (S) {
1089  CHKERR MatZeroEntries(S);
1090  }
1091  MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble Begin";
1092  MoFEMFunctionReturn(0);
1093  };
1094 
1095  post_proc_schur_lhs_ptr->postProcessHook = [this, post_proc_schur_lhs_ptr,
1096  ao_up]() {
1097  MoFEMFunctionBegin;
1098  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1099  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1100  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1101  mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1102  MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble End";
1103  MoFEMFunctionReturn(0);
1104  };
1105 
1106  auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
1107  ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
1108  ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
1109 
1110  MoFEMFunctionReturn(0);
1111 }
1112 
1113 MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
1114  MoFEMFunctionBegin;
1115  auto simple = mField.getInterface<Simple>();
1116  SmartPetscObj<IS> vol_is;
1117  mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1118  simple->getProblemName(), ROW, "U", 0, SPACE_DIM, vol_is, &volEnts);
1119  CHKERR PCFieldSplitSetIS(pc, NULL, vol_is);
1120  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1121  MoFEMFunctionReturn(0);
1122 }
1123 
1124 MoFEMErrorCode SetUpSchurImpl::setDiagonalPC(PC pc) {
1125  MoFEMFunctionBegin;
1126 
1127  auto get_pc = [](auto ksp) {
1128  PC pc_raw;
1129  CHKERR KSPGetPC(ksp, &pc_raw);
1130  return SmartPetscObj<PC>(pc_raw, true); // bump reference
1131  };
1132 
1133  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
1134  auto simple = mField.getInterface<Simple>();
1135  auto A = createDMBlockMat(simple->getDM());
1136  auto P = createDMNestSchurMat(simple->getDM());
1137  CHKERR PCSetOperators(pc, A, P);
1138 
1139  KSP *subksp;
1140  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
1141  CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1142 
1143  auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
1144  MoFEMFunctionBegin;
1145  CHKERR PCSetDM(pc, dm);
1146  PetscBool same = PETSC_FALSE;
1147  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
1148  if (same) {
1149  CHKERR PCMGSetUpViaApproxOrders(
1150  pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
1151  CHKERR PCSetFromOptions(pc);
1152  }
1153  MoFEMFunctionReturn(0);
1154  };
1155 
1156  CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
1157 
1158  CHKERR PetscFree(subksp);
1159  }
1160  MoFEMFunctionReturn(0);
1161 }
1162 
1163 boost::shared_ptr<SetUpSchur>
1164 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
1165  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1166 }
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
OpK
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradSymTensorGrad< BASE_DIM, SPACE_DIM, SPACE_DIM, 0 > OpK
[Define entities]
Definition: elastic.cpp:39
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
ElasticSpring.hpp
Implementation of elastic spring bc.
SetUpSchurImpl
Definition: test_broken_space.cpp:519
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1194
BASE_DIM
constexpr int BASE_DIM
Definition: elastic.cpp:15
Example::addMatBlockOps
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: elastic.cpp:114
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
shear_modulus_G
constexpr double shear_modulus_G
Definition: elastic.cpp:82
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:157
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:640
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Example::doEvalField
PetscBool doEvalField
Definition: elastic.cpp:95
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
LASTBASE
@ LASTBASE
Definition: definitions.h:69
OpK
Definition: simple_elasticity.cpp:16
Example::Example
Example(MoFEM::Interface &m_field)
Definition: elastic.cpp:88
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:532
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
Example::fieldEvalCoords
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: elastic.cpp:96
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:630
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1113
bulk_modulus_K
constexpr double bulk_modulus_K
Definition: elastic.cpp:81
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
PostProcEleByDim< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: dynamic_first_order_con_law.cpp:43
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: dynamic_first_order_con_law.cpp:52
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1225
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::EssentialPostProcRhs< MPCsType >
Specialization for DisplacementCubitBcData.
Definition: EssentialMPCsData.hpp:81
is_plane_strain
PetscBool is_plane_strain
Definition: elastic.cpp:84
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
SetUpSchurImpl::schurDM
SmartPetscObj< DM > schurDM
Definition: contact.cpp:1003
SetUpSchurImpl::createSubDM
MoFEMErrorCode createSubDM()
Definition: contact.cpp:1082
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition: DMMoFEM.hpp:1076
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
SetUpSchurImpl::volEnts
Range volEnts
Definition: elastic.cpp:949
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:462
PostProcEleByDim
Definition: dynamic_first_order_con_law.cpp:38
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
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
CalculateTraction.hpp
Calculate traction for linear problem.
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
SetUpSchurImpl::setEntities
MoFEMErrorCode setEntities()
Definition: elastic.cpp:997
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
SetUpSchurImpl::subEnts
Range subEnts
Definition: elastic.cpp:950
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
Example::vectorFieldPtr
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Definition: elastic.cpp:98
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:216
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
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:528
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:890
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:201
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:884
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
EntData
EntitiesFieldData::EntData EntData
[Define entities]
Definition: elastic.cpp:29
MoFEM::EssentialPostProcLhs< MPCsType >
Specialization for MPCsType.
Definition: EssentialMPCsData.hpp:102
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2580
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
main
int main(int argc, char *argv[])
Definition: elastic.cpp:878
A
constexpr AssemblyType A
[Define dimension]
Definition: elastic.cpp:21
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
help
static char help[]
[Check]
Definition: elastic.cpp:876
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1149
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
I
constexpr IntegrationType I
Definition: elastic.cpp:25
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: dynamic_first_order_con_law.cpp:54
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpTensorTimesSymmetricTensor
Calculate .
Definition: UserDataOperators.hpp:1909
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
BiLinearForm
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
NaturalDomainBC.hpp
Boundary conditions in domain, i.e. body forces.
young_modulus
constexpr double young_modulus
Definition: elastic.cpp:79
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2627
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: elastic.cpp:33
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< BASE_DIM, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:43
Range
DomainEleOp
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::OpSymmetrizeTensor
Definition: UserDataOperators.hpp:1974
FluidLevel.hpp
Natural boundary condition applying pressure from fluid.
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:275
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:535
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
Example::fieldEvalData
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
Definition: elastic.cpp:97
DomainBCs
[OpInternalForce]
Definition: elastic.cpp:45
approx_order
int approx_order
Definition: test_broken_space.cpp:54
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
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:766
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
NaturalBoundaryBC.hpp
Implementation of natural boundary conditions.
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:256
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: elastic.cpp:924
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:44
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
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
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
QUIET
@ QUIET
Definition: definitions.h:221
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
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:589
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< KSP >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:802
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
SetUpSchurImpl::setDiagonalPC
MoFEMErrorCode setDiagonalPC(PC pc)
Definition: contact.cpp:1233
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
ElasticExample::OpCalculateTraction
Definition: CalculateTraction.hpp:12
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
PostProcEleByDim< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: dynamic_first_order_con_law.cpp:49
poisson_ratio
constexpr double poisson_ratio
Definition: elastic.cpp:80
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
BoundaryBCs
Definition: elastic.cpp:46
F
@ F
Definition: free_surface.cpp:396
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
PlasticOps::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
Definition: PlasticOps.hpp:248
Generated by Doxygen 1.8.17 and hosted at University of Glasgow