v0.14.0
contact.cpp
Go to the documentation of this file.
1 /**
2  * \file contact.cpp
3  * \CONTACT contact.cpp
4  *
5  * CONTACT of contact problem
6  *
7  * @copyright Copyright (c) 2023
8  */
9 
10 /* The above code is a preprocessor directive in C++ that checks if the macro
11 "EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it is set
12 to 3" */
13 #ifndef EXECUTABLE_DIMENSION
14 #define EXECUTABLE_DIMENSION 3
15 #endif
16 
17 #ifndef SCHUR_ASSEMBLE
18 #define SCHUR_ASSEMBLE 0
19 #endif
20 
21 #include <MoFEM.hpp>
22 #include <MatrixFunction.hpp>
23 
24 #ifdef PYTHON_SDF
25 #include <boost/python.hpp>
26 #include <boost/python/def.hpp>
27 #include <boost/python/numpy.hpp>
28 namespace bp = boost::python;
29 namespace np = boost::python::numpy;
30 #endif
31 
32 using namespace MoFEM;
33 
34 constexpr AssemblyType AT =
35  (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
36  : AssemblyType::PETSC; //< selected assembly type
37 constexpr IntegrationType IT =
38  IntegrationType::GAUSS; //< selected integration type
39 
40 template <int DIM> struct ElementsAndOps;
41 
42 template <> struct ElementsAndOps<2> : PipelineManager::ElementsAndOpsByDim<2> {
43  static constexpr FieldSpace CONTACT_SPACE = HCURL;
44 };
45 
46 template <> struct ElementsAndOps<3> : PipelineManager::ElementsAndOpsByDim<3> {
47  static constexpr FieldSpace CONTACT_SPACE = HDIV;
48 };
49 
50 constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
51 constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
52 
53 constexpr int SPACE_DIM =
54  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
55 
56 /* The above code is defining an alias `EntData` for the type
57 `EntitiesFieldData::EntData`. This is a C++ syntax for creating a new name for
58 an existing type. */
59 using EntData = EntitiesFieldData::EntData;
60 using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
61 using DomainEleOp = DomainEle::UserDataOperator;
62 using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
63 using BoundaryEleOp = BoundaryEle::UserDataOperator;
64 
65 //! [Specialisation for assembly]
66 
67 constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
68 
69 //! [Operators used for contact]
70 using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::BiLinearForm<
71  IT>::OpMass<1, SPACE_DIM>;
72 using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::LinearForm<
73  IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
74 //! [Operators used for contact]
75 
76 PetscBool is_quasi_static = PETSC_TRUE;
77 
78 int order = 2;
79 int contact_order = 2;
80 int geom_order = 1;
81 double young_modulus = 100;
82 double poisson_ratio = 0.25;
83 double rho = 0.0;
84 double spring_stiffness = 0.0;
85 double vis_spring_stiffness = 0.0;
86 double alpha_damping = 0;
87 
88 double scale = 1.;
89 
90 PetscBool is_axisymmetric = PETSC_FALSE; //< Axisymmetric model
91 
92 // ##define HENCKY_SMALL_STRAIN
93 
94 int atom_test = 0;
95 
96 namespace ContactOps {
97 double cn_contact = 0.1;
98 }; // namespace ContactOps
99 
100 #include <HenckyOps.hpp>
101 using namespace HenckyOps;
102 #include <ContactOps.hpp>
103 #include <PostProcContact.hpp>
104 #include <ContactNaturalBC.hpp>
105 
106 #ifdef WITH_MODULE_MFRONT_INTERFACE
107 #include <MFrontMoFEMInterface.hpp>
108 #endif
109 
110 using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
111 using OpDomainRhsBCs =
112  DomainRhsBCs::OpFlux<ContactOps::DomainBCs, 1, SPACE_DIM>;
113 using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
114 using OpBoundaryRhsBCs =
115  BoundaryRhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
116 using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
117 using OpBoundaryLhsBCs =
118  BoundaryLhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
119 
120 using namespace ContactOps;
121 
122 struct Contact {
123 
124  Contact(MoFEM::Interface &m_field) : mField(m_field) {}
125 
126  MoFEMErrorCode runProblem();
127 
128 private:
129  MoFEM::Interface &mField;
130 
131  MoFEMErrorCode setupProblem();
132  MoFEMErrorCode createCommonData();
133  MoFEMErrorCode bC();
134  MoFEMErrorCode OPs();
135  MoFEMErrorCode tsSolve();
136  MoFEMErrorCode checkResults();
137 
138  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
139  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
140  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
141 
142  boost::shared_ptr<GenericElementInterface> mfrontInterface;
143  boost::shared_ptr<Monitor> monitorPtr;
144 
145 #ifdef PYTHON_SDF
146  boost::shared_ptr<SDFPython> sdfPythonPtr;
147 #endif
148 
149  struct ScaledTimeScale : public MoFEM::TimeScale {
150  using MoFEM::TimeScale::TimeScale;
151  double getScale(const double time) {
152  return scale * MoFEM::TimeScale::getScale(time);
153  };
154  };
155 };
156 
157 //! [Run problem]
158 MoFEMErrorCode Contact::runProblem() {
159  MoFEMFunctionBegin;
160  CHKERR setupProblem();
161  CHKERR createCommonData();
162  CHKERR bC();
163  CHKERR OPs();
164  CHKERR tsSolve();
165  CHKERR checkResults();
166  MoFEMFunctionReturn(0);
167 }
168 //! [Run problem]
169 
170 //! [Set up problem]
171 MoFEMErrorCode Contact::setupProblem() {
172  MoFEMFunctionBegin;
173  Simple *simple = mField.getInterface<Simple>();
174 
175  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
176  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-contact_order", &contact_order,
177  PETSC_NULL);
178  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
179  PETSC_NULL);
180 
181  MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
182  if (contact_order != order)
183  MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
184  MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
185 
186  // Select base
187  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
188  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
189  PetscInt choice_base_value = AINSWORTH;
190  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
191  LASBASETOPT, &choice_base_value, PETSC_NULL);
192 
193  FieldApproximationBase base;
194  switch (choice_base_value) {
195  case AINSWORTH:
196  base = AINSWORTH_LEGENDRE_BASE;
197  MOFEM_LOG("CONTACT", Sev::inform)
198  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
199  break;
200  case DEMKOWICZ:
201  base = DEMKOWICZ_JACOBI_BASE;
202  MOFEM_LOG("CONTACT", Sev::inform)
203  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
204  break;
205  default:
206  base = LASTBASE;
207  break;
208  }
209 
210  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
211  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
212  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
213  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
214  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
215  SPACE_DIM);
216  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
217  SPACE_DIM);
218  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
219 
220  CHKERR simple->setFieldOrder("U", order);
221  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
222 
223  auto get_skin = [&]() {
224  Range body_ents;
225  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
226  Skinner skin(&mField.get_moab());
227  Range skin_ents;
228  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
229  return skin_ents;
230  };
231 
232  auto filter_blocks = [&](auto skin) {
233  bool is_contact_block = false;
234  Range contact_range;
235  for (auto m :
236  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
237 
238  (boost::format("%s(.*)") % "CONTACT").str()
239 
240  ))
241 
242  ) {
243  is_contact_block =
244  true; ///< blocs interation is collective, so that is set irrespective
245  ///< if there are entities in given rank or not in the block
246  MOFEM_LOG("CONTACT", Sev::inform)
247  << "Find contact block set: " << m->getName();
248  auto meshset = m->getMeshset();
249  Range contact_meshset_range;
250  CHKERR mField.get_moab().get_entities_by_dimension(
251  meshset, SPACE_DIM - 1, contact_meshset_range, true);
252 
253  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
254  contact_meshset_range);
255  contact_range.merge(contact_meshset_range);
256  }
257  if (is_contact_block) {
258  MOFEM_LOG("SYNC", Sev::inform)
259  << "Nb entities in contact surface: " << contact_range.size();
260  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
261  skin = intersect(skin, contact_range);
262  }
263  return skin;
264  };
265 
266  auto filter_true_skin = [&](auto skin) {
267  Range boundary_ents;
268  ParallelComm *pcomm =
269  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
270  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
271  PSTATUS_NOT, -1, &boundary_ents);
272  return boundary_ents;
273  };
274 
275  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
276  CHKERR simple->setFieldOrder("SIGMA", 0);
277  int sigma_order = std::max(order, contact_order) - 1;
278  CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
279 
280  if (contact_order > order) {
281  Range ho_ents;
282  if constexpr (SPACE_DIM == 3) {
283  CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
284  moab::Interface::UNION);
285  } else {
286  ho_ents = boundary_ents;
287  }
288  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
289  CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
290  CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
291  }
292 
293  CHKERR simple->setUp();
294 
295  auto project_ho_geometry = [&]() {
296  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
297  return mField.loop_dofs("GEOMETRY", ent_method);
298  };
299 
300  PetscBool project_geometry = PETSC_TRUE;
301  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
302  &project_geometry, PETSC_NULL);
303  if (project_geometry) {
304  CHKERR project_ho_geometry();
305  }
306 
307  MoFEMFunctionReturn(0);
308 } //! [Set up problem]
309 
310 //! [Create common data]
311 MoFEMErrorCode Contact::createCommonData() {
312  MoFEMFunctionBegin;
313 
314  PetscBool use_mfront = PETSC_FALSE;
315  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_mfront", &use_mfront,
316  PETSC_NULL);
317  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_axisymmetric",
318  &is_axisymmetric, PETSC_NULL);
319  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test,
320  PETSC_NULL);
321 
322  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
323  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus", &young_modulus,
324  PETSC_NULL);
325  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio", &poisson_ratio,
326  PETSC_NULL);
327  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
328  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn_contact, PETSC_NULL);
329  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
330  &spring_stiffness, PETSC_NULL);
331  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-vis_spring_stiffness",
332  &vis_spring_stiffness, PETSC_NULL);
333  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping", &alpha_damping,
334  PETSC_NULL);
335 
336  if (!mfrontInterface) {
337  MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
338  MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
339  } else {
340  MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
341  }
342 
343  MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
344  MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
345  MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
346  MOFEM_LOG("CONTACT", Sev::inform)
347  << "Vis spring_stiffness " << vis_spring_stiffness;
348 
349  MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
350 
351  PetscBool use_scale = PETSC_FALSE;
352  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_scale", &use_scale,
353  PETSC_NULL);
354  if (use_scale) {
355  scale /= young_modulus;
356  }
357 
358  MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
359 
360  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_quasi_static",
361  &is_quasi_static, PETSC_NULL);
362  MOFEM_LOG("CONTACT", Sev::inform)
363  << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
364 
365 #ifdef PYTHON_SDF
366  char sdf_file_name[255];
367  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
368  sdf_file_name, 255, PETSC_NULL);
369 
370  sdfPythonPtr = boost::make_shared<SDFPython>();
371  CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
372  sdfPythonWeakPtr = sdfPythonPtr;
373 #endif
374 
375  if (is_axisymmetric) {
376  if (SPACE_DIM == 3) {
377  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
378  "Use executable contact_2d with axisymmetric model");
379  } else {
380  if (!use_mfront) {
381  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
382  "Axisymmetric model is only available with MFront (set "
383  "use_mfront to 1)");
384  } else {
385  MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
386  }
387  }
388  } else {
389  if (SPACE_DIM == 2) {
390  MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
391  }
392  }
393 
394  if (use_mfront) {
395 #ifndef WITH_MODULE_MFRONT_INTERFACE
396  SETERRQ(
397  PETSC_COMM_SELF, MOFEM_NOT_FOUND,
398  "MFrontInterface module was not found while use_mfront was set to 1");
399 #else
400  if (SPACE_DIM == 3) {
401  mfrontInterface =
402  boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
403  mField, "U", "GEOMETRY", true, is_quasi_static);
404  } else if (SPACE_DIM == 2) {
405  if (is_axisymmetric) {
406  mfrontInterface =
407  boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
408  mField, "U", "GEOMETRY", true, is_quasi_static);
409  } else {
410  mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
411  mField, "U", "GEOMETRY", true, is_quasi_static);
412  }
413  }
414 #endif
415  CHKERR mfrontInterface->getCommandLineParameters();
416  }
417 
418  Simple *simple = mField.getInterface<Simple>();
419  auto dm = simple->getDM();
420  monitorPtr =
421  boost::make_shared<Monitor>(dm, scale, mfrontInterface, is_axisymmetric);
422 
423  if (use_mfront) {
424  mfrontInterface->setMonitorPtr(monitorPtr);
425  }
426 
427  MoFEMFunctionReturn(0);
428 }
429 //! [Create common data]
430 
431 //! [Boundary condition]
432 MoFEMErrorCode Contact::bC() {
433  MoFEMFunctionBegin;
434  auto bc_mng = mField.getInterface<BcManager>();
435  auto simple = mField.getInterface<Simple>();
436 
437  for (auto f : {"U", "SIGMA"}) {
438  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
439  "REMOVE_X", f, 0, 0);
440  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
441  "REMOVE_Y", f, 1, 1);
442  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
443  "REMOVE_Z", f, 2, 2);
444  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
445  "REMOVE_ALL", f, 0, 3);
446  }
447 
448  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
449  "SIGMA", 0, 0, false, true);
450  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
451  "SIGMA", 1, 1, false, true);
452  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
453  "SIGMA", 2, 2, false, true);
454  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
455  "SIGMA", 0, 3, false, true);
456  CHKERR bc_mng->removeBlockDOFsOnEntities(
457  simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
458 
459  // Note remove has to be always before push. Then node marking will be
460  // corrupted.
461  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
462  simple->getProblemName(), "U");
463 
464  MoFEMFunctionReturn(0);
465 }
466 //! [Boundary condition]
467 
468 //! [Push operators to pip]
469 MoFEMErrorCode Contact::OPs() {
470  MoFEMFunctionBegin;
471  auto simple = mField.getInterface<Simple>();
472  auto *pip_mng = mField.getInterface<PipelineManager>();
473  auto bc_mng = mField.getInterface<BcManager>();
474  auto time_scale = boost::make_shared<ScaledTimeScale>();
475  auto body_force_time_scale =
476  boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
477 
478  auto integration_rule_vol = [](int, int, int approx_order) {
479  return 2 * approx_order + geom_order - 1;
480  };
481  auto integration_rule_boundary = [](int, int, int approx_order) {
482  return 2 * approx_order + geom_order - 1;
483  };
484 
485  auto add_domain_base_ops = [&](auto &pip) {
486  MoFEMFunctionBegin;
487  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
488  "GEOMETRY");
489  MoFEMFunctionReturn(0);
490  };
491 
492  auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
493  henky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
494  henky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
495 
496  auto add_domain_ops_lhs = [&](auto &pip) {
497  MoFEMFunctionBegin;
498 
499  //! [Only used for dynamics]
500  using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
501  GAUSS>::OpMass<1, SPACE_DIM>;
502  //! [Only used for dynamics]
503  if (is_quasi_static == PETSC_FALSE) {
504 
505  auto *pip_mng = mField.getInterface<PipelineManager>();
506  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
507 
508  auto get_inertia_and_mass_damping =
509  [this, fe_domain_lhs](const double, const double, const double) {
510  return (rho * scale) * fe_domain_lhs->ts_aa +
511  (alpha_damping * scale) * fe_domain_lhs->ts_a;
512  };
513  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
514  } else {
515 
516  auto *pip_mng = mField.getInterface<PipelineManager>();
517  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
518 
519  auto get_inertia_and_mass_damping =
520  [this, fe_domain_lhs](const double, const double, const double) {
521  return (alpha_damping * scale) * fe_domain_lhs->ts_a;
522  };
523  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
524  }
525 
526  if (!mfrontInterface) {
527  CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
528  mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
529  } else {
530  CHKERR mfrontInterface->opFactoryDomainLhs(pip);
531  }
532 
533  MoFEMFunctionReturn(0);
534  };
535 
536  auto add_domain_ops_rhs = [&](auto &pip) {
537  MoFEMFunctionBegin;
538 
539  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
540  pip, mField, "U", {body_force_time_scale}, Sev::inform);
541 
542  //! [Only used for dynamics]
543  using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
544  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
545  //! [Only used for dynamics]
546 
547  // only in case of dynamics
548  if (is_quasi_static == PETSC_FALSE) {
549  auto mat_acceleration = boost::make_shared<MatrixDouble>();
550  pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
551  "U", mat_acceleration));
552  pip.push_back(
553  new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
554  return rho * scale;
555  }));
556  }
557 
558  // only in case of viscosity
559  if (alpha_damping > 0) {
560  auto mat_velocity = boost::make_shared<MatrixDouble>();
561  pip.push_back(
562  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
563  pip.push_back(
564  new OpInertiaForce("U", mat_velocity, [](double, double, double) {
565  return alpha_damping * scale;
566  }));
567  }
568 
569  if (!mfrontInterface) {
570  CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
571  mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
572  } else {
573  CHKERR mfrontInterface->opFactoryDomainRhs(pip);
574  }
575 
576  CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
577  pip, "SIGMA", "U", is_axisymmetric);
578 
579  MoFEMFunctionReturn(0);
580  };
581 
582  auto add_boundary_base_ops = [&](auto &pip) {
583  MoFEMFunctionBegin;
584  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
585  "GEOMETRY");
586  // We have to integrate on curved face geometry, thus integration weight
587  // have to adjusted.
588  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
589  MoFEMFunctionReturn(0);
590  };
591 
592  auto add_boundary_ops_lhs = [&](auto &pip) {
593  MoFEMFunctionBegin;
594 
595  //! [Operators used for contact]
596  using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
597  AT>::BiLinearForm<IT>::OpMass<1, SPACE_DIM>;
598  //! [Operators used for contact]
599 
600  // Add Natural BCs to LHS
601  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
602  pip, mField, "U", Sev::inform);
603 
604  if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
605 
606  auto *pip_mng = mField.getInterface<PipelineManager>();
607  auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
608 
609  pip.push_back(new OpSpringLhs(
610  "U", "U",
611 
612  [this, fe_boundary_lhs](double, double, double) {
613  return spring_stiffness * scale +
614  (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
615  }
616 
617  ));
618  }
619 
620  CHKERR
621  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
622  pip, "SIGMA", "U", is_axisymmetric);
623  CHKERR ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, GAUSS,
624  DomainEle>(
625  mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
626  integration_rule_vol, is_axisymmetric);
627 
628  MoFEMFunctionReturn(0);
629  };
630 
631  auto add_boundary_ops_rhs = [&](auto &pip) {
632  MoFEMFunctionBegin;
633 
634  //! [Operators used for contact]
635  using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
636  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
637  //! [Operators used for contact]
638 
639  // Add Natural BCs to RHS
640  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
641  pip, mField, "U", {time_scale}, Sev::inform);
642 
643  if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
644  auto u_disp = boost::make_shared<MatrixDouble>();
645  auto dot_u_disp = boost::make_shared<MatrixDouble>();
646  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
647  pip.push_back(
648  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
649  pip.push_back(
650  new OpSpringRhs("U", u_disp, [this](double, double, double) {
651  return spring_stiffness * scale;
652  }));
653  pip.push_back(
654  new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
655  return vis_spring_stiffness * scale;
656  }));
657  }
658 
659  CHKERR
660  ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
661  pip, "SIGMA", "U", is_axisymmetric);
662 
663  MoFEMFunctionReturn(0);
664  };
665 
666  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
667  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
668  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
669  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
670 
671  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
672  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
673  CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
674  CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
675 
676  if (mfrontInterface) {
677  CHKERR mfrontInterface->setUpdateElementVariablesOperators();
678  }
679 
680  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
681  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
682  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
683  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
684 
685  MoFEMFunctionReturn(0);
686 }
687 //! [Push operators to pip]
688 
689 //! [Solve]
690 struct SetUpSchur {
691  static boost::shared_ptr<SetUpSchur>
692  createSetUpSchur(MoFEM::Interface &m_field);
693  virtual MoFEMErrorCode setUp(SmartPetscObj<TS> solver) = 0;
694 
695 protected:
696  SetUpSchur() = default;
697 };
698 
699 MoFEMErrorCode Contact::tsSolve() {
700  MoFEMFunctionBegin;
701 
702  Simple *simple = mField.getInterface<Simple>();
703  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
704  ISManager *is_manager = mField.getInterface<ISManager>();
705 
706  auto set_section_monitor = [&](auto solver) {
707  MoFEMFunctionBegin;
708  SNES snes;
709  CHKERR TSGetSNES(solver, &snes);
710  PetscViewerAndFormat *vf;
711  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
712  PETSC_VIEWER_DEFAULT, &vf);
713  CHKERR SNESMonitorSet(
714  snes,
715  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
716  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
717  MoFEMFunctionReturn(0);
718  };
719 
720  auto scatter_create = [&](auto D, auto coeff) {
721  SmartPetscObj<IS> is;
722  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
723  ROW, "U", coeff, coeff, is);
724  int loc_size;
725  CHKERR ISGetLocalSize(is, &loc_size);
726  Vec v;
727  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
728  VecScatter scatter;
729  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
730  return std::make_tuple(SmartPetscObj<Vec>(v),
731  SmartPetscObj<VecScatter>(scatter));
732  };
733 
734  auto set_time_monitor = [&](auto dm, auto solver) {
735  MoFEMFunctionBegin;
736  monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
737  boost::shared_ptr<ForcesAndSourcesCore> null;
738  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
739  monitorPtr, null, null);
740  MoFEMFunctionReturn(0);
741  };
742 
743  auto set_essential_bc = [&]() {
744  MoFEMFunctionBegin;
745  // This is low level pushing finite elements (pipelines) to solver
746  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
747  auto pre_proc_ptr = boost::make_shared<FEMethod>();
748  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
749  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
750 
751  // Add boundary condition scaling
752  auto time_scale = boost::make_shared<TimeScale>();
753 
754  auto get_bc_hook_rhs = [&]() {
755  EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
756  {time_scale}, false);
757  return hook;
758  };
759  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
760 
761  auto get_post_proc_hook_rhs = [&]() {
762  return EssentialPostProcRhs<DisplacementCubitBcData>(
763  mField, post_proc_rhs_ptr, 1.);
764  };
765  auto get_post_proc_hook_lhs = [&]() {
766  return EssentialPostProcLhs<DisplacementCubitBcData>(
767  mField, post_proc_lhs_ptr, 1.);
768  };
769  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
770 
771  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
772  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
773  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
774  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
775  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
776  MoFEMFunctionReturn(0);
777  };
778 
779  // Set up Schur preconditioner
780  auto set_schur_pc = [&](auto solver) {
781  boost::shared_ptr<SetUpSchur> schur_ptr;
782  if (AT == AssemblyType::BLOCK_SCHUR) {
783  // Set up Schur preconditioner
784  schur_ptr = SetUpSchur::createSetUpSchur(mField);
785  CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
786  }
787  return schur_ptr;
788  };
789 
790  auto dm = simple->getDM();
791  auto D = createDMVector(dm);
792 
793  ContactOps::CommonData::createTotalTraction(mField);
794 
795  uXScatter = scatter_create(D, 0);
796  uYScatter = scatter_create(D, 1);
797  if (SPACE_DIM == 3)
798  uZScatter = scatter_create(D, 2);
799 
800  // Add extra finite elements to SNES solver pipelines to resolve essential
801  // boundary conditions
802  CHKERR set_essential_bc();
803 
804  if (is_quasi_static == PETSC_TRUE) {
805  auto solver = pip_mng->createTSIM();
806  CHKERR TSSetFromOptions(solver);
807  auto schur_pc_ptr = set_schur_pc(solver);
808 
809  auto D = createDMVector(dm);
810  CHKERR set_section_monitor(solver);
811  CHKERR set_time_monitor(dm, solver);
812  CHKERR TSSetSolution(solver, D);
813  CHKERR TSSetUp(solver);
814  CHKERR TSSolve(solver, NULL);
815  } else {
816  auto solver = pip_mng->createTSIM2();
817  CHKERR TSSetFromOptions(solver);
818  auto schur_pc_ptr = set_schur_pc(solver);
819 
820  auto dm = simple->getDM();
821  auto D = createDMVector(dm);
822  auto DD = vectorDuplicate(D);
823 
824  CHKERR set_section_monitor(solver);
825  CHKERR set_time_monitor(dm, solver);
826  CHKERR TS2SetSolution(solver, D, DD);
827  CHKERR TSSetUp(solver);
828  CHKERR TSSolve(solver, NULL);
829  }
830 
831  MoFEMFunctionReturn(0);
832 }
833 //! [Solve]
834 
835 //! [Check]
836 MoFEMErrorCode Contact::checkResults() {
837  MoFEMFunctionBegin;
838  if (atom_test && !mField.get_comm_rank()) {
839  const double *t_ptr;
840  CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
841  double hertz_force;
842  double fem_force;
843  double norm = 1e-5;
844  double tol = 1e-3;
845  double tol_norm = 7.5; // change when analytical functions are updated
846  switch (atom_test) {
847  case 1: // plane stress
848  hertz_force = 3.927;
849  fem_force = t_ptr[1];
850  break;
851  case 2: // plane strain
852  hertz_force = 4.675;
853  fem_force = t_ptr[1];
854  norm = monitorPtr->getErrorNorm(1);
855  break;
856  case 3: // 3D
857  hertz_force = 3.968;
858  fem_force = t_ptr[2];
859  case 4: // axisymmetric
860  tol = 5e3;
861  case 5: // axisymmetric
862  hertz_force = 15.873;
863  fem_force = t_ptr[1];
864  norm = monitorPtr->getErrorNorm(1);
865  break;
866  case 6: // wavy 2d
867  hertz_force = 0.374;
868  fem_force = t_ptr[1];
869  break;
870  case 7: // wavy 3d
871  hertz_force = 0.5289;
872  fem_force = t_ptr[2];
873  break;
874  default:
875  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
876  "atom test %d does not exist", atom_test);
877  }
878  if (fabs(fem_force - hertz_force) / hertz_force > tol) {
879  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
880  "atom test %d diverged! %3.4e != %3.4e", atom_test, fem_force,
881  hertz_force);
882  }
883  if (norm > tol_norm) {
884  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
885  "atom test %d diverged! %3.4e > %3.4e", atom_test, norm,
886  tol_norm);
887  }
888  CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
889  }
890 
891  ContactOps::CommonData::totalTraction.reset();
892 
893  MoFEMFunctionReturn(0);
894 }
895 //! [Check]
896 
897 static char help[] = "...\n\n";
898 
899 int main(int argc, char *argv[]) {
900 
901 #ifdef PYTHON_SDF
902  Py_Initialize();
903  np::initialize();
904 #endif
905 
906  // Initialisation of MoFEM/PETSc and MOAB data structures
907  const char param_file[] = "param_file.petsc";
908  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
909 
910  // Add logging channel for CONTACT
911  auto core_log = logging::core::get();
912  core_log->add_sink(
913  LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
914  LogManager::setLog("CONTACT");
915  MOFEM_LOG_TAG("CONTACT", "Indent");
916 
917  try {
918 
919  //! [Register MoFEM discrete manager in PETSc]
920  DMType dm_name = "DMMOFEM";
921  CHKERR DMRegister_MoFEM(dm_name);
922  //! [Register MoFEM discrete manager in PETSc
923 
924  //! [Create MoAB]
925  moab::Core mb_instance; ///< mesh database
926  moab::Interface &moab = mb_instance; ///< mesh database interface
927  //! [Create MoAB]
928 
929  //! [Create MoFEM]
930  MoFEM::Core core(moab); ///< finite element database
931  MoFEM::Interface &m_field = core; ///< finite element database interface
932  //! [Create MoFEM]
933 
934  //! [Load mesh]
935  Simple *simple = m_field.getInterface<Simple>();
936  CHKERR simple->getOptions();
937  CHKERR simple->loadFile("");
938  //! [Load mesh]
939 
940  //! [CONTACT]
941  Contact ex(m_field);
942  CHKERR ex.runProblem();
943  //! [CONTACT]
944  }
945  CATCH_ERRORS;
946 
947  CHKERR MoFEM::Core::Finalize();
948 
949 #ifdef PYTHON_SDF
950  if (Py_FinalizeEx() < 0) {
951  exit(120);
952  }
953 #endif
954 
955  return 0;
956 }
957 
958 struct SetUpSchurImpl : public SetUpSchur {
959 
960  SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {}
961 
962  virtual ~SetUpSchurImpl() {}
963 
964  MoFEMErrorCode setUp(SmartPetscObj<TS> solver);
965 
966 private:
967  MoFEMErrorCode createSubDM();
968  MoFEMErrorCode setOperator();
969  MoFEMErrorCode setPC(PC pc);
970  MoFEMErrorCode setDiagonalPC(PC pc);
971 
972  MoFEM::Interface &mField;
973 
974  SmartPetscObj<Mat> S;
975  SmartPetscObj<DM> schurDM;
976  SmartPetscObj<DM> blockDM;
977 };
978 
979 MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<TS> solver) {
980  MoFEMFunctionBegin;
981  auto simple = mField.getInterface<Simple>();
982  auto pip = mField.getInterface<PipelineManager>();
983 
984  SNES snes;
985  CHKERR TSGetSNES(solver, &snes);
986  KSP ksp;
987  CHKERR SNESGetKSP(snes, &ksp);
988  CHKERR KSPSetFromOptions(ksp);
989 
990  PC pc;
991  CHKERR KSPGetPC(ksp, &pc);
992 
993  PetscBool is_pcfs = PETSC_FALSE;
994  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
995  if (is_pcfs) {
996 
997  MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
998 
999  if (S) {
1000  CHK_THROW_MESSAGE(
1001  MOFEM_DATA_INCONSISTENCY,
1002  "It is expected that Schur matrix is not allocated. This is "
1003  "possible only if PC is set up twice");
1004  }
1005 
1006  CHKERR createSubDM();
1007 
1008  // Add data to DM storage
1009  S = createDMMatrix(schurDM);
1010  CHKERR MatSetBlockSize(S, SPACE_DIM);
1011  // CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
1012 
1013  // Set DM to use shell block matrix
1014  DM solver_dm;
1015  CHKERR TSGetDM(solver, &solver_dm);
1016  CHKERR DMSetMatType(solver_dm, MATSHELL);
1017 
1018  auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1019  auto A = createDMBlockMat(simple->getDM());
1020  auto P = createDMNestSchurMat(simple->getDM());
1021 
1022  if (is_quasi_static == PETSC_TRUE) {
1023  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1024  Mat A, Mat B, void *ctx) {
1025  return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1026  };
1027  CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1028  } else {
1029  auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1030  PetscReal a, PetscReal aa, Mat A, Mat B,
1031  void *ctx) {
1032  return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1033  };
1034  CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1035  }
1036  CHKERR KSPSetOperators(ksp, A, P);
1037 
1038  CHKERR setOperator();
1039  CHKERR setPC(pc);
1040  CHKERR TSSetUp(solver);
1041  CHKERR KSPSetUp(ksp);
1042  CHKERR setDiagonalPC(pc);
1043 
1044  } else {
1045  MOFEM_LOG("CONTACT", Sev::inform) << "No Schur PC";
1046  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1047  pip->getOpBoundaryLhsPipeline().push_back(
1048  createOpSchurAssembleEnd({}, {}, {}, {}, {}, false));
1049  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1050  pip->getOpDomainLhsPipeline().push_back(
1051  createOpSchurAssembleEnd({}, {}, {}, {}, {}, false));
1052  }
1053  MoFEMFunctionReturn(0);
1054 }
1055 
1056 MoFEMErrorCode SetUpSchurImpl::createSubDM() {
1057  MoFEMFunctionBegin;
1058  auto simple = mField.getInterface<Simple>();
1059 
1060  auto create_dm = [&](const char *name, const char *field_name) {
1061  auto dm = createDM(mField.get_comm(), "DMMOFEM");
1062  auto create_dm_imp = [&]() {
1063  MoFEMFunctionBegin;
1064  CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
1065  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
1066  CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
1067  CHKERR DMMoFEMAddSubFieldRow(dm, field_name);
1068  CHKERR DMMoFEMAddSubFieldCol(dm, field_name);
1069  CHKERR DMSetUp(dm);
1070  MoFEMFunctionReturn(0);
1071  };
1072  CHK_THROW_MESSAGE(
1073  create_dm_imp(),
1074  "Error in creating schurDM. It is possible that schurDM is "
1075  "already created");
1076  return dm;
1077  };
1078 
1079  // Note: here we can make block with bubbles of "U" and "SIGMA" fields. See
1080  // vec-0 where bubbles are added.
1081 
1082  schurDM = create_dm("SCHUR", "U");
1083  blockDM = create_dm("BLOCK", "SIGMA");
1084 
1085  if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1086 
1087  auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1088  auto block_mat_data = createBlockMatStructure(
1089  simple->getDM(),
1090 
1091  {{
1092 
1093  simple->getDomainFEName(),
1094 
1095  {
1096 
1097  {"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
1098 
1099  }}}
1100 
1101  );
1102 
1103  return getNestSchurData(
1104 
1105  {schurDM, blockDM}, block_mat_data,
1106 
1107  {"SIGMA"}, {nullptr}, true
1108 
1109  );
1110  };
1111 
1112  auto nested_mat_data = get_nested_mat_data(schurDM, blockDM);
1113  CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1114 
1115  } else {
1116  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1117  "Only BLOCK_SCHUR is implemented");
1118  }
1119 
1120  MoFEMFunctionReturn(0);
1121 }
1122 
1123 MoFEMErrorCode SetUpSchurImpl::setOperator() {
1124  MoFEMFunctionBegin;
1125 
1126  double eps_stab = 1e-4;
1127  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1128  PETSC_NULL);
1129 
1130  using B = FormsIntegrators<BoundaryEleOp>::Assembly<
1131  BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
1132  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1133 
1134  auto simple = mField.getInterface<Simple>();
1135  auto pip = mField.getInterface<PipelineManager>();
1136 
1137  // block data structure
1138  boost::shared_ptr<BlockStructure> block_data;
1139  CHKERR DMMoFEMGetBlocMatData(simple->getDM(), block_data);
1140 
1141  // Boundary
1142  auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
1143  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
1144 
1145  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1146  pip->getOpBoundaryLhsPipeline().push_back(
1147  new OpMassStab("SIGMA", "SIGMA",
1148  [eps_stab](double, double, double) { return eps_stab; }));
1149  pip->getOpBoundaryLhsPipeline().push_back(
1150 
1151  createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, {ao_up}, {S}, {false},
1152  false, block_data)
1153 
1154  );
1155 
1156  // Domain
1157  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1158  pip->getOpDomainLhsPipeline().push_back(
1159 
1160  createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, {ao_up}, {S}, {false},
1161  false, block_data)
1162 
1163  );
1164 
1165  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1166  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1167 
1168  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1169  MoFEMFunctionBegin;
1170  CHKERR MatZeroEntries(S);
1171  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1172  MoFEMFunctionReturn(0);
1173  };
1174 
1175  post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1176  post_proc_schur_lhs_ptr]() {
1177  MoFEMFunctionBegin;
1178  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1179  auto print_mat_norm = [this](auto a, std::string prefix) {
1180  MoFEMFunctionBegin;
1181  double nrm;
1182  CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1183  MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1184  MoFEMFunctionReturn(0);
1185  };
1186  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1187  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1188  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1189  mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1190 #ifndef NDEBUG
1191  CHKERR print_mat_norm(S, "S");
1192 #endif // NDEBUG
1193  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1194  MoFEMFunctionReturn(0);
1195  };
1196 
1197  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1198  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1199  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1200 
1201  MoFEMFunctionReturn(0);
1202 }
1203 
1204 MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
1205  MoFEMFunctionBegin;
1206  auto simple = mField.getInterface<Simple>();
1207  auto block_is = getDMSubData(blockDM)->getSmartRowIs();
1208  CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
1209  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1210  MoFEMFunctionReturn(0);
1211 }
1212 
1213 MoFEMErrorCode SetUpSchurImpl::setDiagonalPC(PC pc) {
1214  MoFEMFunctionBegin;
1215  KSP *subksp;
1216  CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
1217  auto get_pc = [](auto ksp) {
1218  PC pc_raw;
1219  CHKERR KSPGetPC(ksp, &pc_raw);
1220  return SmartPetscObj<PC>(pc_raw, true); // bump reference
1221  };
1222  CHKERR setSchurMatSolvePC(get_pc(subksp[0]));
1223  CHKERR PetscFree(subksp);
1224  MoFEMFunctionReturn(0);
1225 }
1226 
1227 boost::shared_ptr<SetUpSchur>
1228 SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
1229  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1230 }
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::DMMoFEMGetBlocMatData
MoFEMErrorCode DMMoFEMGetBlocMatData(DM dm, boost::shared_ptr< BlockStructure > &)
Get data for block mat.
Definition: DMMoFEM.cpp:1542
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
SetUpSchurImpl
Definition: test_broken_space.cpp:512
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
spring_stiffness
double spring_stiffness
Definition: contact.cpp:84
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Contact::monitorPtr
boost::shared_ptr< Monitor > monitorPtr
Definition: contact.cpp:143
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
Contact::bC
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:432
LASTBASE
@ LASTBASE
Definition: definitions.h:69
ContactOps
Definition: EshelbianContact.hpp:10
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:525
IT
constexpr IntegrationType IT
Definition: contact.cpp:37
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
main
int main(int argc, char *argv[])
Definition: contact.cpp:899
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
order
int order
Definition: contact.cpp:78
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:90
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
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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: contact.cpp:962
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1204
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
Contact::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:138
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::TsSetI2Jacobian
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition: TsCtx.cpp:511
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
rho
double rho
Definition: contact.cpp:83
HenckyOps
Definition: HenckyOps.hpp:11
Contact::mField
MoFEM::Interface & mField
Definition: contact.cpp:129
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
ContactOps::opFactoryBoundaryToDomainLhs
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1144
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
SetUpSchurImpl::schurDM
SmartPetscObj< DM > schurDM
Definition: contact.cpp:975
SetUpSchurImpl::createSubDM
MoFEMErrorCode createSubDM()
Definition: contact.cpp:1056
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: contact.cpp:53
scale
double scale
Definition: contact.cpp:88
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
SetUpSchurImpl::blockDM
SmartPetscObj< DM > blockDM
Definition: contact.cpp:976
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
EntData
EntitiesFieldData::EntData EntData
Definition: contact.cpp:59
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
atom_test
int atom_test
Definition: contact.cpp:94
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, std::vector< SmartPetscObj< AO >> sequence_of_aos, std::vector< SmartPetscObj< Mat >> sequence_of_mats, std::vector< bool > sym_schur, bool symm_op, boost::shared_ptr< BlockStructure > diag_blocks)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2771
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
vis_spring_stiffness
double vis_spring_stiffness
Definition: contact.cpp:85
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
OpSpringRhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:73
alpha_damping
double alpha_damping
Definition: contact.cpp:86
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
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
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:312
poisson_ratio
double poisson_ratio
Definition: contact.cpp:82
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
Contact
Definition: contact.cpp:122
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
Contact::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:140
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:43
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2766
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
contact_order
int contact_order
Definition: contact.cpp:79
MatrixFunction.hpp
Contact::mfrontInterface
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: contact.cpp:142
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
Contact::Contact
Contact(MoFEM::Interface &m_field)
Definition: contact.cpp:124
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:14
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1123
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
Contact::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:311
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
BiLinearForm
Contact::ScaledTimeScale
Definition: contact.cpp:149
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
Contact::ScaledTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: contact.cpp:151
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ContactOps.hpp
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::setSchurMatSolvePC
MoFEMErrorCode setSchurMatSolvePC(SmartPetscObj< PC > pc)
Definition: Schur.cpp:2823
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
ContactOps::cn_contact
double cn_contact
Definition: EshelbianContact.hpp:19
Contact::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:469
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
DomainEleOp
is_quasi_static
PetscBool is_quasi_static
[Operators used for contact]
Definition: contact.cpp:76
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
Contact::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:171
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
young_modulus
double young_modulus
Definition: contact.cpp:81
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
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
approx_order
int approx_order
Definition: test_broken_space.cpp:50
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
Contact::tsSolve
MoFEMErrorCode tsSolve()
Definition: contact.cpp:699
Contact::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: contact.cpp:158
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: contact.cpp:960
geom_order
int geom_order
Definition: contact.cpp:80
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
ContactNaturalBC.hpp
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:40
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1196
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
help
static char help[]
[Check]
Definition: contact.cpp:897
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< VecScatter >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:760
OpSpringLhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:71
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
SetUpSchurImpl::setDiagonalPC
MoFEMErrorCode setDiagonalPC(PC pc)
Definition: contact.cpp:1213
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
AT
constexpr AssemblyType AT
Definition: contact.cpp:34
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
PostProcContact.hpp
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Contact::checkResults
MoFEMErrorCode checkResults()
[Solve]
Definition: contact.cpp:836
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
tol
double tol
Definition: mesh_smoothing.cpp:27
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition: contact.cpp:67
MoFEM::getNestSchurData
boost::shared_ptr< NestSchurData > getNestSchurData(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:2249
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Contact::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:139
Generated by Doxygen 1.8.17 and hosted at University of Glasgow