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  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
229  CHK_THROW_MESSAGE(
230  (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
231  integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
232  "Apply transform");
233  // We have to integrate on curved face geometry, thus integration weight
234  // have to adjusted.
235  integrate_traction->getOpPtrVector().push_back(
236  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
237  integrate_traction->getRuleHook = [](int, int, int approx_order) {
238  return 2 * approx_order + geom_order - 1;
239  };
240 
241  CHK_THROW_MESSAGE(
242  (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
243  integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
244  "push operators to calculate traction");
245 
246  return integrate_traction;
247  };
248 
249  postProcDomainFe = get_post_proc_domain_fe();
250  if constexpr (SPACE_DIM == 2)
251  postProcBdyFe = get_post_proc_bdy_fe();
252  integrateTraction = get_integrate_traction();
253 
254  normsVec = createVectorMPI(
255  m_field_ptr->get_comm(),
256  (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
257  }
258 
259  MoFEMErrorCode preProcess() { return 0; }
260  MoFEMErrorCode operator()() { return 0; }
261 
262  MoFEMErrorCode postProcess() {
263  MoFEMFunctionBegin;
264  MoFEM::Interface *m_field_ptr;
265  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
266 
267  auto post_proc = [&]() {
268  MoFEMFunctionBegin;
269 
270  if (!mfrontInterface) {
271  auto post_proc_begin =
272  boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
273  postProcMesh);
274  auto post_proc_end =
275  boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
276  postProcMesh);
277 
278  CHKERR DMoFEMPreProcessFiniteElements(dM,
279  post_proc_begin->getFEMethod());
280  if (!postProcBdyFe) {
281  postProcDomainFe->copyTs(*this); // this here is a Monitor
282  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
283  } else {
284  postProcDomainFe->copyTs(*this); // this here is a Monitor
285  postProcBdyFe->copyTs(*this);
286  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
287  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
288  }
289  CHKERR DMoFEMPostProcessFiniteElements(dM,
290  post_proc_end->getFEMethod());
291 
292  CHKERR post_proc_end->writeFile(
293  "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
294  } else {
295  CHKERR mfrontInterface->postProcessElement(
296  ts_step, dM,
297  m_field_ptr->getInterface<Simple>()->getDomainFEName());
298  }
299 
300  MoFEMFunctionReturn(0);
301  };
302 
303  auto calculate_traction = [&] {
304  MoFEMFunctionBegin;
305  CHKERR VecZeroEntries(CommonData::totalTraction);
306  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
307  CHKERR VecAssemblyBegin(CommonData::totalTraction);
308  CHKERR VecAssemblyEnd(CommonData::totalTraction);
309  MoFEMFunctionReturn(0);
310  };
311 
312  auto calculate_reactions = [&]() {
313  MoFEMFunctionBegin;
314 
315  auto res = createDMVector(dM);
316 
317  auto assemble_domain = [&]() {
318  MoFEMFunctionBegin;
319  auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
320  auto &pip = fe_rhs->getOpPtrVector();
321  fe_rhs->f = res;
322 
323  auto integration_rule = [](int, int, int approx_order) {
324  return 2 * approx_order + geom_order - 1;
325  };
326  fe_rhs->getRuleHook = integration_rule;
327 
328  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
329  "GEOMETRY");
330  CHKERR
331  ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
332  pip, "SIGMA", "U", is_axisymmetric);
333 
334  if (!mfrontInterface) {
335  CHKERR
336  HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
337  *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
338  } else {
339  CHKERR mfrontInterface->opFactoryDomainRhs(pip);
340  }
341  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
342 
343  CHKERR VecAssemblyBegin(res);
344  CHKERR VecAssemblyEnd(res);
345  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
346  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
347 
348  MoFEMFunctionReturn(0);
349  };
350 
351  auto assemble_boundary = [&]() {
352  MoFEMFunctionBegin;
353  auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
354  auto &pip = fe_rhs->getOpPtrVector();
355  fe_rhs->f = res;
356 
357  auto integration_rule = [](int, int, int approx_order) {
358  return 2 * approx_order + geom_order - 1;
359  };
360  fe_rhs->getRuleHook = integration_rule;
361 
362  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
363  "GEOMETRY");
364  // We have to integrate on curved face geometry, thus integration weight
365  // have to adjusted.
366  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
367 
368  auto u_disp = boost::make_shared<MatrixDouble>();
369  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
370  pip.push_back(
371  new OpSpringRhs("U", u_disp, [this](double, double, double) {
372  return spring_stiffness;
373  }));
374 
375  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
376 
377  MoFEMFunctionReturn(0);
378  };
379 
380  CHKERR assemble_domain();
381  CHKERR assemble_boundary();
382 
383  auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
384  auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
385  m_field_ptr]() {
386  MoFEMFunctionBegin;
387  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
388  *m_field_ptr, fe_post_proc_ptr, res)();
389  MoFEMFunctionReturn(0);
390  };
391  fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
392  CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
393 
394  MoFEMFunctionReturn(0);
395  };
396 
397  auto print_max_min = [&](auto &tuple, const std::string msg) {
398  MoFEMFunctionBegin;
399  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
400  INSERT_VALUES, SCATTER_FORWARD);
401  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
402  INSERT_VALUES, SCATTER_FORWARD);
403  double max, min;
404  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
405  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
406  MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %6.4e min %6.4e max %6.4e",
407  msg.c_str(), ts_t, min, max);
408  MoFEMFunctionReturn(0);
409  };
410 
411  auto print_traction = [&](const std::string msg) {
412  MoFEMFunctionBegin;
413  MoFEM::Interface *m_field_ptr;
414  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
415  if (!m_field_ptr->get_comm_rank()) {
416  const double *t_ptr;
417  CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
418  MOFEM_LOG_C("CONTACT", Sev::inform,
419  "%s time %6.4e %6.16e %6.16e %6.16e", msg.c_str(), ts_t,
420  t_ptr[0], t_ptr[1], t_ptr[2]);
421  CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
422  }
423  MoFEMFunctionReturn(0);
424  };
425 
426  if (mfrontInterface) {
427  CHKERR mfrontInterface->updateElementVariables(
428  dM, m_field_ptr->getInterface<Simple>()->getDomainFEName());
429  }
430 
431  auto calculate_error = [&](MoFEM::ScalarFun &fun) {
432  MoFEMFunctionBegin;
433  struct OpCalcTractions : public BoundaryEleOp {
434  OpCalcTractions(boost::shared_ptr<MatrixDouble> m_ptr,
435  boost::shared_ptr<VectorDouble> p_ptr,
436  boost::shared_ptr<VectorDouble> mag_ptr,
437  boost::shared_ptr<VectorDouble> traction_y_ptr,
438  boost::shared_ptr<MatrixDouble> t_ptr,
439  boost::shared_ptr<MatrixDouble> grad_sdf_ptr)
440  : BoundaryEleOp(NOSPACE, OPSPACE), mPtr(m_ptr), pPtr(p_ptr),
441  magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
442  gradSDFPtr(grad_sdf_ptr) {}
443  MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
444  MoFEMFunctionBegin;
445  FTensor::Index<'i', SPACE_DIM> i;
446  mPtr->resize(SPACE_DIM, pPtr->size());
447  mPtr->clear();
448  magPtr->resize(pPtr->size());
449  magPtr->clear();
450  tyPtr->resize(pPtr->size());
451  tyPtr->clear();
452 
453  auto t_traction = getFTensor1FromMat<SPACE_DIM>(*mPtr);
454  auto t_contact_traction = getFTensor1FromMat<SPACE_DIM>(*tPtr);
455  auto t_p = getFTensor0FromVec(*pPtr);
456  int nb_gauss_pts = pPtr->size();
457  auto t_normal = getFTensor1FromMat<SPACE_DIM>(*gradSDFPtr);
458  auto t_normal_at_gauss = getFTensor1NormalsAtGaussPts();
459  auto t_mag = getFTensor0FromVec(*magPtr);
460  auto t_ty = getFTensor0FromVec(*tyPtr);
461 
462  for (int gg = 0; gg != nb_gauss_pts; gg++) {
463  FTensor::Tensor1<double, SPACE_DIM> normal;
464  t_traction(i) = t_p * (-(t_normal(i) / t_normal.l2()));
465  t_mag = t_contact_traction.l2();
466  t_ty = t_contact_traction(1);
467 
468  ++t_normal;
469  ++t_traction;
470  ++t_p;
471  ++t_mag;
472  ++t_contact_traction;
473  ++t_ty;
474  ++t_normal_at_gauss;
475  }
476  MoFEMFunctionReturn(0);
477  }
478 
479  private:
480  boost::shared_ptr<MatrixDouble> mPtr;
481  boost::shared_ptr<VectorDouble> pPtr;
482  boost::shared_ptr<VectorDouble> magPtr;
483  boost::shared_ptr<MatrixDouble> tPtr;
484  boost::shared_ptr<VectorDouble> tyPtr;
485  boost::shared_ptr<MatrixDouble> gradSDFPtr;
486  };
487 
488  auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
489  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
490  auto simple = m_field_ptr->getInterface<Simple>();
491  Range contact_range;
492  for (auto m :
493  m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
494  std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
495  auto meshset = m->getMeshset();
496  Range contact_meshset_range;
497  CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
498  meshset, SPACE_DIM - 1, contact_meshset_range, true);
499 
500  CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
501  contact_meshset_range);
502  contact_range.merge(contact_meshset_range);
503  }
504  // std::cout << "contact_range.size() " << contact_range.size() <<
505  // std::endl;
506 
507  CHK_THROW_MESSAGE(
508  (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
509  post_proc_norm_fe->getOpPtrVector(), {HDIV}, "GEOMETRY")),
510  "Apply transform");
511  // We have to integrate on curved face geometry, thus integration weight
512  // have to adjusted.
513  post_proc_norm_fe->getOpPtrVector().push_back(
514  new OpSetHOWeightsOnSubDim<SPACE_DIM>());
515  post_proc_norm_fe->getRuleHook = [](int, int, int approx_order) {
516  return 2 * approx_order + geom_order - 1;
517  };
518 
519  post_proc_norm_fe->getOpPtrVector().push_back(
520  new OpCalculateVectorFieldValues<SPACE_DIM>(
521  "U", common_data_ptr->contactDispPtr()));
522  post_proc_norm_fe->getOpPtrVector().push_back(
523  new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
524  "SIGMA", common_data_ptr->contactTractionPtr()));
525  using C = ContactIntegrators<BoundaryEleOp>;
526  post_proc_norm_fe->getOpPtrVector().push_back(
527  new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
528  common_data_ptr));
529 
530  // post_proc_norm_fe->getOpPtrVector().push_back(
531  // new OpScale(common_data_ptr->contactTractionPtr(), scale));
532 
533  auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
534  auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
535  auto mag_traction_ptr = boost::make_shared<VectorDouble>();
536  auto traction_y_ptr = boost::make_shared<VectorDouble>();
537  auto contact_range_ptr = boost::make_shared<Range>(contact_range);
538 
539  post_proc_norm_fe->getOpPtrVector().push_back(
540  new OpGetTensor0fromFunc(analytical_pressure_ptr, fun));
541 
542  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcTractions(
543  analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
544  traction_y_ptr, common_data_ptr->contactTractionPtr(),
545  common_data_ptr->gradSdfPtr()));
546 
547  post_proc_norm_fe->getOpPtrVector().push_back(
548  new OpCalcNormL2Tensor1<SPACE_DIM>(
549  common_data_ptr->contactTractionPtr(), normsVec, TRACTION_NORM_L2,
550  analytical_traction_ptr, contact_range_ptr));
551 
552  // calculate magnitude of traction
553 
554  post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcNormL2Tensor0(
555  mag_traction_ptr, normsVec, MAG_TRACTION_NORM_L2,
556  analytical_pressure_ptr, contact_range_ptr));
557 
558  post_proc_norm_fe->getOpPtrVector().push_back(
559  new OpCalcNormL2Tensor0(traction_y_ptr, normsVec, TRACTION_Y_NORM_L2,
560  analytical_pressure_ptr, contact_range_ptr));
561 
562  CHKERR VecZeroEntries(normsVec);
563  post_proc_norm_fe->copyTs(*this); // set time as is in Monitor
564  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", post_proc_norm_fe);
565  CHKERR VecAssemblyBegin(normsVec);
566  CHKERR VecAssemblyEnd(normsVec);
567 
568  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
569  if (m_field_ptr->get_comm_rank() == 0) {
570  const double *norms;
571  CHKERR VecGetArrayRead(normsVec, &norms);
572  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
573  << "norm_traction: " << std::scientific
574  << std::sqrt(norms[TRACTION_NORM_L2]);
575  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
576  << "norm_mag_traction: " << std::scientific
577  << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
578  MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
579  << "norm_traction_y: " << std::scientific
580  << std::sqrt(norms[TRACTION_Y_NORM_L2]);
581  CHKERR VecRestoreArrayRead(normsVec, &norms);
582  }
583  MoFEMFunctionReturn(0);
584  };
585 
586  int se = 1;
587  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &se, PETSC_NULL);
588 
589  if (!(ts_step % se)) {
590  MOFEM_LOG("CONTACT", Sev::inform)
591  << "Write file at time " << ts_t << " write step " << sTEP;
592  CHKERR post_proc();
593  }
594  CHKERR calculate_traction();
595  CHKERR calculate_reactions();
596  if (atom_test == 1 && sTEP > 0)
597  CHKERR calculate_error(analyticalHertzPressurePlaneStress);
598  if (atom_test == 2 && sTEP > 0)
599  CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
600  if (atom_test == 5 && sTEP > 0)
601  CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
602  if (atom_test == 6 && sTEP > 0)
603  CHKERR calculate_error(analyticalWavy2DPressure);
604 
605  CHKERR print_max_min(uXScatter, "Ux");
606  CHKERR print_max_min(uYScatter, "Uy");
607  if (SPACE_DIM == 3)
608  CHKERR print_max_min(uZScatter, "Uz");
609  CHKERR print_traction("Contact force");
610 
611  ++sTEP;
612 
613  MoFEMFunctionReturn(0);
614  }
615 
616  //! [Analytical function]
617  MoFEM::ScalarFun analyticalWavy2DPressure = [](double x, double y, double z) {
618  // update to atom test values
619  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
620  // delta
621  double delta = 0.0002;
622  // lambda
623  double lambda = 2;
624  // pressure star
625  double p_star = M_PI * E_star * delta / lambda;
626 
627  // Pressure = p_bar + p_star * cos(2 * pi * x / lambda)
628  return p_star + p_star * std::cos(2. * M_PI * x / lambda);
629  };
630 
631  MoFEM::ScalarFun analyticalHertzPressureAxisymmetric = [](double x, double y,
632  double z) {
633  // update to atom test values
634  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
635  // Radius
636  double R = 100.;
637  // Indentation
638  double d = 0.01;
639  // Force
640  double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
641  // Contact area radius
642  double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
643  // Maximum pressure
644  double p_max = (3. * F) / (2. * M_PI * a * a);
645 
646  double r = std::sqrt((x * x) + (y * y));
647 
648  if (r > a) {
649  return 0.;
650  }
651  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
652  return p_max * std::sqrt(1 - ((r * r) / (a * a)));
653  };
654 
655  MoFEM::ScalarFun analyticalHertzPressurePlaneStrain = [](double x, double y,
656  double z) {
657  // update to atom test values
658  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
659  // Radius
660  double R = 100.;
661  // Indentation
662  double d = 0.02745732273553991;
663  // Contact area radius
664  double a = 1;
665  // current radius
666  double r = std::sqrt((x * x) + (y * y));
667 
668  if (r > a) {
669  return 0.;
670  }
671  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
672  return E_star / (2. * R) * std::sqrt(a * a - r * r);
673  };
674  MoFEM::ScalarFun analyticalHertzPressurePlaneStress = [](double x, double y,
675  double z) {
676  // update to atom test values
677  double E_star = young_modulus;
678  // Radius
679  double R = 100.;
680  // Indentation
681  double d = 0.02745732273553991;
682  // Contact area radius
683  double a = 1;
684  // current radius
685  double r = std::sqrt((x * x) + (y * y));
686 
687  if (r > a) {
688  return 0.;
689  }
690  // Pressure = p_max * sqrt(1 - (r^2 / a^2))
691  return E_star / (2. * R) * std::sqrt(a * a - r * r);
692  };
693 
694  // ***DISPLACMENT NOT TESTED***
695  MoFEM::VectorFun<SPACE_DIM> analyticalHertzDisplacement3D = [](double x,
696  double y,
697  double z) {
698  // update to atom test values
699  double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
700  // Radius
701  double R = 100.;
702  // Contact area radius
703  double a = 1;
704  // max pressure
705  double p_0 = (2. * E_star * a) / (M_PI * R);
706  // current radius
707  double r = std::sqrt((x * x) + (y * y));
708  FTensor::Tensor1<double, SPACE_DIM> u;
709  std::vector<double> v_u;
710 
711  double u_z = 0.;
712  double u_r = 0.;
713  // outside contact zone
714  if (r > a) {
715  u_z = (1. - std::pow(poisson_ratio, 2.)) / young_modulus *
716  ((p_0) / (2. * a)) *
717  ((2. * std::pow(a, 2.) - std::pow(r, 2.)) * asin(a / r) +
718  std::pow(r, 2.) * (a / r) *
719  std::pow(1 - (std::pow(a, 2.) / std::pow(r, 2.)), 2.));
720  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
721  (3. * young_modulus) * ((std::pow(a, 2) / r)) * p_0;
722 
723  if (SPACE_DIM == 2)
724  v_u = {u_r, u_z};
725  else
726  v_u = {u_r, u_z, u_r};
727 
728  for (int i = 0; i < SPACE_DIM; ++i)
729  u(i) = v_u[i];
730 
731  return u;
732  }
733 
734  // In contact zone
735  u_z = ((1. - std::pow(poisson_ratio, 2.)) / young_modulus) *
736  ((M_PI * p_0) / 4. * a) * (2. * std::pow(a, 2.) - std::pow(r, 2.));
737  u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
738  (3. * young_modulus) * ((std::pow(a, 2.) / r)) * p_0 *
739  (1 - std::pow(1 - (std::pow(r, 2.) / std::pow(a, 2.)), 1.5));
740 
741  if (SPACE_DIM == 2)
742  v_u = {u_r, u_z};
743  else
744  v_u = {u_r, u_z, u_r};
745 
746  for (int i = 0; i < SPACE_DIM; ++i)
747  u(i) = v_u[i];
748 
749  return u;
750  };
751 
752  MoFEMErrorCode setScatterVectors(
753  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
754  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
755  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter) {
756  MoFEMFunctionBegin;
757  uXScatter = ux_scatter;
758  uYScatter = uy_scatter;
759  uZScatter = uz_scatter;
760  MoFEMFunctionReturn(0);
761  }
762 
763  MoFEMErrorCode getErrorNorm(int normType) {
764  const double *norm;
765  CHKERR VecGetArrayRead(normsVec, &norm);
766  double norm_val = std::sqrt(norm[normType]);
767  CHKERR VecRestoreArrayRead(normsVec, &norm);
768  return norm_val;
769  }
770 
771 private:
772  SmartPetscObj<DM> dM;
773  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
774  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
775  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
776 
777  boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
778 
779  boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
780  boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
781 
782  boost::shared_ptr<BoundaryEle> integrateTraction;
783 
784  enum NORMS {
785  TRACTION_NORM_L2 = 0,
786  MAG_TRACTION_NORM_L2,
787  TRACTION_Y_NORM_L2,
788  LAST_NORM
789  };
790  SmartPetscObj<Vec> normsVec;
791 
792  moab::Core mbVertexPostproc;
793  moab::Interface &moabVertex;
794 
795  double lastTime;
796  double deltaTime;
797  int sTEP;
798 
799  boost::shared_ptr<GenericElementInterface> mfrontInterface;
800 };
801 
802 } // namespace ContactOps
NOSPACE
@ NOSPACE
Definition: definitions.h:83
ContactOps::Monitor::postProcBdyFe
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
Definition: PostProcContact.hpp:780
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:790
ContactOps::PostProcEleByDim
Definition: PostProcContact.hpp:10
spring_stiffness
double spring_stiffness
Definition: contact.cpp:84
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: EshelbianContact.hpp:10
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:90
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:799
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:797
ContactOps::PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: PostProcContact.hpp:24
ContactOps::Monitor::postProcMesh
boost::shared_ptr< moab::Core > postProcMesh
Definition: PostProcContact.hpp:777
ContactOps::Monitor::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: PostProcContact.hpp:773
ContactOps::Monitor::MAG_TRACTION_NORM_L2
@ MAG_TRACTION_NORM_L2
Definition: PostProcContact.hpp:786
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:792
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
atom_test
int atom_test
Definition: contact.cpp:94
ContactOps::Monitor::moabVertex
moab::Interface & moabVertex
Definition: PostProcContact.hpp:793
ContactOps::Monitor::deltaTime
double deltaTime
Definition: PostProcContact.hpp:796
FEMethod
ContactOps::Monitor
Definition: PostProcContact.hpp:28
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
ContactOps::Monitor::getErrorNorm
MoFEMErrorCode getErrorNorm(int normType)
Definition: PostProcContact.hpp:763
ContactOps::PostProcEleByDim< 2 >::PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
Definition: PostProcContact.hpp:14
ContactOps::Monitor::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: PostProcContact.hpp:775
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
ContactOps::Monitor::lastTime
double lastTime
Definition: PostProcContact.hpp:795
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:772
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:73
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:523
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
ContactOps::SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: PostProcContact.hpp:25
ContactOps::Monitor::postProcDomainFe
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
Definition: PostProcContact.hpp:779
ContactOps::Monitor::postProcess
MoFEMErrorCode postProcess()
Definition: PostProcContact.hpp:262
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:787
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', SPACE_DIM >
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:752
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:784
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:1140
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:782
ContactOps::Monitor::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: PostProcContact.hpp:774
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
ContactOps::Monitor::operator()
MoFEMErrorCode operator()()
Definition: PostProcContact.hpp:260
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
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