433 {
436
437 enum { X = 0, Y, Z, MX, MY, MZ, LAST };
438
439 if (
auto fe_ptr =
fePtr.lock()) {
440
441 SmartPetscObj<Vec>
f =
vRhs ?
vRhs : SmartPetscObj<Vec>(fe_ptr->f,
true);
442
443 if (fe_ptr->vecAssembleSwitch) {
444 if ((*fe_ptr->vecAssembleSwitch) && !
vRhs) {
445 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
446 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
447 CHKERR VecAssemblyBegin(f);
449 *fe_ptr->vecAssembleSwitch = false;
450 }
451 }
452
453 auto get_low_hi_uid_by_entities = [](
auto bit_number,
auto f,
auto s) {
456 };
457
458 auto get_low_hi = [fe_ptr](auto lo_uid, auto hi_uid) {
459 auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
460 .lower_bound(lo_uid);
461 auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
462 .upper_bound(hi_uid);
463 return std::make_pair(it, hi_it);
464 };
465
466 auto mpi_array_reduce = [this](auto &array) {
467 std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
468 MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
470 return array_sum;
471 };
472
475
477 const auto problem_name = fe_ptr->problemPtr->getName();
478 const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
479
480 std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
483 }
484
485 for (auto bc : bc_mng->getBcMapByBlockName()) {
486 if (auto disp_bc = bc.second->dispBcPtr) {
487
488 auto &bc_id = bc.first;
489
490 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
491 if (std::regex_match(bc_id, std::regex(regex_str))) {
492
495
497 << "EssentialPreProc<DisplacementCubitBcData>: " << problem_name
500
502 if (auto ext_disp_bc =
503 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
504 disp_bc.get())) {
505 for (
int a = 0;
a != 3; ++
a)
506 t_off(
a) = ext_disp_bc->rotOffset[
a];
507 }
508
509 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
510
514
515 auto get_coords_vec = [&]() {
517 if (verts.size()) {
519 }
520 return coords_vec;
521 };
522
523 auto coords_vec = get_coords_vec();
524 std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
525
526 for (auto pit = verts.const_pair_begin();
527 pit != verts.const_pair_end(); ++pit) {
528 auto [lo_uid, hi_uid] =
529 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
530 auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
531
532 for (; lo != hi; ++lo) {
533 const auto loc_dof = (*lo)->getPetscLocalDofIdx();
534 if (loc_dof < nb_local_dofs) {
535
536 const auto coeff = (*lo)->getDofCoeffIdx();
537
538 if (
539
540 ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
541 coeff == 0) ||
542 ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
543 coeff == 1) ||
544 ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
545 coeff == 2)) {
546
547 const auto ent = (*lo)->getEnt();
548 reactions[coeff] +=
a[loc_dof];
549
550 auto force = [&]() {
552 t_force(coeff) =
a[loc_dof];
553 return t_force;
554 };
555
556 auto coord = [&]() {
557 const auto idx = verts.index(ent);
559 coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
560 coords_vec[3 * idx + Z]};
561 t_coords(
i) -= t_off(
i);
562 return t_coords;
563 };
564
565 auto moment = [&](auto &&t_force, auto &&t_coords) {
568 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords(
k)) *
570 return t_moment;
571 };
572
573 auto t_moment = moment(force(), coord());
574 reactions[MX] += t_moment(X);
575 reactions[MY] += t_moment(Y);
576 reactions[MZ] += t_moment(Z);
577 }
578 }
579 }
580 }
581
583 reactions[Z]};
585 reactions[MZ]};
587 &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
589 &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
590 t_total_force(
i) += t_force(
i);
593 (FTensor::levi_civita<double>(
i,
j,
k) * t_off(
k)) * t_force(
j);
594
595 auto mpi_reactions = mpi_array_reduce(reactions);
597 (*blockReactions)[block_name] = mpi_reactions;
598 }
601 "Block %s Offset: %6.4e %6.4e %6.4e",
602 block_name.c_str(), t_off(X), t_off(Y),
603 t_off(Z));
605 "Block %s Force: %6.4e %6.4e %6.4e",
606 block_name.c_str(), mpi_reactions[X],
607 mpi_reactions[Y], mpi_reactions[Z]);
609 "Block %s Moment: %6.4e %6.4e %6.4e",
610 block_name.c_str(), mpi_reactions[MX],
611 mpi_reactions[MY], mpi_reactions[MZ]);
614 (*reactionPtr)[X] = mpi_reactions[X];
615 (*reactionPtr)[Y] = mpi_reactions[Y];
616 (*reactionPtr)[Z] = mpi_reactions[Z];
617 }
618 } else {
620 "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
621 t_off(Z));
623 "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
624 mpi_reactions[Y], mpi_reactions[Z]);
626 "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
627 mpi_reactions[MY], mpi_reactions[MZ]);
628 }
629 }
630 }
631 }
632
633 CHKERR VecRestoreArrayRead(f, &
a);
634
635 auto mpi_total_reactions = mpi_array_reduce(total_reactions);
637 "WORLD",
sevLevel,
"Essential",
"Total force: %6.4e %6.4e %6.4e",
638 mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
640 "Total moment: %6.4e %6.4e %6.4e",
641 mpi_total_reactions[MX], mpi_total_reactions[MY],
642 mpi_total_reactions[MZ]);
643
644 } else {
646 "Can not lock shared pointer");
647 }
648
650}
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
UBlasVector< double > VectorDouble
constexpr auto field_name
static std::pair< std::string, std::string > extractStringFromBlockId(const std::string block_id, const std::string prb_name)
Extract block name and block name from block id.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
