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