25  if (
auto fe_method_ptr = fePtr.lock()) {
 
   27    auto bc_mng = mField.getInterface<
BcManager>();
 
   29    const auto problem_name = fe_method_ptr->problemPtr->getName();
 
   31    for (
auto bc : bc_mng->getBcMapByBlockName()) {
 
   32      if (
auto disp_bc = bc.second->dispBcPtr) {
 
   34        auto &bc_id = bc.first;
 
   36        auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
 
   37        if (std::regex_match(bc_id, std::regex(regex_str))) {
 
   43            auto field_ptr = mField.get_field_structure(
field_name);
 
   44            return field_ptr->getNbOfCoeffs();
 
   46          const auto nb_field_coeffs = get_field_coeffs(
field_name);
 
   49              << 
"Apply EssentialPreProc<DisplacementCubitBcData>: " 
   50              << problem_name << 
"_" << 
field_name << 
"_" << block_name;
 
   56          if (
auto ext_disp_bc =
 
   59            for (
int a = 0; 
a != 3; ++
a)
 
   60              t_off(
a) = ext_disp_bc->rotOffset[
a];
 
   63          auto scale_value = [&](
const double &
c) {
 
   65            for (
auto s : vecOfTimeScalingMethods) {
 
   66              val *= s->getScale(fe_method_ptr->ts_t);
 
   71          if (disp_bc->data.flag1 == 1)
 
   72            t_vals(0) = scale_value(-disp_bc->data.value1);
 
   73          if (disp_bc->data.flag2 == 1)
 
   74            t_vals(1) = scale_value(-disp_bc->data.value2);
 
   75          if (disp_bc->data.flag3 == 1)
 
   76            t_vals(2) = scale_value(-disp_bc->data.value3);
 
   77          if (disp_bc->data.flag4 == 1)
 
   78            t_angles(0) = scale_value(-disp_bc->data.value4);
 
   79          if (disp_bc->data.flag5 == 1)
 
   80            t_angles(1) = scale_value(-disp_bc->data.value5);
 
   81          if (disp_bc->data.flag6 == 1)
 
   82            t_angles(2) = scale_value(-disp_bc->data.value6);
 
   85          std::array<std::vector<double>, 3> coords;
 
   88          const bool is_rotation =
 
   89              disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
 
   91          auto lambda = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
 
   94            auto v = t_vals(coeff);
 
   97                  coords[0][idx], coords[1][idx], coords[2][idx]);
 
   99                  t_angles, t_coords, t_off)(coeff);
 
  102              v += coords[coeff][idx];
 
  105            field_entity_ptr->getEntFieldData()[coeff] = 
v;
 
  112              [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
 
  114                auto size = field_entity_ptr->getEntFieldData().size();
 
  115                for (
int i = coeff; 
i < size; 
i += nb_field_coeffs)
 
  116                  field_entity_ptr->getEntFieldData()[
i] = 0;
 
  120          auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
 
  121          auto not_verts = subtract(bc.second->bcEnts, verts);
 
  123          if (getCoords || is_rotation) {
 
  124            for (
auto d : {0, 1, 2})
 
  125              coords[d].resize(verts.size());
 
  126            CHKERR mField.get_moab().get_coords(verts, &*coords[0].begin(),
 
  128                                                &*coords[2].begin());
 
  131          if (disp_bc->data.flag1 || disp_bc->data.flag5 ||
 
  132              disp_bc->data.flag6) {
 
  139          if (disp_bc->data.flag2 || disp_bc->data.flag4 ||
 
  140              disp_bc->data.flag6) {
 
  143            if (nb_field_coeffs > 1) {
 
  149          if (disp_bc->data.flag3 || disp_bc->data.flag4 ||
 
  150              disp_bc->data.flag5 || is_rotation) {
 
  153            if (nb_field_coeffs > 2) {
 
  165            "Can not lock shared pointer");
 
 
  180  if (
auto fe_method_ptr = fePtr.lock()) {
 
  182    auto bc_mng = mField.getInterface<
BcManager>();
 
  186    const auto problem_name = fe_method_ptr->problemPtr->getName();
 
  190    for (
auto bc : bc_mng->getBcMapByBlockName()) {
 
  191      if (
auto disp_bc = bc.second->dispBcPtr) {
 
  193        auto &bc_id = bc.first;
 
  195        auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
 
  196        if (std::regex_match(bc_id, std::regex(regex_str))) {
 
  202              << 
"Apply EssentialPreProc<DisplacementCubitBcData>: " 
  203              << problem_name << 
"_" << 
field_name << 
"_" << block_name;
 
  205          const bool is_rotation =
 
  206              disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
 
  208          auto ents = bc.second->bcEnts;
 
  210          std::array<SmartPetscObj<IS>, 3> is_xyz;
 
  211          auto prb_name = fe_method_ptr->problemPtr->getName();
 
  213          if (disp_bc->data.flag1 || is_rotation) {
 
  214            CHKERR is_mng->isCreateProblemFieldAndRankLocal(
 
  217          if (disp_bc->data.flag2 || is_rotation) {
 
  218            CHKERR is_mng->isCreateProblemFieldAndRankLocal(
 
  221          if (disp_bc->data.flag3 || is_rotation) {
 
  222            CHKERR is_mng->isCreateProblemFieldAndRankLocal(
 
  226          auto get_is_sum = [](
auto is1, 
auto is2) {
 
  232          for (
auto &is : is_xyz) {
 
  237                is_sum = get_is_sum(is_sum, is);
 
  246      if (
auto fe_ptr = fePtr.lock()) {
 
  248        auto snes_ctx = fe_ptr->snes_ctx;
 
  249        auto ts_ctx = fe_ptr->ts_ctx;
 
  254        if (fe_ptr->vecAssembleSwitch) {
 
  255          if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
 
  256            CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
 
  257            CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
 
  258            CHKERR VecAssemblyBegin(f);
 
  260            *fe_ptr->vecAssembleSwitch = 
false;
 
  264        const int *index_ptr;
 
  265        CHKERR ISGetIndices(is_sum, &index_ptr);
 
  267        CHKERR ISGetLocalSize(is_sum, &size);
 
  272        CHKERR vec_mng->setLocalGhostVector(problem_name, 
ROW, tmp_x,
 
  273                                            INSERT_VALUES, SCATTER_FORWARD);
 
  275        CHKERR VecGetArrayRead(tmp_x, &u);
 
  281          CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
 
  282          CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
  287          for (
auto i = 0; 
i != size; ++
i) {
 
  288            a[index_ptr[
i]] = vDiag * (
v[index_ptr[
i]] - u[index_ptr[
i]]);
 
  291          CHKERR VecRestoreArrayRead(x, &
v);
 
  294          for (
auto i = 0; 
i != size; ++
i) {
 
  295            a[index_ptr[
i]] = vDiag * u[index_ptr[
i]];
 
  299        CHKERR VecRestoreArrayRead(tmp_x, &u);
 
  301        CHKERR ISRestoreIndices(is_sum, &index_ptr);
 
  307            "Can not lock shared pointer");
 
 
  322  if (
auto fe_method_ptr = fePtr.lock()) {
 
  324    auto bc_mng = mField.getInterface<
BcManager>();
 
  327    const auto problem_name = fe_method_ptr->problemPtr->getName();
 
  331    for (
auto bc : bc_mng->getBcMapByBlockName()) {
 
  332      if (
auto disp_bc = bc.second->dispBcPtr) {
 
  334        auto &bc_id = bc.first;
 
  336        auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
 
  337        if (std::regex_match(bc_id, std::regex(regex_str))) {
 
  343              << 
"Apply EssentialPreProc<DisplacementCubitBcData>: " 
  344              << problem_name << 
"_" << 
field_name << 
"_" << block_name;
 
  346          const bool is_rotation =
 
  347              disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
 
  349          auto ents = bc.second->bcEnts;
 
  351          std::array<SmartPetscObj<IS>, 3> is_xyz;
 
  352          auto prb_name = fe_method_ptr->problemPtr->getName();
 
  354          if (disp_bc->data.flag1 || is_rotation) {
 
  355            CHKERR is_mng->isCreateProblemFieldAndRank(
 
  358          if (disp_bc->data.flag2 || is_rotation) {
 
  359            CHKERR is_mng->isCreateProblemFieldAndRank(
 
  362          if (disp_bc->data.flag3 || is_rotation) {
 
  363            CHKERR is_mng->isCreateProblemFieldAndRank(
 
  367          auto get_is_sum = [](
auto is1, 
auto is2) {
 
  373          for (
auto &is : is_xyz) {
 
  378                is_sum = get_is_sum(is_sum, is);
 
  387      if (
auto fe_ptr = fePtr.lock()) {
 
  391        if ((*fe_ptr->matAssembleSwitch) && !vLhs) {
 
  392          if (*fe_ptr->matAssembleSwitch) {
 
  393            CHKERR MatAssemblyBegin(
B, MAT_FINAL_ASSEMBLY);
 
  394            CHKERR MatAssemblyEnd(
B, MAT_FINAL_ASSEMBLY);
 
  395            *fe_ptr->matAssembleSwitch = 
false;
 
  399          MOFEM_LOG(
"WORLD", Sev::noisy) << 
"Apply AO to IS";
 
  400          CHKERR AOApplicationToPetscIS(vAO, is_sum);
 
  403        CHKERR MatZeroRowsColumnsIS(
B, is_sum, vDiag, PETSC_NULLPTR, PETSC_NULLPTR);
 
  409            "Can not lock shared pointer");
 
 
  435  enum { X = 0, Y, Z, MX, MY, MZ, LAST };
 
  437  if (
auto fe_ptr = fePtr.lock()) {
 
  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;
 
  451    auto get_low_hi_uid_by_entities = [](
auto bit_number, 
auto f, 
auto s) {
 
  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);
 
  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,
 
  474    auto bc_mng = mField.getInterface<
BcManager>();
 
  475    const auto problem_name = fe_ptr->problemPtr->getName();
 
  476    const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
 
  478    std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
 
  480    for (
auto bc : bc_mng->getBcMapByBlockName()) {
 
  481      if (
auto disp_bc = bc.second->dispBcPtr) {
 
  483        auto &bc_id = bc.first;
 
  485        auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
 
  486        if (std::regex_match(bc_id, std::regex(regex_str))) {
 
  492              << 
"EssentialPreProc<DisplacementCubitBcData>: " << problem_name
 
  494          auto bit_number = mField.get_field_bit_number(
field_name);
 
  497          if (
auto ext_disp_bc =
 
  500            for (
int a = 0; 
a != 3; ++
a)
 
  501              t_off(
a) = ext_disp_bc->rotOffset[
a];
 
  504          auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
 
  510          auto get_coords_vec = [&]() {
 
  513              CHKERR mField.get_moab().get_coords(verts, &*coords_vec.begin());
 
  518          auto coords_vec = get_coords_vec();
 
  519          std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
 
  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);
 
  527            for (; lo != hi; ++lo) {
 
  528              const auto loc_dof = (*lo)->getPetscLocalDofIdx();
 
  529              if (loc_dof < nb_local_dofs) {
 
  531                const auto coeff = (*lo)->getDofCoeffIdx();
 
  535                    ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
 
  537                    ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
 
  539                    ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
 
  542                  const auto ent = (*lo)->getEnt();
 
  543                  reactions[coeff] += 
a[loc_dof];
 
  547                    t_force(coeff) = 
a[loc_dof];
 
  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);
 
  560                  auto moment = [&](
auto &&t_force, 
auto &&t_coords) {
 
  563                        (FTensor::levi_civita<double>(
i, 
j, 
k) * t_coords(
k)) *
 
  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);
 
  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);
 
  590          auto mpi_reactions = mpi_array_reduce(reactions);
 
  591          if (printBlockName) {
 
  593                                "Block %s Offset: %6.4e %6.4e %6.4e",
 
  594                                block_name.c_str(), t_off(X), t_off(Y),
 
  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]);
 
  604            if ((reactionPtr != 
nullptr) &&
 
  605                (block_name.compare(reactionBlockName) == 0)) {
 
  606              (*reactionPtr)[X] = mpi_reactions[X];
 
  607              (*reactionPtr)[Y] = mpi_reactions[Y];
 
  608              (*reactionPtr)[Z] = mpi_reactions[Z];
 
  612                                "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
 
  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]);
 
  625    CHKERR VecRestoreArrayRead(f, &
a);
 
  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]);
 
  638            "Can not lock shared pointer");
 
 
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.