431 enum { X = 0, Y, Z, MX, MY, MZ,
LAST };
433 if (
auto fe_ptr =
fePtr.lock()) {
435 SmartPetscObj<Vec>
f =
vRhs ?
vRhs : SmartPetscObj<Vec>(fe_ptr->f,
true);
437 if (fe_ptr->vecAssembleSwitch) {
438 if ((*fe_ptr->vecAssembleSwitch) && !
vRhs) {
439 CHKERR VecGhostUpdateBegin(
f, ADD_VALUES, SCATTER_REVERSE);
440 CHKERR VecGhostUpdateEnd(
f, ADD_VALUES, SCATTER_REVERSE);
443 *fe_ptr->vecAssembleSwitch =
false;
447 auto get_low_hi_uid_by_entities = [](
auto bit_number,
auto f,
auto s) {
452 auto get_low_hi = [fe_ptr](
auto lo_uid,
auto hi_uid) {
453 auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
454 .lower_bound(lo_uid);
455 auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
456 .upper_bound(hi_uid);
457 return std::make_pair(it, hi_it);
460 auto mpi_array_reduce = [
this](
auto &array) {
461 std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
462 MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
471 const auto problem_name = fe_ptr->problemPtr->getName();
472 const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
474 std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
476 for (
auto bc : bc_mng->getBcMapByBlockName()) {
477 if (
auto disp_bc = bc.second->dispBcPtr) {
479 auto &bc_id = bc.first;
481 auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
482 if (std::regex_match(bc_id, std::regex(regex_str))) {
488 <<
"EssentialPreProc<DisplacementCubitBcData>: " << problem_name
493 if (
auto ext_disp_bc =
494 dynamic_cast<DisplacementCubitBcDataWithRotation
const *
>(
496 for (
int a = 0;
a != 3; ++
a)
497 t_off(
a) = ext_disp_bc->rotOffset[
a];
500 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
506 auto get_coords_vec = [&]() {
514 auto coords_vec = get_coords_vec();
515 std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
517 for (
auto pit = verts.const_pair_begin();
518 pit != verts.const_pair_end(); ++pit) {
519 auto [lo_uid, hi_uid] =
520 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
521 auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
523 for (; lo != hi; ++lo) {
524 const auto loc_dof = (*lo)->getPetscLocalDofIdx();
525 if (loc_dof < nb_local_dofs) {
527 const auto coeff = (*lo)->getDofCoeffIdx();
531 ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
533 ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
535 ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
538 const auto ent = (*lo)->getEnt();
539 reactions[coeff] +=
a[loc_dof];
543 t_force(coeff) =
a[loc_dof];
548 const auto idx = verts.index(ent);
550 coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
551 coords_vec[3 * idx + Z]};
552 t_coords(
i) -= t_off(
i);
556 auto moment = [&](
auto &&t_force,
auto &&t_coords) {
559 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords(
k)) *
564 auto t_moment = moment(force(), coord());
565 reactions[MX] += t_moment(X);
566 reactions[MY] += t_moment(Y);
567 reactions[MZ] += t_moment(Z);
578 &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
580 &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
581 t_total_force(
i) += t_force(
i);
584 (FTensor::levi_civita<double>(
i,
j,
k) * t_off(
k)) * t_force(
j);
586 auto mpi_reactions = mpi_array_reduce(reactions);
589 "Block %s Offset: %6.4e %6.4e %6.4e",
590 block_name.c_str(), t_off(X), t_off(Y),
593 "Block %s Force: %6.4e %6.4e %6.4e",
594 block_name.c_str(), mpi_reactions[X],
595 mpi_reactions[Y], mpi_reactions[Z]);
597 "Block %s Moment: %6.4e %6.4e %6.4e",
598 block_name.c_str(), mpi_reactions[MX],
599 mpi_reactions[MY], mpi_reactions[MZ]);
602 "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
605 "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
606 mpi_reactions[Y], mpi_reactions[Z]);
608 "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
609 mpi_reactions[MY], mpi_reactions[MZ]);
617 auto mpi_total_reactions = mpi_array_reduce(total_reactions);
619 "WORLD",
sevLevel,
"Essential",
"Total force: %6.4e %6.4e %6.4e",
620 mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
622 "Total moment: %6.4e %6.4e %6.4e",
623 mpi_total_reactions[MX], mpi_total_reactions[MY],
624 mpi_total_reactions[MZ]);
628 "Can not lock shared pointer");