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::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 pip = mField.getInterface<PipelineManager>();
350  auto simple = mField.getInterface<Simple>();
351  auto bc_mng = mField.getInterface<BcManager>();
352 
353  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
354  "U", 0, 0);
355  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
356  "U", 1, 1);
357  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
358  "U", 2, 2);
359  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
360  "REMOVE_ALL", "U", 0, 3);
361  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
362  simple->getProblemName(), "U");
363 
364  // adding MPCs
365  CHKERR bc_mng->addBlockDOFsToMPCs(simple->getProblemName(), "U");
366 
367  MoFEMFunctionReturn(0);
368 }
369 //! [Boundary condition]
370 
371 //! [Push operators to pipeline]
372 MoFEMErrorCode Example::assembleSystem() {
373  MoFEMFunctionBegin;
374  auto pip = mField.getInterface<PipelineManager>();
375  auto simple = mField.getInterface<Simple>();
376  auto bc_mng = mField.getInterface<BcManager>();
377 
378  //! [Integration rule]
379  auto integration_rule = [](int, int, int approx_order) {
380  return 2 * approx_order + 1;
381  };
382 
383  auto integration_rule_bc = [](int, int, int approx_order) {
384  return 2 * approx_order + 1;
385  };
386 
387  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
388  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
389  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
390  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
391  //! [Integration rule]
392 
393  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
394  pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
395  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
396  pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
397  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
398  pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
399  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
400  pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
401 
402 
403  //! [Push domain stiffness matrix to pipeline]
404  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
405 
406  // Assemble domain stiffness matrix
407  CHKERR addMatBlockOps(pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC",
408  mat_D_ptr, Sev::verbose);
409  pip->getOpDomainLhsPipeline().push_back(new OpK("U", "U", mat_D_ptr));
410  //! [Push domain stiffness matrix to pipeline]
411 
412  //! [Push Internal forces]
413  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
414  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
415  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
416  pip->getOpDomainRhsPipeline().push_back(
417  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
418  mat_grad_ptr));
419  CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
420  mat_D_ptr, Sev::inform);
421  pip->getOpDomainRhsPipeline().push_back(
422  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
423  pip->getOpDomainRhsPipeline().push_back(
424  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
425  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
426 
427  pip->getOpDomainRhsPipeline().push_back(
428  new OpInternalForce("U", mat_stress_ptr,
429  [](double, double, double) constexpr { return -1; }));
430  //! [Push Internal forces]
431 
432  //! [Push Body forces]
433  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
434  pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
435  //! [Push Body forces]
436 
437  //! [Push natural boundary conditions]
438  // Add force boundary condition
439  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
440  pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
441  // Add case for mix type of BCs
442  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
443  pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
444  //! [Push natural boundary conditions]
445  MoFEMFunctionReturn(0);
446 }
447 //! [Push operators to pipeline]
448 
449 
450 struct SetUpSchur {
451  static boost::shared_ptr<SetUpSchur>
452  createSetUpSchur(MoFEM::Interface &m_field);
453  virtual MoFEMErrorCode setUp(SmartPetscObj<KSP> solver) = 0;
454 
455 protected:
456  SetUpSchur() = default;
457 };
458 //! [Solve]
459 MoFEMErrorCode Example::solveSystem() {
460  MoFEMFunctionBegin;
461  auto simple = mField.getInterface<Simple>();
462  auto pip = mField.getInterface<PipelineManager>();
463  auto solver = pip->createKSP();
464  CHKERR KSPSetFromOptions(solver);
465 
466  auto dm = simple->getDM();
467  auto D = createDMVector(dm);
468  auto F = vectorDuplicate(D);
469 
470  auto set_essential_bc = [&]() {
471  MoFEMFunctionBegin;
472  // This is low level pushing finite elements (pipelines) to solver
473  auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
474 
475  auto pre_proc_rhs = boost::make_shared<FEMethod>();
476  auto post_proc_rhs = boost::make_shared<FEMethod>();
477  auto post_proc_lhs = boost::make_shared<FEMethod>();
478 
479  auto get_pre_proc_hook = [&]() {
480  return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
481  {});
482  };
483  pre_proc_rhs->preProcessHook = get_pre_proc_hook();
484 
485  auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
486  MoFEMFunctionBeginHot;
487 
488  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs,
489  1.)();
490  CHKERR EssentialPostProcRhs<MPCsType>(mField, post_proc_rhs, 1.)();
491  MoFEMFunctionReturnHot(0);
492  };
493 
494  auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
495  MoFEMFunctionBeginHot;
496 
497  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField, post_proc_lhs,
498  1.)();
499  CHKERR EssentialPostProcLhs<MPCsType>(mField, post_proc_lhs)();
500  MoFEMFunctionReturnHot(0);
501  };
502 
503  post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
504  post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
505 
506  ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
507  ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
508  ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
509  MoFEMFunctionReturn(0);
510  };
511 
512  auto setup_and_solve = [&]() {
513  MoFEMFunctionBegin;
514  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
515  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
516  CHKERR KSPSetUp(solver);
517  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
518  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
519  CHKERR KSPSolve(solver, F, D);
520  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
521  MoFEMFunctionReturn(0);
522  };
523 
524  MOFEM_LOG_CHANNEL("TIMER");
525  MOFEM_LOG_TAG("TIMER", "timer");
526 
527  CHKERR set_essential_bc();
528 
529  if (A == AssemblyType::SCHUR) {
530  auto schur_ptr = SetUpSchur::createSetUpSchur(mField);
531  CHKERR schur_ptr->setUp(solver);
532  CHKERR setup_and_solve();
533  } else {
534  CHKERR setup_and_solve();
535  }
536 
537  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
538  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
539  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
540 
541  if (doEvalField) {
542  if constexpr (SPACE_DIM == 3) {
543  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
544  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
545  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
546  mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
547  } else {
548  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
549  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
550  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
551  mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
552  }
553 
554  if (vectorFieldPtr->size1()) {
555  auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
556  if constexpr (SPACE_DIM == 2)
557  MOFEM_LOG("FieldEvaluator", Sev::inform)
558  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
559  else
560  MOFEM_LOG("FieldEvaluator", Sev::inform)
561  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
562  << " U_Z: " << t_disp(2);
563  }
564 
565  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
566  }
567  MoFEMFunctionReturn(0);
568 }
569 //! [Solve]
570 
571 //! [Postprocess results]
572 MoFEMErrorCode Example::outputResults() {
573  MoFEMFunctionBegin;
574  auto simple = mField.getInterface<Simple>();
575  auto pip = mField.getInterface<PipelineManager>();
576  auto det_ptr = boost::make_shared<VectorDouble>();
577  auto jac_ptr = boost::make_shared<MatrixDouble>();
578  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
579  //! [Postprocess clean]
580  pip->getDomainRhsFE().reset();
581  pip->getDomainLhsFE().reset();
582  pip->getBoundaryRhsFE().reset();
583  pip->getBoundaryLhsFE().reset();
584  //! [Postprocess clean]
585 
586  //! [Postprocess initialise]
587  auto post_proc_mesh = boost::make_shared<moab::Core>();
588  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
589  mField, post_proc_mesh);
590  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
591  mField, post_proc_mesh);
592  //! [Postprocess initialise]
593 
594  auto calculate_stress_ops = [&](auto &pip) {
595  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
596  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
597  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
598  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
599  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
600  "U", mat_grad_ptr));
601  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
602  CHKERR addMatBlockOps(pip, "U", "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
603  pip.push_back(
604  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
605  pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
606  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
607  auto u_ptr = boost::make_shared<MatrixDouble>();
608  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
609  auto x_ptr = boost::make_shared<MatrixDouble>();
610  pip.push_back(
611  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
612  return boost::make_tuple(u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr);
613  };
614 
615  auto post_proc_domain = [&](auto post_proc_mesh) {
616  auto post_proc_fe =
617  boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
618  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
619 
620  auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
621  calculate_stress_ops(post_proc_fe->getOpPtrVector());
622 
623  post_proc_fe->getOpPtrVector().push_back(
624 
625  new OpPPMap(
626 
627  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
628 
629  {},
630 
631  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
632 
633  {},
634 
635  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
636 
637  )
638 
639  );
640  return post_proc_fe;
641  };
642 
643  auto post_proc_boundary = [&](auto post_proc_mesh) {
644  auto post_proc_fe =
645  boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
646  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
647  post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
648  auto op_loop_side =
649  new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
650  // push ops to side element, through op_loop_side operator
651  auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
652  calculate_stress_ops(op_loop_side->getOpPtrVector());
653  post_proc_fe->getOpPtrVector().push_back(op_loop_side);
654  auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
655  post_proc_fe->getOpPtrVector().push_back(
656  new ElasticExample::OpCalculateTraction(mat_stress_ptr,
657  mat_traction_ptr));
658 
659  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
660 
661  post_proc_fe->getOpPtrVector().push_back(
662 
663  new OpPPMap(
664 
665  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
666 
667  {},
668 
669  {{"U", u_ptr}, {"GEOMETRY", x_ptr}, {"T", mat_traction_ptr}},
670 
671  {},
672 
673  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
674 
675  )
676 
677  );
678  return post_proc_fe;
679  };
680 
681  PetscBool post_proc_skin_only = PETSC_FALSE;
682  if (SPACE_DIM == 3) {
683  post_proc_skin_only = PETSC_TRUE;
684  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin_only",
685  &post_proc_skin_only, PETSC_NULL);
686  }
687  if (post_proc_skin_only == PETSC_FALSE) {
688  pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
689  } else {
690  pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
691  }
692 
693  CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
694  post_proc_begin->getFEMethod());
695  CHKERR pip->loopFiniteElements();
696  CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
697  post_proc_end->getFEMethod());
698 
699  CHKERR post_proc_end->writeFile("out_elastic.h5m");
700  MoFEMFunctionReturn(0);
701 }
702 //! [Postprocess results]
703 
704 //! [Check]
705 MoFEMErrorCode Example::checkResults() {
706  MOFEM_LOG_CHANNEL("WORLD");
707  auto simple = mField.getInterface<Simple>();
708  auto pip = mField.getInterface<PipelineManager>();
709  MoFEMFunctionBegin;
710  pip->getDomainRhsFE().reset();
711  pip->getDomainLhsFE().reset();
712  pip->getBoundaryRhsFE().reset();
713  pip->getBoundaryLhsFE().reset();
714 
715  auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
716  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
717  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule);
718 
719  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
720  pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
721  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
722  pip->getOpBoundaryRhsPipeline(), {}, "GEOMETRY");
723 
724  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
725  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
726  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
727 
728  pip->getOpDomainRhsPipeline().push_back(
729  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
730  mat_grad_ptr));
731  pip->getOpDomainRhsPipeline().push_back(
732  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
733 
734  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
735  CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
736  mat_D_ptr, Sev::verbose);
737  pip->getOpDomainRhsPipeline().push_back(
738  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
739  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
740 
741  pip->getOpDomainRhsPipeline().push_back(
742  new OpInternalForce("U", mat_stress_ptr));
743 
744  pip->getOpBoundaryRhsPipeline().push_back(
745  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
746  mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
747  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
748  pip->getOpDomainRhsPipeline(), mField, "U", Sev::verbose);
749  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
750  pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::verbose);
751 
752  auto dm = simple->getDM();
753  auto res = createDMVector(dm);
754  pip->getDomainRhsFE()->f = res;
755  pip->getBoundaryRhsFE()->f = res;
756 
757  CHKERR VecZeroEntries(res);
758 
759  CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
760  CHKERR pip->loopFiniteElements();
761  CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
762 
763  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
764  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
765  CHKERR VecAssemblyBegin(res);
766  CHKERR VecAssemblyEnd(res);
767 
768  auto zero_residual_at_constrains = [&]() {
769  MoFEMFunctionBegin;
770  auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
771  auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res]() {
772  MoFEMFunctionBegin;
773  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
774  mField, fe_post_proc_ptr, res)();
775  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
776  mField, fe_post_proc_ptr, 0, res)();
777  CHKERR EssentialPostProcRhs<MPCsType>(
778  mField, fe_post_proc_ptr, 0, res)();
779  MoFEMFunctionReturn(0);
780  };
781  fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
782  CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
783  MoFEMFunctionReturn(0);
784  };
785 
786  CHKERR zero_residual_at_constrains();
787 
788  double nrm2;
789  CHKERR VecNorm(res, NORM_2, &nrm2);
790  MOFEM_LOG_CHANNEL("WORLD");
791  MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
792 
793  PetscBool test = PETSC_FALSE;
794  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
795  if (test == PETSC_TRUE) {
796 
797  auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
798  MoFEMFunctionBegin;
799  auto post_proc_fe =
800  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
801  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
802  auto u_vec = boost::make_shared<MatrixDouble>();
803  post_proc_fe->getOpPtrVector().push_back(
804  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
805  post_proc_fe->getOpPtrVector().push_back(
806 
807  new OpPPMap(
808 
809  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
810 
811  {},
812 
813  {{"RES", u_vec}},
814 
815  {}, {})
816 
817  );
818 
819  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
820  post_proc_fe);
821  post_proc_fe->writeFile(out_name);
822  MoFEMFunctionReturn(0);
823  };
824 
825  CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
826 
827  constexpr double eps = 1e-8;
828  if (nrm2 > eps)
829  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
830  "Residual is not zero");
831  }
832 
833  MoFEMFunctionReturn(0);
834 }
835 //! [Check]
836 
837 static char help[] = "...\n\n";
838 
839 int main(int argc, char *argv[]) {
840 
841  // Initialisation of MoFEM/PETSc and MOAB data structures
842  const char param_file[] = "param_file.petsc";
843  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
844 
845  auto core_log = logging::core::get();
846  core_log->add_sink(
847  LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
848 
849  core_log->add_sink(
850  LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
851  LogManager::setLog("FieldEvaluator");
852  MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
853 
854  try {
855 
856  //! [Register MoFEM discrete manager in PETSc]
857  DMType dm_name = "DMMOFEM";
858  CHKERR DMRegister_MoFEM(dm_name);
859  //! [Register MoFEM discrete manager in PETSc
860 
861  //! [Create MoAB]
862  moab::Core mb_instance; ///< mesh database
863  moab::Interface &moab = mb_instance; ///< mesh database interface
864  //! [Create MoAB]
865 
866  //! [Create MoFEM]
867  MoFEM::Core core(moab); ///< finite element database
868  MoFEM::Interface &m_field = core; ///< finite element database interface
869  //! [Create MoFEM]
870 
871  //! [Example]
872  Example ex(m_field);
873  CHKERR ex.runProblem();
874  //! [Example]
875  }
876  CATCH_ERRORS;
877 
878  CHKERR MoFEM::Core::Finalize();
879 }
880 
881 struct SetUpSchurImpl : public SetUpSchur {
882 
883  SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {
884  if (S) {
885  CHK_THROW_MESSAGE(
886  MOFEM_DATA_INCONSISTENCY,
887  "Is expected that schur matrix is not allocated. This is "
888  "possible only is PC is set up twice");
889  }
890  }
891  virtual ~SetUpSchurImpl() { S.reset(); }
892 
893  MoFEMErrorCode setUp(SmartPetscObj<KSP> solver);
894 
895 private:
896  MoFEMErrorCode setEntities();
897  MoFEMErrorCode setUpSubDM();
898  MoFEMErrorCode setOperator();
899  MoFEMErrorCode setPC(PC pc);
900 
901  SmartPetscObj<Mat> S;
902 
903  MoFEM::Interface &mField;
904 
905  SmartPetscObj<DM> subDM;
906  Range volEnts;
907  Range subEnts;
908 };
909 
910 MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<KSP> solver) {
911  MoFEMFunctionBegin;
912  auto pip = mField.getInterface<PipelineManager>();
913  PC pc;
914  CHKERR KSPGetPC(solver, &pc);
915  PetscBool is_pcfs = PETSC_FALSE;
916  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
917  if (is_pcfs) {
918  if (S) {
919  CHK_THROW_MESSAGE(
920  MOFEM_DATA_INCONSISTENCY,
921  "Is expected that schur matrix is not allocated. This is "
922  "possible only is PC is set up twice");
923  }
924  CHKERR setEntities();
925  CHKERR setUpSubDM();
926  S = createDMMatrix(subDM);
927  CHKERR MatSetBlockSize(S, SPACE_DIM);
928  CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
929  CHKERR setOperator();
930  CHKERR setPC(pc);
931  } else {
932  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
933  pip->getOpBoundaryLhsPipeline().push_back(
934  new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
935  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
936  pip->getOpDomainLhsPipeline().push_back(
937  new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
938  }
939  MoFEMFunctionReturn(0);
940 }
941 
942 MoFEMErrorCode SetUpSchurImpl::setEntities() {
943  MoFEMFunctionBegin;
944  auto simple = mField.getInterface<Simple>();
945  CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
946  SPACE_DIM, volEnts);
947  CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
948  subEnts);
949  subEnts = subtract(subEnts, volEnts);
950  MoFEMFunctionReturn(0);
951 };
952 
953 MoFEMErrorCode SetUpSchurImpl::setUpSubDM() {
954  MoFEMFunctionBegin;
955  auto simple = mField.getInterface<Simple>();
956  subDM = createDM(mField.get_comm(), "DMMOFEM");
957  CHKERR DMMoFEMCreateSubDM(subDM, simple->getDM(), "SUB");
958  CHKERR DMMoFEMSetSquareProblem(subDM, PETSC_TRUE);
959  CHKERR DMMoFEMAddElement(subDM, simple->getDomainFEName());
960  auto sub_ents_ptr = boost::make_shared<Range>(subEnts);
961  CHKERR DMMoFEMAddSubFieldRow(subDM, "U", sub_ents_ptr);
962  CHKERR DMMoFEMAddSubFieldCol(subDM, "U", sub_ents_ptr);
963  CHKERR DMSetUp(subDM);
964  MoFEMFunctionReturn(0);
965 }
966 
967 MoFEMErrorCode SetUpSchurImpl::setOperator() {
968  MoFEMFunctionBegin;
969  auto pip = mField.getInterface<PipelineManager>();
970 
971  // Boundary
972  auto dm_is = getDMSubData(subDM)->getSmartRowIs();
973  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
974  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
975  pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
976  {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
977  // Domain
978  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
979  pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
980  {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
981 
982  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
983  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
984 
985  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
986  MoFEMFunctionBegin;
987  if (S) {
988  CHKERR MatZeroEntries(S);
989  }
990  MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble Begin";
991  MoFEMFunctionReturn(0);
992  };
993 
994  post_proc_schur_lhs_ptr->postProcessHook = [this, post_proc_schur_lhs_ptr,
995  ao_up]() {
996  MoFEMFunctionBegin;
997  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
998  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
999  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1000  mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1001  MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble End";
1002  MoFEMFunctionReturn(0);
1003  };
1004 
1005  auto simple = mField.getInterface<Simple>();
1006  auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
1007 
1008  ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
1009  ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
1010 
1011  MoFEMFunctionReturn(0);
1012 }
1013 
1014 MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
1015  MoFEMFunctionBegin;
1016  auto simple = mField.getInterface<Simple>();
1017  SmartPetscObj<IS> vol_is;
1018  mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1019  simple->getProblemName(), ROW, "U", 0, SPACE_DIM, vol_is, &volEnts);
1020  CHKERR PCFieldSplitSetIS(pc, NULL, vol_is);
1021  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1022  MoFEMFunctionReturn(0);
1023 }
1024 
1025 boost::shared_ptr<SetUpSchur>
1026 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
1027  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1028 }
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: adolc_plasticity.cpp:97
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:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce 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:302
ElasticSpring.hpp
Implementation of elastic spring bc.
SetUpSchurImpl
Definition: SetUpSchurImpl.cpp:9
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
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
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(KSP solver)
Definition: SetUpSchurImpl.cpp:39
shear_modulus_G
constexpr double shear_modulus_G
Definition: elastic.cpp:82
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:151
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
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:284
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
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:130
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:460
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1032
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:466
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: elastic.cpp:891
PostProcEleByDim< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: adolc_plasticity.cpp:88
SetUpSchurImpl::setUpSubDM
MoFEMErrorCode setUpSubDM()
Definition: elastic.cpp:953
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1189
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:111
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:104
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:527
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
SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: adolc_plasticity.cpp:98
MoFEM::EssentialPostProcRhs< MPCsType >
Specialization for DisplacementCubitBcData.
Definition: EssentialMPCsData.hpp:80
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:1076
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1003
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
order
constexpr int order
Definition: dg_projection.cpp:18
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:906
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: adolc_plasticity.cpp:99
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
PostProcEleByDim
Definition: adolc_plasticity.cpp:82
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
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::OpSchurAssembleBegin
Clear Schur complement internal data.
Definition: Schur.hpp:22
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:137
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
SetUpSchurImpl::setEntities
MoFEMErrorCode setEntities()
Definition: elastic.cpp:942
MoFEM::OpSchurAssembleEnd< SCHUR_DSYSV >
Definition: Schur.hpp:107
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
SetUpSchurImpl::subEnts
Range subEnts
Definition: elastic.cpp:907
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:228
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:219
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:880
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
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:47
EntData
EntitiesFieldData::EntData EntData
[Define entities]
Definition: elastic.cpp:29
MoFEM::EssentialPostProcLhs< MPCsType >
Specialization for MPCsType.
Definition: EssentialMPCsData.hpp:99
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128
main
int main(int argc, char *argv[])
Definition: elastic.cpp:839
A
constexpr AssemblyType A
[Define dimension]
Definition: elastic.cpp:21
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
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:550
help
static char help[]
[Check]
Definition: elastic.cpp:837
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1101
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:560
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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpTensorTimesSymmetricTensor
Calculate .
Definition: UserDataOperators.hpp:1885
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
MoFEM::PostProcBrokenMeshInMoabBase
Definition: PostProcBrokenMeshInMoabBase.hpp:94
SetUpSchur::SetUpSchur
SetUpSchur()=default
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:128
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:503
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: elastic.cpp:33
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
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:217
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
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:1949
FluidLevel.hpp
Natural boundary condition applying pressure from fluid.
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:276
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
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
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:752
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:264
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: elastic.cpp:883
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
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:575
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< KSP >
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:590
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
ElasticExample::OpCalculateTraction
Definition: CalculateTraction.hpp:12
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:104
PostProcEleByDim< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: adolc_plasticity.cpp:94
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:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
BoundaryBCs
Definition: elastic.cpp:46
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
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