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>;
16 };
17 
18 template <> struct PostProcEleByDim<3> {
19  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
20  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
22 };
23 
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)
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) {
57  pip, {H1, HDIV}, "GEOMETRY")),
58  "Apply base transform");
59  auto hencky_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(hencky_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));
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 [hencky_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", hencky_common_data_ptr->matGradPtr},
130 
131  {"PK1", hencky_common_data_ptr->getMatFirstPiolaStress()}
132 
133  },
134  {}
135 
136  )
137 
138  );
139 
140  if (SPACE_DIM == 3) {
141 
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 
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 [hencky_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);
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 
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 
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);
281  op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
282 
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 
309  MoFEM::Interface *m_field_ptr;
310  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
311 
312  auto post_proc = [&]() {
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 
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  }
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 
346  };
347 
348  auto calculate_force = [&] {
350  CHKERR VecZeroEntries(CommonData::totalTraction);
351  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
352  CHKERR VecAssemblyBegin(CommonData::totalTraction);
353  CHKERR VecAssemblyEnd(CommonData::totalTraction);
355  };
356 
357  auto calculate_area = [&] {
359  integrateArea->copyTs(*this);
360  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateArea);
361  CHKERR VecAssemblyBegin(CommonData::totalTraction);
362  CHKERR VecAssemblyEnd(CommonData::totalTraction);
364  };
365 
366  auto calculate_reactions = [&]() {
368 
369  auto res = createDMVector(dM);
370 
371  auto assemble_domain = [&]() {
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 
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 
403  };
404 
405  auto assemble_boundary = [&]() {
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 
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 
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]() {
441  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
442  *m_field_ptr, fe_post_proc_ptr, res)();
444  };
445  fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
446  CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
447 
449  };
450 
451  auto print_max_min = [&](auto &tuple, const std::string msg) {
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);
463  };
464 
465  auto print_force_and_area = [&]() {
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  }
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) {
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 &) {
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++) {
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  }
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  Range contact_range;
548  for (auto m :
549  m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
550  std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
551  auto meshset = m->getMeshset();
552  Range contact_meshset_range;
553  CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
554  meshset, SPACE_DIM - 1, contact_meshset_range, true);
555 
556  CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
557  contact_meshset_range);
558  contact_range.merge(contact_meshset_range);
559  }
560 
563  post_proc_norm_fe->getOpPtrVector(), {HDIV}, "GEOMETRY")),
564  "Apply transform");
565  // We have to integrate on curved face geometry, thus integration weight
566  // have to adjusted.
567  post_proc_norm_fe->getOpPtrVector().push_back(
568  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
569  post_proc_norm_fe->getRuleHook = [](int, int, int approx_order) {
570  return 2 * approx_order + geom_order - 1;
571  };
572 
573  post_proc_norm_fe->getOpPtrVector().push_back(
574  new OpCalculateVectorFieldValues<SPACE_DIM>(
575  "U", common_data_ptr->contactDispPtr()));
576  post_proc_norm_fe->getOpPtrVector().push_back(
577  new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
578  "SIGMA", common_data_ptr->contactTractionPtr()));
580  post_proc_norm_fe->getOpPtrVector().push_back(
581  new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
582  common_data_ptr));
583 
584  auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
585  auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
586  auto mag_traction_ptr = boost::make_shared<VectorDouble>();
587  auto traction_y_ptr = boost::make_shared<VectorDouble>();
588  auto contact_range_ptr = boost::make_shared<Range>(contact_range);
589 
590  post_proc_norm_fe->getOpPtrVector().push_back(
591  new OpGetTensor0fromFunc(analytical_pressure_ptr, fun));
592 
593  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcTractions(
594  analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
595  traction_y_ptr, common_data_ptr->contactTractionPtr(),
596  common_data_ptr->gradSdfPtr()));
597 
598  post_proc_norm_fe->getOpPtrVector().push_back(
599  new OpCalcNormL2Tensor1<SPACE_DIM>(
600  common_data_ptr->contactTractionPtr(), normsVec, TRACTION_NORM_L2,
601  analytical_traction_ptr, contact_range_ptr));
602 
603  // calculate magnitude of traction
604 
605  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcNormL2Tensor0(
606  mag_traction_ptr, normsVec, MAG_TRACTION_NORM_L2,
607  analytical_pressure_ptr, contact_range_ptr));
608 
609  post_proc_norm_fe->getOpPtrVector().push_back(
610  new OpCalcNormL2Tensor0(traction_y_ptr, normsVec, TRACTION_Y_NORM_L2,
611  analytical_pressure_ptr, contact_range_ptr));
612 
613  CHKERR VecZeroEntries(normsVec);
614  post_proc_norm_fe->copyTs(*this); // set time as is in Monitor
615  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", post_proc_norm_fe);
616  CHKERR VecAssemblyBegin(normsVec);
617  CHKERR VecAssemblyEnd(normsVec);
618 
619  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
620  if (m_field_ptr->get_comm_rank() == 0) {
621  const double *norms;
622  CHKERR VecGetArrayRead(normsVec, &norms);
623  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
624  << "norm_traction: " << std::scientific
625  << std::sqrt(norms[TRACTION_NORM_L2]);
626  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
627  << "norm_mag_traction: " << std::scientific
628  << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
629  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
630  << "norm_traction_y: " << std::scientific
631  << std::sqrt(norms[TRACTION_Y_NORM_L2]);
632  CHKERR VecRestoreArrayRead(normsVec, &norms);
633  }
635  };
636 
637  int se = 1;
638  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &se, PETSC_NULL);
639 
640  if (!(ts_step % se)) {
641  MOFEM_LOG("CONTACT", Sev::inform)
642  << "Write file at time " << ts_t << " write step " << sTEP;
643  CHKERR post_proc();
644  }
645  CHKERR calculate_force();
646  CHKERR calculate_area();
647 
648  CHKERR calculate_reactions();
649 
650  if (atom_test && sTEP) {
651  switch (atom_test) {
652  case 1:
653  CHKERR calculate_error(analyticalHertzPressurePlaneStress);
654  break;
655  case 2:
656  CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
657  break;
658  case 5:
659  CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
660  break;
661  case 6:
662  CHKERR calculate_error(analyticalWavy2DPressure);
663  break;
664  default:
665  break;
666  }
667  }
668 
669  CHKERR print_max_min(uXScatter, "Ux");
670  CHKERR print_max_min(uYScatter, "Uy");
671  if (SPACE_DIM == 3)
672  CHKERR print_max_min(uZScatter, "Uz");
673  CHKERR print_force_and_area();
674  ++sTEP;
675 
677  }
678 
679  //! [Analytical function]
680  MoFEM::ScalarFun analyticalWavy2DPressure = [](double x, double y, double z) {
681  // update to atom test values
682  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
683  // delta
684  double delta = 0.0002;
685  // lambda
686  double lambda = 2;
687  // pressure star
688  double p_star = M_PI * E_star * delta / lambda;
689 
690  // Pressure = p_bar + p_star * cos(2 * pi * x / lambda)
691  return p_star + p_star * std::cos(2. * M_PI * x / lambda);
692  };
693 
694  MoFEM::ScalarFun analyticalHertzPressureAxisymmetric = [](double x, double y,
695  double z) {
696  // update to atom test values
697  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
698  // Radius
699  double R = 100.;
700  // Indentation
701  double d = 0.01;
702  // Force
703  double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
704  // Contact area radius
705  double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
706  // Maximum pressure
707  double p_max = (3. * F) / (2. * M_PI * a * a);
708 
709  double r = std::sqrt((x * x) + (y * y));
710 
711  if (r > a) {
712  return 0.;
713  }
714  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
715  return p_max * std::sqrt(1 - ((r * r) / (a * a)));
716  };
717 
718  MoFEM::ScalarFun analyticalHertzPressurePlaneStrain = [](double x, double y,
719  double z) {
720  // update to atom test values
721  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
722  // Radius
723  double R = 100.;
724  // Indentation
725  double d = 0.02745732273553991;
726  // Contact area radius
727  double a = 1;
728  // current radius
729  double r = std::sqrt((x * x) + (y * y));
730 
731  if (r > a) {
732  return 0.;
733  }
734  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
735  return E_star / (2. * R) * std::sqrt(a * a - r * r);
736  };
737  MoFEM::ScalarFun analyticalHertzPressurePlaneStress = [](double x, double y,
738  double z) {
739  // update to atom test values
740  double E_star = young_modulus;
741  // Radius
742  double R = 100.;
743  // Indentation
744  double d = 0.02745732273553991;
745  // Contact area radius
746  double a = 1;
747  // current radius
748  double r = std::sqrt((x * x) + (y * y));
749 
750  if (r > a) {
751  return 0.;
752  }
753  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
754  return E_star / (2. * R) * std::sqrt(a * a - r * r);
755  };
756 
757  // ***DISPLACMENT NOT TESTED***
758  MoFEM::VectorFun<SPACE_DIM> analyticalHertzDisplacement3D = [](double x,
759  double y,
760  double z) {
761  // update to atom test values
762  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
763  // Radius
764  double R = 100.;
765  // Contact area radius
766  double a = 1;
767  // max pressure
768  double p_0 = (2. * E_star * a) / (M_PI * R);
769  // current radius
770  double r = std::sqrt((x * x) + (y * y));
772  std::vector<double> v_u;
773 
774  double u_z = 0.;
775  double u_r = 0.;
776  // outside contact zone
777  if (r > a) {
778  u_z = (1. - std::pow(poisson_ratio, 2.)) / young_modulus *
779  ((p_0) / (2. * a)) *
780  ((2. * std::pow(a, 2.) - std::pow(r, 2.)) * asin(a / r) +
781  std::pow(r, 2.) * (a / r) *
782  std::pow(1 - (std::pow(a, 2.) / std::pow(r, 2.)), 2.));
783  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
784  (3. * young_modulus) * ((std::pow(a, 2) / r)) * p_0;
785 
786  if (SPACE_DIM == 2)
787  v_u = {u_r, u_z};
788  else
789  v_u = {u_r, u_z, u_r};
790 
791  for (int i = 0; i < SPACE_DIM; ++i)
792  u(i) = v_u[i];
793 
794  return u;
795  }
796 
797  // In contact zone
798  u_z = ((1. - std::pow(poisson_ratio, 2.)) / young_modulus) *
799  ((M_PI * p_0) / 4. * a) * (2. * std::pow(a, 2.) - std::pow(r, 2.));
800  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
801  (3. * young_modulus) * ((std::pow(a, 2.) / r)) * p_0 *
802  (1 - std::pow(1 - (std::pow(r, 2.) / std::pow(a, 2.)), 1.5));
803 
804  if (SPACE_DIM == 2)
805  v_u = {u_r, u_z};
806  else
807  v_u = {u_r, u_z, u_r};
808 
809  for (int i = 0; i < SPACE_DIM; ++i)
810  u(i) = v_u[i];
811 
812  return u;
813  };
814 
816  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
817  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
818  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter) {
820  uXScatter = ux_scatter;
821  uYScatter = uy_scatter;
822  uZScatter = uz_scatter;
824  }
825 
826  MoFEMErrorCode getErrorNorm(int normType) {
827  const double *norm;
828  CHKERR VecGetArrayRead(normsVec, &norm);
829  double norm_val = std::sqrt(norm[normType]);
830  CHKERR VecRestoreArrayRead(normsVec, &norm);
831  return norm_val;
832  }
833 
834 private:
835  SmartPetscObj<DM> dM;
836  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
837  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
838  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
839 
840  boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
841 
842  boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
843  boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
844 
845  boost::shared_ptr<BoundaryEle> integrateTraction;
846  boost::shared_ptr<BoundaryEle> integrateArea;
847 
848  enum NORMS {
849  TRACTION_NORM_L2 = 0,
852  LAST_NORM
853  };
854  SmartPetscObj<Vec> normsVec;
855 
858 
859  double lastTime;
860  double deltaTime;
861  int sTEP;
862 
863  boost::shared_ptr<GenericElementInterface> mfrontInterface;
864 };
865 
866 } // namespace ContactOps
NOSPACE
@ NOSPACE
Definition: definitions.h:83
ContactOps::Monitor::postProcBdyFe
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
Definition: PostProcContact.hpp:843
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:854
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:125
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:863
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:861
ContactOps::PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: PostProcContact.hpp:24
ContactOps::Monitor::postProcMesh
boost::shared_ptr< moab::Core > postProcMesh
Definition: PostProcContact.hpp:840
ContactOps::Monitor::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: PostProcContact.hpp:836
AddHOOps
Definition: EshelbianOperators.hpp:599
scale
double scale
Definition: plastic.cpp:123
ContactOps::Monitor::MAG_TRACTION_NORM_L2
@ MAG_TRACTION_NORM_L2
Definition: PostProcContact.hpp:850
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:141
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:856
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
ContactOps::Monitor::moabVertex
moab::Interface & moabVertex
Definition: PostProcContact.hpp:857
ContactOps::Monitor::deltaTime
double deltaTime
Definition: PostProcContact.hpp:860
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:826
ContactOps::PostProcEleByDim< 2 >::PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
Definition: PostProcContact.hpp:14
ContactOps::Monitor::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: PostProcContact.hpp:838
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
ContactOps::Monitor::lastTime
double lastTime
Definition: PostProcContact.hpp:859
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:835
a
constexpr double a
Definition: approx_sphere.cpp:30
BoundaryEleOp
R
@ R
Definition: free_surface.cpp:396
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
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:126
ContactOps::Monitor::postProcDomainFe
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
Definition: PostProcContact.hpp:842
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:851
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:187
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:815
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:848
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:54
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:845
ContactOps::Monitor::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: PostProcContact.hpp:837
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
atom_test
int atom_test
Atom test.
Definition: plastic.cpp:121
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:846
F
@ F
Definition: free_surface.cpp:396
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709