31 boost::shared_ptr<MFrontInterface> mfront_interface_ptr =
nullptr,
38 CHKERR DMoFEMGetInterfacePtr(
dM, &m_field_ptr);
40 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
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 &) {
52 boost::shared_ptr<MatrixDouble> mPtr;
56 auto push_domain_ops = [&](
auto &pip) {
58 pip, {
H1,
HDIV},
"GEOMETRY")),
59 "Apply base transform");
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));
70 auto hooke_common_data_ptr =
71 HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
72 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
74 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
75 return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
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),
86 auto push_bdy_ops = [&](
auto &pip) {
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>(
97 return common_data_ptr;
100 auto get_domain_pip = [&](
auto &pip)
101 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
103 auto op_loop_side =
new OpLoopSide<SideEle>(
104 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
106 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
107 pip.push_back(op_loop_side);
108 return op_loop_side->getOpPtrVector();
114 auto get_post_proc_domain_fe = [&]() {
116 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
117 auto &pip = post_proc_fe->getOpPtrVector();
119 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
120 push_domain_ops(get_domain_pip(pip));
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>();
126 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
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>>) {
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}},
148 {{
"STRESS", common_data_ptr->getMatCauchyStress()},
149 {
"STRAIN", common_data_ptr->getMatStrain()}}));
152 else if constexpr (std::is_same_v<
153 std::decay_t<
decltype(common_data_ptr)>,
154 boost::shared_ptr<HenckyOps::CommonData>>) {
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()}
165 {{
"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
168 hencky_or_hooke_common_data_ptr);
173 pip, {
HDIV},
"GEOMETRY")),
175 auto common_data_ptr = push_bdy_ops(pip);
181 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
183 {{
"SDF", common_data_ptr->sdfPtr()},
184 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
188 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
189 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
195 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
205 auto get_post_proc_bdy_fe = [&]() {
207 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
208 auto &pip = post_proc_fe->getOpPtrVector();
211 pip, {
HDIV},
"GEOMETRY")),
213 auto common_data_ptr = push_bdy_ops(pip);
216 auto op_loop_side =
new OpLoopSide<SideEle>(
217 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
219 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
220 pip.push_back(op_loop_side);
222 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
223 push_domain_ops(op_loop_side->getOpPtrVector());
225 auto X_ptr = boost::make_shared<MatrixDouble>();
227 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
233 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
235 {{
"SDF", common_data_ptr->sdfPtr()},
236 {
"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
238 {{
"U", common_data_ptr->contactDispPtr()},
240 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
241 {
"GRAD_SDF", common_data_ptr->gradSdfPtr()}
247 {{
"HESS_SDF", common_data_ptr->hessSdfPtr()}}
256 auto get_integrate_traction = [&]() {
257 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
260 integrate_traction->getOpPtrVector(), {HDIV},
"GEOMETRY")),
264 integrate_traction->getOpPtrVector().push_back(
265 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
271 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
273 "push operators to calculate traction");
275 return integrate_traction;
278 auto get_integrate_area = [&]() {
279 auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
283 integrate_area->getOpPtrVector(), {HDIV},
"GEOMETRY")),
287 integrate_area->getOpPtrVector().push_back(
288 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
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);
301 CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
302 contact_meshset_range);
303 contact_range.merge(contact_meshset_range);
306 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
308 auto op_loop_side =
new OpLoopSide<SideEle>(
309 *m_field_ptr, m_field_ptr->getInterface<Simple>()->getDomainFEName(),
312 op_loop_side->getOpPtrVector(), {H1},
"GEOMETRY");
315 (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
316 integrate_area->getOpPtrVector(), op_loop_side,
"SIGMA",
"U",
318 "push operators to calculate area");
320 return integrate_area;
323 postProcDomainFe = get_post_proc_domain_fe();
325 postProcBdyFe = get_post_proc_bdy_fe();
327 integrateTraction = get_integrate_traction();
328 integrateArea = get_integrate_area();
331 m_field_ptr->get_comm(),
332 (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
341 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
343 auto post_proc = [&]() {
346 if (!mfrontInterfacePtr) {
347 auto post_proc_begin =
348 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
351 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
354 CHKERR DMoFEMPreProcessFiniteElements(dM,
355 post_proc_begin->getFEMethod());
356 if (!postProcBdyFe) {
357 postProcDomainFe->copyTs(*
this);
358 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcDomainFe);
360 postProcDomainFe->copyTs(*
this);
361 postProcBdyFe->copyTs(*
this);
362 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", postProcDomainFe);
363 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcBdyFe);
365 CHKERR DMoFEMPostProcessFiniteElements(dM,
366 post_proc_end->getFEMethod());
368 CHKERR post_proc_end->writeFile(
369 "out_contact_" + boost::lexical_cast<std::string>(sTEP) +
".h5m");
371 CHKERR mfrontInterfacePtr->postProcess(
373 m_field_ptr->
getInterface<Simple>()->getDomainFEName());
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);
388 auto calculate_area = [&] {
390 integrateArea->copyTs(*
this);
391 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", integrateArea);
392 CHKERR VecAssemblyBegin(CommonData::totalTraction);
393 CHKERR VecAssemblyEnd(CommonData::totalTraction);
397 auto calculate_reactions = [&]() {
400 auto res = createDMVector(dM);
402 auto assemble_domain = [&]() {
404 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
405 auto &pip = fe_rhs->getOpPtrVector();
416 ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
419 if (!mfrontInterfacePtr) {
421 HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
422 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
424 CHKERR mfrontInterfacePtr->opFactoryDomainRhs(pip,
"U");
426 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", fe_rhs);
428 CHKERR VecAssemblyBegin(res);
429 CHKERR VecAssemblyEnd(res);
430 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
431 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
436 auto assemble_boundary = [&]() {
438 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
439 auto &pip = fe_rhs->getOpPtrVector();
451 pip.push_back(
new OpSetHOWeightsOnSubDim<SPACE_DIM>());
453 auto u_disp = boost::make_shared<MatrixDouble>();
454 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_disp));
456 new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
460 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", fe_rhs);
466 CHKERR assemble_boundary();
468 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
469 auto get_post_proc_hook_rhs = [
this, fe_post_proc_ptr, res,
472 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
473 *m_field_ptr, fe_post_proc_ptr, res)();
476 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
477 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
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);
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);
496 auto print_force_and_area = [&]() {
499 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
502 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
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]);
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);
514 if (mfrontInterfacePtr) {
515 CHKERR mfrontInterfacePtr->updateInternalVariables(
516 dM, m_field_ptr->
getInterface<Simple>()->getDomainFEName());
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)
529 magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
530 gradSDFPtr(grad_sdf_ptr) {}
531 MoFEMErrorCode doWork(
int, EntityType, EntitiesFieldData::EntData &) {
536 magPtr->resize(pPtr->size());
538 tyPtr->resize(pPtr->size());
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);
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);
560 ++t_contact_traction;
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;
576 auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
577 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
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;
585 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
588 contact_meshset_range);
589 contact_range.merge(contact_meshset_range);
594 post_proc_norm_fe->getOpPtrVector(), {HDIV},
"GEOMETRY")),
598 post_proc_norm_fe->getOpPtrVector().push_back(
599 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
600 post_proc_norm_fe->getRuleHook = [](int, int,
int approx_order) {
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>(
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);
621 post_proc_norm_fe->getOpPtrVector().push_back(
622 new OpGetTensor0fromFunc(analytical_pressure_ptr,
fun));
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()));
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));
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));
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));
644 CHKERR VecZeroEntries(normsVec);
645 post_proc_norm_fe->copyTs(*
this);
646 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", post_proc_norm_fe);
647 CHKERR VecAssemblyBegin(normsVec);
648 CHKERR VecAssemblyEnd(normsVec);
653 CHKERR VecGetArrayRead(normsVec, &norms);
655 <<
"norm_traction: " << std::scientific
656 << std::sqrt(norms[TRACTION_NORM_L2]);
658 <<
"norm_mag_traction: " << std::scientific
659 << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
661 <<
"norm_traction_y: " << std::scientific
662 << std::sqrt(norms[TRACTION_Y_NORM_L2]);
663 CHKERR VecRestoreArrayRead(normsVec, &norms);
669 CHKERR PetscOptionsGetInt(PETSC_NULLPTR,
"",
"-save_every", &se,
672 if (!(ts_step % se)) {
674 <<
"Write file at time " << ts_t <<
" write step " << sTEP;
680 CHKERR calculate_reactions();
685 CHKERR calculate_error(analyticalHertzPressurePlaneStress);
688 CHKERR calculate_error(analyticalHertzPressurePlaneStrain);
691 CHKERR calculate_error(analyticalHertzPressureAxisymmetric);
694 CHKERR calculate_error(analyticalWavy2DPressure);
701 CHKERR print_max_min(uXScatter,
"Ux");
702 CHKERR print_max_min(uYScatter,
"Uy");
704 CHKERR print_max_min(uZScatter,
"Uz");
705 CHKERR print_force_and_area();