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 using SetPtsData = FieldEvaluatorInterface::SetPtsData;
63 
64 inline double iso_hardening_exp(double tau, double b_iso) {
65  return std::exp(
66  std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
67  -b_iso * tau));
68 }
69 
70 /**
71  * Isotropic hardening
72  */
73 inline double iso_hardening(double tau, double H, double Qinf, double b_iso,
74  double sigmaY) {
75  return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;
76 }
77 
78 inline double iso_hardening_dtau(double tau, double H, double Qinf,
79  double b_iso) {
80  auto r = [&](auto tau) {
81  return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);
82  };
83  constexpr double eps = 1e-12;
84  return std::max(r(tau), eps * r(0));
85 }
86 
87 /**
88  * Kinematic hardening
89  */
90 template <typename T, int DIM>
91 inline auto
92 kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,
93  double C1_k) {
94  FTensor::Index<'i', DIM> i;
95  FTensor::Index<'j', DIM> j;
96  FTensor::Tensor2_symmetric<double, DIM> t_alpha;
97  if (C1_k < std::numeric_limits<double>::epsilon()) {
98  t_alpha(i, j) = 0;
99  return t_alpha;
100  }
101  t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
102  return t_alpha;
103 }
104 
105 template <int DIM>
106 inline auto kinematic_hardening_dplastic_strain(double C1_k) {
107  FTensor::Index<'i', DIM> i;
108  FTensor::Index<'j', DIM> j;
109  FTensor::Index<'k', DIM> k;
110  FTensor::Index<'l', DIM> l;
111  FTensor::Ddg<double, DIM, DIM> t_diff;
112  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
113  t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
114  return t_diff;
115 }
116 
117 PetscBool is_large_strains = PETSC_TRUE; ///< Large strains
118 PetscBool set_timer = PETSC_FALSE; ///< Set timer
119 PetscBool do_eval_field = PETSC_FALSE; ///< Evaluate field
120 
121 int atom_test = 0; ///< Atom test
122 
123 double scale = 1.;
124 
125 double young_modulus = 206913; ///< Young modulus
126 double poisson_ratio = 0.29; ///< Poisson ratio
127 double sigmaY = 450; ///< Yield stress
128 double H = 129; ///< Hardening
129 double visH = 0; ///< Viscous hardening
130 double zeta = 5e-2; ///< Viscous hardening
131 double Qinf = 265; ///< Saturation yield stress
132 double b_iso = 16.93; ///< Saturation exponent
133 double C1_k = 0; ///< Kinematic hardening
134 
135 double cn0 = 1;
136 double cn1 = 1;
137 
138 int order = 2; ///< Order displacement
139 int tau_order = order - 2; ///< Order of tau files
140 int ep_order = order - 1; ///< Order of ep files
141 int geom_order = 2; ///< Order if fixed.
142 
143 PetscBool is_quasi_static = PETSC_TRUE;
144 double rho = 0.0;
145 double alpha_damping = 0;
146 
147 #include <HenckyOps.hpp>
148 #include <PlasticOps.hpp>
149 #include <PlasticNaturalBCs.hpp>
150 
151 #ifdef ADD_CONTACT
152  #ifdef PYTHON_SDF
153  #include <boost/python.hpp>
154  #include <boost/python/def.hpp>
155  #include <boost/python/numpy.hpp>
156 namespace bp = boost::python;
157 namespace np = boost::python::numpy;
158  #endif
159 
160 namespace ContactOps {
161 
162 double cn_contact = 0.1;
163 
164 }; // namespace ContactOps
165 
166  #include <ContactOps.hpp>
167 #endif // ADD_CONTACT
168 
169 using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
170 using OpDomainRhsBCs =
171  DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
172 using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
173 using OpBoundaryRhsBCs =
174  BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
175 using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
176 using OpBoundaryLhsBCs =
177  BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
178 
179 using namespace PlasticOps;
180 using namespace HenckyOps;
181 
182 namespace PlasticOps {
183 
184 template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
185 
186 template <> struct AddHOOps<2, 3, 3> {
187  AddHOOps() = delete;
188  static MoFEMErrorCode
189  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
190  std::vector<FieldSpace> space, std::string geom_field_name);
191 };
192 
193 template <> struct AddHOOps<1, 2, 2> {
194  AddHOOps() = delete;
195  static MoFEMErrorCode
196  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
197  std::vector<FieldSpace> space, std::string geom_field_name);
198 };
199 
200 template <> struct AddHOOps<3, 3, 3> {
201  AddHOOps() = delete;
202  static MoFEMErrorCode
203  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
204  std::vector<FieldSpace> space, std::string geom_field_name);
205 };
206 
207 template <> struct AddHOOps<2, 2, 2> {
208  AddHOOps() = delete;
209  static MoFEMErrorCode
210  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
211  std::vector<FieldSpace> space, std::string geom_field_name);
212 };
213 
214 } // namespace PlasticOps
215 
216 struct Example {
217 
218  Example(MoFEM::Interface &m_field) : mField(m_field) {}
219 
220  MoFEMErrorCode runProblem();
221 
222  enum { VOL, COUNT };
223  static inline std::array<double, 2> meshVolumeAndCount = {0, 0};
224 
225 private:
226  MoFEM::Interface &mField;
227 
228  MoFEMErrorCode setupProblem();
229  MoFEMErrorCode createCommonData();
230  MoFEMErrorCode bC();
231  MoFEMErrorCode OPs();
232  MoFEMErrorCode tsSolve();
233  MoFEMErrorCode testOperators();
234 
235  boost::shared_ptr<DomainEle> reactionFe;
236 
237  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
238  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
239  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
240 
241  struct ScaledTimeScale : public MoFEM::TimeScale {
242  using MoFEM::TimeScale::TimeScale;
243  double getScale(const double time) {
244  return scale * MoFEM::TimeScale::getScale(time);
245  };
246  };
247 
248 #ifdef ADD_CONTACT
249  #ifdef PYTHON_SDF
250  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
251  #endif
252 #endif // ADD_CONTACT
253 };
254 
255 //! [Run problem]
256 MoFEMErrorCode Example::runProblem() {
257  MoFEMFunctionBegin;
258  CHKERR createCommonData();
259  CHKERR setupProblem();
260  CHKERR bC();
261  CHKERR OPs();
262  PetscBool test_ops = PETSC_FALSE;
263  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_operators", &test_ops,
264  PETSC_NULL);
265  if (test_ops == PETSC_FALSE) {
266  CHKERR tsSolve();
267  } else {
268  CHKERR testOperators();
269  }
270  MoFEMFunctionReturn(0);
271 }
272 //! [Run problem]
273 
274 //! [Set up problem]
275 MoFEMErrorCode Example::setupProblem() {
276  MoFEMFunctionBegin;
277  Simple *simple = mField.getInterface<Simple>();
278 
279  Range domain_ents;
280  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
281  true);
282  auto get_ents_by_dim = [&](const auto dim) {
283  if (dim == SPACE_DIM) {
284  return domain_ents;
285  } else {
286  Range ents;
287  if (dim == 0)
288  CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
289  else
290  CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
291  return ents;
292  }
293  };
294 
295  auto get_base = [&]() {
296  auto domain_ents = get_ents_by_dim(SPACE_DIM);
297  if (domain_ents.empty())
298  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
299  const auto type = type_from_handle(domain_ents[0]);
300  switch (type) {
301  case MBQUAD:
302  return DEMKOWICZ_JACOBI_BASE;
303  case MBHEX:
304  return DEMKOWICZ_JACOBI_BASE;
305  case MBTRI:
306  return AINSWORTH_LEGENDRE_BASE;
307  case MBTET:
308  return AINSWORTH_LEGENDRE_BASE;
309  default:
310  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
311  }
312  return NOBASE;
313  };
314 
315  const auto base = get_base();
316  MOFEM_LOG("PLASTICITY", Sev::inform)
317  << "Base " << ApproximationBaseNames[base];
318 
319  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
320  CHKERR simple->addDomainField("EP", L2, base, size_symm);
321  CHKERR simple->addDomainField("TAU", L2, base, 1);
322  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
323 
324  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
325 
326  PetscBool order_edge = PETSC_FALSE;
327  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,
328  PETSC_NULL);
329  PetscBool order_face = PETSC_FALSE;
330  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
331  PETSC_NULL);
332  PetscBool order_volume = PETSC_FALSE;
333  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
334  PETSC_NULL);
335 
336  if (order_edge || order_face || order_volume) {
337 
338  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
339  ? "true"
340  : "false";
341  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
342  ? "true"
343  : "false";
344  MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
345  ? "true"
346  : "false";
347 
348  auto ents = get_ents_by_dim(0);
349  if (order_edge)
350  ents.merge(get_ents_by_dim(1));
351  if (order_face)
352  ents.merge(get_ents_by_dim(2));
353  if (order_volume)
354  ents.merge(get_ents_by_dim(3));
355  CHKERR simple->setFieldOrder("U", order, &ents);
356  } else {
357  CHKERR simple->setFieldOrder("U", order);
358  }
359  CHKERR simple->setFieldOrder("EP", ep_order);
360  CHKERR simple->setFieldOrder("TAU", tau_order);
361 
362  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
363 
364 #ifdef ADD_CONTACT
365  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
366  SPACE_DIM);
367  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
368  SPACE_DIM);
369 
370  auto get_skin = [&]() {
371  Range body_ents;
372  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
373  Skinner skin(&mField.get_moab());
374  Range skin_ents;
375  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
376  return skin_ents;
377  };
378 
379  auto filter_blocks = [&](auto skin) {
380  bool is_contact_block = true;
381  Range contact_range;
382  for (auto m :
383  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
384 
385  (boost::format("%s(.*)") % "CONTACT").str()
386 
387  ))
388 
389  ) {
390  is_contact_block =
391  true; ///< blocs interation is collective, so that is set irrespective
392  ///< if there are entities in given rank or not in the block
393  MOFEM_LOG("CONTACT", Sev::inform)
394  << "Find contact block set: " << m->getName();
395  auto meshset = m->getMeshset();
396  Range contact_meshset_range;
397  CHKERR mField.get_moab().get_entities_by_dimension(
398  meshset, SPACE_DIM - 1, contact_meshset_range, true);
399 
400  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
401  contact_meshset_range);
402  contact_range.merge(contact_meshset_range);
403  }
404  if (is_contact_block) {
405  MOFEM_LOG("SYNC", Sev::inform)
406  << "Nb entities in contact surface: " << contact_range.size();
407  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
408  skin = intersect(skin, contact_range);
409  }
410  return skin;
411  };
412 
413  auto filter_true_skin = [&](auto skin) {
414  Range boundary_ents;
415  ParallelComm *pcomm =
416  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
417  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
418  PSTATUS_NOT, -1, &boundary_ents);
419  return boundary_ents;
420  };
421 
422  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
423  CHKERR simple->setFieldOrder("SIGMA", 0);
424  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
425 #endif
426 
427  CHKERR simple->setUp();
428  CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
429 
430  auto project_ho_geometry = [&]() {
431  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
432  return mField.loop_dofs("GEOMETRY", ent_method);
433  };
434  PetscBool project_geometry = PETSC_TRUE;
435  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
436  &project_geometry, PETSC_NULL);
437  if (project_geometry) {
438  CHKERR project_ho_geometry();
439  }
440 
441  auto get_volume = [&]() {
442  using VolOp = DomainEle::UserDataOperator;
443  auto *op_ptr = new VolOp(NOSPACE, VolOp::OPSPACE);
444  std::array<double, 2> volume_and_count;
445  op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
446  EntityType type,
447  EntitiesFieldData::EntData &data) {
448  MoFEMFunctionBegin;
449  auto op_ptr = static_cast<VolOp *>(base_op_ptr);
450  volume_and_count[VOL] += op_ptr->getMeasure();
451  volume_and_count[COUNT] += 1;
452  // in necessary at integration over Gauss points.
453  MoFEMFunctionReturn(0);
454  };
455  volume_and_count = {0, 0};
456  auto fe = boost::make_shared<DomainEle>(mField);
457  fe->getOpPtrVector().push_back(op_ptr);
458 
459  auto dm = simple->getDM();
460  CHK_THROW_MESSAGE(
461  DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe),
462  "cac volume");
463  std::array<double, 2> tot_volume_and_count;
464  MPI_Allreduce(volume_and_count.data(), tot_volume_and_count.data(),
465  volume_and_count.size(), MPI_DOUBLE, MPI_SUM,
466  mField.get_comm());
467  return tot_volume_and_count;
468  };
469 
470  meshVolumeAndCount = get_volume();
471  MOFEM_LOG("PLASTICITY", Sev::inform)
472  << "Mesh volume " << meshVolumeAndCount[VOL] << " nb. of elements "
473  << meshVolumeAndCount[COUNT];
474 
475  MoFEMFunctionReturn(0);
476 }
477 //! [Set up problem]
478 
479 //! [Create common data]
480 MoFEMErrorCode Example::createCommonData() {
481  MoFEMFunctionBegin;
482 
483  auto get_command_line_parameters = [&]() {
484  MoFEMFunctionBegin;
485 
486  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
487  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
488  &young_modulus, PETSC_NULL);
489  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
490  &poisson_ratio, PETSC_NULL);
491  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
492  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
493  PETSC_NULL);
494  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
495  PETSC_NULL);
496  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
497  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
498  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
499  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
500  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
501  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);
502  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
503  &is_large_strains, PETSC_NULL);
504  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
505  PETSC_NULL);
506  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test,
507  PETSC_NULL);
508 
509  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
510  PetscBool tau_order_is_set; ///< true if tau order is set
511  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
512  &tau_order_is_set);
513  PetscBool ep_order_is_set; ///< true if tau order is set
514  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
515  &ep_order_is_set);
516  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
517  PETSC_NULL);
518 
519  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
520  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
521  &alpha_damping, PETSC_NULL);
522 
523  MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
524  MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
525  MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
526  MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
527  MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
528  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
529  MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
530  MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
531  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
532  MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
533  MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
534 
535  if (tau_order_is_set == PETSC_FALSE)
536  tau_order = order - 2;
537  if (ep_order_is_set == PETSC_FALSE)
538  ep_order = order - 1;
539 
540  MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
541  MOFEM_LOG("PLASTICITY", Sev::inform)
542  << "Ep approximation order " << ep_order;
543  MOFEM_LOG("PLASTICITY", Sev::inform)
544  << "Tau approximation order " << tau_order;
545  MOFEM_LOG("PLASTICITY", Sev::inform)
546  << "Geometry approximation order " << geom_order;
547 
548  MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
549  MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
550 
551  PetscBool is_scale = PETSC_TRUE;
552  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
553  PETSC_NULL);
554  if (is_scale) {
555  scale /= young_modulus;
556  }
557 
558  MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
559 
560 #ifdef ADD_CONTACT
561  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
562  &ContactOps::cn_contact, PETSC_NULL);
563  MOFEM_LOG("CONTACT", Sev::inform)
564  << "cn_contact " << ContactOps::cn_contact;
565 #endif // ADD_CONTACT
566 
567  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
568  &is_quasi_static, PETSC_NULL);
569  MOFEM_LOG("PLASTICITY", Sev::inform)
570  << "Is quasi static: " << (is_quasi_static ? "true" : "false");
571 
572  MoFEMFunctionReturn(0);
573  };
574 
575  CHKERR get_command_line_parameters();
576 
577 #ifdef ADD_CONTACT
578  #ifdef PYTHON_SDF
579  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
580  CHKERR sdfPythonPtr->sdfInit("sdf.py");
581  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
582  #endif
583 #endif // ADD_CONTACT
584 
585  MoFEMFunctionReturn(0);
586 }
587 //! [Create common data]
588 
589 //! [Boundary condition]
590 MoFEMErrorCode Example::bC() {
591  MoFEMFunctionBegin;
592 
593  auto simple = mField.getInterface<Simple>();
594  auto bc_mng = mField.getInterface<BcManager>();
595 
596  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
597  "U", 0, 0);
598  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
599  "U", 1, 1);
600  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
601  "U", 2, 2);
602  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
603  "REMOVE_ALL", "U", 0, 3);
604 
605 #ifdef ADD_CONTACT
606  for (auto b : {"FIX_X", "REMOVE_X"})
607  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
608  "SIGMA", 0, 0, false, true);
609  for (auto b : {"FIX_Y", "REMOVE_Y"})
610  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
611  "SIGMA", 1, 1, false, true);
612  for (auto b : {"FIX_Z", "REMOVE_Z"})
613  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
614  "SIGMA", 2, 2, false, true);
615  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
616  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
617  "SIGMA", 0, 3, false, true);
618  CHKERR bc_mng->removeBlockDOFsOnEntities(
619  simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
620 #endif
621 
622  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
623  simple->getProblemName(), "U");
624 
625  auto &bc_map = bc_mng->getBcMapByBlockName();
626  for (auto bc : bc_map)
627  MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
628 
629  MoFEMFunctionReturn(0);
630 }
631 //! [Boundary condition]
632 
633 //! [Push operators to pipeline]
634 MoFEMErrorCode Example::OPs() {
635  MoFEMFunctionBegin;
636  auto pip_mng = mField.getInterface<PipelineManager>();
637 
638  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
639 
640  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
641 
642  auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
643  MoFEMFunctionBegin;
644 
645  CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
646  pip, {HDIV}, "GEOMETRY");
647  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
648 
649  // Add Natural BCs to LHS
650  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
651  pip, mField, "U", Sev::inform);
652 
653 #ifdef ADD_CONTACT
654  auto simple = mField.getInterface<Simple>();
655  CHKERR
656  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
657  pip, "SIGMA", "U");
658  CHKERR
659  ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
660  mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
661  vol_rule);
662 #endif // ADD_CONTACT
663 
664  MoFEMFunctionReturn(0);
665  };
666 
667  auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
668  MoFEMFunctionBegin;
669 
670  CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
671  pip, {HDIV}, "GEOMETRY");
672  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
673 
674  // Add Natural BCs to RHS
675  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
676  pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
677 
678 #ifdef ADD_CONTACT
679  CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
680  pip, "SIGMA", "U");
681 #endif // ADD_CONTACT
682 
683  MoFEMFunctionReturn(0);
684  };
685 
686  auto add_domain_ops_lhs = [this](auto &pip) {
687  MoFEMFunctionBegin;
688  CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
689  pip, {H1, HDIV}, "GEOMETRY");
690 
691  if (is_quasi_static == PETSC_FALSE) {
692 
693  //! [Only used for dynamics]
694  using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
695  GAUSS>::OpMass<1, SPACE_DIM>;
696  //! [Only used for dynamics]
697 
698  auto get_inertia_and_mass_damping = [this](const double, const double,
699  const double) {
700  auto *pip = mField.getInterface<PipelineManager>();
701  auto &fe_domain_lhs = pip->getDomainLhsFE();
702  return (rho / scale) * fe_domain_lhs->ts_aa +
703  (alpha_damping / scale) * fe_domain_lhs->ts_a;
704  };
705  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
706  }
707 
708  CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
709  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
710 
711  MoFEMFunctionReturn(0);
712  };
713 
714  auto add_domain_ops_rhs = [this](auto &pip) {
715  MoFEMFunctionBegin;
716 
717  CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
718  pip, {H1, HDIV}, "GEOMETRY");
719 
720  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
721  pip, mField, "U",
722  {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
723  Sev::inform);
724 
725  // only in case of dynamics
726  if (is_quasi_static == PETSC_FALSE) {
727 
728  //! [Only used for dynamics]
729  using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
730  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
731  //! [Only used for dynamics]
732 
733  auto mat_acceleration = boost::make_shared<MatrixDouble>();
734  pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
735  "U", mat_acceleration));
736  pip.push_back(
737  new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
738  return rho / scale;
739  }));
740  if (alpha_damping > 0) {
741  auto mat_velocity = boost::make_shared<MatrixDouble>();
742  pip.push_back(
743  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
744  pip.push_back(
745  new OpInertiaForce("U", mat_velocity, [](double, double, double) {
746  return alpha_damping / scale;
747  }));
748  }
749  }
750 
751  CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
752  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
753 
754 #ifdef ADD_CONTACT
755  CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
756  pip, "SIGMA", "U");
757 #endif // ADD_CONTACT
758 
759  MoFEMFunctionReturn(0);
760  };
761 
762  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
763  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
764 
765  // Boundary
766  CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
767  CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
768 
769  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
770  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
771 
772  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
773  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
774 
775  auto create_reaction_pipeline = [&](auto &pip) {
776  MoFEMFunctionBegin;
777  CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
778  pip, {H1}, "GEOMETRY");
779  CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
780  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
781  MoFEMFunctionReturn(0);
782  };
783 
784  reactionFe = boost::make_shared<DomainEle>(mField);
785  reactionFe->getRuleHook = vol_rule;
786  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
787  reactionFe->postProcessHook =
788  EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
789 
790  MoFEMFunctionReturn(0);
791 }
792 //! [Push operators to pipeline]
793 
794 //! [Solve]
795 struct SetUpSchur {
796 
797  /**
798  * @brief Create data structure for handling Schur complement
799  *
800  * @param m_field
801  * @param sub_dm Schur complement sub dm
802  * @param field_split_it IS of Schur block
803  * @param ao_map AO map from sub dm to main problem
804  * @return boost::shared_ptr<SetUpSchur>
805  */
806  static boost::shared_ptr<SetUpSchur> createSetUpSchur(
807 
808  MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
809  SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
810 
811  );
812  virtual MoFEMErrorCode setUp(TS solver) = 0;
813 
814 protected:
815  SetUpSchur() = default;
816 };
817 
818 MoFEMErrorCode Example::tsSolve() {
819  MoFEMFunctionBegin;
820 
821  Simple *simple = mField.getInterface<Simple>();
822  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
823  ISManager *is_manager = mField.getInterface<ISManager>();
824 
825  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
826 
827  auto set_section_monitor = [&](auto solver) {
828  MoFEMFunctionBegin;
829  SNES snes;
830  CHKERR TSGetSNES(solver, &snes);
831  CHKERR SNESMonitorSet(snes,
832  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
833  void *))MoFEMSNESMonitorFields,
834  (void *)(snes_ctx_ptr.get()), nullptr);
835  MoFEMFunctionReturn(0);
836  };
837 
838  auto create_post_process_elements = [&]() {
839  auto push_vol_ops = [this](auto &pip) {
840  CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
841  pip, {H1, HDIV}, "GEOMETRY");
842 
843  auto [common_plastic_ptr, common_hencky_ptr] =
844  PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
845  mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
846 
847  if (common_hencky_ptr) {
848  if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
849  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
850  }
851 
852  return std::make_pair(common_plastic_ptr, common_hencky_ptr);
853  };
854 
855  auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
856  MoFEMFunctionBegin;
857 
858  auto &pip = pp_fe->getOpPtrVector();
859 
860  auto [common_plastic_ptr, common_hencky_ptr] = p;
861 
862  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
863 
864  auto x_ptr = boost::make_shared<MatrixDouble>();
865  pip.push_back(
866  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
867  auto u_ptr = boost::make_shared<MatrixDouble>();
868  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
869 
870  if (is_large_strains) {
871 
872  pip.push_back(
873 
874  new OpPPMap(
875 
876  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
877 
878  {{"PLASTIC_SURFACE",
879  common_plastic_ptr->getPlasticSurfacePtr()},
880  {"PLASTIC_MULTIPLIER",
881  common_plastic_ptr->getPlasticTauPtr()}},
882 
883  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
884 
885  {{"GRAD", common_hencky_ptr->matGradPtr},
886  {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
887 
888  {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
889  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
890  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
891 
892  )
893 
894  );
895 
896  } else {
897 
898  pip.push_back(
899 
900  new OpPPMap(
901 
902  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
903 
904  {{"PLASTIC_SURFACE",
905  common_plastic_ptr->getPlasticSurfacePtr()},
906  {"PLASTIC_MULTIPLIER",
907  common_plastic_ptr->getPlasticTauPtr()}},
908 
909  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
910 
911  {},
912 
913  {{"STRAIN", common_plastic_ptr->mStrainPtr},
914  {"STRESS", common_plastic_ptr->mStressPtr},
915  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
916  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
917 
918  )
919 
920  );
921  }
922 
923  MoFEMFunctionReturn(0);
924  };
925 
926  PetscBool post_proc_vol;
927  PetscBool post_proc_skin;
928 
929  if constexpr (SPACE_DIM == 2) {
930  post_proc_vol = PETSC_TRUE;
931  post_proc_skin = PETSC_FALSE;
932  } else {
933  post_proc_vol = PETSC_FALSE;
934  post_proc_skin = PETSC_TRUE;
935  }
936  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol", &post_proc_vol,
937  PETSC_NULL);
938  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
939  &post_proc_skin, PETSC_NULL);
940 
941  auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
942  post_proc_vol]() {
943  if (post_proc_vol == PETSC_FALSE)
944  return boost::shared_ptr<PostProcEle>();
945  auto pp_fe = boost::make_shared<PostProcEle>(mField);
946  CHK_MOAB_THROW(
947  push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
948  "push_vol_post_proc_ops");
949  return pp_fe;
950  };
951 
952  auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
953  post_proc_skin]() {
954  if (post_proc_skin == PETSC_FALSE)
955  return boost::shared_ptr<SkinPostProcEle>();
956 
957  auto simple = mField.getInterface<Simple>();
958  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
959  auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
960  SPACE_DIM, Sev::verbose);
961  pp_fe->getOpPtrVector().push_back(op_side);
962  CHK_MOAB_THROW(push_vol_post_proc_ops(
963  pp_fe, push_vol_ops(op_side->getOpPtrVector())),
964  "push_vol_post_proc_ops");
965  return pp_fe;
966  };
967 
968  return std::make_pair(vol_post_proc(), skin_post_proc());
969  };
970 
971  auto scatter_create = [&](auto D, auto coeff) {
972  SmartPetscObj<IS> is;
973  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
974  ROW, "U", coeff, coeff, is);
975  int loc_size;
976  CHKERR ISGetLocalSize(is, &loc_size);
977  Vec v;
978  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
979  VecScatter scatter;
980  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
981  return std::make_tuple(SmartPetscObj<Vec>(v),
982  SmartPetscObj<VecScatter>(scatter));
983  };
984 
985  boost::shared_ptr<SetPtsData> field_eval_data;
986  boost::shared_ptr<MatrixDouble> u_field_ptr;
987 
988  std::array<double, SPACE_DIM> field_eval_coords;
989  int coords_dim = SPACE_DIM;
990  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
991  field_eval_coords.data(), &coords_dim,
992  &do_eval_field);
993 
994  boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
995  scalar_field_ptrs = boost::make_shared<
996  std::map<std::string, boost::shared_ptr<VectorDouble>>>();
997  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
998  vector_field_ptrs = boost::make_shared<
999  std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1000  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1001  sym_tensor_field_ptrs = boost::make_shared<
1002  std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1003  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1004  tensor_field_ptrs = boost::make_shared<
1005  std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1006 
1007  if (do_eval_field) {
1008  auto u_field_ptr = boost::make_shared<MatrixDouble>();
1009  field_eval_data =
1010  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1011 
1012  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
1013  field_eval_data, simple->getDomainFEName());
1014 
1015  field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1016  auto no_rule = [](int, int, int) { return -1; };
1017  auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1018  field_eval_fe_ptr->getRuleHook = no_rule;
1019 
1020  CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1021  field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1022 
1023  auto [common_plastic_ptr, common_hencky_ptr] =
1024  PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1025  mField, "MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U",
1026  "EP", "TAU", 1., Sev::inform);
1027 
1028  field_eval_fe_ptr->getOpPtrVector().push_back(
1029  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_field_ptr));
1030 
1031  if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1032  if (is_large_strains) {
1033  scalar_field_ptrs->insert(
1034  {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1035  scalar_field_ptrs->insert(
1036  {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1037  vector_field_ptrs->insert({"U", u_field_ptr});
1038  sym_tensor_field_ptrs->insert(
1039  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1040  sym_tensor_field_ptrs->insert(
1041  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1042  sym_tensor_field_ptrs->insert(
1043  {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1044  tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1045  tensor_field_ptrs->insert(
1046  {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1047  } else {
1048  scalar_field_ptrs->insert(
1049  {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1050  scalar_field_ptrs->insert(
1051  {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1052  vector_field_ptrs->insert({"U", u_field_ptr});
1053  sym_tensor_field_ptrs->insert(
1054  {"STRAIN", common_plastic_ptr->mStrainPtr});
1055  sym_tensor_field_ptrs->insert(
1056  {"STRESS", common_plastic_ptr->mStressPtr});
1057  sym_tensor_field_ptrs->insert(
1058  {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1059  sym_tensor_field_ptrs->insert(
1060  {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1061  }
1062  }
1063  }
1064 
1065  auto test_monitor_ptr = boost::make_shared<FEMethod>();
1066 
1067  auto set_time_monitor = [&](auto dm, auto solver) {
1068  MoFEMFunctionBegin;
1069  boost::shared_ptr<Monitor<SPACE_DIM>> monitor_ptr(new Monitor<SPACE_DIM>(
1070  dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
1071  uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
1072  vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1073  boost::shared_ptr<ForcesAndSourcesCore> null;
1074 
1075  test_monitor_ptr->postProcessHook = [&]() {
1076  MoFEMFunctionBegin;
1077 
1078  if (atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1079  test_monitor_ptr->ts_step == 25) {
1080 
1081  if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1082  auto t_tau =
1083  getFTensor0FromVec(*scalar_field_ptrs->at("PLASTIC_MULTIPLIER"));
1084  MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point tau: " << t_tau;
1085 
1086  if (atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1087  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1088  "atom test %d failed: wrong plastic multiplier value",
1089  atom_test);
1090  }
1091  }
1092 
1093  if (vector_field_ptrs->at("U")->size1()) {
1094  FTensor::Index<'i', SPACE_DIM> i;
1095  auto t_disp =
1096  getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1097  MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point U: " << t_disp;
1098 
1099  if (atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1100  fabs(t_disp(1) + 0.0526736) > 1e-5 || fabs(t_disp(2)) > 1e-5) {
1101  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1102  "atom test %d failed: wrong displacement value",
1103  atom_test);
1104  }
1105  }
1106 
1107  if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1108  auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1109  *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1110  MOFEM_LOG("PlasticSync", Sev::inform)
1111  << "Eval point EP: " << t_plastic_strain;
1112 
1113  if (atom_test == 1 &&
1114  fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1115  fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1116  fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1117  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1118  "atom test %d failed: wrong plastic strain value",
1119  atom_test);
1120  }
1121  }
1122 
1123  if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1124  auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1125  *tensor_field_ptrs->at("FIRST_PIOLA"));
1126  MOFEM_LOG("PlasticSync", Sev::inform)
1127  << "Eval point Piola stress: " << t_piola_stress;
1128 
1129  if (atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1130  t_piola_stress(0, 0)) > 1e-5 ||
1131  fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1132  fabs(t_piola_stress(1, 1)) >
1133  1e-5) {
1134  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1135  "atom test %d failed: wrong Piola stress value",
1136  atom_test);
1137  }
1138  }
1139  }
1140 
1141  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
1142  MoFEMFunctionReturn(0);
1143  };
1144 
1145  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
1146  monitor_ptr, null, test_monitor_ptr);
1147 
1148  MoFEMFunctionReturn(0);
1149  };
1150 
1151  auto set_schur_pc = [&](auto solver,
1152  boost::shared_ptr<SetUpSchur> &schur_ptr) {
1153  MoFEMFunctionBeginHot;
1154 
1155  auto name_prb = simple->getProblemName();
1156 
1157  // create sub dm for Schur complement
1158  auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
1159  SmartPetscObj<DM> &dm_sub) {
1160  MoFEMFunctionBegin;
1161  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1162  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
1163  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1164  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1165  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1166  for (auto f : {"U"}) {
1167  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1168  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1169  }
1170  CHKERR DMSetUp(dm_sub);
1171 
1172  MoFEMFunctionReturn(0);
1173  };
1174 
1175  auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
1176  SmartPetscObj<DM> &dm_sub) {
1177  MoFEMFunctionBegin;
1178  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1179  CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
1180  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1181  CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1182  CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1183 #ifdef ADD_CONTACT
1184  for (auto f : {"SIGMA", "EP", "TAU"}) {
1185  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1186  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1187  }
1188 #else
1189  for (auto f : {"EP", "TAU"}) {
1190  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1191  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1192  }
1193 #endif
1194  CHKERR DMSetUp(dm_sub);
1195  MoFEMFunctionReturn(0);
1196  };
1197 
1198  // Create nested (sub BC) Schur DM
1199  if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1200 
1201  SmartPetscObj<DM> dm_schur;
1202  CHKERR create_schur_dm(simple->getDM(), dm_schur);
1203  SmartPetscObj<DM> dm_block;
1204  CHKERR create_block_dm(simple->getDM(), dm_block);
1205 
1206 #ifdef ADD_CONTACT
1207 
1208  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1209  auto block_mat_data = createBlockMatStructure(
1210  simple->getDM(),
1211 
1212  {
1213 
1214  {simple->getDomainFEName(),
1215 
1216  {{"U", "U"},
1217  {"SIGMA", "SIGMA"},
1218  {"U", "SIGMA"},
1219  {"SIGMA", "U"},
1220  {"EP", "EP"},
1221  {"TAU", "TAU"},
1222  {"U", "EP"},
1223  {"EP", "U"},
1224  {"EP", "TAU"},
1225  {"TAU", "EP"},
1226  {"TAU", "U"}
1227 
1228  }},
1229 
1230  {simple->getBoundaryFEName(),
1231 
1232  {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1233 
1234  }}
1235 
1236  }
1237 
1238  );
1239 
1240  return createSchurNestedMatrixStruture(
1241 
1242  {dm_schur, dm_block}, block_mat_data,
1243 
1244  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1245 
1246  );
1247  };
1248 
1249 #else
1250 
1251  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1252  auto block_mat_data =
1253  createBlockMatStructure(simple->getDM(),
1254 
1255  {{simple->getDomainFEName(),
1256 
1257  {{"U", "U"},
1258  {"EP", "EP"},
1259  {"TAU", "TAU"},
1260  {"U", "EP"},
1261  {"EP", "U"},
1262  {"EP", "TAU"},
1263  {"TAU", "U"},
1264  {"TAU", "EP"}
1265 
1266  }}}
1267 
1268  );
1269 
1270  return createSchurNestedMatrixStruture(
1271 
1272  {dm_schur, dm_block}, block_mat_data,
1273 
1274  {"EP", "TAU"}, {nullptr, nullptr}, false
1275 
1276  );
1277  };
1278 
1279 #endif
1280 
1281  auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1282  CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1283 
1284  auto block_is = getDMSubData(dm_block)->getSmartRowIs();
1285  auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
1286 
1287  // Indices has to be map fro very to level, while assembling Schur
1288  // complement.
1289  schur_ptr =
1290  SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
1291  CHKERR schur_ptr->setUp(solver);
1292  }
1293 
1294  MoFEMFunctionReturnHot(0);
1295  };
1296 
1297  auto dm = simple->getDM();
1298  auto D = createDMVector(dm);
1299  auto DD = vectorDuplicate(D);
1300  uXScatter = scatter_create(D, 0);
1301  uYScatter = scatter_create(D, 1);
1302  if constexpr (SPACE_DIM == 3)
1303  uZScatter = scatter_create(D, 2);
1304 
1305  auto create_solver = [pip_mng]() {
1306  if (is_quasi_static == PETSC_TRUE)
1307  return pip_mng->createTSIM();
1308  else
1309  return pip_mng->createTSIM2();
1310  };
1311 
1312  auto solver = create_solver();
1313 
1314  auto active_pre_lhs = []() {
1315  MoFEMFunctionBegin;
1316  std::fill(PlasticOps::CommonData::activityData.begin(),
1317  PlasticOps::CommonData::activityData.end(), 0);
1318  MoFEMFunctionReturn(0);
1319  };
1320 
1321  auto active_post_lhs = [&]() {
1322  MoFEMFunctionBegin;
1323  auto get_iter = [&]() {
1324  SNES snes;
1325  CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1326  int iter;
1327  CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1328  "Can not get iter");
1329  return iter;
1330  };
1331 
1332  auto iter = get_iter();
1333  if (iter >= 0) {
1334 
1335  std::array<int, 5> activity_data;
1336  std::fill(activity_data.begin(), activity_data.end(), 0);
1337  MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1338  activity_data.data(), activity_data.size(), MPI_INT,
1339  MPI_SUM, mField.get_comm());
1340 
1341  int &active_points = activity_data[0];
1342  int &avtive_full_elems = activity_data[1];
1343  int &avtive_elems = activity_data[2];
1344  int &nb_points = activity_data[3];
1345  int &nb_elements = activity_data[4];
1346 
1347  if (nb_points) {
1348 
1349  double proc_nb_points =
1350  100 * static_cast<double>(active_points) / nb_points;
1351  double proc_nb_active =
1352  100 * static_cast<double>(avtive_elems) / nb_elements;
1353  double proc_nb_full_active = 100;
1354  if (avtive_elems)
1355  proc_nb_full_active =
1356  100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1357 
1358  MOFEM_LOG_C("PLASTICITY", Sev::inform,
1359  "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1360  "elements %d "
1361  "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1362  iter, nb_points, active_points, proc_nb_points,
1363  avtive_elems, proc_nb_active, avtive_full_elems,
1364  proc_nb_full_active, iter);
1365  }
1366  }
1367 
1368  MoFEMFunctionReturn(0);
1369  };
1370 
1371  auto add_active_dofs_elem = [&](auto dm) {
1372  MoFEMFunctionBegin;
1373  auto fe_pre_proc = boost::make_shared<FEMethod>();
1374  fe_pre_proc->preProcessHook = active_pre_lhs;
1375  auto fe_post_proc = boost::make_shared<FEMethod>();
1376  fe_post_proc->postProcessHook = active_post_lhs;
1377  auto ts_ctx_ptr = getDMTsCtx(dm);
1378  ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1379  ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1380  MoFEMFunctionReturn(0);
1381  };
1382 
1383  auto set_essential_bc = [&](auto dm, auto solver) {
1384  MoFEMFunctionBegin;
1385  // This is low level pushing finite elements (pipelines) to solver
1386 
1387  auto pre_proc_ptr = boost::make_shared<FEMethod>();
1388  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1389  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1390 
1391  // Add boundary condition scaling
1392  auto disp_time_scale = boost::make_shared<TimeScale>();
1393 
1394  auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1395  EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1396  {disp_time_scale}, false);
1397  return hook;
1398  };
1399  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1400 
1401  auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1402  MoFEMFunctionBegin;
1403  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
1404  mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1405  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
1406  mField, post_proc_rhs_ptr, 1.)();
1407  MoFEMFunctionReturn(0);
1408  };
1409  auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1410  return EssentialPostProcLhs<DisplacementCubitBcData>(
1411  mField, post_proc_lhs_ptr, 1.);
1412  };
1413  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1414 
1415  auto ts_ctx_ptr = getDMTsCtx(dm);
1416  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1417  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1418  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1419  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1420  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1421 
1422  MoFEMFunctionReturn(0);
1423  };
1424 
1425  if (is_quasi_static == PETSC_TRUE) {
1426  CHKERR TSSetSolution(solver, D);
1427  } else {
1428  CHKERR TS2SetSolution(solver, D, DD);
1429  }
1430 
1431  CHKERR set_section_monitor(solver);
1432  CHKERR set_time_monitor(dm, solver);
1433  CHKERR TSSetFromOptions(solver);
1434  CHKERR TSSetUp(solver);
1435 
1436  CHKERR add_active_dofs_elem(dm);
1437  boost::shared_ptr<SetUpSchur> schur_ptr;
1438  CHKERR set_schur_pc(solver, schur_ptr);
1439  CHKERR set_essential_bc(dm, solver);
1440 
1441  MOFEM_LOG_CHANNEL("TIMER");
1442  MOFEM_LOG_TAG("TIMER", "timer");
1443  if (set_timer)
1444  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1445  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1446  CHKERR TSSetUp(solver);
1447  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1448  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1449  CHKERR TSSolve(solver, NULL);
1450  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1451 
1452  MoFEMFunctionReturn(0);
1453 }
1454 //! [Solve]
1455 
1456 //! [TestOperators]
1457 MoFEMErrorCode Example::testOperators() {
1458  MoFEMFunctionBegin;
1459 
1460  // get operators tester
1461  auto simple = mField.getInterface<Simple>();
1462  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1463  // OperatorsTester
1464  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1465  // pipeline manager
1466 
1467  constexpr double eps = 1e-9;
1468 
1469  auto x = opt->setRandomFields(simple->getDM(), {
1470 
1471  {"U", {-1e-6, 1e-6}},
1472 
1473  {"EP", {-1e-6, 1e-6}},
1474 
1475  {"TAU", {0, 1e-4}}
1476 
1477  });
1478 
1479  auto dot_x_plastic_active =
1480  opt->setRandomFields(simple->getDM(), {
1481 
1482  {"U", {-1, 1}},
1483 
1484  {"EP", {-1, 1}},
1485 
1486  {"TAU", {0.1, 1}}
1487 
1488  });
1489  auto diff_x_plastic_active =
1490  opt->setRandomFields(simple->getDM(), {
1491 
1492  {"U", {-1, 1}},
1493 
1494  {"EP", {-1, 1}},
1495 
1496  {"TAU", {-1, 1}}
1497 
1498  });
1499 
1500  auto dot_x_elastic =
1501  opt->setRandomFields(simple->getDM(), {
1502 
1503  {"U", {-1, 1}},
1504 
1505  {"EP", {-1, 1}},
1506 
1507  {"TAU", {-1, -0.1}}
1508 
1509  });
1510  auto diff_x_elastic =
1511  opt->setRandomFields(simple->getDM(), {
1512 
1513  {"U", {-1, 1}},
1514 
1515  {"EP", {-1, 1}},
1516 
1517  {"TAU", {-1, 1}}
1518 
1519  });
1520 
1521  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
1522  auto dot_x, auto diff_x) {
1523  MoFEMFunctionBegin;
1524 
1525  auto diff_res = opt->checkCentralFiniteDifference(
1526  simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1527  SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
1528 
1529  // Calculate norm of difference between directional derivative calculated
1530  // from finite difference, and tangent matrix.
1531  double fnorm;
1532  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1533  MOFEM_LOG_C("PLASTICITY", Sev::inform,
1534  "Test consistency of tangent matrix %3.4e", fnorm);
1535 
1536  constexpr double err = 1e-5;
1537  if (fnorm > err)
1538  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1539  "Norm of directional derivative too large err = %3.4e", fnorm);
1540 
1541  MoFEMFunctionReturn(0);
1542  };
1543 
1544  MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
1545  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1546  pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
1547 
1548  MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
1549  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1550  pip->getDomainRhsFE(), dot_x_plastic_active,
1551  diff_x_plastic_active);
1552 
1553  MoFEMFunctionReturn(0);
1554 };
1555 
1556 //! [TestOperators]
1557 
1558 static char help[] = "...\n\n";
1559 
1560 int main(int argc, char *argv[]) {
1561 
1562 #ifdef ADD_CONTACT
1563  #ifdef PYTHON_SDF
1564  Py_Initialize();
1565  np::initialize();
1566  #endif
1567 #endif // ADD_CONTACT
1568 
1569  // Initialisation of MoFEM/PETSc and MOAB data structures
1570  const char param_file[] = "param_file.petsc";
1571  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1572 
1573  // Add logging channel for example
1574  auto core_log = logging::core::get();
1575  core_log->add_sink(
1576  LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
1577  core_log->add_sink(
1578  LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
1579  LogManager::setLog("PLASTICITY");
1580  MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
1581 
1582 #ifdef ADD_CONTACT
1583  core_log->add_sink(
1584  LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
1585  LogManager::setLog("CONTACT");
1586  MOFEM_LOG_TAG("CONTACT", "Contact");
1587 #endif // ADD_CONTACT
1588 
1589  core_log->add_sink(
1590  LogManager::createSink(LogManager::getStrmSync(), "PlasticSync"));
1591  LogManager::setLog("PlasticSync");
1592  MOFEM_LOG_TAG("PlasticSync", "PlasticSync");
1593 
1594  try {
1595 
1596  //! [Register MoFEM discrete manager in PETSc]
1597  DMType dm_name = "DMMOFEM";
1598  CHKERR DMRegister_MoFEM(dm_name);
1599  //! [Register MoFEM discrete manager in PETSc
1600 
1601  //! [Create MoAB]
1602  moab::Core mb_instance; ///< mesh database
1603  moab::Interface &moab = mb_instance; ///< mesh database interface
1604  //! [Create MoAB]
1605 
1606  //! [Create MoFEM]
1607  MoFEM::Core core(moab); ///< finite element database
1608  MoFEM::Interface &m_field = core; ///< finite element database interface
1609  //! [Create MoFEM]
1610 
1611  //! [Load mesh]
1612  Simple *simple = m_field.getInterface<Simple>();
1613  CHKERR simple->getOptions();
1614  CHKERR simple->loadFile();
1615  //! [Load mesh]
1616 
1617  //! [Example]
1618  Example ex(m_field);
1619  CHKERR ex.runProblem();
1620  //! [Example]
1621  }
1622  CATCH_ERRORS;
1623 
1624  CHKERR MoFEM::Core::Finalize();
1625 
1626 #ifdef ADD_CONTACT
1627  #ifdef PYTHON_SDF
1628  if (Py_FinalizeEx() < 0) {
1629  exit(120);
1630  }
1631  #endif
1632 #endif // ADD_CONTACT
1633 
1634  return 0;
1635 }
1636 
1637 struct SetUpSchurImpl : public SetUpSchur {
1638 
1639  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
1640  SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
1641  : SetUpSchur(), mField(m_field), subDM(sub_dm),
1642  fieldSplitIS(field_split_is), aoSchur(ao_up) {
1643  if (S) {
1644  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1645  "Is expected that schur matrix is not "
1646  "allocated. This is "
1647  "possible only is PC is set up twice");
1648  }
1649  }
1650  virtual ~SetUpSchurImpl() { S.reset(); }
1651 
1652  MoFEMErrorCode setUp(TS solver);
1653  MoFEMErrorCode preProc();
1654  MoFEMErrorCode postProc();
1655 
1656 private:
1657  SmartPetscObj<Mat> S;
1658 
1659  MoFEM::Interface &mField;
1660  SmartPetscObj<DM> subDM; ///< field split sub dm
1661  SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
1662  SmartPetscObj<AO> aoSchur; ///> main DM to subDM
1663 };
1664 
1665 MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
1666  MoFEMFunctionBegin;
1667  auto simple = mField.getInterface<Simple>();
1668  auto pip_mng = mField.getInterface<PipelineManager>();
1669 
1670  SNES snes;
1671  CHKERR TSGetSNES(solver, &snes);
1672  KSP ksp;
1673  CHKERR SNESGetKSP(snes, &ksp);
1674  CHKERR KSPSetFromOptions(ksp);
1675 
1676  PC pc;
1677  CHKERR KSPGetPC(ksp, &pc);
1678  PetscBool is_pcfs = PETSC_FALSE;
1679  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1680  if (is_pcfs) {
1681  if (S) {
1682  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1683  "Is expected that schur matrix is not "
1684  "allocated. This is "
1685  "possible only is PC is set up twice");
1686  }
1687 
1688  S = createDMMatrix(subDM);
1689  CHKERR MatSetBlockSize(S, SPACE_DIM);
1690 
1691  // Set DM to use shell block matrix
1692  DM solver_dm;
1693  CHKERR TSGetDM(solver, &solver_dm);
1694  CHKERR DMSetMatType(solver_dm, MATSHELL);
1695 
1696  auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1697  auto A = createDMBlockMat(simple->getDM());
1698  auto P = createDMNestSchurMat(simple->getDM());
1699 
1700  if (is_quasi_static == PETSC_TRUE) {
1701  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1702  Mat A, Mat B, void *ctx) {
1703  return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1704  };
1705  CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1706  } else {
1707  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1708  PetscReal a, PetscReal aa, Mat A, Mat B,
1709  void *ctx) {
1710  return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1711  };
1712  CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1713  }
1714  CHKERR KSPSetOperators(ksp, A, P);
1715 
1716  auto set_ops = [&]() {
1717  MoFEMFunctionBegin;
1718  auto pip_mng = mField.getInterface<PipelineManager>();
1719 
1720 #ifndef ADD_CONTACT
1721  // Boundary
1722  pip_mng->getOpBoundaryLhsPipeline().push_front(
1723  createOpSchurAssembleBegin());
1724  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1725 
1726  {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1727 
1728  ));
1729  // Domain
1730  pip_mng->getOpDomainLhsPipeline().push_front(
1731  createOpSchurAssembleBegin());
1732  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1733 
1734  {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1735 
1736  ));
1737 #else
1738 
1739  double eps_stab = 1e-4;
1740  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1741  PETSC_NULL);
1742 
1743  using B = FormsIntegrators<BoundaryEleOp>::Assembly<
1744  BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
1745  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1746 
1747  // Boundary
1748  pip_mng->getOpBoundaryLhsPipeline().push_front(
1749  createOpSchurAssembleBegin());
1750  pip_mng->getOpBoundaryLhsPipeline().push_back(
1751  new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1752  return eps_stab;
1753  }));
1754  pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1755 
1756  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1757  false, false
1758 
1759  ));
1760  // Domain
1761  pip_mng->getOpDomainLhsPipeline().push_front(
1762  createOpSchurAssembleBegin());
1763  pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1764 
1765  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1766  false, false
1767 
1768  ));
1769 #endif // ADD_CONTACT
1770  MoFEMFunctionReturn(0);
1771  };
1772 
1773  auto set_assemble_elems = [&]() {
1774  MoFEMFunctionBegin;
1775  auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1776  schur_asmb_pre_proc->preProcessHook = [this]() {
1777  MoFEMFunctionBegin;
1778  CHKERR MatZeroEntries(S);
1779  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1780  MoFEMFunctionReturn(0);
1781  };
1782  auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1783 
1784  schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1785  MoFEMFunctionBegin;
1786  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1787 
1788  // Apply essential constrains to Schur complement
1789  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1790  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1791  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1792  mField, schur_asmb_post_proc, 1, S, aoSchur)();
1793 
1794  MoFEMFunctionReturn(0);
1795  };
1796  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1797  ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1798  ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1799  MoFEMFunctionReturn(0);
1800  };
1801 
1802  auto set_pc = [&]() {
1803  MoFEMFunctionBegin;
1804  CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1805  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1806  MoFEMFunctionReturn(0);
1807  };
1808 
1809  auto set_diagonal_pc = [&]() {
1810  MoFEMFunctionBegin;
1811  KSP *subksp;
1812  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
1813  auto get_pc = [](auto ksp) {
1814  PC pc_raw;
1815  CHKERR KSPGetPC(ksp, &pc_raw);
1816  return SmartPetscObj<PC>(pc_raw,
1817  true); // bump reference
1818  };
1819  CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1820  CHKERR PetscFree(subksp);
1821  MoFEMFunctionReturn(0);
1822  };
1823 
1824  CHKERR set_ops();
1825  CHKERR set_pc();
1826  CHKERR set_assemble_elems();
1827 
1828  CHKERR TSSetUp(solver);
1829  CHKERR KSPSetUp(ksp);
1830  CHKERR set_diagonal_pc();
1831 
1832  } else {
1833  pip_mng->getOpBoundaryLhsPipeline().push_front(
1834  createOpSchurAssembleBegin());
1835  pip_mng->getOpBoundaryLhsPipeline().push_back(
1836  createOpSchurAssembleEnd({}, {}));
1837  pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1838  pip_mng->getOpDomainLhsPipeline().push_back(
1839  createOpSchurAssembleEnd({}, {}));
1840  }
1841 
1842  // fieldSplitIS.reset();
1843  // aoSchur.reset();
1844  MoFEMFunctionReturn(0);
1845 }
1846 
1847 boost::shared_ptr<SetUpSchur>
1848 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
1849  SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1850  SmartPetscObj<AO> ao_up) {
1851  return boost::shared_ptr<SetUpSchur>(
1852  new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
1853 }
1854 
1855 namespace PlasticOps {
1856 
1857 MoFEMErrorCode AddHOOps<2, 3, 3>::add(
1858  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1859  std::vector<FieldSpace> spaces, std::string geom_field_name) {
1860  MoFEMFunctionBegin;
1861  CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
1862  MoFEMFunctionReturn(0);
1863 }
1864 
1865 MoFEMErrorCode AddHOOps<1, 2, 2>::add(
1866  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1867  std::vector<FieldSpace> spaces, std::string geom_field_name) {
1868  MoFEMFunctionBegin;
1869  CHKERR MoFEM::AddHOOps<1, 2, 2>::add(pipeline, spaces, geom_field_name);
1870  MoFEMFunctionReturn(0);
1871 }
1872 
1873 template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM>
1874 MoFEMErrorCode
1875 scaleL2(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1876  std::string geom_field_name) {
1877  MoFEMFunctionBegin;
1878 
1879  auto jac_ptr = boost::make_shared<MatrixDouble>();
1880  auto det_ptr = boost::make_shared<VectorDouble>();
1881  pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
1882  geom_field_name, jac_ptr));
1883  pipeline.push_back(new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, nullptr));
1884 
1885  auto scale_ptr = boost::make_shared<double>(1.);
1886  auto scale = Example::meshVolumeAndCount[0] /
1887  Example::meshVolumeAndCount[1]; // average volume of elements
1888  using OP = ForcesAndSourcesCore::UserDataOperator;
1889  auto op_scale = new OP(NOSPACE, OP::OPSPACE);
1890  op_scale->doWorkRhsHook = [scale_ptr, det_ptr,
1891  scale](DataOperator *base_op_ptr, int, EntityType,
1892  EntitiesFieldData::EntData &) {
1893  *scale_ptr = scale / det_ptr->size(); // distribute average element size
1894  // over integration points
1895  return 0;
1896  };
1897  pipeline.push_back(op_scale);
1898 
1899  for (auto base : {AINSWORTH_LEGENDRE_BASE, AINSWORTH_LOBATTO_BASE,
1900  DEMKOWICZ_JACOBI_BASE}) {
1901  pipeline.push_back(
1902  new OpScaleBaseBySpaceInverseOfMeasure(L2, base, det_ptr, scale_ptr));
1903  }
1904 
1905  MoFEMFunctionReturn(0);
1906 }
1907 
1908 MoFEMErrorCode AddHOOps<3, 3, 3>::add(
1909  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1910  std::vector<FieldSpace> spaces, std::string geom_field_name) {
1911  MoFEMFunctionBegin;
1912  constexpr bool scale_l2 = false;
1913  if (scale_l2) {
1914  CHKERR scaleL2<3, 3, 3>(pipeline, geom_field_name);
1915  }
1916  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name,
1917  nullptr, nullptr, nullptr);
1918  MoFEMFunctionReturn(0);
1919 }
1920 
1921 MoFEMErrorCode AddHOOps<2, 2, 2>::add(
1922  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1923  std::vector<FieldSpace> spaces, std::string geom_field_name) {
1924  MoFEMFunctionBegin;
1925  constexpr bool scale_l2 = false;
1926  if (scale_l2) {
1927  CHKERR scaleL2<2, 2, 2>(pipeline, geom_field_name);
1928  }
1929  CHKERR MoFEM::AddHOOps<2, 2, 2>::add(pipeline, spaces, geom_field_name,
1930  nullptr, nullptr, nullptr);
1931  MoFEMFunctionReturn(0);
1932 }
1933 
1934 } // namespace PlasticOps
NOSPACE
@ NOSPACE
Definition: definitions.h:83
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:818
SetUpSchurImpl
Definition: test_broken_space.cpp:519
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:118
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
PlasticOps::AddHOOps
Definition: plastic.cpp:184
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
PlasticOps::scaleL2
MoFEMErrorCode scaleL2(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string geom_field_name)
Definition: plastic.cpp:1875
Example::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:634
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
SetUpSchurImpl::aoSchur
SmartPetscObj< AO > aoSchur
Definition: plastic.cpp:1662
Example::Example
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:218
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPSPACE
@ OPSPACE
operator do Work is execute on space data
Definition: ForcesAndSourcesCore.hpp:570
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:125
ContactOps
Definition: contact.cpp:99
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:532
rho
double rho
Definition: plastic.cpp:144
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:127
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
do_eval_field
PetscBool do_eval_field
Evaluate field.
Definition: plastic.cpp:119
Example::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:239
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: plastic.cpp:1650
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
filter_true_skin
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Definition: EshelbianPlasticity.cpp:142
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:519
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:78
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:243
HenckyOps
Definition: HenckyOps.hpp:11
MoFEM::OpScaleBaseBySpaceInverseOfMeasure
Scale base functions by inverses of measure of element.
Definition: HODataOperators.hpp:390
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
scale
double scale
Definition: plastic.cpp:123
HenckyOps.hpp
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
int order
Order displacement.
Definition: plastic.cpp:138
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:130
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:141
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1660
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::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1276
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:597
iso_hardening_exp
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:64
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:62
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
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:2585
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: test_broken_space.cpp:529
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:1560
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::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:1661
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:92
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:237
is_quasi_static
PetscBool is_quasi_static
Definition: plastic.cpp:143
MoFEM::FieldEvaluatorInterface::SetPtsData
Definition: FieldEvaluator.hpp:29
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:216
help
static char help[]
[TestOperators]
Definition: plastic.cpp:1558
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:528
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
cn1
double cn1
Definition: plastic.cpp:136
Example::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:238
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:201
C1_k
double C1_k
Kinematic hardening.
Definition: plastic.cpp:133
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:106
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:126
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2580
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
cn0
double cn0
Definition: plastic.cpp:135
MatrixFunction.hpp
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
get_skin
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
Definition: EshelbianPlasticity.cpp:163
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:169
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
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:129
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
VolOp
VolEle::UserDataOperator VolOp
Definition: test_cache_on_entities.cpp:23
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:109
SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: plastic.cpp:62
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:137
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:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
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:2627
b_iso
double b_iso
Saturation exponent.
Definition: plastic.cpp:132
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:451
DomainEleOp
PlasticOps
Definition: plastic.cpp:182
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
Example::testOperators
MoFEMErrorCode testOperators()
[Solve]
Definition: plastic.cpp:1457
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:590
IntegrationRules.hpp
iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:73
Example::VOL
@ VOL
Definition: plastic.cpp:222
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:117
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:275
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:2343
Example::ScaledTimeScale
Definition: plastic.cpp:241
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:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ep_order
int ep_order
Order of ep files.
Definition: plastic.cpp:140
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3684
Qinf
double Qinf
Saturation yield stress.
Definition: plastic.cpp:131
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:239
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:256
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:226
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:44
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
Example::meshVolumeAndCount
static std::array< double, 2 > meshVolumeAndCount
Definition: plastic.cpp:223
atom_test
int atom_test
Atom test.
Definition: plastic.cpp:121
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:235
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:1082
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:480
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:589
MoFEM::SmartPetscObj< VecScatter >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:802
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
OP
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:128
tau_order
int tau_order
Order of tau files.
Definition: plastic.cpp:139
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
Definition: plastic.cpp:1639
alpha_damping
double alpha_damping
Definition: plastic.cpp:145
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow