31 boost::shared_ptr<GenericElementInterface> mfront_interface =
nullptr,
37 CHKERR DMoFEMGetInterfacePtr(
dM, &m_field_ptr);
39 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
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 &) {
51 boost::shared_ptr<MatrixDouble> mPtr;
55 auto push_domain_ops = [&](
auto &pip) {
57 pip, {
H1,
HDIV},
"GEOMETRY")),
58 "Apply base transform");
61 using CommonDataVariant =
62 std::variant<boost::shared_ptr<HookeOps::CommonData>,
63 boost::shared_ptr<HenckyOps::CommonData>>;
64 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
65 pip.push_back(
new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
66 "SIGMA", contact_stress_ptr));
69 auto hooke_common_data_ptr =
71 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
73 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
74 return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
77 auto hencky_common_data_ptr =
79 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
80 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
81 return std::make_tuple(CommonDataVariant(hencky_common_data_ptr),
85 auto push_bdy_ops = [&](
auto &pip) {
87 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
88 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
89 "U", common_data_ptr->contactDispPtr()));
90 pip.push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
91 "SIGMA", common_data_ptr->contactTractionPtr()));
92 pip.push_back(
new OpScale(common_data_ptr->contactTractionPtr(),
scale));
94 pip.push_back(
new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
96 return common_data_ptr;
99 auto get_domain_pip = [&](
auto &pip)
100 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
102 auto op_loop_side =
new OpLoopSide<SideEle>(
103 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
105 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
106 pip.push_back(op_loop_side);
107 return op_loop_side->getOpPtrVector();
115 auto get_post_proc_domain_fe = [&]() {
117 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
118 auto &pip = post_proc_fe->getOpPtrVector();
120 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
121 push_domain_ops(get_domain_pip(pip));
123 auto u_ptr = boost::make_shared<MatrixDouble>();
124 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_ptr));
125 auto X_ptr = boost::make_shared<MatrixDouble>();
127 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
136 [&](
auto &common_data_ptr) {
137 if constexpr (std::is_same_v<
138 std::decay_t<
decltype(common_data_ptr)>,
139 boost::shared_ptr<HookeOps::CommonData>>) {
141 post_proc_fe->getPostProcMesh(),
142 post_proc_fe->getMapGaussPts(), {},
143 {{
"U", u_ptr}, {
"GEOMETRY", X_ptr}},
144 {{
"SIGMA", contact_stress_ptr},
145 {
"G", common_data_ptr->matGradPtr}},
146 {{
"STRESS", common_data_ptr->getMatCauchyStress()},
147 {
"STRAIN", common_data_ptr->getMatStrain()}}));
150 else if constexpr (std::is_same_v<
151 std::decay_t<
decltype(common_data_ptr)>,
156 post_proc_fe->getPostProcMesh(),
157 post_proc_fe->getMapGaussPts(), {},
158 {{
"U", u_ptr}, {
"GEOMETRY", X_ptr}},
159 {{
"SIGMA", contact_stress_ptr},
160 {
"G", common_data_ptr->matGradPtr},
161 {
"PK1", common_data_ptr->getMatFirstPiolaStress()}
164 {{
"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
167 hencky_or_hooke_common_data_ptr);
172 pip, {
HDIV},
"GEOMETRY")),
174 auto common_data_ptr = push_bdy_ops(pip);
180 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
182 {{
"SDF", common_data_ptr->sdfPtr()},
183 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
187 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
188 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
194 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
204 auto get_post_proc_bdy_fe = [&]() {
206 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
207 auto &pip = post_proc_fe->getOpPtrVector();
210 pip, {
HDIV},
"GEOMETRY")),
212 auto common_data_ptr = push_bdy_ops(pip);
215 auto op_loop_side =
new OpLoopSide<SideEle>(
216 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
218 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
219 pip.push_back(op_loop_side);
221 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
222 push_domain_ops(op_loop_side->getOpPtrVector());
224 auto X_ptr = boost::make_shared<MatrixDouble>();
226 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
232 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
234 {{
"SDF", common_data_ptr->sdfPtr()},
235 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
237 {{
"U", common_data_ptr->contactDispPtr()},
239 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
240 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
246 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
255 auto get_integrate_traction = [&]() {
256 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
259 integrate_traction->getOpPtrVector(), {HDIV},
"GEOMETRY")),
262 integrate_traction->getOpPtrVector().push_back(
263 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
264 integrate_traction->getRuleHook = [](int, int,
int approx_order) {
269 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
271 "push operators to calculate traction");
273 return integrate_traction;
276 auto get_integrate_area = [&]() {
277 auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
281 integrate_area->getOpPtrVector(), {HDIV},
"GEOMETRY")),
285 integrate_area->getOpPtrVector().push_back(
286 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
287 integrate_area->getRuleHook = [](int, int,
int approx_order) {
292 m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
293 std::regex((boost::format(
"%s(.*)") %
"CONTACT").str()))) {
294 auto meshset =
m->getMeshset();
295 Range contact_meshset_range;
296 CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
297 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
299 CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
300 contact_meshset_range);
301 contact_range.merge(contact_meshset_range);
304 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
306 auto op_loop_side =
new OpLoopSide<SideEle>(
307 *m_field_ptr, m_field_ptr->getInterface<Simple>()->getDomainFEName(),
310 op_loop_side->getOpPtrVector(), {H1},
"GEOMETRY");
313 (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
314 integrate_area->getOpPtrVector(), op_loop_side,
"SIGMA",
"U",
316 "push operators to calculate area");
318 return integrate_area;
321 postProcDomainFe = get_post_proc_domain_fe();
323 postProcBdyFe = get_post_proc_bdy_fe();
325 integrateTraction = get_integrate_traction();
326 integrateArea = get_integrate_area();
329 m_field_ptr->get_comm(),
330 (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
339 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
341 auto post_proc = [&]() {
344 if (!mfrontInterface) {
345 auto post_proc_begin =
346 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
349 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
352 CHKERR DMoFEMPreProcessFiniteElements(dM,
353 post_proc_begin->getFEMethod());
354 if (!postProcBdyFe) {
355 postProcDomainFe->copyTs(*
this);
356 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcDomainFe);
358 postProcDomainFe->copyTs(*
this);
359 postProcBdyFe->copyTs(*
this);
360 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", postProcDomainFe);
361 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcBdyFe);
363 CHKERR DMoFEMPostProcessFiniteElements(dM,
364 post_proc_end->getFEMethod());
366 CHKERR post_proc_end->writeFile(
367 "out_contact_" + boost::lexical_cast<std::string>(sTEP) +
".h5m");
369 CHKERR mfrontInterface->postProcessElement(
371 m_field_ptr->
getInterface<Simple>()->getDomainFEName());
377 auto calculate_force = [&] {
379 CHKERR VecZeroEntries(CommonData::totalTraction);
380 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", integrateTraction);
381 CHKERR VecAssemblyBegin(CommonData::totalTraction);
382 CHKERR VecAssemblyEnd(CommonData::totalTraction);
386 auto calculate_area = [&] {
388 integrateArea->copyTs(*
this);
389 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", integrateArea);
390 CHKERR VecAssemblyBegin(CommonData::totalTraction);
391 CHKERR VecAssemblyEnd(CommonData::totalTraction);
395 auto calculate_reactions = [&]() {
398 auto res = createDMVector(dM);
400 auto assemble_domain = [&]() {
402 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
403 auto &pip = fe_rhs->getOpPtrVector();
417 if (!mfrontInterface) {
420 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
422 CHKERR mfrontInterface->opFactoryDomainRhs(pip);
424 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", fe_rhs);
426 CHKERR VecAssemblyBegin(res);
427 CHKERR VecAssemblyEnd(res);
428 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
429 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
434 auto assemble_boundary = [&]() {
436 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
437 auto &pip = fe_rhs->getOpPtrVector();
449 pip.push_back(
new OpSetHOWeightsOnSubDim<SPACE_DIM>());
451 auto u_disp = boost::make_shared<MatrixDouble>();
452 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_disp));
454 new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
458 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", fe_rhs);
464 CHKERR assemble_boundary();
466 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
467 auto get_post_proc_hook_rhs = [
this, fe_post_proc_ptr, res,
470 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
471 *m_field_ptr, fe_post_proc_ptr, res)();
474 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
475 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
480 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
482 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
483 INSERT_VALUES, SCATTER_FORWARD);
484 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
485 INSERT_VALUES, SCATTER_FORWARD);
487 CHKERR VecMax(std::get<0>(tuple), PETSC_NULLPTR, &max);
488 CHKERR VecMin(std::get<0>(tuple), PETSC_NULLPTR, &min);
489 MOFEM_LOG_C(
"CONTACT", Sev::inform,
"%s time %6.4e min %6.4e max %6.4e",
490 msg.c_str(), ts_t, min, max);
494 auto print_force_and_area = [&]() {
497 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
500 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
502 "Contact force: time %6.3e Fx: %6.6e Fy: %6.6e Fz: %6.6e",
503 ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
505 "Contact area: time %6.3e Active: %6.6e Potential: %6.6e",
506 ts_t, t_ptr[3], t_ptr[4]);
507 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
512 if (mfrontInterface) {
513 CHKERR mfrontInterface->updateElementVariables(
514 dM, m_field_ptr->
getInterface<Simple>()->getDomainFEName());
520 OpCalcTractions(boost::shared_ptr<MatrixDouble> m_ptr,
521 boost::shared_ptr<VectorDouble> p_ptr,
522 boost::shared_ptr<VectorDouble> mag_ptr,
523 boost::shared_ptr<VectorDouble> traction_y_ptr,
524 boost::shared_ptr<MatrixDouble> t_ptr,
525 boost::shared_ptr<MatrixDouble> grad_sdf_ptr)
527 magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
528 gradSDFPtr(grad_sdf_ptr) {}
529 MoFEMErrorCode doWork(
int, EntityType, EntitiesFieldData::EntData &) {
534 magPtr->resize(pPtr->size());
536 tyPtr->resize(pPtr->size());
539 auto t_traction = getFTensor1FromMat<SPACE_DIM>(*mPtr);
540 auto t_contact_traction = getFTensor1FromMat<SPACE_DIM>(*tPtr);
541 auto t_p = getFTensor0FromVec(*pPtr);
542 int nb_gauss_pts = pPtr->size();
543 auto t_normal = getFTensor1FromMat<SPACE_DIM>(*gradSDFPtr);
544 auto t_normal_at_gauss = getFTensor1NormalsAtGaussPts();
545 auto t_mag = getFTensor0FromVec(*magPtr);
546 auto t_ty = getFTensor0FromVec(*tyPtr);
548 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
550 t_traction(
i) = t_p * (-(t_normal(
i) / t_normal.l2()));
551 t_mag = t_contact_traction.
l2();
552 t_ty = t_contact_traction(1);
558 ++t_contact_traction;
566 boost::shared_ptr<MatrixDouble> mPtr;
567 boost::shared_ptr<VectorDouble> pPtr;
568 boost::shared_ptr<VectorDouble> magPtr;
569 boost::shared_ptr<MatrixDouble> tPtr;
570 boost::shared_ptr<VectorDouble> tyPtr;
571 boost::shared_ptr<MatrixDouble> gradSDFPtr;
574 auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
575 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
578 m_field_ptr->
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
579 std::regex((boost::format(
"%s(.*)") %
"CONTACT").str()))) {
580 auto meshset =
m->getMeshset();
581 Range contact_meshset_range;
583 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
586 contact_meshset_range);
587 contact_range.merge(contact_meshset_range);
592 post_proc_norm_fe->getOpPtrVector(), {HDIV},
"GEOMETRY")),
596 post_proc_norm_fe->getOpPtrVector().push_back(
597 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
598 post_proc_norm_fe->getRuleHook = [](int, int,
int approx_order) {
602 post_proc_norm_fe->getOpPtrVector().push_back(
603 new OpCalculateVectorFieldValues<SPACE_DIM>(
604 "U", common_data_ptr->contactDispPtr()));
605 post_proc_norm_fe->getOpPtrVector().push_back(
606 new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
607 "SIGMA", common_data_ptr->contactTractionPtr()));
609 post_proc_norm_fe->getOpPtrVector().push_back(
610 new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
613 auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
614 auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
615 auto mag_traction_ptr = boost::make_shared<VectorDouble>();
616 auto traction_y_ptr = boost::make_shared<VectorDouble>();
617 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
619 post_proc_norm_fe->getOpPtrVector().push_back(
620 new OpGetTensor0fromFunc(analytical_pressure_ptr,
fun));
622 post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcTractions(
623 analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
624 traction_y_ptr, common_data_ptr->contactTractionPtr(),
625 common_data_ptr->gradSdfPtr()));
627 post_proc_norm_fe->getOpPtrVector().push_back(
628 new OpCalcNormL2Tensor1<SPACE_DIM>(
629 common_data_ptr->contactTractionPtr(), normsVec, TRACTION_NORM_L2,
630 analytical_traction_ptr, contact_range_ptr));
634 post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor0(
635 mag_traction_ptr, normsVec, MAG_TRACTION_NORM_L2,
636 analytical_pressure_ptr, contact_range_ptr));
638 post_proc_norm_fe->getOpPtrVector().push_back(
639 new OpCalcNormL2Tensor0(traction_y_ptr, normsVec, TRACTION_Y_NORM_L2,
640 analytical_pressure_ptr, contact_range_ptr));
642 CHKERR VecZeroEntries(normsVec);
643 post_proc_norm_fe->copyTs(*
this);
644 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", post_proc_norm_fe);
645 CHKERR VecAssemblyBegin(normsVec);
646 CHKERR VecAssemblyEnd(normsVec);
651 CHKERR VecGetArrayRead(normsVec, &norms);
653 <<
"norm_traction: " << std::scientific
654 << std::sqrt(norms[TRACTION_NORM_L2]);
656 <<
"norm_mag_traction: " << std::scientific
657 << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
659 <<
"norm_traction_y: " << std::scientific
660 << std::sqrt(norms[TRACTION_Y_NORM_L2]);
661 CHKERR VecRestoreArrayRead(normsVec, &norms);
667 CHKERR PetscOptionsGetInt(PETSC_NULLPTR,
"",
"-save_every", &se, PETSC_NULLPTR);
669 if (!(ts_step % se)) {
671 <<
"Write file at time " << ts_t <<
" write step " << sTEP;
677 CHKERR calculate_reactions();
682 CHKERR calculate_error(analyticalHertzPressurePlaneStress);
685 CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
688 CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
691 CHKERR calculate_error(analyticalWavy2DPressure);
698 CHKERR print_max_min(uXScatter,
"Ux");
699 CHKERR print_max_min(uYScatter,
"Uy");
701 CHKERR print_max_min(uZScatter,
"Uz");
702 CHKERR print_force_and_area();