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