431                                                                             {
  434 
  435  enum { X = 0, Y, Z, MX, MY, MZ, LAST };
  436 
  437  if (
auto fe_ptr = 
fePtr.lock()) {
 
  438 
  439    SmartPetscObj<Vec> 
f = 
vRhs ? 
vRhs : SmartPetscObj<Vec>(fe_ptr->f, 
true);
 
  440 
  441    if (fe_ptr->vecAssembleSwitch) {
  442      if ((*fe_ptr->vecAssembleSwitch) && !
vRhs) {
 
  443        CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
 
  444        CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
 
  445        CHKERR VecAssemblyBegin(f);
 
  447        *fe_ptr->vecAssembleSwitch = false;
  448      }
  449    }
  450 
  451    auto get_low_hi_uid_by_entities = [](
auto bit_number, 
auto f, 
auto s) {
 
  454    };
  455 
  456    auto get_low_hi = [fe_ptr](auto lo_uid, auto hi_uid) {
  457      auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
  458                    .lower_bound(lo_uid);
  459      auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
  460                       .upper_bound(hi_uid);
  461      return std::make_pair(it, hi_it);
  462    };
  463 
  464    auto mpi_array_reduce = [this](auto &array) {
  465      std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
  466      MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
  468      return array_sum;
  469    };
  470 
  473 
  475    const auto problem_name = fe_ptr->problemPtr->getName();
  476    const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
  477 
  478    std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
  479 
  480    for (auto bc : bc_mng->getBcMapByBlockName()) {
  481      if (auto disp_bc = bc.second->dispBcPtr) {
  482 
  483        auto &bc_id = bc.first;
  484 
  485        auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
  486        if (std::regex_match(bc_id, std::regex(regex_str))) {
  487 
  490 
  492              << "EssentialPreProc<DisplacementCubitBcData>: " << problem_name
  495 
  497          if (auto ext_disp_bc =
  498                  dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
  499                      disp_bc.get())) {
  500            for (
int a = 0; 
a != 3; ++
a)
 
  501              t_off(
a) = ext_disp_bc->rotOffset[
a];
 
  502          }
  503 
  504          auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
  505 
  509 
  510          auto get_coords_vec = [&]() {
  512            if (verts.size()) {
  514            }
  515            return coords_vec;
  516          };
  517 
  518          auto coords_vec = get_coords_vec();
  519          std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
  520 
  521          for (auto pit = verts.const_pair_begin();
  522               pit != verts.const_pair_end(); ++pit) {
  523            auto [lo_uid, hi_uid] =
  524                get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
  525            auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
  526 
  527            for (; lo != hi; ++lo) {
  528              const auto loc_dof = (*lo)->getPetscLocalDofIdx();
  529              if (loc_dof < nb_local_dofs) {
  530 
  531                const auto coeff = (*lo)->getDofCoeffIdx();
  532 
  533                if (
  534 
  535                    ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
  536                     coeff == 0) ||
  537                    ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
  538                     coeff == 1) ||
  539                    ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
  540                     coeff == 2)) {
  541 
  542                  const auto ent = (*lo)->getEnt();
  543                  reactions[coeff] += 
a[loc_dof];
 
  544 
  545                  auto force = [&]() {
  547                    t_force(coeff) = 
a[loc_dof];
 
  548                    return t_force;
  549                  };
  550 
  551                  auto coord = [&]() {
  552                    const auto idx = verts.index(ent);
  554                        coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
  555                        coords_vec[3 * idx + Z]};
  556                    t_coords(
i) -= t_off(
i);
 
  557                    return t_coords;
  558                  };
  559 
  560                  auto moment = [&](auto &&t_force, auto &&t_coords) {
  563                        (FTensor::levi_civita<double>(
i, 
j, 
k) * t_coords(
k)) *
 
  565                    return t_moment;
  566                  };
  567 
  568                  auto t_moment = moment(force(), coord());
  569                  reactions[MX] += t_moment(X);
  570                  reactions[MY] += t_moment(Y);
  571                  reactions[MZ] += t_moment(Z);
  572                }
  573              }
  574            }
  575          }
  576 
  578                                              reactions[Z]};
  580                                               reactions[MZ]};
  582              &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
  584              &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
  585          t_total_force(
i) += t_force(
i);
 
  588              (FTensor::levi_civita<double>(
i, 
j, 
k) * t_off(
k)) * t_force(
j);
 
  589 
  590          auto mpi_reactions = mpi_array_reduce(reactions);
  593                                "Block %s Offset: %6.4e %6.4e %6.4e",
  594                                block_name.c_str(), t_off(X), t_off(Y),
  595                                t_off(Z));
  597                                "Block %s Force: %6.4e %6.4e %6.4e",
  598                                block_name.c_str(), mpi_reactions[X],
  599                                mpi_reactions[Y], mpi_reactions[Z]);
  601                                "Block %s Moment: %6.4e %6.4e %6.4e",
  602                                block_name.c_str(), mpi_reactions[MX],
  603                                mpi_reactions[MY], mpi_reactions[MZ]);
  606              (*reactionPtr)[X] = mpi_reactions[X];
  607              (*reactionPtr)[Y] = mpi_reactions[Y];
  608              (*reactionPtr)[Z] = mpi_reactions[Z];
  609            }
  610          } else {
  612                                "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
  613                                t_off(Z));
  615                                "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
  616                                mpi_reactions[Y], mpi_reactions[Z]);
  618                                "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
  619                                mpi_reactions[MY], mpi_reactions[MZ]);
  620          }
  621        }
  622      }
  623    }
  624 
  625    CHKERR VecRestoreArrayRead(f, &
a);
 
  626 
  627    auto mpi_total_reactions = mpi_array_reduce(total_reactions);
  629        "WORLD", 
sevLevel, 
"Essential", 
"Total force: %6.4e %6.4e %6.4e",
 
  630        mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
  632                        "Total moment: %6.4e %6.4e %6.4e",
  633                        mpi_total_reactions[MX], mpi_total_reactions[MY],
  634                        mpi_total_reactions[MZ]);
  635 
  636  } else {
  638            "Can not lock shared pointer");
  639  }
  640 
  642}
#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.