312 auto post_proc = [&]() {
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());
335 post_proc_end->getFEMethod());
337 CHKERR post_proc_end->writeFile(
338 "out_contact_" + boost::lexical_cast<std::string>(
sTEP) +
".h5m");
342 m_field_ptr->
getInterface<Simple>()->getDomainFEName());
348 auto calculate_force = [&] {
357 auto calculate_area = [&] {
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();
382 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {
H1},
385 ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
390 HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
391 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
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();
416 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
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 = [&]() {
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]);
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>();
550 m_field_ptr->
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
551 std::regex((boost::format(
"%s(.*)") %
"CONTACT").str()))) {
552 auto meshset =
m->getMeshset();
553 Range contact_meshset_range;
555 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
558 contact_meshset_range);
559 contact_range.merge(contact_meshset_range);
563 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
564 post_proc_norm_fe->getOpPtrVector(), {HDIV},
"GEOMETRY")),
568 post_proc_norm_fe->getOpPtrVector().push_back(
569 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
574 post_proc_norm_fe->getOpPtrVector().push_back(
575 new OpCalculateVectorFieldValues<SPACE_DIM>(
576 "U", common_data_ptr->contactDispPtr()));
577 post_proc_norm_fe->getOpPtrVector().push_back(
578 new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
579 "SIGMA", common_data_ptr->contactTractionPtr()));
580 using C = ContactIntegrators<BoundaryEleOp>;
581 post_proc_norm_fe->getOpPtrVector().push_back(
582 new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
585 auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
586 auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
587 auto mag_traction_ptr = boost::make_shared<VectorDouble>();
588 auto traction_y_ptr = boost::make_shared<VectorDouble>();
589 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
591 post_proc_norm_fe->getOpPtrVector().push_back(
592 new OpGetTensor0fromFunc(analytical_pressure_ptr,
fun));
594 post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcTractions(
595 analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
596 traction_y_ptr, common_data_ptr->contactTractionPtr(),
597 common_data_ptr->gradSdfPtr()));
599 post_proc_norm_fe->getOpPtrVector().push_back(
600 new OpCalcNormL2Tensor1<SPACE_DIM>(
602 analytical_traction_ptr, contact_range_ptr));
606 post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor0(
608 analytical_pressure_ptr, contact_range_ptr));
610 post_proc_norm_fe->getOpPtrVector().push_back(
612 analytical_pressure_ptr, contact_range_ptr));
615 post_proc_norm_fe->copyTs(*
this);
625 <<
"norm_traction: " << std::scientific
628 <<
"norm_mag_traction: " << std::scientific
631 <<
"norm_traction_y: " << std::scientific
641 if (!(ts_step % se)) {
643 <<
"Write file at time " << ts_t <<
" write step " <<
sTEP;
649 CHKERR calculate_reactions();
674 CHKERR print_force_and_area();