v0.14.0
incompressible_elasticity.cpp
Go to the documentation of this file.
1 /**
2  * @file incompressible_elasticity.cpp
3  * @brief Incompressible elasticity problem
4  */
5 
6 #include <MoFEM.hpp>
7 
8 using namespace MoFEM;
9 
10 constexpr int SPACE_DIM =
11  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
12 
13 constexpr AssemblyType AT =
15  : AssemblyType::PETSC; //< selected assembly type
16 
17 constexpr IntegrationType IT =
18  IntegrationType::GAUSS; //< selected integration type
20 
24 using BoundaryEle =
30 
31 // struct DomainBCs {};
32 // struct BoundaryBCs {};
33 
34 // using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
35 // using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
36 // using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
37 // using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
38 // using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
39 // using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
40 
41 PetscBool doEvalField;
43 
46  std::pair<boost::shared_ptr<PostProcEle>,
47  boost::shared_ptr<SkinPostProcEle>>
48  pair_post_proc_fe,
49  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
50  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
51  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter,
52  std::array<double, SPACE_DIM> pass_field_eval_coords,
53  boost::shared_ptr<SetPtsData> pass_field_eval_data,
54  boost::shared_ptr<MatrixDouble> vec_field_ptr)
55  : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
56  uZScatter(uz_scatter), fieldEvalCoords(pass_field_eval_coords),
57  fieldEvalData(pass_field_eval_data), vecFieldPtr(vec_field_ptr) {
58  postProcFe = pair_post_proc_fe.first;
59  skinPostProcFe = pair_post_proc_fe.second;
60  };
61 
62  MoFEMErrorCode preProcess() { return 0; }
63  MoFEMErrorCode operator()() { return 0; }
64 
67 
68  MoFEM::Interface *m_field_ptr;
69  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
70  auto *simple = m_field_ptr->getInterface<Simple>();
71 
72  if (doEvalField) {
73  if (SPACE_DIM == 3) {
75  ->evalFEAtThePoint3D(
76  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
77  simple->getDomainFEName(), fieldEvalData,
78  m_field_ptr->get_comm_rank(), m_field_ptr->get_comm_rank(),
79  getCacheWeakPtr().lock(), MF_EXIST, QUIET);
80  } else {
82  ->evalFEAtThePoint2D(
83  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
84  simple->getDomainFEName(), fieldEvalData,
85  m_field_ptr->get_comm_rank(), m_field_ptr->get_comm_rank(),
86  getCacheWeakPtr().lock(), MF_EXIST, QUIET);
87  }
88 
89  if (vecFieldPtr->size1()) {
90  auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vecFieldPtr);
91  if constexpr (SPACE_DIM == 2)
92  MOFEM_LOG("SYNC", Sev::inform)
93  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
94  else
95  MOFEM_LOG("SYNC", Sev::inform)
96  << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
97  << " U_Z: " << t_disp(2);
98  }
99  }
100  MOFEM_LOG_SEVERITY_SYNC(m_field_ptr->get_comm(), Sev::inform);
101 
102  auto make_vtk = [&]() {
104  if (postProcFe) {
105  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcFe,
106  getCacheWeakPtr());
107  CHKERR postProcFe->writeFile("out_incomp_elasticity" +
108  boost::lexical_cast<std::string>(ts_step) +
109  ".h5m");
110  }
111  if (skinPostProcFe) {
112  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", skinPostProcFe,
113  getCacheWeakPtr());
114  CHKERR skinPostProcFe->writeFile(
115  "out_skin_incomp_elasticity_" +
116  boost::lexical_cast<std::string>(ts_step) + ".h5m");
117  }
119  };
120 
121  auto print_max_min = [&](auto &tuple, const std::string msg) {
123  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
124  INSERT_VALUES, SCATTER_FORWARD);
125  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
126  INSERT_VALUES, SCATTER_FORWARD);
127  double max, min;
128  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
129  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
130  MOFEM_LOG_C("INCOMP_ELASTICITY", Sev::inform,
131  "%s time %3.4e min %3.4e max %3.4e", msg.c_str(), ts_t, min,
132  max);
134  };
135 
136  CHKERR make_vtk();
137  CHKERR print_max_min(uXScatter, "Ux");
138  CHKERR print_max_min(uYScatter, "Uy");
139  if constexpr (SPACE_DIM == 3)
140  CHKERR print_max_min(uZScatter, "Uz");
141 
143  }
144 
145 private:
147  boost::shared_ptr<PostProcEle> postProcFe;
148  boost::shared_ptr<SkinPostProcEle> skinPostProcFe;
149  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
150  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
151  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
152  std::array<double, SPACE_DIM> fieldEvalCoords;
153  boost::shared_ptr<SetPtsData> fieldEvalData;
154  boost::shared_ptr<MatrixDouble> vecFieldPtr;
155 };
156 
157 // Assemble to A matrix, by default, however, some terms are assembled only to
158 // preconditioning.
159 
160 template <>
164  const EntitiesFieldData::EntData &row_data,
165  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
166  return MatSetValues<AssemblyTypeSelector<AT>>(
167  op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
168  };
169 
170 /**
171  * @brief Element used to specialise assembly
172  *
173  */
174 struct DomainEleOpStab : public DomainEleOp {
176 };
177 
178 /**
179  * @brief Specialise assembly for Stabilised matrix
180  *
181  * @tparam
182  */
183 template <>
187  const EntitiesFieldData::EntData &row_data,
188  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
189  return MatSetValues<AssemblyTypeSelector<AT>>(
190  op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
191  };
192 //! [Specialisation for assembly]
193 
194 int order = 2;
195 int geom_order = 1;
196 inline static double young_modulus = 100;
197 inline static double poisson_ratio = 0.25;
198 inline static double mu;
199 inline static double lambda;
200 
201 PetscBool isDiscontinuousPressure = PETSC_FALSE;
202 
204 
205  Incompressible(MoFEM::Interface &m_field) : mField(m_field) {}
206 
207  MoFEMErrorCode runProblem();
208 
209 private:
211 
212  MoFEMErrorCode setupProblem();
213  MoFEMErrorCode createCommonData();
214  MoFEMErrorCode bC();
215  MoFEMErrorCode OPs();
216  MoFEMErrorCode tsSolve();
217  MoFEMErrorCode checkResults();
218 
219  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
220  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
221  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
222 };
223 
224 template <int DIM>
226  OpCalculateLameStress(double m_u, boost::shared_ptr<MatrixDouble> stress_ptr,
227  boost::shared_ptr<MatrixDouble> strain_ptr,
228  boost::shared_ptr<VectorDouble> pressure_ptr)
229  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST), mU(m_u),
230  stressPtr(stress_ptr), strainPtr(strain_ptr),
231  pressurePtr(pressure_ptr) {}
232 
233  MoFEMErrorCode doWork(int side, EntityType type,
236  // Define Indicies
239 
240  // Define Kronecker Delta
242 
243  // Number of Gauss points
244  const size_t nb_gauss_pts = getGaussPts().size2();
245 
246  stressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
247  auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(stressPtr));
248  auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(strainPtr));
249  auto t_pressure = getFTensor0FromVec(*(pressurePtr));
250 
251  const double l_mu = mU;
252  // Extract matrix from data matrix
253  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
254 
255  t_stress(i, j) = t_pressure * t_kd(i, j) + 2. * l_mu * t_strain(i, j);
256 
257  ++t_strain;
258  ++t_stress;
259  ++t_pressure;
260  }
261 
263  }
264 
265 private:
266  double mU;
267  boost::shared_ptr<MatrixDouble> stressPtr;
268  boost::shared_ptr<MatrixDouble> strainPtr;
269  boost::shared_ptr<VectorDouble> pressurePtr;
270 };
271 
272 //! [Run problem]
275  CHKERR setupProblem();
276  CHKERR createCommonData();
277  CHKERR bC();
278  CHKERR OPs();
279  CHKERR tsSolve();
280  CHKERR checkResults();
282 }
283 //! [Run problem]
284 
285 //! [Set up problem]
288  Simple *simple = mField.getInterface<Simple>();
289 
290  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
291  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
292  PETSC_NULL);
293  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_discontinuous_pressure",
294  &isDiscontinuousPressure, PETSC_NULL);
295 
296  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
297  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
298 
299  // Select base
300  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
301  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
302  PetscInt choice_base_value = AINSWORTH;
303  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
304  LASBASETOPT, &choice_base_value, PETSC_NULL);
305 
307  switch (choice_base_value) {
308  case AINSWORTH:
310  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
311  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
312  break;
313  case DEMKOWICZ:
314  base = DEMKOWICZ_JACOBI_BASE;
315  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
316  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
317  break;
318  default:
319  base = LASTBASE;
320  break;
321  }
322 
323  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
324  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
325  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
326  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
327 
328  // Adding fields related to incompressible elasticity
329  // Add displacement domain and boundary fields
330  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
331  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
332  CHKERR simple->setFieldOrder("U", order);
333 
334  // Add pressure domain and boundary fields
335  // Choose either Crouzeix-Raviart element:
337  CHKERR simple->addDomainField("P", L2, base, 1);
338  CHKERR simple->setFieldOrder("P", order - 2);
339  } else {
340  // ... or Taylor-Hood element:
341  CHKERR simple->addDomainField("P", H1, base, 1);
342  CHKERR simple->setFieldOrder("P", order - 1);
343  }
344 
345  // Add geometry data field
346  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
347  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
348 
349  CHKERR simple->setUp();
350 
351  auto project_ho_geometry = [&]() {
352  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
353  return mField.loop_dofs("GEOMETRY", ent_method);
354  };
355  CHKERR project_ho_geometry();
356 
358 } //! [Set up problem]
359 
360 //! [Create common data]
363 
364  auto get_options = [&]() {
366  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
367  &young_modulus, PETSC_NULL);
368  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
369  &poisson_ratio, PETSC_NULL);
370 
371  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
372  << "Young modulus " << young_modulus;
373  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
374  << "Poisson_ratio " << poisson_ratio;
375 
376  mu = young_modulus / (2. + 2. * poisson_ratio);
377  const double lambda_denom =
378  (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
379  lambda = young_modulus * poisson_ratio / lambda_denom;
380 
382  };
383 
384  CHKERR get_options();
385 
387 }
388 //! [Create common data]
389 
390 //! [Boundary condition]
393 
394  auto pipeline_mng = mField.getInterface<PipelineManager>();
395  auto simple = mField.getInterface<Simple>();
396  auto bc_mng = mField.getInterface<BcManager>();
397  auto time_scale = boost::make_shared<TimeScale>();
398 
399  auto integration_rule = [](int, int, int approx_order) {
400  return 2 * (approx_order - 1);
401  };
402 
403  using DomainNaturalBC =
404  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
405  using OpBodyForce =
407  using BoundaryNaturalBC =
408  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
410 
411  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
413  pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
414  //! [Natural boundary condition]
416  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
417  "FORCE", Sev::inform);
418 
419  //! [Define gravity vector]
421  pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
422  "BODY_FORCE", Sev::inform);
423 
424  // Essential BC
425  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
426  "U", 0, 0);
427  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
428  "U", 1, 1);
429  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
430  "U", 2, 2);
431  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
432  simple->getProblemName(), "U");
433 
435 }
436 //! [Boundary condition]
437 
438 //! [Push operators to pip]
441  auto pip_mng = mField.getInterface<PipelineManager>();
442 
443  auto integration_rule_vol = [](int, int, int approx_order) {
444  return 2 * approx_order + geom_order - 1;
445  };
446 
447  auto add_domain_base_ops = [&](auto &pip) {
450  "GEOMETRY");
452  };
453 
454  auto add_domain_ops_lhs = [&](auto &pip) {
456 
457  // This assemble A-matrix
458  using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
460  // This assemble B-matrix (preconditioned)
461  using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
463  //! [Operators used for RHS incompressible elasticity]
464 
465  //! [Operators used for incompressible elasticity]
466  using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
470  //! [Operators used for incompressible elasticity]
471 
472  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
473  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
474  mat_D_ptr->resize(size_symm * size_symm, 1);
475 
480 
482  auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
483  t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
484 
485  pip.push_back(new OpMixScalarTimesDiv(
486  "P", "U",
487  [](const double, const double, const double) constexpr { return -1.; },
488  true, false));
489  pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
490 
491  auto get_lambda_reciprocal = [](const double, const double, const double) {
492  return 1. / lambda;
493  };
494  if (poisson_ratio < 0.5)
495  pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
496 
497  double eps_stab = 0;
498  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
499  PETSC_NULL);
500  if (eps_stab > 0)
501  pip.push_back(new OpMassPressureStab(
502  "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
503 
505  };
506 
507  auto add_domain_ops_rhs = [&](auto &pip) {
509 
510  //! [Operators used for RHS incompressible elasticity]
511  using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
512  AT>::LinearForm<GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
513 
514  using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
515  AT>::LinearForm<GAUSS>::OpMixDivTimesU<1, SPACE_DIM, SPACE_DIM>;
516 
517  using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
518  AT>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
519 
520  //! [Operators used for RHS incompressible elasticity]
521 
522  auto pressure_ptr = boost::make_shared<VectorDouble>();
523  pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
524 
525  auto div_u_ptr = boost::make_shared<VectorDouble>();
526  pip.push_back(
528 
529  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
531  "U", grad_u_ptr));
532 
533  auto strain_ptr = boost::make_shared<MatrixDouble>();
534  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
535 
536  auto get_four_mu = [](const double, const double, const double) {
537  return -2. * mu;
538  };
539  auto minus_one = [](const double, const double, const double) constexpr {
540  return -1.;
541  };
542 
543  pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
544 
545  pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
546 
547  pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
548 
549  auto get_lambda_reciprocal = [](const double, const double, const double) {
550  return 1. / lambda;
551  };
552  if (poisson_ratio < 0.5) {
553  pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
554  get_lambda_reciprocal));
555  }
556 
558  };
559 
560  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
561  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
562  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
563  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
564 
565  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
566  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
567 
569 }
570 //! [Push operators to pip]
571 
572 //! [Solve]
573 struct SetUpSchur {
574  static boost::shared_ptr<SetUpSchur>
575  createSetUpSchur(MoFEM::Interface &m_field);
576  virtual MoFEMErrorCode setUp(SmartPetscObj<TS> solver) = 0;
577 
578 protected:
579  SetUpSchur() = default;
580 };
581 
584 
585  auto simple = mField.getInterface<Simple>();
586  auto pip_mng = mField.getInterface<PipelineManager>();
587  auto is_manager = mField.getInterface<ISManager>();
588 
589  auto set_section_monitor = [&](auto solver) {
591  SNES snes;
592  CHKERR TSGetSNES(solver, &snes);
593  PetscViewerAndFormat *vf;
594  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
595  PETSC_VIEWER_DEFAULT, &vf);
596  CHKERR SNESMonitorSet(
597  snes,
598  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
599  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
601  };
602 
603  auto scatter_create = [&](auto D, auto coeff) {
605  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
606  ROW, "U", coeff, coeff, is);
607  int loc_size;
608  CHKERR ISGetLocalSize(is, &loc_size);
609  Vec v;
610  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
611  VecScatter scatter;
612  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
613  return std::make_tuple(SmartPetscObj<Vec>(v),
614  SmartPetscObj<VecScatter>(scatter));
615  };
616 
617  auto create_post_process_elements = [&]() {
618  auto pp_fe = boost::make_shared<PostProcEle>(mField);
619 
620  auto push_vol_ops = [this](auto &pip) {
623  "GEOMETRY");
624 
626  };
627 
628  auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
630 
631  auto &pip = pp_fe->getOpPtrVector();
632 
634 
635  auto x_ptr = boost::make_shared<MatrixDouble>();
636  pip.push_back(
637  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
638  auto u_ptr = boost::make_shared<MatrixDouble>();
639  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
640 
641  auto pressure_ptr = boost::make_shared<VectorDouble>();
642  pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
643 
644  auto div_u_ptr = boost::make_shared<VectorDouble>();
646  "U", div_u_ptr));
647 
648  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
650  "U", grad_u_ptr));
651 
652  auto strain_ptr = boost::make_shared<MatrixDouble>();
653  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
654 
655  auto stress_ptr = boost::make_shared<MatrixDouble>();
656  pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
657  mu, stress_ptr, strain_ptr, pressure_ptr));
658 
659  pip.push_back(
660 
661  new OpPPMap(
662 
663  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
664 
665  {{"P", pressure_ptr}},
666 
667  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
668 
669  {},
670 
671  {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
672 
673  )
674 
675  );
676 
678  };
679 
680  auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
681  PetscBool post_proc_vol = PETSC_FALSE;
682  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
683  &post_proc_vol, PETSC_NULL);
684  if (post_proc_vol == PETSC_FALSE)
685  return boost::shared_ptr<PostProcEle>();
686  auto pp_fe = boost::make_shared<PostProcEle>(mField);
688  push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
689  "push_vol_post_proc_ops");
690  return pp_fe;
691  };
692 
693  auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
694  PetscBool post_proc_skin = PETSC_TRUE;
695  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
696  &post_proc_skin, PETSC_NULL);
697  if (post_proc_skin == PETSC_FALSE)
698  return boost::shared_ptr<SkinPostProcEle>();
699 
700  auto simple = mField.getInterface<Simple>();
701  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
702  auto op_side = new OpLoopSide<DomainEle>(
703  mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
704  pp_fe->getOpPtrVector().push_back(op_side);
705  CHK_MOAB_THROW(push_vol_post_proc_ops(
706  pp_fe, push_vol_ops(op_side->getOpPtrVector())),
707  "push_vol_post_proc_ops");
708  return pp_fe;
709  };
710 
711  return std::make_pair(vol_post_proc(), skin_post_proc());
712  };
713 
714  boost::shared_ptr<SetPtsData> field_eval_data;
715  boost::shared_ptr<MatrixDouble> vector_field_ptr;
716 
717  std::array<double, SPACE_DIM> field_eval_coords;
718  int coords_dim = SPACE_DIM;
719  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
720  field_eval_coords.data(), &coords_dim,
721  &doEvalField);
722 
723  if (doEvalField) {
724  vector_field_ptr = boost::make_shared<MatrixDouble>();
725  field_eval_data =
726  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
727  if constexpr (SPACE_DIM == 3) {
728  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
729  field_eval_data, simple->getDomainFEName());
730  } else {
731  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
732  field_eval_data, simple->getDomainFEName());
733  }
734 
735  field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
736  auto no_rule = [](int, int, int) { return -1; };
737  auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
738  field_eval_fe_ptr->getRuleHook = no_rule;
739  field_eval_fe_ptr->getOpPtrVector().push_back(
740  new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
741  }
742 
743  auto set_time_monitor = [&](auto dm, auto solver) {
745  boost::shared_ptr<MonitorIncompressible> monitor_ptr(
747  create_post_process_elements(), uXScatter,
748  uYScatter, uZScatter, field_eval_coords,
749  field_eval_data, vector_field_ptr));
750  boost::shared_ptr<ForcesAndSourcesCore> null;
751  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
752  monitor_ptr, null, null);
754  };
755 
756  auto set_essential_bc = [&]() {
758  // This is low level pushing finite elements (pipelines) to solver
759  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
760  auto pre_proc_ptr = boost::make_shared<FEMethod>();
761  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
762 
763  // Add boundary condition scaling
764  auto time_scale = boost::make_shared<TimeScale>();
765 
766  pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
767  mField, pre_proc_ptr, {time_scale}, false);
768  post_proc_rhs_ptr->postProcessHook =
769  EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs_ptr,
770  1.);
771  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
772  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
773  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
775  };
776 
777  auto set_schur_pc = [&](auto solver) {
778  SNES snes;
779  CHKERR TSGetSNES(solver, &snes);
780  KSP ksp;
781  CHKERR SNESGetKSP(snes, &ksp);
782  PC pc;
783  CHKERR KSPGetPC(ksp, &pc);
784  PetscBool is_pcfs = PETSC_FALSE;
785  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
786  boost::shared_ptr<SetUpSchur> schur_ptr;
787  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
788 
789  auto A = createDMMatrix(simple->getDM());
790  auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
791 
793  TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
794  "set operators");
795  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
796  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
797  pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
799  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
800  CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
801  CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
803  };
804  post_proc_schur_lhs_ptr->postProcessHook = [this,
805  post_proc_schur_lhs_ptr]() {
807  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
808  *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
809  CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
810  CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
812  mField, post_proc_schur_lhs_ptr, 1.,
813  SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
814 
815  CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
816  CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
817  CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
818  SAME_NONZERO_PATTERN);
819  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
821  };
822  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
823  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
824 
825  if (is_pcfs == PETSC_TRUE) {
826  if (AT == AssemblyType::SCHUR) {
827  schur_ptr = SetUpSchur::createSetUpSchur(mField);
828  CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
829  } else {
830  auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
832  auto name_prb = simple->getProblemName();
833  SmartPetscObj<IS> is_p;
834  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
835  name_prb, ROW, "P", 0, 1, is_p);
836  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_p);
838  };
839  CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
840  "set pcfieldsplit preconditioned");
841  }
842  return boost::make_tuple(schur_ptr, A, B);
843  }
844 
845  return boost::make_tuple(schur_ptr, A, B);
846  };
847 
848  auto dm = simple->getDM();
849  auto D = createDMVector(dm);
850 
851  uXScatter = scatter_create(D, 0);
852  uYScatter = scatter_create(D, 1);
853  if (SPACE_DIM == 3)
854  uZScatter = scatter_create(D, 2);
855 
856  // Add extra finite elements to SNES solver pipelines to resolve essential
857  // boundary conditions
858  CHKERR set_essential_bc();
859 
860  auto solver = pip_mng->createTSIM();
861  CHKERR TSSetFromOptions(solver);
862 
863  CHKERR set_section_monitor(solver);
864  CHKERR set_time_monitor(dm, solver);
865  auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
866 
867  CHKERR TSSetSolution(solver, D);
868  CHKERR TSSetUp(solver);
869  CHKERR TSSolve(solver, NULL);
870 
872 }
873 //! [Solve]
874 
875 //! [Check]
877 //! [Check]
878 
879 static char help[] = "...\n\n";
880 
881 int main(int argc, char *argv[]) {
882 
883  // Initialisation of MoFEM/PETSc and MOAB data structures
884  const char param_file[] = "param_file.petsc";
885  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
886 
887  // Add logging channel for CONTACT
888  auto core_log = logging::core::get();
889  core_log->add_sink(
890  LogManager::createSink(LogManager::getStrmWorld(), "INCOMP_ELASTICITY"));
891  LogManager::setLog("INCOMP_ELASTICITY");
892  MOFEM_LOG_TAG("INCOMP_ELASTICITY", "Indent");
893 
894  try {
895 
896  //! [Register MoFEM discrete manager in PETSc]
897  DMType dm_name = "DMMOFEM";
898  CHKERR DMRegister_MoFEM(dm_name);
899  //! [Register MoFEM discrete manager in PETSc
900 
901  //! [Create MoAB]
902  moab::Core mb_instance; ///< mesh database
903  moab::Interface &moab = mb_instance; ///< mesh database interface
904  //! [Create MoAB]
905 
906  //! [Create MoFEM]
907  MoFEM::Core core(moab); ///< finite element database
908  MoFEM::Interface &m_field = core; ///< finite element database interface
909  //! [Create MoFEM]
910 
911  //! [Load mesh]
912  Simple *simple = m_field.getInterface<Simple>();
913  CHKERR simple->getOptions();
914  CHKERR simple->loadFile("");
915  //! [Load mesh]
916 
917  //! [CONTACT]
918  Incompressible ex(m_field);
919  CHKERR ex.runProblem();
920  //! [CONTACT]
921  }
922  CATCH_ERRORS;
923 
925 
926  return 0;
927 }
928 
929 struct SetUpSchurImpl : public SetUpSchur {
930 
932  virtual ~SetUpSchurImpl() { S.reset(); }
934 
935 private:
937  MoFEMErrorCode setPC(PC pc);
938 
940 
942 
944 
946 };
947 
950  auto simple = mField.getInterface<Simple>();
951  auto pip = mField.getInterface<PipelineManager>();
952  auto dm = simple->getDM();
953 
954  SNES snes;
955  CHKERR TSGetSNES(solver, &snes);
956  KSP ksp;
957  CHKERR SNESGetKSP(snes, &ksp);
958  PC pc;
959  CHKERR KSPGetPC(ksp, &pc);
960 
961  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Setup Schur pc";
962 
963  if (S) {
965  "Is expected that schur matrix is not allocated. This is "
966  "possible only is PC is set up twice");
967  }
968 
969  auto create_sub_dm = [&]() {
970  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
971  auto set_up = [&]() {
973  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
974  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
975  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
976  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
977  CHKERR DMSetUp(sub_dm);
979  };
980  CHK_THROW_MESSAGE(set_up(), "sey up dm");
981  return sub_dm;
982  };
983 
984  subDM = create_sub_dm();
986  CHKERR MatSetBlockSize(S, SPACE_DIM);
987 
988  auto dm_is = getDMSubData(subDM)->getSmartRowIs();
989  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
990  // Domain
991  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
992  pip->getOpDomainLhsPipeline().push_back(
993  createOpSchurAssembleEnd({"P"}, {nullptr}, ao_up, S, false, false));
994 
995  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
996  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
997 
998  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1000  CHKERR MatZeroEntries(S);
1002  };
1003 
1004  post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1005  post_proc_schur_lhs_ptr]() {
1007 
1008  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1009  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1011  mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1012 
1013 #ifndef NDEBUG
1014  auto print_mat_norm = [this](auto a, std::string prefix) {
1016  double nrm;
1017  CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1018  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose)
1019  << prefix << " norm = " << nrm;
1021  };
1022  CHKERR print_mat_norm(S, "S");
1023 #endif // NDEBUG
1025  };
1026 
1027  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1028  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1029  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1030 
1031  SmartPetscObj<IS> is_p;
1032  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1033  simple->getProblemName(), ROW, "P", 0, 1, is_p);
1034  CHKERR PCFieldSplitSetIS(pc, NULL, is_p);
1035  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1036 
1038 }
1039 
1040 boost::shared_ptr<SetUpSchur>
1042  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1043 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
Incompressible::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: incompressible_elasticity.cpp:219
SetUpSchurImpl
Definition: test_broken_space.cpp:519
young_modulus
static double young_modulus
Definition: incompressible_elasticity.cpp:196
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MonitorIncompressible::dM
SmartPetscObj< DM > dM
Definition: incompressible_elasticity.cpp:146
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
OpCalculateLameStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: incompressible_elasticity.cpp:233
MonitorIncompressible
Definition: incompressible_elasticity.cpp:42
poisson_ratio
static double poisson_ratio
Definition: incompressible_elasticity.cpp:197
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MonitorIncompressible::postProcFe
boost::shared_ptr< PostProcEle > postProcFe
Definition: incompressible_elasticity.cpp:147
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
LASTBASE
@ LASTBASE
Definition: definitions.h:69
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
Incompressible::bC
MoFEMErrorCode bC()
[Create common data]
Definition: incompressible_elasticity.cpp:391
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:488
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:532
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
OpCalculateLameStress::OpCalculateLameStress
OpCalculateLameStress(double m_u, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > pressure_ptr)
Definition: incompressible_elasticity.cpp:226
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: incompressible_elasticity.cpp:23
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MonitorIncompressible::MonitorIncompressible
MonitorIncompressible(SmartPetscObj< DM > dm, std::pair< boost::shared_ptr< PostProcEle >, boost::shared_ptr< SkinPostProcEle >> pair_post_proc_fe, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uz_scatter, std::array< double, SPACE_DIM > pass_field_eval_coords, boost::shared_ptr< SetPtsData > pass_field_eval_data, boost::shared_ptr< MatrixDouble > vec_field_ptr)
Definition: incompressible_elasticity.cpp:44
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
isDiscontinuousPressure
PetscBool isDiscontinuousPressure
Definition: incompressible_elasticity.cpp:201
SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: incompressible_elasticity.cpp:29
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
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: incompressible_elasticity.cpp:932
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1225
mu
static double mu
Definition: incompressible_elasticity.cpp:198
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
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MonitorIncompressible::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: incompressible_elasticity.cpp:150
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
MonitorIncompressible::fieldEvalCoords
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: incompressible_elasticity.cpp:152
CoordinateTypes
CoordinateTypes
Coodinate system.
Definition: definitions.h:127
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::createSubDM
MoFEMErrorCode createSubDM()
Definition: contact.cpp:1082
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1660
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
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:127
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
Incompressible::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: incompressible_elasticity.cpp:439
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPA
Mat getKSPA() const
Definition: ForcesAndSourcesCore.hpp:1097
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
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
OpCalculateLameStress
Definition: incompressible_elasticity.cpp:225
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:529
OpCalculateLameStress::mU
double mU
Definition: incompressible_elasticity.cpp:266
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
EntData
EntitiesFieldData::EntData EntData
Definition: incompressible_elasticity.cpp:21
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
geom_order
int geom_order
Definition: incompressible_elasticity.cpp:195
Incompressible::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: incompressible_elasticity.cpp:221
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Incompressible::mField
MoFEM::Interface & mField
Definition: incompressible_elasticity.cpp:210
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
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
Incompressible::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: incompressible_elasticity.cpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::FieldEvaluatorInterface::SetPtsData
Definition: FieldEvaluator.hpp:29
MonitorIncompressible::fieldEvalData
boost::shared_ptr< SetPtsData > fieldEvalData
Definition: incompressible_elasticity.cpp:153
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
OpCalculateLameStress::stressPtr
boost::shared_ptr< MatrixDouble > stressPtr
Definition: incompressible_elasticity.cpp:267
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
double
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:528
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
SPACE_DIM
constexpr int SPACE_DIM
Definition: incompressible_elasticity.cpp:10
Incompressible::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: incompressible_elasticity.cpp:361
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
Incompressible::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: incompressible_elasticity.cpp:286
OpCalculateLameStress::strainPtr
boost::shared_ptr< MatrixDouble > strainPtr
Definition: incompressible_elasticity.cpp:268
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2580
Incompressible::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: incompressible_elasticity.cpp:220
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MonitorIncompressible::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: incompressible_elasticity.cpp:62
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:169
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1149
Incompressible::Incompressible
Incompressible(MoFEM::Interface &m_field)
Definition: incompressible_elasticity.cpp:205
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
help
static char help[]
[Check]
Definition: incompressible_elasticity.cpp:879
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:234
order
int order
[Specialisation for assembly]
Definition: incompressible_elasticity.cpp:194
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpBaseImpl::MatSetValuesHook
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Definition: FormsIntegrators.hpp:218
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
BiLinearForm
MonitorIncompressible::vecFieldPtr
boost::shared_ptr< MatrixDouble > vecFieldPtr
Definition: incompressible_elasticity.cpp:154
FTensor::Index< 'i', SPACE_DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
Incompressible::tsSolve
MoFEMErrorCode tsSolve()
Definition: incompressible_elasticity.cpp:582
AT
constexpr AssemblyType AT
Definition: incompressible_elasticity.cpp:13
OpCalculateLameStress::pressurePtr
boost::shared_ptr< VectorDouble > pressurePtr
Definition: incompressible_elasticity.cpp:269
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
DomainEleOpStab
Element used to specialise assembly.
Definition: incompressible_elasticity.cpp:174
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MonitorIncompressible::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: incompressible_elasticity.cpp:65
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
main
int main(int argc, char *argv[])
Definition: incompressible_elasticity.cpp:881
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
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
approx_order
int approx_order
Definition: test_broken_space.cpp:54
MonitorIncompressible::skinPostProcFe
boost::shared_ptr< SkinPostProcEle > skinPostProcFe
Definition: incompressible_elasticity.cpp:148
coord_type
constexpr CoordinateTypes coord_type
Definition: incompressible_elasticity.cpp:19
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
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
MonitorIncompressible::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: incompressible_elasticity.cpp:63
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: incompressible_elasticity.cpp:931
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
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
Incompressible
Definition: incompressible_elasticity.cpp:203
QUIET
@ QUIET
Definition: definitions.h:221
MonitorIncompressible::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: incompressible_elasticity.cpp:149
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
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
IT
constexpr IntegrationType IT
Definition: incompressible_elasticity.cpp:17
MoFEM::SmartPetscObj< DM >
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1105
Incompressible::checkResults
MoFEMErrorCode checkResults()
[Solve]
Definition: incompressible_elasticity.cpp:876
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
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
MonitorIncompressible::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: incompressible_elasticity.cpp:151
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