v0.14.0
plastic.cpp
Go to the documentation of this file.
1 /**
2  * \file plastic.cpp
3  * \example plastic.cpp
4  *
5  * Plasticity in 2d and 3d
6  *
7  */
8 
9 /* The above code is a preprocessor directive in C++ that checks if the macro
10 "EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it replaces
11 the " */
12 #ifndef EXECUTABLE_DIMENSION
13 #define EXECUTABLE_DIMENSION 3
14 #endif
15 
16 // #undef ADD_CONTACT
17 
18 #include <MoFEM.hpp>
19 #include <MatrixFunction.hpp>
20 #include <IntegrationRules.hpp>
21 
22 using namespace MoFEM;
23 
24 template <int DIM> struct ElementsAndOps;
25 
26 template <> struct ElementsAndOps<2> {
27  using DomainEle = PipelineManager::FaceEle;
28  using BoundaryEle = PipelineManager::EdgeEle;
29  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
30  static constexpr FieldSpace CONTACT_SPACE = HCURL;
31 };
32 
33 template <> struct ElementsAndOps<3> {
34  using DomainEle = PipelineManager::VolEle;
35  using BoundaryEle = PipelineManager::FaceEle;
36  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
37  static constexpr FieldSpace CONTACT_SPACE = HDIV;
38 };
39 
40 constexpr int SPACE_DIM =
41  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
42 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
43 
44 constexpr AssemblyType AT =
45  (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
46  : AssemblyType::PETSC; //< selected assembly type
47 constexpr IntegrationType IT =
48  IntegrationType::GAUSS; //< selected integration type
49 
50 constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
51 constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
52 constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
53 
54 using EntData = EntitiesFieldData::EntData;
55 using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
56 using DomainEleOp = DomainEle::UserDataOperator;
57 using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
58 using BoundaryEleOp = BoundaryEle::UserDataOperator;
59 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
60 using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
61 using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
62 
63 inline double iso_hardening_exp(double tau, double b_iso) {
64  return std::exp(
65  std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
66  -b_iso * tau));
67 }
68 
69 /**
70  * Isotropic hardening
71  */
72 inline double iso_hardening(double tau, double H, double Qinf, double b_iso,
73  double sigmaY) {
74  return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;
75 }
76 
77 inline double iso_hardening_dtau(double tau, double H, double Qinf,
78  double b_iso) {
79  auto r = [&](auto tau) {
80  return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);
81  };
82  constexpr double eps = 1e-12;
83  return std::max(r(tau), eps * r(0));
84 }
85 
86 /**
87  * Kinematic hardening
88  */
89 template <typename T, int DIM>
90 inline auto
91 kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,
92  double C1_k) {
93  FTensor::Index<'i', DIM> i;
94  FTensor::Index<'j', DIM> j;
95  FTensor::Tensor2_symmetric<double, DIM> t_alpha;
96  if (C1_k < std::numeric_limits<double>::epsilon()) {
97  t_alpha(i, j) = 0;
98  return t_alpha;
99  }
100  t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
101  return t_alpha;
102 }
103 
104 template <int DIM>
105 inline auto kinematic_hardening_dplastic_strain(double C1_k) {
106  FTensor::Index<'i', DIM> i;
107  FTensor::Index<'j', DIM> j;
108  FTensor::Index<'k', DIM> k;
109  FTensor::Index<'l', DIM> l;
110  FTensor::Ddg<double, DIM, DIM> t_diff;
111  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
112  t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
113  return t_diff;
114 }
115 
116 PetscBool is_large_strains = PETSC_TRUE; ///< Large strains
117 PetscBool set_timer = PETSC_FALSE; ///< Set timer
118 
119 double scale = 1.;
120 
121 double young_modulus = 206913; ///< Young modulus
122 double poisson_ratio = 0.29; ///< Poisson ratio
123 double sigmaY = 450; ///< Yield stress
124 double H = 129; ///< Hardening
125 double visH = 0; ///< Viscous hardening
126 double zeta = 5e-2; ///< Viscous hardening
127 double Qinf = 265; ///< Saturation yield stress
128 double b_iso = 16.93; ///< Saturation exponent
129 double C1_k = 0; ///< Kinematic hardening
130 
131 double cn0 = 1;
132 double cn1 = 1;
133 
134 int order = 2; ///< Order displacement
135 int tau_order = order - 2; ///< Order of tau files
136 int ep_order = order - 1; ///< Order of ep files
137 int geom_order = 2; ///< Order if fixed.
138 
139 PetscBool is_quasi_static = PETSC_TRUE;
140 double rho = 0.0;
141 double alpha_damping = 0;
142 
143 #include <HenckyOps.hpp>
144 #include <PlasticOps.hpp>
145 #include <PlasticNaturalBCs.hpp>
146 
147 #ifdef ADD_CONTACT
148 #ifdef PYTHON_SDF
149 #include <boost/python.hpp>
150 #include <boost/python/def.hpp>
151 #include <boost/python/numpy.hpp>
152 namespace bp = boost::python;
153 namespace np = boost::python::numpy;
154 #endif
155 
156 namespace ContactOps {
157 
158 double cn_contact = 0.1;
159 
160 }; // namespace ContactOps
161 
162 #include <ContactOps.hpp>
163 #endif // ADD_CONTACT
164 
165 using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
166 using OpDomainRhsBCs =
167  DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
168 using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
169 using OpBoundaryRhsBCs =
170  BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
171 using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
172 using OpBoundaryLhsBCs =
173  BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
174 
175 using namespace PlasticOps;
176 using namespace HenckyOps;
177 struct Example {
178 
179  Example(MoFEM::Interface &m_field) : mField(m_field) {}
180 
181  MoFEMErrorCode runProblem();
182 
183 private:
184  MoFEM::Interface &mField;
185 
186  MoFEMErrorCode setupProblem();
187  MoFEMErrorCode createCommonData();
188  MoFEMErrorCode bC();
189  MoFEMErrorCode OPs();
190  MoFEMErrorCode tsSolve();
191 
192  boost::shared_ptr<DomainEle> reactionFe;
193 
194  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
195  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
196  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
197 
198  struct ScaledTimeScale : public MoFEM::TimeScale {
199  using MoFEM::TimeScale::TimeScale;
200  double getScale(const double time) {
201  return scale * MoFEM::TimeScale::getScale(time);
202  };
203  };
204 
205 #ifdef ADD_CONTACT
206 #ifdef PYTHON_SDF
207  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
208 #endif
209 #endif // ADD_CONTACT
210 };
211 
212 //! [Run problem]
213 MoFEMErrorCode Example::runProblem() {
214  MoFEMFunctionBegin;
215  CHKERR createCommonData();
216  CHKERR setupProblem();
217  CHKERR bC();
218  CHKERR OPs();
219  CHKERR tsSolve();
220  MoFEMFunctionReturn(0);
221 }
222 //! [Run problem]
223 
224 //! [Set up problem]
225 MoFEMErrorCode Example::setupProblem() {
226  MoFEMFunctionBegin;
227  Simple *simple = mField.getInterface<Simple>();
228 
229  Range domain_ents;
230  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
231  true);
232  auto get_ents_by_dim = [&](const auto dim) {
233  if (dim == SPACE_DIM) {
234  return domain_ents;
235  } else {
236  Range ents;
237  if (dim == 0)
238  CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
239  else
240  CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
241  return ents;
242  }
243  };
244 
245  auto get_base = [&]() {
246  auto domain_ents = get_ents_by_dim(SPACE_DIM);
247  if (domain_ents.empty())
248  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
249  const auto type = type_from_handle(domain_ents[0]);
250  switch (type) {
251  case MBQUAD:
252  return DEMKOWICZ_JACOBI_BASE;
253  case MBHEX:
254  return DEMKOWICZ_JACOBI_BASE;
255  case MBTRI:
256  return AINSWORTH_LEGENDRE_BASE;
257  case MBTET:
258  return AINSWORTH_LEGENDRE_BASE;
259  default:
260  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
261  }
262  return NOBASE;
263  };
264 
265  const auto base = get_base();
266  MOFEM_LOG("PLASTICITY", Sev::inform)
267  << "Base " << ApproximationBaseNames[base];
268 
269  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
270  CHKERR simple->addDomainField("EP", L2, base, size_symm);
271  CHKERR simple->addDomainField("TAU", L2, base, 1);
272  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
273 
274  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
275 
276  PetscBool order_edge = PETSC_FALSE;
277  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,
278  PETSC_NULL);
279  PetscBool order_face = PETSC_FALSE;
280  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
281  PETSC_NULL);
282  PetscBool order_volume = PETSC_FALSE;
283  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
284  PETSC_NULL);
285 
286  if (order_edge || order_face || order_volume) {
287 
288  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
289  ? "true"
290  : "false";
291  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
292  ? "true"
293  : "false";
294  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
295  ? "true"
296  : "false";
297 
298  auto ents = get_ents_by_dim(0);
299  if (order_edge)
300  ents.merge(get_ents_by_dim(1));
301  if (order_face)
302  ents.merge(get_ents_by_dim(2));
303  if (order_volume)
304  ents.merge(get_ents_by_dim(3));
305  CHKERR simple->setFieldOrder("U", order, &ents);
306  } else {
307  CHKERR simple->setFieldOrder("U", order);
308  }
309  CHKERR simple->setFieldOrder("EP", ep_order);
310  CHKERR simple->setFieldOrder("TAU", tau_order);
311 
312  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
313 
314 #ifdef ADD_CONTACT
315  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
316  SPACE_DIM);
317  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
318  SPACE_DIM);
319 
320  auto get_skin = [&]() {
321  Range body_ents;
322  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
323  Skinner skin(&mField.get_moab());
324  Range skin_ents;
325  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
326  return skin_ents;
327  };
328 
329  auto filter_blocks = [&](auto skin) {
330  bool is_contact_block = false;
331  Range contact_range;
332  for (auto m :
333  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
334 
335  (boost::format("%s(.*)") % "CONTACT").str()
336 
337  ))
338 
339  ) {
340  is_contact_block =
341  true; ///< blocs interation is collective, so that is set irrespective
342  ///< if there are entities in given rank or not in the block
343  MOFEM_LOG("CONTACT", Sev::inform)
344  << "Find contact block set: " << m->getName();
345  auto meshset = m->getMeshset();
346  Range contact_meshset_range;
347  CHKERR mField.get_moab().get_entities_by_dimension(
348  meshset, SPACE_DIM - 1, contact_meshset_range, true);
349 
350  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
351  contact_meshset_range);
352  contact_range.merge(contact_meshset_range);
353  }
354  if (is_contact_block) {
355  MOFEM_LOG("SYNC", Sev::inform)
356  << "Nb entities in contact surface: " << contact_range.size();
357  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
358  skin = intersect(skin, contact_range);
359  }
360  return skin;
361  };
362 
363  auto filter_true_skin = [&](auto skin) {
364  Range boundary_ents;
365  ParallelComm *pcomm =
366  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
367  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
368  PSTATUS_NOT, -1, &boundary_ents);
369  return boundary_ents;
370  };
371 
372  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
373  CHKERR simple->setFieldOrder("SIGMA", 0);
374  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
375 #endif
376 
377  CHKERR simple->setUp();
378  CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
379 
380  auto project_ho_geometry = [&]() {
381  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
382  return mField.loop_dofs("GEOMETRY", ent_method);
383  };
384  PetscBool project_geometry = PETSC_TRUE;
385  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
386  &project_geometry, PETSC_NULL);
387  if (project_geometry) {
388  CHKERR project_ho_geometry();
389  }
390 
391  MoFEMFunctionReturn(0);
392 }
393 //! [Set up problem]
394 
395 //! [Create common data]
396 MoFEMErrorCode Example::createCommonData() {
397  MoFEMFunctionBegin;
398 
399  auto get_command_line_parameters = [&]() {
400  MoFEMFunctionBegin;
401 
402  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
403  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
404  &young_modulus, PETSC_NULL);
405  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
406  &poisson_ratio, PETSC_NULL);
407  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
408  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
409  PETSC_NULL);
410  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
411  PETSC_NULL);
412  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
413  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
414  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
415  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
416  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
417  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);
418  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
419  &is_large_strains, PETSC_NULL);
420  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
421  PETSC_NULL);
422 
423  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
424  PetscBool tau_order_is_set; ///< true if tau order is set
425  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
426  &tau_order_is_set);
427  PetscBool ep_order_is_set; ///< true if tau order is set
428  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
429  &ep_order_is_set);
430  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
431  PETSC_NULL);
432 
433  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
434  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
435  &alpha_damping, PETSC_NULL);
436 
437  MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
438  MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
439  MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
440  MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
441  MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
442  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
443  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
444  MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
445  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
446  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
447  MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
448 
449  if (tau_order_is_set == PETSC_FALSE)
450  tau_order = order - 2;
451  if (ep_order_is_set == PETSC_FALSE)
452  ep_order = order - 1;
453 
454  MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
455  MOFEM_LOG("PLASTICITY", Sev::inform)
456  << "Ep approximation order " << ep_order;
457  MOFEM_LOG("PLASTICITY", Sev::inform)
458  << "Tau approximation order " << tau_order;
459  MOFEM_LOG("PLASTICITY", Sev::inform)
460  << "Geometry approximation order " << geom_order;
461 
462  MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
463  MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
464 
465  PetscBool is_scale = PETSC_TRUE;
466  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
467  PETSC_NULL);
468  if (is_scale) {
469  scale /= young_modulus;
470  }
471 
472  MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
473 
474 #ifdef ADD_CONTACT
475  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
476  &ContactOps::cn_contact, PETSC_NULL);
477  MOFEM_LOG("CONTACT", Sev::inform)
478  << "cn_contact " << ContactOps::cn_contact;
479 #endif // ADD_CONTACT
480 
481  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
482  &is_quasi_static, PETSC_NULL);
483  MOFEM_LOG("PLASTICITY", Sev::inform)
484  << "Is quasi static: " << (is_quasi_static ? "true" : "false");
485 
486  MoFEMFunctionReturn(0);
487  };
488 
489  CHKERR get_command_line_parameters();
490 
491 #ifdef ADD_CONTACT
492 #ifdef PYTHON_SDF
493  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
494  CHKERR sdfPythonPtr->sdfInit("sdf.py");
495  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
496 #endif
497 #endif // ADD_CONTACT
498 
499  MoFEMFunctionReturn(0);
500 }
501 //! [Create common data]
502 
503 //! [Boundary condition]
504 MoFEMErrorCode Example::bC() {
505  MoFEMFunctionBegin;
506 
507  auto simple = mField.getInterface<Simple>();
508  auto bc_mng = mField.getInterface<BcManager>();
509  auto prb_mng = mField.getInterface<ProblemsManager>();
510 
511  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
512  "U", 0, 0);
513  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
514  "U", 1, 1);
515  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
516  "U", 2, 2);
517  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
518  "REMOVE_ALL", "U", 0, 3);
519 
520 #ifdef ADD_CONTACT
521  for (auto b : {"FIX_X", "REMOVE_X"})
522  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
523  "SIGMA", 0, 0, false, true);
524  for (auto b : {"FIX_Y", "REMOVE_Y"})
525  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
526  "SIGMA", 1, 1, false, true);
527  for (auto b : {"FIX_Z", "REMOVE_Z"})
528  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
529  "SIGMA", 2, 2, false, true);
530  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
531  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
532  "SIGMA", 0, 3, false, true);
533  CHKERR bc_mng->removeBlockDOFsOnEntities(
534  simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
535 #endif
536 
537  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
538  simple->getProblemName(), "U");
539 
540  auto &bc_map = bc_mng->getBcMapByBlockName();
541  for (auto bc : bc_map)
542  MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
543 
544  MoFEMFunctionReturn(0);
545 }
546 //! [Boundary condition]
547 
548 //! [Push operators to pipeline]
549 MoFEMErrorCode Example::OPs() {
550  MoFEMFunctionBegin;
551  auto pip_mng = mField.getInterface<PipelineManager>();
552  auto simple = mField.getInterface<Simple>();
553  auto bc_mng = mField.getInterface<BcManager>();
554 
555  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
556 
557  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
558 
559  auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
560  MoFEMFunctionBegin;
561 
562  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
563  "GEOMETRY");
564  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
565 
566  // Add Natural BCs to LHS
567  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
568  pip, mField, "U", Sev::inform);
569 
570 #ifdef ADD_CONTACT
571  CHKERR
572  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
573  pip, "SIGMA", "U");
574  CHKERR
575  ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
576  mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
577  vol_rule);
578 #endif // ADD_CONTACT
579 
580  MoFEMFunctionReturn(0);
581  };
582 
583  auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
584  MoFEMFunctionBegin;
585 
586  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
587  "GEOMETRY");
588  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
589 
590  // Add Natural BCs to RHS
591  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
592  pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
593 
594 #ifdef ADD_CONTACT
595  CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
596  pip, "SIGMA", "U");
597 #endif // ADD_CONTACT
598 
599  MoFEMFunctionReturn(0);
600  };
601 
602  auto add_domain_ops_lhs = [this](auto &pip) {
603  MoFEMFunctionBegin;
604  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
605  "GEOMETRY");
606 
607  if (is_quasi_static == PETSC_FALSE) {
608 
609  //! [Only used for dynamics]
610  using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
611  GAUSS>::OpMass<1, SPACE_DIM>;
612  //! [Only used for dynamics]
613 
614  auto get_inertia_and_mass_damping = [this](const double, const double,
615  const double) {
616  auto *pip = mField.getInterface<PipelineManager>();
617  auto &fe_domain_lhs = pip->getDomainLhsFE();
618  return (rho / scale) * fe_domain_lhs->ts_aa +
619  (alpha_damping / scale) * fe_domain_lhs->ts_a;
620  };
621  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
622  }
623 
624  CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
625  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
626 
627  MoFEMFunctionReturn(0);
628  };
629 
630  auto add_domain_ops_rhs = [this](auto &pip) {
631  MoFEMFunctionBegin;
632 
633  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
634  "GEOMETRY");
635 
636  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
637  pip, mField, "U",
638  {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
639  Sev::inform);
640 
641  // only in case of dynamics
642  if (is_quasi_static == PETSC_FALSE) {
643 
644  //! [Only used for dynamics]
645  using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
646  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
647  //! [Only used for dynamics]
648 
649  auto mat_acceleration = boost::make_shared<MatrixDouble>();
650  pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
651  "U", mat_acceleration));
652  pip.push_back(
653  new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
654  return rho / scale;
655  }));
656  if (alpha_damping > 0) {
657  auto mat_velocity = boost::make_shared<MatrixDouble>();
658  pip.push_back(
659  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
660  pip.push_back(
661  new OpInertiaForce("U", mat_velocity, [](double, double, double) {
662  return alpha_damping / scale;
663  }));
664  }
665  }
666 
667  CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
668  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
669 
670 #ifdef ADD_CONTACT
671  CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
672  pip, "SIGMA", "U");
673 #endif // ADD_CONTACT
674 
675  MoFEMFunctionReturn(0);
676  };
677 
678  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
679  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
680 
681  // Boundary
682  CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
683  CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
684 
685  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
686  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
687 
688  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
689  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
690 
691  auto create_reaction_pipeline = [&](auto &pip) {
692  MoFEMFunctionBegin;
693  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
694  "GEOMETRY");
695  CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
696  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
697  MoFEMFunctionReturn(0);
698  };
699 
700  reactionFe = boost::make_shared<DomainEle>(mField);
701  reactionFe->getRuleHook = vol_rule;
702  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
703  reactionFe->postProcessHook =
704  EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
705 
706  MoFEMFunctionReturn(0);
707 }
708 //! [Push operators to pipeline]
709 
710 //! [Solve]
711 struct SetUpSchur {
712 
713  /**
714  * @brief Create data structure for handling Schur complement
715  *
716  * @param m_field
717  * @param sub_dm Schur complement sub dm
718  * @param field_split_it IS of Schur block
719  * @param ao_map AO map from sub dm to main problem
720  * @return boost::shared_ptr<SetUpSchur>
721  */
722  static boost::shared_ptr<SetUpSchur> createSetUpSchur(
723 
724  MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
725  SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
726 
727  );
728  virtual MoFEMErrorCode setUp(TS solver) = 0;
729 
730 protected:
731  SetUpSchur() = default;
732 };
733 
734 MoFEMErrorCode Example::tsSolve() {
735  MoFEMFunctionBegin;
736 
737  Simple *simple = mField.getInterface<Simple>();
738  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
739  ISManager *is_manager = mField.getInterface<ISManager>();
740 
741  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
742 
743  auto set_section_monitor = [&](auto solver) {
744  MoFEMFunctionBegin;
745  SNES snes;
746  CHKERR TSGetSNES(solver, &snes);
747  CHKERR SNESMonitorSet(snes,
748  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
749  void *))MoFEMSNESMonitorFields,
750  (void *)(snes_ctx_ptr.get()), nullptr);
751  MoFEMFunctionReturn(0);
752  };
753 
754  auto create_post_process_elements = [&]() {
755  auto pp_fe = boost::make_shared<PostProcEle>(mField);
756  auto &pip = pp_fe->getOpPtrVector();
757 
758  auto push_vol_ops = [this](auto &pip) {
759  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
760  "GEOMETRY");
761 
762  auto [common_plastic_ptr, common_henky_ptr] =
763  PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
764  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
765 
766  if (common_henky_ptr) {
767  if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
768  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
769  }
770 
771  return std::make_pair(common_plastic_ptr, common_henky_ptr);
772  };
773 
774  auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
775  MoFEMFunctionBegin;
776 
777  auto &pip = pp_fe->getOpPtrVector();
778 
779  auto [common_plastic_ptr, common_henky_ptr] = p;
780 
781  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
782 
783  auto x_ptr = boost::make_shared<MatrixDouble>();
784  pip.push_back(
785  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
786  auto u_ptr = boost::make_shared<MatrixDouble>();
787  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
788 
789  if (is_large_strains) {
790 
791  pip.push_back(
792 
793  new OpPPMap(
794 
795  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
796 
797  {{"PLASTIC_SURFACE",
798  common_plastic_ptr->getPlasticSurfacePtr()},
799  {"PLASTIC_MULTIPLIER",
800  common_plastic_ptr->getPlasticTauPtr()}},
801 
802  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
803 
804  {{"GRAD", common_henky_ptr->matGradPtr},
805  {"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
806 
807  {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
808  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
809 
810  )
811 
812  );
813 
814  } else {
815 
816  pip.push_back(
817 
818  new OpPPMap(
819 
820  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
821 
822  {{"PLASTIC_SURFACE",
823  common_plastic_ptr->getPlasticSurfacePtr()},
824  {"PLASTIC_MULTIPLIER",
825  common_plastic_ptr->getPlasticTauPtr()}},
826 
827  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
828 
829  {},
830 
831  {{"STRAIN", common_plastic_ptr->mStrainPtr},
832  {"STRESS", common_plastic_ptr->mStressPtr},
833  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
834  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
835 
836  )
837 
838  );
839  }
840 
841  MoFEMFunctionReturn(0);
842  };
843 
844  auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
845  PetscBool post_proc_vol = PETSC_FALSE;
846  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
847  &post_proc_vol, PETSC_NULL);
848  if (post_proc_vol == PETSC_FALSE)
849  return boost::shared_ptr<PostProcEle>();
850  auto pp_fe = boost::make_shared<PostProcEle>(mField);
851  CHK_MOAB_THROW(
852  push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
853  "push_vol_post_proc_ops");
854  return pp_fe;
855  };
856 
857  auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
858  PetscBool post_proc_skin = PETSC_TRUE;
859  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
860  &post_proc_skin, PETSC_NULL);
861  if (post_proc_skin == PETSC_FALSE)
862  return boost::shared_ptr<SkinPostProcEle>();
863 
864  auto simple = mField.getInterface<Simple>();
865  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
866  auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
867  SPACE_DIM, Sev::verbose);
868  pp_fe->getOpPtrVector().push_back(op_side);
869  CHK_MOAB_THROW(push_vol_post_proc_ops(
870  pp_fe, push_vol_ops(op_side->getOpPtrVector())),
871  "push_vol_post_proc_ops");
872  return pp_fe;
873  };
874 
875  return std::make_pair(vol_post_proc(), skin_post_proc());
876  };
877 
878  auto scatter_create = [&](auto D, auto coeff) {
879  SmartPetscObj<IS> is;
880  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
881  ROW, "U", coeff, coeff, is);
882  int loc_size;
883  CHKERR ISGetLocalSize(is, &loc_size);
884  Vec v;
885  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
886  VecScatter scatter;
887  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
888  return std::make_tuple(SmartPetscObj<Vec>(v),
889  SmartPetscObj<VecScatter>(scatter));
890  };
891 
892  auto set_time_monitor = [&](auto dm, auto solver) {
893  MoFEMFunctionBegin;
894  boost::shared_ptr<Monitor> monitor_ptr(
895  new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
896  uYScatter, uZScatter));
897  boost::shared_ptr<ForcesAndSourcesCore> null;
898  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
899  monitor_ptr, null, null);
900  MoFEMFunctionReturn(0);
901  };
902 
903  auto set_schur_pc = [&](auto solver,
904  boost::shared_ptr<SetUpSchur> &schur_ptr) {
905  MoFEMFunctionBeginHot;
906 
907  auto bc_mng = mField.getInterface<BcManager>();
908  auto name_prb = simple->getProblemName();
909 
910  // create sub dm for Schur complement
911  auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
912  SmartPetscObj<DM> &dm_sub) {
913  MoFEMFunctionBegin;
914  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
915  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
916  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
917  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
918  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
919  for (auto f : {"U"}) {
920  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
921  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
922  }
923  CHKERR DMSetUp(dm_sub);
924 
925  MoFEMFunctionReturn(0);
926  };
927 
928  auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
929  SmartPetscObj<DM> &dm_sub) {
930  MoFEMFunctionBegin;
931  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
932  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
933  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
934  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
935  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
936 #ifdef ADD_CONTACT
937  for (auto f : {"SIGMA", "EP", "TAU"}) {
938  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
939  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
940  }
941 #else
942  for (auto f : {"EP", "TAU"}) {
943  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
944  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
945  }
946 #endif
947  CHKERR DMSetUp(dm_sub);
948  MoFEMFunctionReturn(0);
949  };
950 
951  // Create nested (sub BC) Schur DM
952  if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
953 
954  SmartPetscObj<DM> dm_schur;
955  CHKERR create_schur_dm(simple->getDM(), dm_schur);
956  SmartPetscObj<DM> dm_block;
957  CHKERR create_block_dm(simple->getDM(), dm_block);
958 
959 #ifdef ADD_CONTACT
960 
961  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
962  auto block_mat_data = createBlockMatStructure(
963  simple->getDM(),
964 
965  {
966 
967  {simple->getDomainFEName(),
968 
969  {{"U", "U"},
970  {"SIGMA", "SIGMA"},
971  {"U", "SIGMA"},
972  {"SIGMA", "U"},
973  {"EP", "EP"},
974  {"TAU", "TAU"},
975  {"U", "EP"},
976  {"EP", "U"},
977  {"EP", "TAU"},
978  {"TAU", "EP"},
979  {"TAU", "U"}
980 
981  }},
982 
983  {simple->getBoundaryFEName(),
984 
985  {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
986 
987  }}
988 
989  }
990 
991  );
992 
993  return getNestSchurData(
994 
995  {dm_schur, dm_block}, block_mat_data,
996 
997  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
998 
999  );
1000  };
1001 
1002 #else
1003 
1004  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1005  auto block_mat_data =
1006  createBlockMatStructure(simple->getDM(),
1007 
1008  {{simple->getDomainFEName(),
1009 
1010  {{"U", "U"},
1011  {"EP", "EP"},
1012  {"TAU", "TAU"},
1013  {"U", "EP"},
1014  {"EP", "U"},
1015  {"EP", "TAU"},
1016  {"TAU", "U"},
1017  {"TAU", "EP"}
1018 
1019  }}}
1020 
1021  );
1022 
1023  return getNestSchurData(
1024 
1025  {dm_schur, dm_block}, block_mat_data,
1026 
1027  {"EP", "TAU"}, {nullptr, nullptr}, false
1028 
1029  );
1030  };
1031 
1032 #endif
1033 
1034  auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1035  CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1036 
1037  auto block_is = getDMSubData(dm_block)->getSmartRowIs();
1038  auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
1039 
1040  // Indices has to be map fro very to level, while assembling Schur
1041  // complement.
1042  schur_ptr =
1043  SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
1044  CHKERR schur_ptr->setUp(solver);
1045  }
1046 
1047  MoFEMFunctionReturnHot(0);
1048  };
1049 
1050  auto dm = simple->getDM();
1051  auto D = createDMVector(dm);
1052  auto DD = vectorDuplicate(D);
1053  uXScatter = scatter_create(D, 0);
1054  uYScatter = scatter_create(D, 1);
1055  if constexpr (SPACE_DIM == 3)
1056  uZScatter = scatter_create(D, 2);
1057 
1058  auto create_solver = [pip_mng]() {
1059  if (is_quasi_static == PETSC_TRUE)
1060  return pip_mng->createTSIM();
1061  else
1062  return pip_mng->createTSIM2();
1063  };
1064 
1065  auto solver = create_solver();
1066 
1067  auto active_pre_lhs = []() {
1068  MoFEMFunctionBegin;
1069  std::fill(PlasticOps::CommonData::activityData.begin(),
1070  PlasticOps::CommonData::activityData.end(), 0);
1071  MoFEMFunctionReturn(0);
1072  };
1073 
1074  auto active_post_lhs = [&]() {
1075  MoFEMFunctionBegin;
1076  auto get_iter = [&]() {
1077  SNES snes;
1078  CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1079  int iter;
1080  CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1081  "Can not get iter");
1082  return iter;
1083  };
1084 
1085  auto iter = get_iter();
1086  if (iter >= 0) {
1087 
1088  std::array<int, 5> activity_data;
1089  std::fill(activity_data.begin(), activity_data.end(), 0);
1090  MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1091  activity_data.data(), activity_data.size(), MPI_INT,
1092  MPI_SUM, mField.get_comm());
1093 
1094  int &active_points = activity_data[0];
1095  int &avtive_full_elems = activity_data[1];
1096  int &avtive_elems = activity_data[2];
1097  int &nb_points = activity_data[3];
1098  int &nb_elements = activity_data[4];
1099 
1100  if (nb_points) {
1101 
1102  double proc_nb_points =
1103  100 * static_cast<double>(active_points) / nb_points;
1104  double proc_nb_active =
1105  100 * static_cast<double>(avtive_elems) / nb_elements;
1106  double proc_nb_full_active = 100;
1107  if (avtive_elems)
1108  proc_nb_full_active =
1109  100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1110 
1111  MOFEM_LOG_C("PLASTICITY", Sev::inform,
1112  "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "
1113  "elements %d "
1114  "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1115  iter, nb_points, active_points, proc_nb_points,
1116  avtive_elems, proc_nb_active, avtive_full_elems,
1117  proc_nb_full_active, iter);
1118  }
1119  }
1120 
1121  MoFEMFunctionReturn(0);
1122  };
1123 
1124  auto add_active_dofs_elem = [&](auto dm) {
1125  MoFEMFunctionBegin;
1126  auto fe_pre_proc = boost::make_shared<FEMethod>();
1127  fe_pre_proc->preProcessHook = active_pre_lhs;
1128  auto fe_post_proc = boost::make_shared<FEMethod>();
1129  fe_post_proc->postProcessHook = active_post_lhs;
1130  auto ts_ctx_ptr = getDMTsCtx(dm);
1131  ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1132  ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1133  MoFEMFunctionReturn(0);
1134  };
1135 
1136  auto set_essential_bc = [&](auto dm, auto solver) {
1137  MoFEMFunctionBegin;
1138  // This is low level pushing finite elements (pipelines) to solver
1139 
1140  auto pre_proc_ptr = boost::make_shared<FEMethod>();
1141  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1142  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1143 
1144  // Add boundary condition scaling
1145  auto disp_time_scale = boost::make_shared<TimeScale>();
1146 
1147  auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1148  EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1149  {disp_time_scale}, false);
1150  return hook;
1151  };
1152  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1153 
1154  auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1155  MoFEMFunctionBegin;
1156  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
1157  mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1158  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
1159  mField, post_proc_rhs_ptr, 1.)();
1160  MoFEMFunctionReturn(0);
1161  };
1162  auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1163  return EssentialPostProcLhs<DisplacementCubitBcData>(
1164  mField, post_proc_lhs_ptr, 1.);
1165  };
1166  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1167 
1168  auto ts_ctx_ptr = getDMTsCtx(dm);
1169  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1170  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1171  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1172  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1173  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1174 
1175  MoFEMFunctionReturn(0);
1176  };
1177 
1178  if (is_quasi_static == PETSC_TRUE) {
1179  CHKERR TSSetSolution(solver, D);
1180  } else {
1181  CHKERR TS2SetSolution(solver, D, DD);
1182  }
1183 
1184  CHKERR set_section_monitor(solver);
1185  CHKERR set_time_monitor(dm, solver);
1186  CHKERR TSSetFromOptions(solver);
1187  CHKERR TSSetUp(solver);
1188 
1189  CHKERR add_active_dofs_elem(dm);
1190  boost::shared_ptr<SetUpSchur> schur_ptr;
1191  CHKERR set_schur_pc(solver, schur_ptr);
1192  CHKERR set_essential_bc(dm, solver);
1193 
1194  MOFEM_LOG_CHANNEL("TIMER");
1195  MOFEM_LOG_TAG("TIMER", "timer");
1196  if (set_timer)
1197  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1198  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1199  CHKERR TSSetUp(solver);
1200  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1201  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1202  CHKERR TSSolve(solver, NULL);
1203  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1204 
1205  MoFEMFunctionReturn(0);
1206 }
1207 //! [Solve]
1208 
1209 static char help[] = "...\n\n";
1210 
1211 int main(int argc, char *argv[]) {
1212 
1213 #ifdef ADD_CONTACT
1214 #ifdef PYTHON_SDF
1215  Py_Initialize();
1216  np::initialize();
1217 #endif
1218 #endif // ADD_CONTACT
1219 
1220  // Initialisation of MoFEM/PETSc and MOAB data structures
1221  const char param_file[] = "param_file.petsc";
1222  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1223 
1224  // Add logging channel for example
1225  auto core_log = logging::core::get();
1226  core_log->add_sink(
1227  LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
1228  core_log->add_sink(
1229  LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
1230  LogManager::setLog("PLASTICITY");
1231  MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
1232 
1233 #ifdef ADD_CONTACT
1234  core_log->add_sink(
1235  LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
1236  LogManager::setLog("CONTACT");
1237  MOFEM_LOG_TAG("CONTACT", "Contact");
1238 #endif // ADD_CONTACT
1239 
1240  try {
1241 
1242  //! [Register MoFEM discrete manager in PETSc]
1243  DMType dm_name = "DMMOFEM";
1244  CHKERR DMRegister_MoFEM(dm_name);
1245  //! [Register MoFEM discrete manager in PETSc
1246 
1247  //! [Create MoAB]
1248  moab::Core mb_instance; ///< mesh database
1249  moab::Interface &moab = mb_instance; ///< mesh database interface
1250  //! [Create MoAB]
1251 
1252  //! [Create MoFEM]
1253  MoFEM::Core core(moab); ///< finite element database
1254  MoFEM::Interface &m_field = core; ///< finite element database interface
1255  //! [Create MoFEM]
1256 
1257  //! [Load mesh]
1258  Simple *simple = m_field.getInterface<Simple>();
1259  CHKERR simple->getOptions();
1260  CHKERR simple->loadFile();
1261  //! [Load mesh]
1262 
1263  //! [Example]
1264  Example ex(m_field);
1265  CHKERR ex.runProblem();
1266  //! [Example]
1267  }
1268  CATCH_ERRORS;
1269 
1270  CHKERR MoFEM::Core::Finalize();
1271 
1272 #ifdef ADD_CONTACT
1273 #ifdef PYTHON_SDF
1274  if (Py_FinalizeEx() < 0) {
1275  exit(120);
1276  }
1277 #endif
1278 #endif // ADD_CONTACT
1279 
1280  return 0;
1281 }
1282 
1283 struct SetUpSchurImpl : public SetUpSchur {
1284 
1285  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
1286  SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
1287  : SetUpSchur(), mField(m_field), subDM(sub_dm),
1288  fieldSplitIS(field_split_is), aoSchur(ao_up) {
1289  if (S) {
1290  CHK_THROW_MESSAGE(
1291  MOFEM_DATA_INCONSISTENCY,
1292  "Is expected that schur matrix is not allocated. This is "
1293  "possible only is PC is set up twice");
1294  }
1295  }
1296  virtual ~SetUpSchurImpl() { S.reset(); }
1297 
1298  MoFEMErrorCode setUp(TS solver);
1299  MoFEMErrorCode preProc();
1300  MoFEMErrorCode postProc();
1301 
1302 private:
1303  SmartPetscObj<Mat> S;
1304 
1305  MoFEM::Interface &mField;
1306  SmartPetscObj<DM> subDM; ///< field split sub dm
1307  SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
1308  SmartPetscObj<AO> aoSchur; ///> main DM to subDM
1309 };
1310 
1311 MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
1312  MoFEMFunctionBegin;
1313  auto simple = mField.getInterface<Simple>();
1314  auto pip_mng = mField.getInterface<PipelineManager>();
1315 
1316  SNES snes;
1317  CHKERR TSGetSNES(solver, &snes);
1318  KSP ksp;
1319  CHKERR SNESGetKSP(snes, &ksp);
1320  CHKERR KSPSetFromOptions(ksp);
1321 
1322  PC pc;
1323  CHKERR KSPSetFromOptions(ksp);
1324  CHKERR KSPGetPC(ksp, &pc);
1325  PetscBool is_pcfs = PETSC_FALSE;
1326  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1327  if (is_pcfs) {
1328  if (S) {
1329  CHK_THROW_MESSAGE(
1330  MOFEM_DATA_INCONSISTENCY,
1331  "Is expected that schur matrix is not allocated. This is "
1332  "possible only is PC is set up twice");
1333  }
1334 
1335  S = createDMMatrix(subDM);
1336  CHKERR MatSetBlockSize(S, SPACE_DIM);
1337 
1338  // Set DM to use shell block matrix
1339  DM solver_dm;
1340  CHKERR TSGetDM(solver, &solver_dm);
1341  CHKERR DMSetMatType(solver_dm, MATSHELL);
1342 
1343  auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1344  auto A = createDMBlockMat(simple->getDM());
1345  auto P = createDMNestSchurMat(simple->getDM());
1346 
1347  if (is_quasi_static == PETSC_TRUE) {
1348  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1349  Mat A, Mat B, void *ctx) {
1350  return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1351  };
1352  CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1353  } else {
1354  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1355  PetscReal a, PetscReal aa, Mat A, Mat B,
1356  void *ctx) {
1357  return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1358  };
1359  CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1360  }
1361  CHKERR KSPSetOperators(ksp, A, P);
1362 
1363  auto set_ops = [&]() {
1364  MoFEMFunctionBegin;
1365  auto pip_mng = mField.getInterface<PipelineManager>();
1366 
1367  boost::shared_ptr<BlockStructure> block_data;
1368  CHKERR DMMoFEMGetBlocMatData(simple->getDM(), block_data);
1369 
1370 #ifndef ADD_CONTACT
1371  // Boundary
1372  pip_mng->getOpBoundaryLhsPipeline().push_front(
1373  createOpSchurAssembleBegin());
1374  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1375 
1376  {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoSchur},
1377  {SmartPetscObj<Mat>(), S}, {false, false}, false, block_data
1378 
1379  ));
1380  // Domain
1381  pip_mng->getOpDomainLhsPipeline().push_front(
1382  createOpSchurAssembleBegin());
1383  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1384 
1385  {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoSchur},
1386  {SmartPetscObj<Mat>(), S}, {false, false}, false, block_data
1387 
1388  ));
1389 #else
1390 
1391  double eps_stab = 1e-4;
1392  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1393  PETSC_NULL);
1394 
1395  using B = FormsIntegrators<BoundaryEleOp>::Assembly<
1396  BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
1397  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1398 
1399  // Boundary
1400  pip_mng->getOpBoundaryLhsPipeline().push_front(
1401  createOpSchurAssembleBegin());
1402  pip_mng->getOpBoundaryLhsPipeline().push_back(
1403  new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1404  return eps_stab;
1405  }));
1406  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1407 
1408  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1409  {SmartPetscObj<AO>(), SmartPetscObj<AO>(), aoSchur},
1410  {SmartPetscObj<Mat>(), SmartPetscObj<Mat>(), S},
1411  {false, false, false}, false, block_data
1412 
1413  ));
1414  // Domain
1415  pip_mng->getOpDomainLhsPipeline().push_front(
1416  createOpSchurAssembleBegin());
1417  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1418 
1419  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1420  {SmartPetscObj<AO>(), SmartPetscObj<AO>(), aoSchur},
1421  {SmartPetscObj<Mat>(), SmartPetscObj<Mat>(), S},
1422  {false, false, false}, false, block_data
1423 
1424  ));
1425 #endif // ADD_CONTACT
1426  MoFEMFunctionReturn(0);
1427  };
1428 
1429  auto set_assemble_elems = [&]() {
1430  MoFEMFunctionBegin;
1431  auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1432  schur_asmb_pre_proc->preProcessHook = [this]() {
1433  MoFEMFunctionBegin;
1434  CHKERR MatZeroEntries(S);
1435  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1436  MoFEMFunctionReturn(0);
1437  };
1438  auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1439 
1440  schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1441  MoFEMFunctionBegin;
1442  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1443 
1444  // Apply essential constrains to Schur complement
1445  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1446  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1447  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1448  mField, schur_asmb_post_proc, 1, S, aoSchur)();
1449 
1450  MoFEMFunctionReturn(0);
1451  };
1452  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1453  ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1454  ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1455  MoFEMFunctionReturn(0);
1456  };
1457 
1458  auto set_pc = [&]() {
1459  MoFEMFunctionBegin;
1460  CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1461  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1462  MoFEMFunctionReturn(0);
1463  };
1464 
1465  auto set_diagonal_pc = [&]() {
1466  MoFEMFunctionBegin;
1467  KSP *subksp;
1468  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
1469  auto get_pc = [](auto ksp) {
1470  PC pc_raw;
1471  CHKERR KSPGetPC(ksp, &pc_raw);
1472  return SmartPetscObj<PC>(pc_raw, true); // bump reference
1473  };
1474  CHKERR setSchurMatSolvePC(get_pc(subksp[0]));
1475  CHKERR PetscFree(subksp);
1476  MoFEMFunctionReturn(0);
1477  };
1478 
1479  CHKERR set_ops();
1480  CHKERR set_pc();
1481  CHKERR set_assemble_elems();
1482 
1483  CHKERR TSSetUp(solver);
1484  CHKERR KSPSetUp(ksp);
1485  CHKERR set_diagonal_pc();
1486 
1487  } else {
1488  pip_mng->getOpBoundaryLhsPipeline().push_front(
1489  createOpSchurAssembleBegin());
1490  pip_mng->getOpBoundaryLhsPipeline().push_back(
1491  createOpSchurAssembleEnd({}, {}, {}, {}, {}, false));
1492  pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1493  pip_mng->getOpDomainLhsPipeline().push_back(
1494  createOpSchurAssembleEnd({}, {}, {}, {}, {}, false));
1495  }
1496 
1497  // fieldSplitIS.reset();
1498  // aoSchur.reset();
1499  MoFEMFunctionReturn(0);
1500 }
1501 
1502 boost::shared_ptr<SetUpSchur>
1503 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
1504  SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1505  SmartPetscObj<AO> ao_up) {
1506  return boost::shared_ptr<SetUpSchur>(
1507  new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
1508 }
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::DMMoFEMGetBlocMatData
MoFEMErrorCode DMMoFEMGetBlocMatData(DM dm, boost::shared_ptr< BlockStructure > &)
Get data for block mat.
Definition: DMMoFEM.cpp:1542
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1307
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
Example::tsSolve
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:734
SetUpSchurImpl
Definition: test_broken_space.cpp:512
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
set_timer
PetscBool set_timer
Set timer.
Definition: plastic.cpp:117
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:157
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Example::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:549
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
SetUpSchurImpl::aoSchur
SmartPetscObj< AO > aoSchur
Definition: plastic.cpp:1308
Example::Example
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:179
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
ContactOps
Definition: EshelbianContact.hpp:10
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:525
rho
double rho
Definition: plastic.cpp:140
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:123
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::BLOCK_PRECONDITIONER_SCHUR
@ BLOCK_PRECONDITIONER_SCHUR
Definition: FormsIntegrators.hpp:109
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
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: plastic.cpp:1296
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
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::TsSetI2Jacobian
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition: TsCtx.cpp:511
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
iso_hardening_dtau
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:77
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
Example::ScaledTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:200
HenckyOps
Definition: HenckyOps.hpp:11
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
scale
double scale
Definition: plastic.cpp:119
HenckyOps.hpp
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
int order
Order displacement.
Definition: plastic.cpp:134
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:137
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1306
sdf.r
int r
Definition: sdf.py:8
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
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition: DMMoFEM.hpp:1076
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
iso_hardening_exp
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:63
MoFEM::PipelineManager::EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: PipelineManager.hpp:36
ElementsAndOps< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: plastic.cpp:29
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
PlasticOps.hpp
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
Example::reactionFe
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:192
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:522
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, std::vector< SmartPetscObj< AO >> sequence_of_aos, std::vector< SmartPetscObj< Mat >> sequence_of_mats, std::vector< bool > sym_schur, bool symm_op, boost::shared_ptr< BlockStructure > diag_blocks)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2771
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
main
int main(int argc, char *argv[])
Definition: plastic.cpp:1211
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
AT
constexpr AssemblyType AT
Definition: plastic.cpp:44
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::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
kinematic_hardening
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:91
SPACE_DIM
constexpr int SPACE_DIM
Definition: plastic.cpp:40
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
Definition: SetUpSchurImpl.cpp:202
is_quasi_static
PetscBool is_quasi_static
Definition: plastic.cpp:139
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:177
help
static char help[]
[Solve]
Definition: plastic.cpp:1209
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::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:312
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:521
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
cn1
double cn1
Definition: plastic.cpp:132
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
C1_k
double C1_k
Kinematic hardening.
Definition: plastic.cpp:129
ElementsAndOps< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: plastic.cpp:36
kinematic_hardening_dplastic_strain
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:105
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2766
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
cn0
double cn0
Definition: plastic.cpp:131
MatrixFunction.hpp
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
Example::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:196
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1876
visH
double visH
Viscous hardening.
Definition: plastic.cpp:125
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
EntData
EntitiesFieldData::EntData EntData
Definition: plastic.cpp:54
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
IT
constexpr IntegrationType IT
Definition: plastic.cpp:47
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
ContactOps.hpp
FTensor::Index< 'i', DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::setSchurMatSolvePC
MoFEMErrorCode setSchurMatSolvePC(SmartPetscObj< PC > pc)
Definition: Schur.cpp:2823
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
ContactOps::cn_contact
double cn_contact
Definition: EshelbianContact.hpp:19
b_iso
double b_iso
Saturation exponent.
Definition: plastic.cpp:128
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
MoFEM::PipelineManager::VolEle
MoFEM::VolumeElementForcesAndSourcesCore VolEle
Definition: PipelineManager.hpp:34
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
PlasticOps
Definition: PlasticNaturalBCs.hpp:13
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
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
Example::bC
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:504
IntegrationRules.hpp
iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:72
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
is_large_strains
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:116
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
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
FTensor::Ddg
Definition: Ddg_value.hpp:7
Example::ScaledTimeScale
Definition: plastic.cpp:198
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
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ep_order
int ep_order
Order of ep files.
Definition: plastic.cpp:136
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
Qinf
double Qinf
Saturation yield stress.
Definition: plastic.cpp:127
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
Example::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:195
Example::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:194
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:184
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:40
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1196
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:396
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:578
MoFEM::SmartPetscObj< VecScatter >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:760
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
Definition: SetUpSchurImpl.cpp:213
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
PlasticNaturalBCs.hpp
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
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
H
double H
Hardening.
Definition: plastic.cpp:124
tau_order
int tau_order
Order of tau files.
Definition: plastic.cpp:135
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
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
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
Definition: plastic.cpp:1285
alpha_damping
double alpha_damping
Definition: plastic.cpp:141
MoFEM::getNestSchurData
boost::shared_ptr< NestSchurData > getNestSchurData(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:2249
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
SetUpSchurImpl::P
SmartPetscObj< Mat > P
Definition: SetUpSchurImpl.cpp:29
Generated by Doxygen 1.8.17 and hosted at University of Glasgow