v0.14.0
PostProcContact.hpp
Go to the documentation of this file.
1 /**
2  * \file PostProcContact.hpp
3  *
4  *
5  * @copyright Copyright (c) 2023
6  */
7 
8 namespace ContactOps {
9 
10 template <int DIM> struct PostProcEleByDim;
11 
12 template <> struct PostProcEleByDim<2> {
13  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
14  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
15  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
16 };
17 
18 template <> struct PostProcEleByDim<3> {
19  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
20  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
21  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
22 };
23 
24 using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
25 using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
26 using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
27 
28 struct Monitor : public FEMethod {
29 
30  Monitor(SmartPetscObj<DM> &dm, double scale,
31  boost::shared_ptr<GenericElementInterface> mfront_interface = nullptr,
32  bool is_axisymmetric = false)
33  : dM(dm), moabVertex(mbVertexPostproc), sTEP(0),
34  mfrontInterface(mfront_interface) {
35 
36  MoFEM::Interface *m_field_ptr;
37  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
38 
39  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
40 
41  struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
42  OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
43  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
44  mPtr(m_ptr), scale(s) {}
45  MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
46  *mPtr *= 1./scale;
47  return 0;
48  }
49 
50  private:
51  boost::shared_ptr<MatrixDouble> mPtr;
52  double scale;
53  };
54 
55  auto push_domain_ops = [&](auto &pip) {
56  CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
57  pip, {H1, HDIV}, "GEOMETRY")),
58  "Apply base transform");
59  auto henky_common_data_ptr =
60  commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
61  *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
62  auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
63  pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
64  "SIGMA", contact_stress_ptr));
65  pip.push_back(new OpScale(contact_stress_ptr, scale));
66  return std::make_tuple(henky_common_data_ptr, contact_stress_ptr);
67  };
68 
69  auto push_bdy_ops = [&](auto &pip) {
70  // evaluate traction
71  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
72  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
73  "U", common_data_ptr->contactDispPtr()));
74  pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
75  "SIGMA", common_data_ptr->contactTractionPtr()));
76  pip.push_back(new OpScale(common_data_ptr->contactTractionPtr(), scale));
77  using C = ContactIntegrators<BoundaryEleOp>;
78  pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
79  common_data_ptr));
80  return common_data_ptr;
81  };
82 
83  auto get_domain_pip = [&](auto &pip)
84  -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
85  if constexpr (SPACE_DIM == 3) {
86  auto op_loop_side = new OpLoopSide<SideEle>(
87  *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
88  boost::make_shared<
89  ForcesAndSourcesCore::UserDataOperator::AdjCache>());
90  pip.push_back(op_loop_side);
91  return op_loop_side->getOpPtrVector();
92  } else {
93  return pip;
94  }
95  };
96 
97  auto get_post_proc_domain_fe = [&]() {
98  auto post_proc_fe =
99  boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
100  auto &pip = post_proc_fe->getOpPtrVector();
101 
102  auto [henky_common_data_ptr, contact_stress_ptr] =
103  push_domain_ops(get_domain_pip(pip));
104 
105  auto u_ptr = boost::make_shared<MatrixDouble>();
106  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
107  auto X_ptr = boost::make_shared<MatrixDouble>();
108  pip.push_back(
109  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
110 
111 
112 
113  pip.push_back(
114 
115  new OpPPMap(
116 
117  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
118 
119  {},
120  {
121 
122  {"U", u_ptr}, {"GEOMETRY", X_ptr}
123 
124  },
125  {
126 
127  {"SIGMA", contact_stress_ptr},
128 
129  {"G", henky_common_data_ptr->matGradPtr},
130 
131  {"PK1", henky_common_data_ptr->getMatFirstPiolaStress()}
132 
133  },
134  {}
135 
136  )
137 
138  );
139 
140  if (SPACE_DIM == 3) {
141 
142  CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
143  pip, {HDIV}, "GEOMETRY")),
144  "Apply transform");
145  auto common_data_ptr = push_bdy_ops(pip);
146 
147  pip.push_back(
148 
149  new OpPPMap(
150 
151  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
152 
153  {{"SDF", common_data_ptr->sdfPtr()},
154  {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
155 
156  {
157 
158  {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
159  {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
160 
161  },
162 
163  {},
164 
165  {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
166 
167  )
168 
169  );
170  }
171 
172  return post_proc_fe;
173  };
174 
175  auto get_post_proc_bdy_fe = [&]() {
176  auto post_proc_fe =
177  boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
178  auto &pip = post_proc_fe->getOpPtrVector();
179 
180  CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
181  pip, {HDIV}, "GEOMETRY")),
182  "Apply transform");
183  auto common_data_ptr = push_bdy_ops(pip);
184 
185  // create OP which run element on side
186  auto op_loop_side = new OpLoopSide<SideEle>(
187  *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
188  boost::make_shared<
189  ForcesAndSourcesCore::UserDataOperator::AdjCache>());
190  pip.push_back(op_loop_side);
191 
192  auto [henky_common_data_ptr, contact_stress_ptr] =
193  push_domain_ops(op_loop_side->getOpPtrVector());
194 
195  auto X_ptr = boost::make_shared<MatrixDouble>();
196  pip.push_back(
197  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
198 
199  pip.push_back(
200 
201  new OpPPMap(
202 
203  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
204 
205  {{"SDF", common_data_ptr->sdfPtr()},
206  {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
207 
208  {{"U", common_data_ptr->contactDispPtr()},
209  {"GEOMETRY", X_ptr},
210  {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
211  {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
212 
213  },
214 
215  {},
216 
217  {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
218 
219  )
220 
221  );
222 
223  return post_proc_fe;
224  };
225 
226  auto get_integrate_traction = [&]() {
227  auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
228  CHK_THROW_MESSAGE(
229  (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
230  integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
231  "Apply transform");
232  // We have to integrate on curved face geometry, thus integration weight
233  // have to adjusted.
234  integrate_traction->getOpPtrVector().push_back(
235  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
236  integrate_traction->getRuleHook = [](int, int, int approx_order) {
237  return 2 * approx_order + geom_order - 1;
238  };
239 
240  CHK_THROW_MESSAGE(
241  (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
242  integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
243  "push operators to calculate traction");
244 
245  return integrate_traction;
246  };
247 
248  auto get_integrate_area = [&]() {
249  auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
250 
251  CHK_THROW_MESSAGE(
252  (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
253  integrate_area->getOpPtrVector(), {HDIV}, "GEOMETRY")),
254  "Apply transform");
255  // We have to integrate on curved face geometry, thus integration weight have to adjusted.
256  integrate_area->getOpPtrVector().push_back(
257  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
258  integrate_area->getRuleHook = [](int, int, int approx_order) {
259  return 2 * approx_order + geom_order - 1;
260  };
261  Range contact_range;
262  for (auto m :
263  m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
264  std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
265  auto meshset = m->getMeshset();
266  Range contact_meshset_range;
267  CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
268  meshset, SPACE_DIM - 1, contact_meshset_range, true);
269 
270  CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
271  contact_meshset_range);
272  contact_range.merge(contact_meshset_range);
273  }
274 
275  auto contact_range_ptr = boost::make_shared<Range>(contact_range);
276 
277  auto op_loop_side = new OpLoopSide<SideEle>(
278  *m_field_ptr, m_field_ptr->getInterface<Simple>()->getDomainFEName(),
279  SPACE_DIM);
280  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
281  op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
282 
283  CHK_THROW_MESSAGE(
284  (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
285  integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
286  is_axisymmetric, contact_range_ptr)),
287  "push operators to calculate area");
288 
289  return integrate_area;
290  };
291 
292  postProcDomainFe = get_post_proc_domain_fe();
293  if constexpr (SPACE_DIM == 2)
294  postProcBdyFe = get_post_proc_bdy_fe();
295 
296  integrateTraction = get_integrate_traction();
297  integrateArea = get_integrate_area();
298 
299  normsVec = createVectorMPI(
300  m_field_ptr->get_comm(),
301  (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
302  }
303 
304  MoFEMErrorCode preProcess() { return 0; }
305  MoFEMErrorCode operator()() { return 0; }
306 
307  MoFEMErrorCode postProcess() {
308  MoFEMFunctionBegin;
309  MoFEM::Interface *m_field_ptr;
310  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
311 
312  auto post_proc = [&]() {
313  MoFEMFunctionBegin;
314 
315  if (!mfrontInterface) {
316  auto post_proc_begin =
317  boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
318  postProcMesh);
319  auto post_proc_end =
320  boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
321  postProcMesh);
322 
323  CHKERR DMoFEMPreProcessFiniteElements(dM,
324  post_proc_begin->getFEMethod());
325  if (!postProcBdyFe) {
326  postProcDomainFe->copyTs(*this); // this here is a Monitor
327  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
328  } else {
329  postProcDomainFe->copyTs(*this); // this here is a Monitor
330  postProcBdyFe->copyTs(*this);
331  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
332  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
333  }
334  CHKERR DMoFEMPostProcessFiniteElements(dM,
335  post_proc_end->getFEMethod());
336 
337  CHKERR post_proc_end->writeFile(
338  "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
339  } else {
340  CHKERR mfrontInterface->postProcessElement(
341  ts_step, dM,
342  m_field_ptr->getInterface<Simple>()->getDomainFEName());
343  }
344 
345  MoFEMFunctionReturn(0);
346  };
347 
348  auto calculate_force = [&] {
349  MoFEMFunctionBegin;
350  CHKERR VecZeroEntries(CommonData::totalTraction);
351  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
352  CHKERR VecAssemblyBegin(CommonData::totalTraction);
353  CHKERR VecAssemblyEnd(CommonData::totalTraction);
354  MoFEMFunctionReturn(0);
355  };
356 
357  auto calculate_area = [&] {
358  MoFEMFunctionBegin;
359  integrateArea->copyTs(*this);
360  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateArea);
361  CHKERR VecAssemblyBegin(CommonData::totalTraction);
362  CHKERR VecAssemblyEnd(CommonData::totalTraction);
363  MoFEMFunctionReturn(0);
364  };
365 
366  auto calculate_reactions = [&]() {
367  MoFEMFunctionBegin;
368 
369  auto res = createDMVector(dM);
370 
371  auto assemble_domain = [&]() {
372  MoFEMFunctionBegin;
373  auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
374  auto &pip = fe_rhs->getOpPtrVector();
375  fe_rhs->f = res;
376 
377  auto integration_rule = [](int, int, int approx_order) {
378  return 2 * approx_order + geom_order - 1;
379  };
380  fe_rhs->getRuleHook = integration_rule;
381 
382  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
383  "GEOMETRY");
384  CHKERR
385  ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
386  pip, "SIGMA", "U", is_axisymmetric);
387 
388  if (!mfrontInterface) {
389  CHKERR
390  HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
391  *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
392  } else {
393  CHKERR mfrontInterface->opFactoryDomainRhs(pip);
394  }
395  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
396 
397  CHKERR VecAssemblyBegin(res);
398  CHKERR VecAssemblyEnd(res);
399  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
400  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
401 
402  MoFEMFunctionReturn(0);
403  };
404 
405  auto assemble_boundary = [&]() {
406  MoFEMFunctionBegin;
407  auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
408  auto &pip = fe_rhs->getOpPtrVector();
409  fe_rhs->f = res;
410 
411  auto integration_rule = [](int, int, int approx_order) {
412  return 2 * approx_order + geom_order - 1;
413  };
414  fe_rhs->getRuleHook = integration_rule;
415 
416  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
417  "GEOMETRY");
418  // We have to integrate on curved face geometry, thus integration weight
419  // have to adjusted.
420  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
421 
422  auto u_disp = boost::make_shared<MatrixDouble>();
423  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
424  pip.push_back(
425  new OpSpringRhs("U", u_disp, [this](double, double, double) {
426  return spring_stiffness;
427  }));
428 
429  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
430 
431  MoFEMFunctionReturn(0);
432  };
433 
434  CHKERR assemble_domain();
435  CHKERR assemble_boundary();
436 
437  auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
438  auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
439  m_field_ptr]() {
440  MoFEMFunctionBegin;
441  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
442  *m_field_ptr, fe_post_proc_ptr, res)();
443  MoFEMFunctionReturn(0);
444  };
445  fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
446  CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
447 
448  MoFEMFunctionReturn(0);
449  };
450 
451  auto print_max_min = [&](auto &tuple, const std::string msg) {
452  MoFEMFunctionBegin;
453  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
454  INSERT_VALUES, SCATTER_FORWARD);
455  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
456  INSERT_VALUES, SCATTER_FORWARD);
457  double max, min;
458  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
459  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
460  MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %6.4e min %6.4e max %6.4e",
461  msg.c_str(), ts_t, min, max);
462  MoFEMFunctionReturn(0);
463  };
464 
465  auto print_force_and_area = [&]() {
466  MoFEMFunctionBegin;
467  MoFEM::Interface *m_field_ptr;
468  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
469  if (!m_field_ptr->get_comm_rank()) {
470  const double *t_ptr;
471  CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
472  MOFEM_LOG_C("CONTACT", Sev::inform,
473  "Contact force: time %6.3e Fx: %6.6e Fy: %6.6e Fz: %6.6e",
474  ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
475  MOFEM_LOG_C("CONTACT", Sev::inform,
476  "Contact area: time %6.3e Active: %6.6e Potential: %6.6e",
477  ts_t, t_ptr[3], t_ptr[4]);
478  CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
479  }
480  MoFEMFunctionReturn(0);
481  };
482 
483  if (mfrontInterface) {
484  CHKERR mfrontInterface->updateElementVariables(
485  dM, m_field_ptr->getInterface<Simple>()->getDomainFEName());
486  }
487 
488  auto calculate_error = [&](MoFEM::ScalarFun &fun) {
489  MoFEMFunctionBegin;
490  struct OpCalcTractions : public BoundaryEleOp {
491  OpCalcTractions(boost::shared_ptr<MatrixDouble> m_ptr,
492  boost::shared_ptr<VectorDouble> p_ptr,
493  boost::shared_ptr<VectorDouble> mag_ptr,
494  boost::shared_ptr<VectorDouble> traction_y_ptr,
495  boost::shared_ptr<MatrixDouble> t_ptr,
496  boost::shared_ptr<MatrixDouble> grad_sdf_ptr)
497  : BoundaryEleOp(NOSPACE, OPSPACE), mPtr(m_ptr), pPtr(p_ptr),
498  magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
499  gradSDFPtr(grad_sdf_ptr) {}
500  MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
501  MoFEMFunctionBegin;
502  FTensor::Index<'i', SPACE_DIM> i;
503  mPtr->resize(SPACE_DIM, pPtr->size());
504  mPtr->clear();
505  magPtr->resize(pPtr->size());
506  magPtr->clear();
507  tyPtr->resize(pPtr->size());
508  tyPtr->clear();
509 
510  auto t_traction = getFTensor1FromMat<SPACE_DIM>(*mPtr);
511  auto t_contact_traction = getFTensor1FromMat<SPACE_DIM>(*tPtr);
512  auto t_p = getFTensor0FromVec(*pPtr);
513  int nb_gauss_pts = pPtr->size();
514  auto t_normal = getFTensor1FromMat<SPACE_DIM>(*gradSDFPtr);
515  auto t_normal_at_gauss = getFTensor1NormalsAtGaussPts();
516  auto t_mag = getFTensor0FromVec(*magPtr);
517  auto t_ty = getFTensor0FromVec(*tyPtr);
518 
519  for (int gg = 0; gg != nb_gauss_pts; gg++) {
520  FTensor::Tensor1<double, SPACE_DIM> normal;
521  t_traction(i) = t_p * (-(t_normal(i) / t_normal.l2()));
522  t_mag = t_contact_traction.l2();
523  t_ty = t_contact_traction(1);
524 
525  ++t_normal;
526  ++t_traction;
527  ++t_p;
528  ++t_mag;
529  ++t_contact_traction;
530  ++t_ty;
531  ++t_normal_at_gauss;
532  }
533  MoFEMFunctionReturn(0);
534  }
535 
536  private:
537  boost::shared_ptr<MatrixDouble> mPtr;
538  boost::shared_ptr<VectorDouble> pPtr;
539  boost::shared_ptr<VectorDouble> magPtr;
540  boost::shared_ptr<MatrixDouble> tPtr;
541  boost::shared_ptr<VectorDouble> tyPtr;
542  boost::shared_ptr<MatrixDouble> gradSDFPtr;
543  };
544 
545  auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
546  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
547  auto simple = m_field_ptr->getInterface<Simple>();
548  Range contact_range;
549  for (auto m :
550  m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
551  std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
552  auto meshset = m->getMeshset();
553  Range contact_meshset_range;
554  CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
555  meshset, SPACE_DIM - 1, contact_meshset_range, true);
556 
557  CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
558  contact_meshset_range);
559  contact_range.merge(contact_meshset_range);
560  }
561 
562  CHK_THROW_MESSAGE(
563  (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
564  post_proc_norm_fe->getOpPtrVector(), {HDIV}, "GEOMETRY")),
565  "Apply transform");
566  // We have to integrate on curved face geometry, thus integration weight
567  // have to adjusted.
568  post_proc_norm_fe->getOpPtrVector().push_back(
569  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
570  post_proc_norm_fe->getRuleHook = [](int, int, int approx_order) {
571  return 2 * approx_order + geom_order - 1;
572  };
573 
574  post_proc_norm_fe->getOpPtrVector().push_back(
575  new OpCalculateVectorFieldValues<SPACE_DIM>(
576  "U", common_data_ptr->contactDispPtr()));
577  post_proc_norm_fe->getOpPtrVector().push_back(
578  new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
579  "SIGMA", common_data_ptr->contactTractionPtr()));
580  using C = ContactIntegrators<BoundaryEleOp>;
581  post_proc_norm_fe->getOpPtrVector().push_back(
582  new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
583  common_data_ptr));
584 
585  auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
586  auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
587  auto mag_traction_ptr = boost::make_shared<VectorDouble>();
588  auto traction_y_ptr = boost::make_shared<VectorDouble>();
589  auto contact_range_ptr = boost::make_shared<Range>(contact_range);
590 
591  post_proc_norm_fe->getOpPtrVector().push_back(
592  new OpGetTensor0fromFunc(analytical_pressure_ptr, fun));
593 
594  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcTractions(
595  analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
596  traction_y_ptr, common_data_ptr->contactTractionPtr(),
597  common_data_ptr->gradSdfPtr()));
598 
599  post_proc_norm_fe->getOpPtrVector().push_back(
600  new OpCalcNormL2Tensor1<SPACE_DIM>(
601  common_data_ptr->contactTractionPtr(), normsVec, TRACTION_NORM_L2,
602  analytical_traction_ptr, contact_range_ptr));
603 
604  // calculate magnitude of traction
605 
606  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcNormL2Tensor0(
607  mag_traction_ptr, normsVec, MAG_TRACTION_NORM_L2,
608  analytical_pressure_ptr, contact_range_ptr));
609 
610  post_proc_norm_fe->getOpPtrVector().push_back(
611  new OpCalcNormL2Tensor0(traction_y_ptr, normsVec, TRACTION_Y_NORM_L2,
612  analytical_pressure_ptr, contact_range_ptr));
613 
614  CHKERR VecZeroEntries(normsVec);
615  post_proc_norm_fe->copyTs(*this); // set time as is in Monitor
616  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", post_proc_norm_fe);
617  CHKERR VecAssemblyBegin(normsVec);
618  CHKERR VecAssemblyEnd(normsVec);
619 
620  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
621  if (m_field_ptr->get_comm_rank() == 0) {
622  const double *norms;
623  CHKERR VecGetArrayRead(normsVec, &norms);
624  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
625  << "norm_traction: " << std::scientific
626  << std::sqrt(norms[TRACTION_NORM_L2]);
627  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
628  << "norm_mag_traction: " << std::scientific
629  << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
630  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
631  << "norm_traction_y: " << std::scientific
632  << std::sqrt(norms[TRACTION_Y_NORM_L2]);
633  CHKERR VecRestoreArrayRead(normsVec, &norms);
634  }
635  MoFEMFunctionReturn(0);
636  };
637 
638  int se = 1;
639  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &se, PETSC_NULL);
640 
641  if (!(ts_step % se)) {
642  MOFEM_LOG("CONTACT", Sev::inform)
643  << "Write file at time " << ts_t << " write step " << sTEP;
644  CHKERR post_proc();
645  }
646  CHKERR calculate_force();
647  CHKERR calculate_area();
648 
649  CHKERR calculate_reactions();
650 
651  if (atom_test && sTEP) {
652  switch (atom_test) {
653  case 1:
654  CHKERR calculate_error(analyticalHertzPressurePlaneStress);
655  break;
656  case 2:
657  CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
658  break;
659  case 5:
660  CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
661  break;
662  case 6:
663  CHKERR calculate_error(analyticalWavy2DPressure);
664  break;
665  default:
666  break;
667  }
668  }
669 
670  CHKERR print_max_min(uXScatter, "Ux");
671  CHKERR print_max_min(uYScatter, "Uy");
672  if (SPACE_DIM == 3)
673  CHKERR print_max_min(uZScatter, "Uz");
674  CHKERR print_force_and_area();
675  ++sTEP;
676 
677  MoFEMFunctionReturn(0);
678  }
679 
680  //! [Analytical function]
681  MoFEM::ScalarFun analyticalWavy2DPressure = [](double x, double y, double z) {
682  // update to atom test values
683  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
684  // delta
685  double delta = 0.0002;
686  // lambda
687  double lambda = 2;
688  // pressure star
689  double p_star = M_PI * E_star * delta / lambda;
690 
691  // Pressure = p_bar + p_star * cos(2 * pi * x / lambda)
692  return p_star + p_star * std::cos(2. * M_PI * x / lambda);
693  };
694 
695  MoFEM::ScalarFun analyticalHertzPressureAxisymmetric = [](double x, double y,
696  double z) {
697  // update to atom test values
698  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
699  // Radius
700  double R = 100.;
701  // Indentation
702  double d = 0.01;
703  // Force
704  double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
705  // Contact area radius
706  double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
707  // Maximum pressure
708  double p_max = (3. * F) / (2. * M_PI * a * a);
709 
710  double r = std::sqrt((x * x) + (y * y));
711 
712  if (r > a) {
713  return 0.;
714  }
715  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
716  return p_max * std::sqrt(1 - ((r * r) / (a * a)));
717  };
718 
719  MoFEM::ScalarFun analyticalHertzPressurePlaneStrain = [](double x, double y,
720  double z) {
721  // update to atom test values
722  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
723  // Radius
724  double R = 100.;
725  // Indentation
726  double d = 0.02745732273553991;
727  // Contact area radius
728  double a = 1;
729  // current radius
730  double r = std::sqrt((x * x) + (y * y));
731 
732  if (r > a) {
733  return 0.;
734  }
735  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
736  return E_star / (2. * R) * std::sqrt(a * a - r * r);
737  };
738  MoFEM::ScalarFun analyticalHertzPressurePlaneStress = [](double x, double y,
739  double z) {
740  // update to atom test values
741  double E_star = young_modulus;
742  // Radius
743  double R = 100.;
744  // Indentation
745  double d = 0.02745732273553991;
746  // Contact area radius
747  double a = 1;
748  // current radius
749  double r = std::sqrt((x * x) + (y * y));
750 
751  if (r > a) {
752  return 0.;
753  }
754  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
755  return E_star / (2. * R) * std::sqrt(a * a - r * r);
756  };
757 
758  // ***DISPLACMENT NOT TESTED***
759  MoFEM::VectorFun<SPACE_DIM> analyticalHertzDisplacement3D = [](double x,
760  double y,
761  double z) {
762  // update to atom test values
763  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
764  // Radius
765  double R = 100.;
766  // Contact area radius
767  double a = 1;
768  // max pressure
769  double p_0 = (2. * E_star * a) / (M_PI * R);
770  // current radius
771  double r = std::sqrt((x * x) + (y * y));
772  FTensor::Tensor1<double, SPACE_DIM> u;
773  std::vector<double> v_u;
774 
775  double u_z = 0.;
776  double u_r = 0.;
777  // outside contact zone
778  if (r > a) {
779  u_z = (1. - std::pow(poisson_ratio, 2.)) / young_modulus *
780  ((p_0) / (2. * a)) *
781  ((2. * std::pow(a, 2.) - std::pow(r, 2.)) * asin(a / r) +
782  std::pow(r, 2.) * (a / r) *
783  std::pow(1 - (std::pow(a, 2.) / std::pow(r, 2.)), 2.));
784  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
785  (3. * young_modulus) * ((std::pow(a, 2) / r)) * p_0;
786 
787  if (SPACE_DIM == 2)
788  v_u = {u_r, u_z};
789  else
790  v_u = {u_r, u_z, u_r};
791 
792  for (int i = 0; i < SPACE_DIM; ++i)
793  u(i) = v_u[i];
794 
795  return u;
796  }
797 
798  // In contact zone
799  u_z = ((1. - std::pow(poisson_ratio, 2.)) / young_modulus) *
800  ((M_PI * p_0) / 4. * a) * (2. * std::pow(a, 2.) - std::pow(r, 2.));
801  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
802  (3. * young_modulus) * ((std::pow(a, 2.) / r)) * p_0 *
803  (1 - std::pow(1 - (std::pow(r, 2.) / std::pow(a, 2.)), 1.5));
804 
805  if (SPACE_DIM == 2)
806  v_u = {u_r, u_z};
807  else
808  v_u = {u_r, u_z, u_r};
809 
810  for (int i = 0; i < SPACE_DIM; ++i)
811  u(i) = v_u[i];
812 
813  return u;
814  };
815 
816  MoFEMErrorCode setScatterVectors(
817  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
818  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
819  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter) {
820  MoFEMFunctionBegin;
821  uXScatter = ux_scatter;
822  uYScatter = uy_scatter;
823  uZScatter = uz_scatter;
824  MoFEMFunctionReturn(0);
825  }
826 
827  MoFEMErrorCode getErrorNorm(int normType) {
828  const double *norm;
829  CHKERR VecGetArrayRead(normsVec, &norm);
830  double norm_val = std::sqrt(norm[normType]);
831  CHKERR VecRestoreArrayRead(normsVec, &norm);
832  return norm_val;
833  }
834 
835 private:
836  SmartPetscObj<DM> dM;
837  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
838  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
839  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
840 
841  boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
842 
843  boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
844  boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
845 
846  boost::shared_ptr<BoundaryEle> integrateTraction;
847  boost::shared_ptr<BoundaryEle> integrateArea;
848 
849  enum NORMS {
850  TRACTION_NORM_L2 = 0,
851  MAG_TRACTION_NORM_L2,
852  TRACTION_Y_NORM_L2,
853  LAST_NORM
854  };
855  SmartPetscObj<Vec> normsVec;
856 
857  moab::Core mbVertexPostproc;
858  moab::Interface &moabVertex;
859 
860  double lastTime;
861  double deltaTime;
862  int sTEP;
863 
864  boost::shared_ptr<GenericElementInterface> mfrontInterface;
865 };
866 
867 } // namespace ContactOps
NOSPACE
@ NOSPACE
Definition: definitions.h:83
ContactOps::Monitor::postProcBdyFe
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
Definition: PostProcContact.hpp:844
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
ContactOps::Monitor::normsVec
SmartPetscObj< Vec > normsVec
Definition: PostProcContact.hpp:855
ContactOps::PostProcEleByDim
Definition: PostProcContact.hpp:10
spring_stiffness
double spring_stiffness
Definition: contact.cpp:87
H1
@ H1
continuous field
Definition: definitions.h:85
ContactOps::PostProcEleByDim< 3 >::PostProcEleDomain
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleDomain
Definition: PostProcContact.hpp:19
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
ContactOps::Monitor::Monitor
Monitor(SmartPetscObj< DM > &dm, double scale, boost::shared_ptr< GenericElementInterface > mfront_interface=nullptr, bool is_axisymmetric=false)
Definition: PostProcContact.hpp:30
ContactOps::PostProcEleByDim< 2 >::PostProcEleDomain
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomain
Definition: PostProcContact.hpp:13
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
ContactOps
Definition: contact.cpp:99
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:93
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ContactOps::Monitor::mfrontInterface
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: PostProcContact.hpp:864
FaceSideEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::FaceSideEle FaceSideEle
Definition: level_set.cpp:41
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
ContactOps::PostProcEleByDim< 3 >::PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
Definition: PostProcContact.hpp:20
ContactOps::Monitor::sTEP
int sTEP
Definition: PostProcContact.hpp:862
ContactOps::PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: PostProcContact.hpp:24
ContactOps::Monitor::postProcMesh
boost::shared_ptr< moab::Core > postProcMesh
Definition: PostProcContact.hpp:841
ContactOps::Monitor::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: PostProcContact.hpp:837
scale
double scale
Definition: plastic.cpp:119
ContactOps::Monitor::MAG_TRACTION_NORM_L2
@ MAG_TRACTION_NORM_L2
Definition: PostProcContact.hpp:851
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:137
sdf.r
int r
Definition: sdf.py:8
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ContactOps::Monitor::mbVertexPostproc
moab::Core mbVertexPostproc
Definition: PostProcContact.hpp:857
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
atom_test
int atom_test
Definition: contact.cpp:97
ContactOps::Monitor::moabVertex
moab::Interface & moabVertex
Definition: PostProcContact.hpp:858
ContactOps::Monitor::deltaTime
double deltaTime
Definition: PostProcContact.hpp:861
FEMethod
ContactOps::Monitor
Definition: PostProcContact.hpp:28
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:23
ContactOps::Monitor::getErrorNorm
MoFEMErrorCode getErrorNorm(int normType)
Definition: PostProcContact.hpp:827
ContactOps::PostProcEleByDim< 2 >::PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
Definition: PostProcContact.hpp:14
ContactOps::Monitor::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: PostProcContact.hpp:839
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
ContactOps::Monitor::lastTime
double lastTime
Definition: PostProcContact.hpp:860
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::VectorFun
boost::function< FTensor::Tensor1< double, DIM >(const double, const double, const double)> VectorFun
Vector function type.
Definition: FormsIntegrators.hpp:176
ContactOps::Monitor::dM
SmartPetscObj< DM > dM
Definition: PostProcContact.hpp:836
a
constexpr double a
Definition: approx_sphere.cpp:30
BoundaryEleOp
R
@ R
Definition: free_surface.cpp:394
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
OpSpringRhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:75
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: child_and_parent.cpp:40
ContactOps::ContactIntegrators
Definition: ContactOps.hpp:549
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
ContactOps::Monitor::postProcDomainFe
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
Definition: PostProcContact.hpp:843
ContactOps::Monitor::postProcess
MoFEMErrorCode postProcess()
Definition: PostProcContact.hpp:307
MoFEM::ScalarFun
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
Definition: FormsIntegrators.hpp:144
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ContactOps::PostProcEleByDim< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: PostProcContact.hpp:21
ContactOps::Monitor::TRACTION_Y_NORM_L2
@ TRACTION_Y_NORM_L2
Definition: PostProcContact.hpp:852
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', SPACE_DIM >
ContactOps::SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: PostProcContact.hpp:25
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
ContactOps::Monitor::setScatterVectors
MoFEMErrorCode setScatterVectors(std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uz_scatter)
Definition: PostProcContact.hpp:816
Range
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
ContactOps::Monitor::NORMS
NORMS
Definition: PostProcContact.hpp:849
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
ContactOps::PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: PostProcContact.hpp:26
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
ContactOps::PostProcEleByDim< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: PostProcContact.hpp:15
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
approx_order
int approx_order
Definition: test_broken_space.cpp:50
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
ContactOps::Monitor::integrateTraction
boost::shared_ptr< BoundaryEle > integrateTraction
Definition: PostProcContact.hpp:846
ContactOps::Monitor::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: PostProcContact.hpp:838
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
ContactOps::Monitor::operator()
MoFEMErrorCode operator()()
Definition: PostProcContact.hpp:305
fun
auto fun
Function to approximate.
Definition: dg_projection.cpp:36
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
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
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
ContactOps::Monitor::integrateArea
boost::shared_ptr< BoundaryEle > integrateArea
Definition: PostProcContact.hpp:847
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
Generated by Doxygen 1.8.17 and hosted at University of Glasgow