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  MoFEMErrorCode testOperators();
192 
193  boost::shared_ptr<DomainEle> reactionFe;
194 
195  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
196  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
197  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
198 
199  struct ScaledTimeScale : public MoFEM::TimeScale {
200  using MoFEM::TimeScale::TimeScale;
201  double getScale(const double time) {
202  return scale * MoFEM::TimeScale::getScale(time);
203  };
204  };
205 
206 #ifdef ADD_CONTACT
207  #ifdef PYTHON_SDF
208  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
209  #endif
210 #endif // ADD_CONTACT
211 };
212 
213 //! [Run problem]
214 MoFEMErrorCode Example::runProblem() {
215  MoFEMFunctionBegin;
216  CHKERR createCommonData();
217  CHKERR setupProblem();
218  CHKERR bC();
219  CHKERR OPs();
220  PetscBool test_ops = PETSC_FALSE;
221  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_operators", &test_ops,
222  PETSC_NULL);
223  if (test_ops == PETSC_FALSE) {
224  CHKERR tsSolve();
225  } else {
226  CHKERR testOperators();
227  }
228  MoFEMFunctionReturn(0);
229 }
230 //! [Run problem]
231 
232 //! [Set up problem]
233 MoFEMErrorCode Example::setupProblem() {
234  MoFEMFunctionBegin;
235  Simple *simple = mField.getInterface<Simple>();
236 
237  Range domain_ents;
238  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
239  true);
240  auto get_ents_by_dim = [&](const auto dim) {
241  if (dim == SPACE_DIM) {
242  return domain_ents;
243  } else {
244  Range ents;
245  if (dim == 0)
246  CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
247  else
248  CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
249  return ents;
250  }
251  };
252 
253  auto get_base = [&]() {
254  auto domain_ents = get_ents_by_dim(SPACE_DIM);
255  if (domain_ents.empty())
256  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
257  const auto type = type_from_handle(domain_ents[0]);
258  switch (type) {
259  case MBQUAD:
260  return DEMKOWICZ_JACOBI_BASE;
261  case MBHEX:
262  return DEMKOWICZ_JACOBI_BASE;
263  case MBTRI:
264  return AINSWORTH_LEGENDRE_BASE;
265  case MBTET:
266  return AINSWORTH_LEGENDRE_BASE;
267  default:
268  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
269  }
270  return NOBASE;
271  };
272 
273  const auto base = get_base();
274  MOFEM_LOG("PLASTICITY", Sev::inform)
275  << "Base " << ApproximationBaseNames[base];
276 
277  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
278  CHKERR simple->addDomainField("EP", L2, base, size_symm);
279  CHKERR simple->addDomainField("TAU", L2, base, 1);
280  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
281 
282  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
283 
284  PetscBool order_edge = PETSC_FALSE;
285  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,
286  PETSC_NULL);
287  PetscBool order_face = PETSC_FALSE;
288  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
289  PETSC_NULL);
290  PetscBool order_volume = PETSC_FALSE;
291  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
292  PETSC_NULL);
293 
294  if (order_edge || order_face || order_volume) {
295 
296  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
297  ? "true"
298  : "false";
299  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
300  ? "true"
301  : "false";
302  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
303  ? "true"
304  : "false";
305 
306  auto ents = get_ents_by_dim(0);
307  if (order_edge)
308  ents.merge(get_ents_by_dim(1));
309  if (order_face)
310  ents.merge(get_ents_by_dim(2));
311  if (order_volume)
312  ents.merge(get_ents_by_dim(3));
313  CHKERR simple->setFieldOrder("U", order, &ents);
314  } else {
315  CHKERR simple->setFieldOrder("U", order);
316  }
317  CHKERR simple->setFieldOrder("EP", ep_order);
318  CHKERR simple->setFieldOrder("TAU", tau_order);
319 
320  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
321 
322 #ifdef ADD_CONTACT
323  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
324  SPACE_DIM);
325  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
326  SPACE_DIM);
327 
328  auto get_skin = [&]() {
329  Range body_ents;
330  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
331  Skinner skin(&mField.get_moab());
332  Range skin_ents;
333  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
334  return skin_ents;
335  };
336 
337  auto filter_blocks = [&](auto skin) {
338  bool is_contact_block = true;
339  Range contact_range;
340  for (auto m :
341  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
342 
343  (boost::format("%s(.*)") % "CONTACT").str()
344 
345  ))
346 
347  ) {
348  is_contact_block =
349  true; ///< blocs interation is collective, so that is set irrespective
350  ///< if there are entities in given rank or not in the block
351  MOFEM_LOG("CONTACT", Sev::inform)
352  << "Find contact block set: " << m->getName();
353  auto meshset = m->getMeshset();
354  Range contact_meshset_range;
355  CHKERR mField.get_moab().get_entities_by_dimension(
356  meshset, SPACE_DIM - 1, contact_meshset_range, true);
357 
358  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
359  contact_meshset_range);
360  contact_range.merge(contact_meshset_range);
361  }
362  if (is_contact_block) {
363  MOFEM_LOG("SYNC", Sev::inform)
364  << "Nb entities in contact surface: " << contact_range.size();
365  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
366  skin = intersect(skin, contact_range);
367  }
368  return skin;
369  };
370 
371  auto filter_true_skin = [&](auto skin) {
372  Range boundary_ents;
373  ParallelComm *pcomm =
374  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
375  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
376  PSTATUS_NOT, -1, &boundary_ents);
377  return boundary_ents;
378  };
379 
380  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
381  CHKERR simple->setFieldOrder("SIGMA", 0);
382  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
383 #endif
384 
385  CHKERR simple->setUp();
386  CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
387 
388  auto project_ho_geometry = [&]() {
389  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
390  return mField.loop_dofs("GEOMETRY", ent_method);
391  };
392  PetscBool project_geometry = PETSC_TRUE;
393  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
394  &project_geometry, PETSC_NULL);
395  if (project_geometry) {
396  CHKERR project_ho_geometry();
397  }
398 
399  MoFEMFunctionReturn(0);
400 }
401 //! [Set up problem]
402 
403 //! [Create common data]
404 MoFEMErrorCode Example::createCommonData() {
405  MoFEMFunctionBegin;
406 
407  auto get_command_line_parameters = [&]() {
408  MoFEMFunctionBegin;
409 
410  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
411  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
412  &young_modulus, PETSC_NULL);
413  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
414  &poisson_ratio, PETSC_NULL);
415  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
416  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
417  PETSC_NULL);
418  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
419  PETSC_NULL);
420  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
421  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
422  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
423  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
424  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
425  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);
426  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
427  &is_large_strains, PETSC_NULL);
428  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
429  PETSC_NULL);
430 
431  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
432  PetscBool tau_order_is_set; ///< true if tau order is set
433  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
434  &tau_order_is_set);
435  PetscBool ep_order_is_set; ///< true if tau order is set
436  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
437  &ep_order_is_set);
438  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
439  PETSC_NULL);
440 
441  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
442  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
443  &alpha_damping, PETSC_NULL);
444 
445  MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
446  MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
447  MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
448  MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
449  MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
450  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
451  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
452  MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
453  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
454  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
455  MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
456 
457  if (tau_order_is_set == PETSC_FALSE)
458  tau_order = order - 2;
459  if (ep_order_is_set == PETSC_FALSE)
460  ep_order = order - 1;
461 
462  MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
463  MOFEM_LOG("PLASTICITY", Sev::inform)
464  << "Ep approximation order " << ep_order;
465  MOFEM_LOG("PLASTICITY", Sev::inform)
466  << "Tau approximation order " << tau_order;
467  MOFEM_LOG("PLASTICITY", Sev::inform)
468  << "Geometry approximation order " << geom_order;
469 
470  MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
471  MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
472 
473  PetscBool is_scale = PETSC_TRUE;
474  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
475  PETSC_NULL);
476  if (is_scale) {
477  scale /= young_modulus;
478  }
479 
480  MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
481 
482 #ifdef ADD_CONTACT
483  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
484  &ContactOps::cn_contact, PETSC_NULL);
485  MOFEM_LOG("CONTACT", Sev::inform)
486  << "cn_contact " << ContactOps::cn_contact;
487 #endif // ADD_CONTACT
488 
489  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
490  &is_quasi_static, PETSC_NULL);
491  MOFEM_LOG("PLASTICITY", Sev::inform)
492  << "Is quasi static: " << (is_quasi_static ? "true" : "false");
493 
494  MoFEMFunctionReturn(0);
495  };
496 
497  CHKERR get_command_line_parameters();
498 
499 #ifdef ADD_CONTACT
500  #ifdef PYTHON_SDF
501  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
502  CHKERR sdfPythonPtr->sdfInit("sdf.py");
503  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
504  #endif
505 #endif // ADD_CONTACT
506 
507  MoFEMFunctionReturn(0);
508 }
509 //! [Create common data]
510 
511 //! [Boundary condition]
512 MoFEMErrorCode Example::bC() {
513  MoFEMFunctionBegin;
514 
515  auto simple = mField.getInterface<Simple>();
516  auto bc_mng = mField.getInterface<BcManager>();
517 
518  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
519  "U", 0, 0);
520  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
521  "U", 1, 1);
522  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
523  "U", 2, 2);
524  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
525  "REMOVE_ALL", "U", 0, 3);
526 
527 #ifdef ADD_CONTACT
528  for (auto b : {"FIX_X", "REMOVE_X"})
529  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
530  "SIGMA", 0, 0, false, true);
531  for (auto b : {"FIX_Y", "REMOVE_Y"})
532  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
533  "SIGMA", 1, 1, false, true);
534  for (auto b : {"FIX_Z", "REMOVE_Z"})
535  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
536  "SIGMA", 2, 2, false, true);
537  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
538  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
539  "SIGMA", 0, 3, false, true);
540  CHKERR bc_mng->removeBlockDOFsOnEntities(
541  simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
542 #endif
543 
544  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
545  simple->getProblemName(), "U");
546 
547  auto &bc_map = bc_mng->getBcMapByBlockName();
548  for (auto bc : bc_map)
549  MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
550 
551  MoFEMFunctionReturn(0);
552 }
553 //! [Boundary condition]
554 
555 //! [Push operators to pipeline]
556 MoFEMErrorCode Example::OPs() {
557  MoFEMFunctionBegin;
558  auto pip_mng = mField.getInterface<PipelineManager>();
559 
560  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
561 
562  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
563 
564  auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
565  MoFEMFunctionBegin;
566 
567  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
568  "GEOMETRY");
569  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
570 
571  // Add Natural BCs to LHS
572  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
573  pip, mField, "U", Sev::inform);
574 
575 #ifdef ADD_CONTACT
576  auto simple = mField.getInterface<Simple>();
577  CHKERR
578  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
579  pip, "SIGMA", "U");
580  CHKERR
581  ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
582  mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
583  vol_rule);
584 #endif // ADD_CONTACT
585 
586  MoFEMFunctionReturn(0);
587  };
588 
589  auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
590  MoFEMFunctionBegin;
591 
592  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
593  "GEOMETRY");
594  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
595 
596  // Add Natural BCs to RHS
597  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
598  pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
599 
600 #ifdef ADD_CONTACT
601  CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
602  pip, "SIGMA", "U");
603 #endif // ADD_CONTACT
604 
605  MoFEMFunctionReturn(0);
606  };
607 
608  auto add_domain_ops_lhs = [this](auto &pip) {
609  MoFEMFunctionBegin;
610  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
611  "GEOMETRY");
612 
613  if (is_quasi_static == PETSC_FALSE) {
614 
615  //! [Only used for dynamics]
616  using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
617  GAUSS>::OpMass<1, SPACE_DIM>;
618  //! [Only used for dynamics]
619 
620  auto get_inertia_and_mass_damping = [this](const double, const double,
621  const double) {
622  auto *pip = mField.getInterface<PipelineManager>();
623  auto &fe_domain_lhs = pip->getDomainLhsFE();
624  return (rho / scale) * fe_domain_lhs->ts_aa +
625  (alpha_damping / scale) * fe_domain_lhs->ts_a;
626  };
627  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
628  }
629 
630  CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
631  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
632 
633  MoFEMFunctionReturn(0);
634  };
635 
636  auto add_domain_ops_rhs = [this](auto &pip) {
637  MoFEMFunctionBegin;
638 
639  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
640  "GEOMETRY");
641 
642  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
643  pip, mField, "U",
644  {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
645  Sev::inform);
646 
647  // only in case of dynamics
648  if (is_quasi_static == PETSC_FALSE) {
649 
650  //! [Only used for dynamics]
651  using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
652  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
653  //! [Only used for dynamics]
654 
655  auto mat_acceleration = boost::make_shared<MatrixDouble>();
656  pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
657  "U", mat_acceleration));
658  pip.push_back(
659  new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
660  return rho / scale;
661  }));
662  if (alpha_damping > 0) {
663  auto mat_velocity = boost::make_shared<MatrixDouble>();
664  pip.push_back(
665  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
666  pip.push_back(
667  new OpInertiaForce("U", mat_velocity, [](double, double, double) {
668  return alpha_damping / scale;
669  }));
670  }
671  }
672 
673  CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
674  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
675 
676 #ifdef ADD_CONTACT
677  CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
678  pip, "SIGMA", "U");
679 #endif // ADD_CONTACT
680 
681  MoFEMFunctionReturn(0);
682  };
683 
684  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
685  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
686 
687  // Boundary
688  CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
689  CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
690 
691  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
692  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
693 
694  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
695  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
696 
697  auto create_reaction_pipeline = [&](auto &pip) {
698  MoFEMFunctionBegin;
699  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
700  "GEOMETRY");
701  CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
702  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
703  MoFEMFunctionReturn(0);
704  };
705 
706  reactionFe = boost::make_shared<DomainEle>(mField);
707  reactionFe->getRuleHook = vol_rule;
708  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
709  reactionFe->postProcessHook =
710  EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
711 
712  MoFEMFunctionReturn(0);
713 }
714 //! [Push operators to pipeline]
715 
716 //! [Solve]
717 struct SetUpSchur {
718 
719  /**
720  * @brief Create data structure for handling Schur complement
721  *
722  * @param m_field
723  * @param sub_dm Schur complement sub dm
724  * @param field_split_it IS of Schur block
725  * @param ao_map AO map from sub dm to main problem
726  * @return boost::shared_ptr<SetUpSchur>
727  */
728  static boost::shared_ptr<SetUpSchur> createSetUpSchur(
729 
730  MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
731  SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
732 
733  );
734  virtual MoFEMErrorCode setUp(TS solver) = 0;
735 
736 protected:
737  SetUpSchur() = default;
738 };
739 
740 MoFEMErrorCode Example::tsSolve() {
741  MoFEMFunctionBegin;
742 
743  Simple *simple = mField.getInterface<Simple>();
744  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
745  ISManager *is_manager = mField.getInterface<ISManager>();
746 
747  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
748 
749  auto set_section_monitor = [&](auto solver) {
750  MoFEMFunctionBegin;
751  SNES snes;
752  CHKERR TSGetSNES(solver, &snes);
753  CHKERR SNESMonitorSet(snes,
754  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
755  void *))MoFEMSNESMonitorFields,
756  (void *)(snes_ctx_ptr.get()), nullptr);
757  MoFEMFunctionReturn(0);
758  };
759 
760  auto create_post_process_elements = [&]() {
761  auto pp_fe = boost::make_shared<PostProcEle>(mField);
762 
763  auto push_vol_ops = [this](auto &pip) {
764  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
765  "GEOMETRY");
766 
767  auto [common_plastic_ptr, common_hencky_ptr] =
768  PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
769  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
770 
771  if (common_hencky_ptr) {
772  if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
773  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
774  }
775 
776  return std::make_pair(common_plastic_ptr, common_hencky_ptr);
777  };
778 
779  auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
780  MoFEMFunctionBegin;
781 
782  auto &pip = pp_fe->getOpPtrVector();
783 
784  auto [common_plastic_ptr, common_hencky_ptr] = p;
785 
786  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
787 
788  auto x_ptr = boost::make_shared<MatrixDouble>();
789  pip.push_back(
790  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
791  auto u_ptr = boost::make_shared<MatrixDouble>();
792  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
793 
794  if (is_large_strains) {
795 
796  pip.push_back(
797 
798  new OpPPMap(
799 
800  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
801 
802  {{"PLASTIC_SURFACE",
803  common_plastic_ptr->getPlasticSurfacePtr()},
804  {"PLASTIC_MULTIPLIER",
805  common_plastic_ptr->getPlasticTauPtr()}},
806 
807  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
808 
809  {{"GRAD", common_hencky_ptr->matGradPtr},
810  {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
811 
812  {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
813  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
814 
815  )
816 
817  );
818 
819  } else {
820 
821  pip.push_back(
822 
823  new OpPPMap(
824 
825  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
826 
827  {{"PLASTIC_SURFACE",
828  common_plastic_ptr->getPlasticSurfacePtr()},
829  {"PLASTIC_MULTIPLIER",
830  common_plastic_ptr->getPlasticTauPtr()}},
831 
832  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
833 
834  {},
835 
836  {{"STRAIN", common_plastic_ptr->mStrainPtr},
837  {"STRESS", common_plastic_ptr->mStressPtr},
838  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
839  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
840 
841  )
842 
843  );
844  }
845 
846  MoFEMFunctionReturn(0);
847  };
848 
849  auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
850  PetscBool post_proc_vol = PETSC_FALSE;
851  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
852  &post_proc_vol, PETSC_NULL);
853  if (post_proc_vol == PETSC_FALSE)
854  return boost::shared_ptr<PostProcEle>();
855  auto pp_fe = boost::make_shared<PostProcEle>(mField);
856  CHK_MOAB_THROW(
857  push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
858  "push_vol_post_proc_ops");
859  return pp_fe;
860  };
861 
862  auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
863  PetscBool post_proc_skin = PETSC_TRUE;
864  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
865  &post_proc_skin, PETSC_NULL);
866  if (post_proc_skin == PETSC_FALSE)
867  return boost::shared_ptr<SkinPostProcEle>();
868 
869  auto simple = mField.getInterface<Simple>();
870  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
871  auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
872  SPACE_DIM, Sev::verbose);
873  pp_fe->getOpPtrVector().push_back(op_side);
874  CHK_MOAB_THROW(push_vol_post_proc_ops(
875  pp_fe, push_vol_ops(op_side->getOpPtrVector())),
876  "push_vol_post_proc_ops");
877  return pp_fe;
878  };
879 
880  return std::make_pair(vol_post_proc(), skin_post_proc());
881  };
882 
883  auto scatter_create = [&](auto D, auto coeff) {
884  SmartPetscObj<IS> is;
885  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
886  ROW, "U", coeff, coeff, is);
887  int loc_size;
888  CHKERR ISGetLocalSize(is, &loc_size);
889  Vec v;
890  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
891  VecScatter scatter;
892  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
893  return std::make_tuple(SmartPetscObj<Vec>(v),
894  SmartPetscObj<VecScatter>(scatter));
895  };
896 
897  auto set_time_monitor = [&](auto dm, auto solver) {
898  MoFEMFunctionBegin;
899  boost::shared_ptr<Monitor> monitor_ptr(
900  new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
901  uYScatter, uZScatter));
902  boost::shared_ptr<ForcesAndSourcesCore> null;
903  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
904  monitor_ptr, null, null);
905  MoFEMFunctionReturn(0);
906  };
907 
908  auto set_schur_pc = [&](auto solver,
909  boost::shared_ptr<SetUpSchur> &schur_ptr) {
910  MoFEMFunctionBeginHot;
911 
912  auto name_prb = simple->getProblemName();
913 
914  // create sub dm for Schur complement
915  auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
916  SmartPetscObj<DM> &dm_sub) {
917  MoFEMFunctionBegin;
918  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
919  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
920  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
921  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
922  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
923  for (auto f : {"U"}) {
924  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
925  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
926  }
927  CHKERR DMSetUp(dm_sub);
928 
929  MoFEMFunctionReturn(0);
930  };
931 
932  auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
933  SmartPetscObj<DM> &dm_sub) {
934  MoFEMFunctionBegin;
935  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
936  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
937  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
938  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
939  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
940 #ifdef ADD_CONTACT
941  for (auto f : {"SIGMA", "EP", "TAU"}) {
942  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
943  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
944  }
945 #else
946  for (auto f : {"EP", "TAU"}) {
947  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
948  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
949  }
950 #endif
951  CHKERR DMSetUp(dm_sub);
952  MoFEMFunctionReturn(0);
953  };
954 
955  // Create nested (sub BC) Schur DM
956  if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
957 
958  SmartPetscObj<DM> dm_schur;
959  CHKERR create_schur_dm(simple->getDM(), dm_schur);
960  SmartPetscObj<DM> dm_block;
961  CHKERR create_block_dm(simple->getDM(), dm_block);
962 
963 #ifdef ADD_CONTACT
964 
965  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
966  auto block_mat_data = createBlockMatStructure(
967  simple->getDM(),
968 
969  {
970 
971  {simple->getDomainFEName(),
972 
973  {{"U", "U"},
974  {"SIGMA", "SIGMA"},
975  {"U", "SIGMA"},
976  {"SIGMA", "U"},
977  {"EP", "EP"},
978  {"TAU", "TAU"},
979  {"U", "EP"},
980  {"EP", "U"},
981  {"EP", "TAU"},
982  {"TAU", "EP"},
983  {"TAU", "U"}
984 
985  }},
986 
987  {simple->getBoundaryFEName(),
988 
989  {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
990 
991  }}
992 
993  }
994 
995  );
996 
997  return createSchurNestedMatrixStruture(
998 
999  {dm_schur, dm_block}, block_mat_data,
1000 
1001  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1002 
1003  );
1004  };
1005 
1006 #else
1007 
1008  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1009  auto block_mat_data =
1010  createBlockMatStructure(simple->getDM(),
1011 
1012  {{simple->getDomainFEName(),
1013 
1014  {{"U", "U"},
1015  {"EP", "EP"},
1016  {"TAU", "TAU"},
1017  {"U", "EP"},
1018  {"EP", "U"},
1019  {"EP", "TAU"},
1020  {"TAU", "U"},
1021  {"TAU", "EP"}
1022 
1023  }}}
1024 
1025  );
1026 
1027  return createSchurNestedMatrixStruture(
1028 
1029  {dm_schur, dm_block}, block_mat_data,
1030 
1031  {"EP", "TAU"}, {nullptr, nullptr}, false
1032 
1033  );
1034  };
1035 
1036 #endif
1037 
1038  auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1039  CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1040 
1041  auto block_is = getDMSubData(dm_block)->getSmartRowIs();
1042  auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
1043 
1044  // Indices has to be map fro very to level, while assembling Schur
1045  // complement.
1046  schur_ptr =
1047  SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
1048  CHKERR schur_ptr->setUp(solver);
1049  }
1050 
1051  MoFEMFunctionReturnHot(0);
1052  };
1053 
1054  auto dm = simple->getDM();
1055  auto D = createDMVector(dm);
1056  auto DD = vectorDuplicate(D);
1057  uXScatter = scatter_create(D, 0);
1058  uYScatter = scatter_create(D, 1);
1059  if constexpr (SPACE_DIM == 3)
1060  uZScatter = scatter_create(D, 2);
1061 
1062  auto create_solver = [pip_mng]() {
1063  if (is_quasi_static == PETSC_TRUE)
1064  return pip_mng->createTSIM();
1065  else
1066  return pip_mng->createTSIM2();
1067  };
1068 
1069  auto solver = create_solver();
1070 
1071  auto active_pre_lhs = []() {
1072  MoFEMFunctionBegin;
1073  std::fill(PlasticOps::CommonData::activityData.begin(),
1074  PlasticOps::CommonData::activityData.end(), 0);
1075  MoFEMFunctionReturn(0);
1076  };
1077 
1078  auto active_post_lhs = [&]() {
1079  MoFEMFunctionBegin;
1080  auto get_iter = [&]() {
1081  SNES snes;
1082  CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1083  int iter;
1084  CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1085  "Can not get iter");
1086  return iter;
1087  };
1088 
1089  auto iter = get_iter();
1090  if (iter >= 0) {
1091 
1092  std::array<int, 5> activity_data;
1093  std::fill(activity_data.begin(), activity_data.end(), 0);
1094  MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1095  activity_data.data(), activity_data.size(), MPI_INT,
1096  MPI_SUM, mField.get_comm());
1097 
1098  int &active_points = activity_data[0];
1099  int &avtive_full_elems = activity_data[1];
1100  int &avtive_elems = activity_data[2];
1101  int &nb_points = activity_data[3];
1102  int &nb_elements = activity_data[4];
1103 
1104  if (nb_points) {
1105 
1106  double proc_nb_points =
1107  100 * static_cast<double>(active_points) / nb_points;
1108  double proc_nb_active =
1109  100 * static_cast<double>(avtive_elems) / nb_elements;
1110  double proc_nb_full_active = 100;
1111  if (avtive_elems)
1112  proc_nb_full_active =
1113  100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1114 
1115  MOFEM_LOG_C("PLASTICITY", Sev::inform,
1116  "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1117  "elements %d "
1118  "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1119  iter, nb_points, active_points, proc_nb_points,
1120  avtive_elems, proc_nb_active, avtive_full_elems,
1121  proc_nb_full_active, iter);
1122  }
1123  }
1124 
1125  MoFEMFunctionReturn(0);
1126  };
1127 
1128  auto add_active_dofs_elem = [&](auto dm) {
1129  MoFEMFunctionBegin;
1130  auto fe_pre_proc = boost::make_shared<FEMethod>();
1131  fe_pre_proc->preProcessHook = active_pre_lhs;
1132  auto fe_post_proc = boost::make_shared<FEMethod>();
1133  fe_post_proc->postProcessHook = active_post_lhs;
1134  auto ts_ctx_ptr = getDMTsCtx(dm);
1135  ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1136  ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1137  MoFEMFunctionReturn(0);
1138  };
1139 
1140  auto set_essential_bc = [&](auto dm, auto solver) {
1141  MoFEMFunctionBegin;
1142  // This is low level pushing finite elements (pipelines) to solver
1143 
1144  auto pre_proc_ptr = boost::make_shared<FEMethod>();
1145  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1146  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1147 
1148  // Add boundary condition scaling
1149  auto disp_time_scale = boost::make_shared<TimeScale>();
1150 
1151  auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1152  EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1153  {disp_time_scale}, false);
1154  return hook;
1155  };
1156  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1157 
1158  auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1159  MoFEMFunctionBegin;
1160  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
1161  mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1162  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
1163  mField, post_proc_rhs_ptr, 1.)();
1164  MoFEMFunctionReturn(0);
1165  };
1166  auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1167  return EssentialPostProcLhs<DisplacementCubitBcData>(
1168  mField, post_proc_lhs_ptr, 1.);
1169  };
1170  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1171 
1172  auto ts_ctx_ptr = getDMTsCtx(dm);
1173  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1174  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1175  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1176  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1177  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1178 
1179  MoFEMFunctionReturn(0);
1180  };
1181 
1182  if (is_quasi_static == PETSC_TRUE) {
1183  CHKERR TSSetSolution(solver, D);
1184  } else {
1185  CHKERR TS2SetSolution(solver, D, DD);
1186  }
1187 
1188  CHKERR set_section_monitor(solver);
1189  CHKERR set_time_monitor(dm, solver);
1190  CHKERR TSSetFromOptions(solver);
1191  CHKERR TSSetUp(solver);
1192 
1193  CHKERR add_active_dofs_elem(dm);
1194  boost::shared_ptr<SetUpSchur> schur_ptr;
1195  CHKERR set_schur_pc(solver, schur_ptr);
1196  CHKERR set_essential_bc(dm, solver);
1197 
1198  MOFEM_LOG_CHANNEL("TIMER");
1199  MOFEM_LOG_TAG("TIMER", "timer");
1200  if (set_timer)
1201  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1202  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1203  CHKERR TSSetUp(solver);
1204  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1205  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1206  CHKERR TSSolve(solver, NULL);
1207  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1208 
1209  MoFEMFunctionReturn(0);
1210 }
1211 //! [Solve]
1212 
1213 //! [TestOperators]
1214 MoFEMErrorCode Example::testOperators() {
1215  MoFEMFunctionBegin;
1216 
1217  // get operators tester
1218  auto simple = mField.getInterface<Simple>();
1219  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1220  // OperatorsTester
1221  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1222  // pipeline manager
1223 
1224  constexpr double eps = 1e-9;
1225 
1226  auto x = opt->setRandomFields(simple->getDM(), {
1227 
1228  {"U", {-1e-6, 1e-6}},
1229 
1230  {"EP", {-1e-6, 1e-6}},
1231 
1232  {"TAU", {0, 1e-4}}
1233 
1234  });
1235 
1236  auto dot_x_plastic_active =
1237  opt->setRandomFields(simple->getDM(), {
1238 
1239  {"U", {-1, 1}},
1240 
1241  {"EP", {-1, 1}},
1242 
1243  {"TAU", {0.1, 1}}
1244 
1245  });
1246  auto diff_x_plastic_active =
1247  opt->setRandomFields(simple->getDM(), {
1248 
1249  {"U", {-1, 1}},
1250 
1251  {"EP", {-1, 1}},
1252 
1253  {"TAU", {-1, 1}}
1254 
1255  });
1256 
1257  auto dot_x_elastic =
1258  opt->setRandomFields(simple->getDM(), {
1259 
1260  {"U", {-1, 1}},
1261 
1262  {"EP", {-1, 1}},
1263 
1264  {"TAU", {-1, -0.1}}
1265 
1266  });
1267  auto diff_x_elastic =
1268  opt->setRandomFields(simple->getDM(), {
1269 
1270  {"U", {-1, 1}},
1271 
1272  {"EP", {-1, 1}},
1273 
1274  {"TAU", {-1, 1}}
1275 
1276  });
1277 
1278  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
1279  auto dot_x, auto diff_x) {
1280  MoFEMFunctionBegin;
1281 
1282  auto diff_res = opt->checkCentralFiniteDifference(
1283  simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1284  SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
1285 
1286  // Calculate norm of difference between directional derivative calculated
1287  // from finite difference, and tangent matrix.
1288  double fnorm;
1289  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1290  MOFEM_LOG_C("PLASTICITY", Sev::inform,
1291  "Test consistency of tangent matrix %3.4e", fnorm);
1292 
1293  constexpr double err = 1e-5;
1294  if (fnorm > err)
1295  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1296  "Norm of directional derivative too large err = %3.4e", fnorm);
1297 
1298  MoFEMFunctionReturn(0);
1299  };
1300 
1301  MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
1302  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1303  pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
1304 
1305  MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
1306  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1307  pip->getDomainRhsFE(), dot_x_plastic_active,
1308  diff_x_plastic_active);
1309 
1310 
1311 
1312  MoFEMFunctionReturn(0);
1313 };
1314 
1315 //! [TestOperators]
1316 
1317 static char help[] = "...\n\n";
1318 
1319 int main(int argc, char *argv[]) {
1320 
1321 #ifdef ADD_CONTACT
1322  #ifdef PYTHON_SDF
1323  Py_Initialize();
1324  np::initialize();
1325  #endif
1326 #endif // ADD_CONTACT
1327 
1328  // Initialisation of MoFEM/PETSc and MOAB data structures
1329  const char param_file[] = "param_file.petsc";
1330  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1331 
1332  // Add logging channel for example
1333  auto core_log = logging::core::get();
1334  core_log->add_sink(
1335  LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
1336  core_log->add_sink(
1337  LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
1338  LogManager::setLog("PLASTICITY");
1339  MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
1340 
1341 #ifdef ADD_CONTACT
1342  core_log->add_sink(
1343  LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
1344  LogManager::setLog("CONTACT");
1345  MOFEM_LOG_TAG("CONTACT", "Contact");
1346 #endif // ADD_CONTACT
1347 
1348  try {
1349 
1350  //! [Register MoFEM discrete manager in PETSc]
1351  DMType dm_name = "DMMOFEM";
1352  CHKERR DMRegister_MoFEM(dm_name);
1353  //! [Register MoFEM discrete manager in PETSc
1354 
1355  //! [Create MoAB]
1356  moab::Core mb_instance; ///< mesh database
1357  moab::Interface &moab = mb_instance; ///< mesh database interface
1358  //! [Create MoAB]
1359 
1360  //! [Create MoFEM]
1361  MoFEM::Core core(moab); ///< finite element database
1362  MoFEM::Interface &m_field = core; ///< finite element database interface
1363  //! [Create MoFEM]
1364 
1365  //! [Load mesh]
1366  Simple *simple = m_field.getInterface<Simple>();
1367  CHKERR simple->getOptions();
1368  CHKERR simple->loadFile();
1369  //! [Load mesh]
1370 
1371  //! [Example]
1372  Example ex(m_field);
1373  CHKERR ex.runProblem();
1374  //! [Example]
1375  }
1376  CATCH_ERRORS;
1377 
1378  CHKERR MoFEM::Core::Finalize();
1379 
1380 #ifdef ADD_CONTACT
1381  #ifdef PYTHON_SDF
1382  if (Py_FinalizeEx() < 0) {
1383  exit(120);
1384  }
1385  #endif
1386 #endif // ADD_CONTACT
1387 
1388  return 0;
1389 }
1390 
1391 struct SetUpSchurImpl : public SetUpSchur {
1392 
1393  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
1394  SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
1395  : SetUpSchur(), mField(m_field), subDM(sub_dm),
1396  fieldSplitIS(field_split_is), aoSchur(ao_up) {
1397  if (S) {
1398  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1399  "Is expected that schur matrix is not "
1400  "allocated. This is "
1401  "possible only is PC is set up twice");
1402  }
1403  }
1404  virtual ~SetUpSchurImpl() { S.reset(); }
1405 
1406  MoFEMErrorCode setUp(TS solver);
1407  MoFEMErrorCode preProc();
1408  MoFEMErrorCode postProc();
1409 
1410 private:
1411  SmartPetscObj<Mat> S;
1412 
1413  MoFEM::Interface &mField;
1414  SmartPetscObj<DM> subDM; ///< field split sub dm
1415  SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
1416  SmartPetscObj<AO> aoSchur; ///> main DM to subDM
1417 };
1418 
1419 MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
1420  MoFEMFunctionBegin;
1421  auto simple = mField.getInterface<Simple>();
1422  auto pip_mng = mField.getInterface<PipelineManager>();
1423 
1424  SNES snes;
1425  CHKERR TSGetSNES(solver, &snes);
1426  KSP ksp;
1427  CHKERR SNESGetKSP(snes, &ksp);
1428  CHKERR KSPSetFromOptions(ksp);
1429 
1430  PC pc;
1431  CHKERR KSPGetPC(ksp, &pc);
1432  PetscBool is_pcfs = PETSC_FALSE;
1433  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1434  if (is_pcfs) {
1435  if (S) {
1436  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1437  "Is expected that schur matrix is not "
1438  "allocated. This is "
1439  "possible only is PC is set up twice");
1440  }
1441 
1442  S = createDMMatrix(subDM);
1443  CHKERR MatSetBlockSize(S, SPACE_DIM);
1444 
1445  // Set DM to use shell block matrix
1446  DM solver_dm;
1447  CHKERR TSGetDM(solver, &solver_dm);
1448  CHKERR DMSetMatType(solver_dm, MATSHELL);
1449 
1450  auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1451  auto A = createDMBlockMat(simple->getDM());
1452  auto P = createDMNestSchurMat(simple->getDM());
1453 
1454  if (is_quasi_static == PETSC_TRUE) {
1455  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1456  Mat A, Mat B, void *ctx) {
1457  return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1458  };
1459  CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1460  } else {
1461  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1462  PetscReal a, PetscReal aa, Mat A, Mat B,
1463  void *ctx) {
1464  return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1465  };
1466  CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1467  }
1468  CHKERR KSPSetOperators(ksp, A, P);
1469 
1470  auto set_ops = [&]() {
1471  MoFEMFunctionBegin;
1472  auto pip_mng = mField.getInterface<PipelineManager>();
1473 
1474 #ifndef ADD_CONTACT
1475  // Boundary
1476  pip_mng->getOpBoundaryLhsPipeline().push_front(
1477  createOpSchurAssembleBegin());
1478  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1479 
1480  {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1481 
1482  ));
1483  // Domain
1484  pip_mng->getOpDomainLhsPipeline().push_front(
1485  createOpSchurAssembleBegin());
1486  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1487 
1488  {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1489 
1490  ));
1491 #else
1492 
1493  double eps_stab = 1e-4;
1494  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1495  PETSC_NULL);
1496 
1497  using B = FormsIntegrators<BoundaryEleOp>::Assembly<
1498  BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
1499  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1500 
1501  // Boundary
1502  pip_mng->getOpBoundaryLhsPipeline().push_front(
1503  createOpSchurAssembleBegin());
1504  pip_mng->getOpBoundaryLhsPipeline().push_back(
1505  new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1506  return eps_stab;
1507  }));
1508  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1509 
1510  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1511  false, false
1512 
1513  ));
1514  // Domain
1515  pip_mng->getOpDomainLhsPipeline().push_front(
1516  createOpSchurAssembleBegin());
1517  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1518 
1519  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1520  false, false
1521 
1522  ));
1523 #endif // ADD_CONTACT
1524  MoFEMFunctionReturn(0);
1525  };
1526 
1527  auto set_assemble_elems = [&]() {
1528  MoFEMFunctionBegin;
1529  auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1530  schur_asmb_pre_proc->preProcessHook = [this]() {
1531  MoFEMFunctionBegin;
1532  CHKERR MatZeroEntries(S);
1533  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1534  MoFEMFunctionReturn(0);
1535  };
1536  auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1537 
1538  schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1539  MoFEMFunctionBegin;
1540  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1541 
1542  // Apply essential constrains to Schur complement
1543  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1544  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1545  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1546  mField, schur_asmb_post_proc, 1, S, aoSchur)();
1547 
1548  MoFEMFunctionReturn(0);
1549  };
1550  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1551  ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1552  ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1553  MoFEMFunctionReturn(0);
1554  };
1555 
1556  auto set_pc = [&]() {
1557  MoFEMFunctionBegin;
1558  CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1559  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1560  MoFEMFunctionReturn(0);
1561  };
1562 
1563  auto set_diagonal_pc = [&]() {
1564  MoFEMFunctionBegin;
1565  KSP *subksp;
1566  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
1567  auto get_pc = [](auto ksp) {
1568  PC pc_raw;
1569  CHKERR KSPGetPC(ksp, &pc_raw);
1570  return SmartPetscObj<PC>(pc_raw,
1571  true); // bump reference
1572  };
1573  CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1574  CHKERR PetscFree(subksp);
1575  MoFEMFunctionReturn(0);
1576  };
1577 
1578  CHKERR set_ops();
1579  CHKERR set_pc();
1580  CHKERR set_assemble_elems();
1581 
1582  CHKERR TSSetUp(solver);
1583  CHKERR KSPSetUp(ksp);
1584  CHKERR set_diagonal_pc();
1585 
1586  } else {
1587  pip_mng->getOpBoundaryLhsPipeline().push_front(
1588  createOpSchurAssembleBegin());
1589  pip_mng->getOpBoundaryLhsPipeline().push_back(
1590  createOpSchurAssembleEnd({}, {}));
1591  pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1592  pip_mng->getOpDomainLhsPipeline().push_back(
1593  createOpSchurAssembleEnd({}, {}));
1594  }
1595 
1596  // fieldSplitIS.reset();
1597  // aoSchur.reset();
1598  MoFEMFunctionReturn(0);
1599 }
1600 
1601 boost::shared_ptr<SetUpSchur>
1602 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
1603  SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1604  SmartPetscObj<AO> ao_up) {
1605  return boost::shared_ptr<SetUpSchur>(
1606  new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
1607 }
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
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
Example::tsSolve
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:740
SetUpSchurImpl
Definition: test_broken_space.cpp:515
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:556
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
SetUpSchurImpl::aoSchur
SmartPetscObj< AO > aoSchur
Definition: plastic.cpp:1416
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
ContactOps
Definition: contact.cpp:99
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:528
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
Example::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:197
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:1404
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
Example::ScaledTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:201
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:1414
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:2010
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
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:2186
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:525
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:1319
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
SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1415
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()
Example::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:195
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[]
[TestOperators]
Definition: plastic.cpp:1317
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:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: test_broken_space.cpp:524
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
Example::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:196
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
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
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2181
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
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:1898
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:58
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::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2223
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
Example::testOperators
MoFEMErrorCode testOperators()
[Solve]
Definition: plastic.cpp:1214
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:512
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:1148
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:233
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
MoFEM::createSchurNestedMatrixStruture
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:1944
Example::ScaledTimeScale
Definition: plastic.cpp:199
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:214
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_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
Example::reactionFe
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:193
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:1009
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:404
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:768
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
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:1393
alpha_damping
double alpha_damping
Definition: plastic.cpp:141
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow