31 boost::shared_ptr<GenericElementInterface> mfront_interface =
nullptr,
39 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
42 OpScale(boost::shared_ptr<MatrixDouble> m_ptr,
double s)
44 mPtr(m_ptr),
scale(s) {}
51 boost::shared_ptr<MatrixDouble> mPtr;
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);
69 auto push_bdy_ops = [&](
auto &pip) {
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>(
80 return common_data_ptr;
83 auto get_domain_pip = [&](
auto &pip)
84 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
86 auto op_loop_side =
new OpLoopSide<SideEle>(
87 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
89 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
90 pip.push_back(op_loop_side);
91 return op_loop_side->getOpPtrVector();
97 auto get_post_proc_domain_fe = [&]() {
99 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
100 auto &pip = post_proc_fe->getOpPtrVector();
102 auto [hencky_common_data_ptr, contact_stress_ptr] =
103 push_domain_ops(get_domain_pip(pip));
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>();
109 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
117 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
122 {
"U", u_ptr}, {
"GEOMETRY", X_ptr}
127 {
"SIGMA", contact_stress_ptr},
129 {
"G", hencky_common_data_ptr->matGradPtr},
131 {
"PK1", hencky_common_data_ptr->getMatFirstPiolaStress()}
143 pip, {
HDIV},
"GEOMETRY")),
145 auto common_data_ptr = push_bdy_ops(pip);
151 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
153 {{
"SDF", common_data_ptr->sdfPtr()},
154 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
158 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
159 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
165 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
175 auto get_post_proc_bdy_fe = [&]() {
177 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
178 auto &pip = post_proc_fe->getOpPtrVector();
181 pip, {
HDIV},
"GEOMETRY")),
183 auto common_data_ptr = push_bdy_ops(pip);
186 auto op_loop_side =
new OpLoopSide<SideEle>(
187 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
189 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
190 pip.push_back(op_loop_side);
192 auto [hencky_common_data_ptr, contact_stress_ptr] =
193 push_domain_ops(op_loop_side->getOpPtrVector());
195 auto X_ptr = boost::make_shared<MatrixDouble>();
197 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
203 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
205 {{
"SDF", common_data_ptr->sdfPtr()},
206 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
208 {{
"U", common_data_ptr->contactDispPtr()},
210 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
211 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
217 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
226 auto get_integrate_traction = [&]() {
227 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
230 integrate_traction->getOpPtrVector(), {HDIV},
"GEOMETRY")),
234 integrate_traction->getOpPtrVector().push_back(
235 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
241 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
243 "push operators to calculate traction");
245 return integrate_traction;
248 auto get_integrate_area = [&]() {
249 auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
253 integrate_area->getOpPtrVector(), {HDIV},
"GEOMETRY")),
256 integrate_area->getOpPtrVector().push_back(
257 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
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);
270 CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
271 contact_meshset_range);
272 contact_range.merge(contact_meshset_range);
275 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
277 auto op_loop_side =
new OpLoopSide<SideEle>(
278 *m_field_ptr, m_field_ptr->getInterface<Simple>()->getDomainFEName(),
281 op_loop_side->getOpPtrVector(), {H1},
"GEOMETRY");
284 (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
285 integrate_area->getOpPtrVector(), op_loop_side,
"SIGMA",
"U",
287 "push operators to calculate area");
289 return integrate_area;
292 postProcDomainFe = get_post_proc_domain_fe();
294 postProcBdyFe = get_post_proc_bdy_fe();
296 integrateTraction = get_integrate_traction();
297 integrateArea = get_integrate_area();
300 m_field_ptr->get_comm(),
301 (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
312 auto post_proc = [&]() {
315 if (!mfrontInterface) {
316 auto post_proc_begin =
317 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
320 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
324 post_proc_begin->getFEMethod());
325 if (!postProcBdyFe) {
326 postProcDomainFe->copyTs(*
this);
329 postProcDomainFe->copyTs(*
this);
330 postProcBdyFe->copyTs(*
this);
335 post_proc_end->getFEMethod());
337 CHKERR post_proc_end->writeFile(
338 "out_contact_" + boost::lexical_cast<std::string>(sTEP) +
".h5m");
340 CHKERR mfrontInterface->postProcessElement(
342 m_field_ptr->
getInterface<Simple>()->getDomainFEName());
348 auto calculate_force = [&] {
350 CHKERR VecZeroEntries(CommonData::totalTraction);
352 CHKERR VecAssemblyBegin(CommonData::totalTraction);
353 CHKERR VecAssemblyEnd(CommonData::totalTraction);
357 auto calculate_area = [&] {
359 integrateArea->copyTs(*
this);
361 CHKERR VecAssemblyBegin(CommonData::totalTraction);
362 CHKERR VecAssemblyEnd(CommonData::totalTraction);
366 auto calculate_reactions = [&]() {
371 auto assemble_domain = [&]() {
373 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
374 auto &pip = fe_rhs->getOpPtrVector();
385 ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
388 if (!mfrontInterface) {
390 HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
391 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
393 CHKERR mfrontInterface->opFactoryDomainRhs(pip);
397 CHKERR VecAssemblyBegin(res);
398 CHKERR VecAssemblyEnd(res);
399 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
400 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
405 auto assemble_boundary = [&]() {
407 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
408 auto &pip = fe_rhs->getOpPtrVector();
420 pip.push_back(
new OpSetHOWeightsOnSubDim<SPACE_DIM>());
422 auto u_disp = boost::make_shared<MatrixDouble>();
423 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_disp));
425 new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
435 CHKERR assemble_boundary();
437 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
438 auto get_post_proc_hook_rhs = [
this, fe_post_proc_ptr, res,
441 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
442 *m_field_ptr, fe_post_proc_ptr, res)();
445 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
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);
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);
465 auto print_force_and_area = [&]() {
471 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
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]);
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);
483 if (mfrontInterface) {
484 CHKERR mfrontInterface->updateElementVariables(
485 dM, m_field_ptr->
getInterface<Simple>()->getDomainFEName());
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)
498 magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
499 gradSDFPtr(grad_sdf_ptr) {}
505 magPtr->resize(pPtr->size());
507 tyPtr->resize(pPtr->size());
510 auto t_traction = getFTensor1FromMat<SPACE_DIM>(*mPtr);
511 auto t_contact_traction = getFTensor1FromMat<SPACE_DIM>(*tPtr);
513 int nb_gauss_pts = pPtr->size();
514 auto t_normal = getFTensor1FromMat<SPACE_DIM>(*gradSDFPtr);
515 auto t_normal_at_gauss = getFTensor1NormalsAtGaussPts();
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);
529 ++t_contact_traction;
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;
545 auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
546 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
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;
554 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
557 contact_meshset_range);
558 contact_range.merge(contact_meshset_range);
563 post_proc_norm_fe->getOpPtrVector(), {HDIV},
"GEOMETRY")),
567 post_proc_norm_fe->getOpPtrVector().push_back(
568 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
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>(
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);
590 post_proc_norm_fe->getOpPtrVector().push_back(
591 new OpGetTensor0fromFunc(analytical_pressure_ptr,
fun));
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()));
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));
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));
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));
613 CHKERR VecZeroEntries(normsVec);
614 post_proc_norm_fe->copyTs(*
this);
616 CHKERR VecAssemblyBegin(normsVec);
617 CHKERR VecAssemblyEnd(normsVec);
622 CHKERR VecGetArrayRead(normsVec, &norms);
624 <<
"norm_traction: " << std::scientific
625 << std::sqrt(norms[TRACTION_NORM_L2]);
627 <<
"norm_mag_traction: " << std::scientific
628 << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
630 <<
"norm_traction_y: " << std::scientific
631 << std::sqrt(norms[TRACTION_Y_NORM_L2]);
632 CHKERR VecRestoreArrayRead(normsVec, &norms);
640 if (!(ts_step % se)) {
642 <<
"Write file at time " << ts_t <<
" write step " << sTEP;
648 CHKERR calculate_reactions();
653 CHKERR calculate_error(analyticalHertzPressurePlaneStress);
656 CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
659 CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
662 CHKERR calculate_error(analyticalWavy2DPressure);
669 CHKERR print_max_min(uXScatter,
"Ux");
670 CHKERR print_max_min(uYScatter,
"Uy");
672 CHKERR print_max_min(uZScatter,
"Uz");
673 CHKERR print_force_and_area();
684 double delta = 0.0002;
691 return p_star + p_star * std::cos(2. * M_PI * x /
lambda);
703 double F = (4. / 3.) * E_star * std::sqrt(
R) * std::pow(
d, 1.5);
705 double a = std::pow((3. *
F *
R) / (4. * E_star), 1. / 3.);
707 double p_max = (3. *
F) / (2. * M_PI *
a *
a);
709 double r = std::sqrt((x * x) + (y * y));
715 return p_max * std::sqrt(1 - ((
r *
r) / (
a *
a)));
725 double d = 0.02745732273553991;
729 double r = std::sqrt((x * x) + (y * y));
735 return E_star / (2. *
R) * std::sqrt(
a *
a -
r *
r);
744 double d = 0.02745732273553991;
748 double r = std::sqrt((x * x) + (y * y));
754 return E_star / (2. *
R) * std::sqrt(
a *
a -
r *
r);
768 double p_0 = (2. * E_star *
a) / (M_PI *
R);
770 double r = std::sqrt((x * x) + (y * y));
772 std::vector<double> v_u;
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.));
789 v_u = {u_r, u_z, u_r};
799 ((M_PI * p_0) / 4. *
a) * (2. * std::pow(
a, 2.) - std::pow(
r, 2.));
802 (1 - std::pow(1 - (std::pow(
r, 2.) / std::pow(
a, 2.)), 1.5));
807 v_u = {u_r, u_z, u_r};
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;
828 CHKERR VecGetArrayRead(normsVec, &norm);
829 double norm_val = std::sqrt(norm[normType]);
830 CHKERR VecRestoreArrayRead(normsVec, &norm);
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;
840 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
849 TRACTION_NORM_L2 = 0,