15                   {
   16  IdxDataType(const UId uid, const int dof) {
   17    bcopy(&uid, dAta, 4 * sizeof(int));
   18    dAta[4] = dof;
   19  }
   20 
   21private:
   22  int dAta[5];
   23};
   24 
   25struct IdxDataTypePtr {
   26  IdxDataTypePtr(const int *ptr) : pTr(ptr) {}
   27  inline int getDofIdx() const {
   28    int global_dof = pTr[4];
   29    return global_dof;
   30  }
   31  inline UId getUId()
 const {
 
   32    unsigned int b0, b1, b2, b3;
   33    bcopy(&pTr[0], &b0, sizeof(int));
   34    bcopy(&pTr[1], &b1, sizeof(int));
   35    bcopy(&pTr[2], &b2, sizeof(int));
   36    bcopy(&pTr[3], &b3, sizeof(int));
   37    UId uid = 
static_cast<UId>(b0) | 
static_cast<UId>(b1) << 8 * 
sizeof(
int) |
 
   38              static_cast<UId>(b2) << 16 * 
sizeof(
int) |
 
   39              static_cast<UId>(b3) << 24 * 
sizeof(
int);
 
   40    return uid;
   41  }
   42 
   43private:
   44  const int *pTr;
   45};
   46 
   48ProblemsManager::query_interface(boost::typeindex::type_index type_index,
   49                                 UnknownInterface **iface) const {
   50  *iface = const_cast<ProblemsManager *>(this);
   51  return 0;
   52}
   53 
   54ProblemsManager::ProblemsManager(
const MoFEM::Core &core)
 
   56      buildProblemFromFields(PETSC_FALSE),
   57      synchroniseProblemEntities(PETSC_FALSE) {
   58  PetscLogEventRegister("ProblemsManager", 0, &MOFEM_EVENT_ProblemsManager);
   59}
   60 
   64  PetscOptionsBegin(m_field.
get_comm(), 
"", 
"Problem manager", 
"none");
 
   65  {
   67        "-problem_build_from_fields",
   68        "Add DOFs to problem directly from fields not through DOFs on elements",
   69        "", buildProblemFromFields, &buildProblemFromFields, NULL);
   70  }
   71  PetscOptionsEnd();
   73}
   74 
   76    const Range &ents, 
const int dim, 
const int adj_dim, 
const int n_parts,
 
   77    Tag *th_vertex_weights, 
Tag *th_edge_weights, 
Tag *th_part_weights,
 
   78    int verb, 
const bool debug) {
 
   80      .getInterface<CommInterface>()
   81      ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
   82                      th_edge_weights, th_part_weights, verb, 
debug);
 
   83}
   84 
   85MoFEMErrorCode ProblemsManager::buildProblem(
const std::string name,
 
   86                                             const bool square_matrix,
   87                                             int verb) {
   88 
   91  if (!(cOre.getBuildMoFEM() & (1 << 0)))
   93  if (!(cOre.getBuildMoFEM() & (1 << 1)))
   95  if (!(cOre.getBuildMoFEM() & (1 << 2)))
   97  const Problem *problem_ptr;
   99  CHKERR buildProblem(
const_cast<Problem *
>(problem_ptr), square_matrix, verb);
 
  100  cOre.getBuildMoFEM() |= 1 << 3; 
  101                                  
  103}
  104 
  106                                             const bool square_matrix,
  107                                             int verb) {
  109  auto dofs_field_ptr = m_field.
get_dofs();
 
  113  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
  114 
  115  
  116  
  117 
  118  if (problem_ptr->getBitRefLevel().none()) {
  119    SETERRQ(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
  120             problem_ptr->getName().c_str());
  121  }
  123 
  124  
  125  problem_ptr->numeredFiniteElementsPtr->clear();
  126 
  128  auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
  130    for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
  131 
  132      const auto uid = (*it)->getLocalUniqueId();
  133 
  134      auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
 
  135      for (
auto lo = 
r.first; lo != 
r.second; ++lo) {
 
  136 
  137        if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
  138          std::array<bool, 2> row_col = {false, false};
  139 
  140          const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
 
  141          const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
 
  142          const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
 
  143 
  144          
  145          if ((fe_bit & prb_mask) != fe_bit)
  146            continue;
  147          if ((fe_bit & prb_bit).none())
  148            continue;
  149 
  150          auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
  151            auto dit = nb_dofs->lower_bound(uid);
  152            decltype(dit) hi_dit;
  153            if (dit != nb_dofs->end()) {
  154              hi_dit = nb_dofs->upper_bound(
  156              view.insert(dit, hi_dit);
  157              row_col[rc] = true;
  158            }
  159          };
  160 
  161          if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
  162            add_to_view(dofs_field_ptr, dofs_rows, 
ROW);
 
  163 
  164          if (!square_matrix)
  165            if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
  166              add_to_view(dofs_field_ptr, dofs_cols, 
COL);
 
  167 
  168          if (square_matrix && row_col[
ROW])
 
  169            break;
  170          else if (row_col[
ROW] && row_col[
COL])
 
  171            break;
  172        }
  173      }
  174    }
  176  };
  177 
  178  CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
 
  179 
  180  
  181  {
  182    
  183    problem_ptr->nbDofsRow = 0;
  184    problem_ptr->nbLocDofsRow = 0;
  185    problem_ptr->nbGhostDofsRow = 0;
  186 
  187    
  188    DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
  189        hi_miit;
  190    hi_miit = dofs_rows.get<0>().end();
  191 
  192    int count_dofs = dofs_rows.get<1>().count(true);
  193    boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
  194        boost::shared_ptr<std::vector<NumeredDofEntity>>(
  195            new std::vector<NumeredDofEntity>());
  196    problem_ptr->getRowDofsSequence()->push_back(dofs_array);
  197    dofs_array->reserve(count_dofs);
  198    miit = dofs_rows.get<0>().begin();
  199    for (; miit != hi_miit; miit++) {
  200      if ((*miit)->getActive()) {
  201        dofs_array->emplace_back(*miit);
  202        dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
  203      }
  204    }
  205    auto hint = problem_ptr->numeredRowDofsPtr->end();
  206    for (
auto &
v : *dofs_array) {
 
  207      hint = problem_ptr->numeredRowDofsPtr->emplace_hint(hint, dofs_array, &
v);
 
  208    }
  209  }
  210 
  211  
  212  if (!square_matrix) {
  213    
  214    problem_ptr->nbDofsCol = 0;
  215    problem_ptr->nbLocDofsCol = 0;
  216    problem_ptr->nbGhostDofsCol = 0;
  217 
  218    
  219    DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
  220        hi_miit;
  221    hi_miit = dofs_cols.get<0>().end();
  222 
  223    int count_dofs = 0;
  224    miit = dofs_cols.get<0>().begin();
  225    for (; miit != hi_miit; miit++) {
  226      if (!(*miit)->getActive()) {
  227        continue;
  228      }
  229      count_dofs++;
  230    }
  231 
  232    boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
  233        boost::shared_ptr<std::vector<NumeredDofEntity>>(
  234            new std::vector<NumeredDofEntity>());
  235    problem_ptr->getColDofsSequence()->push_back(dofs_array);
  236    dofs_array->reserve(count_dofs);
  237    miit = dofs_cols.get<0>().begin();
  238    for (; miit != hi_miit; miit++) {
  239      if (!(*miit)->getActive()) {
  240        continue;
  241      }
  242      dofs_array->emplace_back(*miit);
  243      dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
  244    }
  245    auto hint = problem_ptr->numeredColDofsPtr->end();
  246    for (
auto &
v : *dofs_array) {
 
  247      hint = problem_ptr->numeredColDofsPtr->emplace_hint(hint, dofs_array, &
v);
 
  248    }
  249  } else {
  250    problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
  251    problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
  252    problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
  253  }
  254 
  255  
  258        << problem_ptr->getName() << " Nb. local dofs "
  259        << problem_ptr->numeredRowDofsPtr->size() << " by "
  260        << problem_ptr->numeredColDofsPtr->size();
  262  }
  263 
  266        << "FEs row dofs " << *problem_ptr << " Nb. row dof "
  267        << problem_ptr->getNbDofsRow();
  268    for (auto &miit : *problem_ptr->numeredRowDofsPtr)
  270 
  272        << "FEs col dofs " << *problem_ptr << " Nb. col dof "
  273        << problem_ptr->getNbDofsCol();
  274    for (auto &miit : *problem_ptr->numeredColDofsPtr)
  277  }
  278 
  279  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; 
  280                                               
  281                                               
  282 
  283  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
  284 
  286}
  287 
  289    const std::string name, const bool square_matrix, int verb) {
  292 
  293  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
  295  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
  297  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
  299 
  300  const Problem *problem_ptr;
  302  CHKERR buildProblemOnDistributedMesh(
const_cast<Problem *
>(problem_ptr),
 
  303                                       square_matrix, verb);
  304 
  305  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM;
  306  cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
  307 
  309}
  310 
  312    Problem *problem_ptr, const bool square_matrix, int verb) {
  316  auto dofs_field_ptr = m_field.
get_dofs();
 
  318  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
  319 
  320  
  322 
  324 
  325  if (problem_ptr->getBitRefLevel().none())
  327             "problem <%s> refinement level not set",
  328             problem_ptr->getName().c_str());
  329 
  330  int loop_size = 2;
  331  if (square_matrix) {
  332    loop_size = 1;
  333    problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
  334  } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
  335    problem_ptr->numeredColDofsPtr =
  336        boost::shared_ptr<NumeredDofEntity_multiIndex>(
  338  }
  339 
  340  const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
 
  341  const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
 
  342 
  343  
  345 
  346  
  347  
  348  if (buildProblemFromFields == PETSC_FALSE) {
  349 
  350    auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
  354      for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
  355           ++it) {
  356 
  357        const auto uid = (*it)->getLocalUniqueId();
  358 
  359        auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
 
  360        for (
auto lo = 
r.first; lo != 
r.second; ++lo) {
 
  361 
  362          if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
  363            std::array<bool, 2> row_col = {false, false};
  364 
  365            const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
 
  366 
  367            
  368            if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
  369 
  370              auto add_to_view = [&](auto dofs, auto &view, auto rc) {
  371                auto dit = dofs->lower_bound(uid);
  372                decltype(dit) hi_dit;
  373                if (dit != dofs->end()) {
  374                  hi_dit = dofs->upper_bound(
  376                  view.insert(dit, hi_dit);
  377                  row_col[rc] = true;
  378                }
  379              };
  380 
  381              if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
  382                add_to_view(dofs_field_ptr, dofs_rows, 
ROW);
 
  383 
  384              if (!square_matrix)
  385                if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
  386                  add_to_view(dofs_field_ptr, dofs_cols, 
COL);
 
  387 
  388              if (square_matrix && row_col[
ROW])
 
  389                break;
  390              else if (row_col[
ROW] && row_col[
COL])
 
  391                break;
  392            }
  393          }
  394        }
  395      }
  397    };
  398 
  399    CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
 
  400  }
  401 
  402  
  403  
  404  if (buildProblemFromFields == PETSC_TRUE) {
  405    
  407    for (auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
  408      if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
  409        fields_ids_row |= fit->get()->getBitFieldIdRow();
  410        fields_ids_col |= fit->get()->getBitFieldIdCol();
  411      }
  412    }
  413    
  414    for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
  415      if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
  416 
  417        auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
  418            FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
  419        auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
  420            FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
  421 
  422        for (; dit != hi_dit; dit++) {
  423 
  424          const int owner_proc = dit->get()->getOwnerProc();
  426            const unsigned char pstatus = dit->get()->getPStatus();
  427            if (pstatus == 0) {
  428              continue;
  429            }
  430          }
  431 
  432          const auto &dof_bit = (*dit)->getBitRefLevel();
  433          
  434          if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
  435 
  436            if ((fit->get()->getId() & fields_ids_row).any()) {
  437              dofs_rows.insert(*dit);
  438            }
  439            if (!square_matrix) {
  440              if ((fit->get()->getId() & fields_ids_col).any()) {
  441                dofs_cols.insert(*dit);
  442              }
  443            }
  444          }
  445        }
  446      }
  447    }
  448  }
  449 
  450  if (synchroniseProblemEntities) {
  451    
  453    BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
 
  454    for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
  455         fit != fe_ptr->end(); fit++) {
  456      if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
  457        fields_ids_row |= fit->get()->getBitFieldIdRow();
  458        fields_ids_col |= fit->get()->getBitFieldIdCol();
  459      }
  460    }
  461 
  463    for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
  464      DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
  465          hi_miit;
  466      miit = dofs_ptr[ss]->get<1>().lower_bound(1);
  467      hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
  468      Range ents_to_synchronise;
 
  469      for (; miit != hi_miit; ++miit) {
  470        if (miit->get()->getEntDofIdx() != 0)
  471          continue;
  472        ents_to_synchronise.insert(miit->get()->getEnt());
  473      }
  474      Range tmp_ents = ents_to_synchronise;
 
  476          ents_to_synchronise, nullptr, verb);
  477      ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
  478      for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
  479        if ((fit->get()->getId() & *fields_ids[ss]).any()) {
  480          const auto bit_number = (*fit)->getBitNumber();
  481          for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
  482               pit != ents_to_synchronise.pair_end(); ++pit) {
  483            const auto f = pit->first;
 
  484            const auto s = pit->second;
  485            const auto lo_uid =
  486                FieldEntity::getLocalUniqueIdCalculate(bit_number, f);
  487            const auto hi_uid =
  488                FieldEntity::getLocalUniqueIdCalculate(bit_number, s);
  489 
  490            auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
  491            auto hi_dit =
  492                dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
  493 
  494            dofs_ptr[ss]->insert(dit, hi_dit);
  495          }
  496        }
  497      }
  498    }
  499  }
  500 
  501  
  503  boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
  504      problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
  505  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
  506  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
  507                            &problem_ptr->nbLocDofsCol};
  508  int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
  509                            &problem_ptr->nbGhostDofsCol};
  510  for (int ss = 0; ss < 2; ss++) {
  511    *(nbdof_ptr[ss]) = 0;
  512    *(local_nbdof_ptr[ss]) = 0;
  513    *(ghost_nbdof_ptr[ss]) = 0;
  514  }
  515 
  516  
  517  
  518  int nb_local_dofs[] = {0, 0};
  519  for (int ss = 0; ss < loop_size; ss++) {
  520    DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
  521        hi_miit;
  522    miit = dofs_ptr[ss]->get<1>().lower_bound(1);
  523    hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
  524    for (; miit != hi_miit; miit++) {
  525      int owner_proc = (*miit)->getOwnerProc();
  527        nb_local_dofs[ss]++;
  528      }
  529    }
  530  }
  531 
  532  
  533  int start_ranges[2], end_ranges[2];
  534  for (int ss = 0; ss != loop_size; ss++) {
  535    PetscLayout layout;
  537    CHKERR PetscLayoutSetBlockSize(layout, 1);
 
  538    CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
 
  539    CHKERR PetscLayoutSetUp(layout);
 
  540    CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
 
  541    CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
 
  542    CHKERR PetscLayoutDestroy(&layout);
 
  543  }
  544  if (square_matrix) {
  545    nbdof_ptr[1] = nbdof_ptr[0];
  546    nb_local_dofs[1] = nb_local_dofs[0];
  547    start_ranges[1] = start_ranges[0];
  548    end_ranges[1] = end_ranges[0];
  549  }
  550 
  551  
  552  const size_t idx_data_type_size = sizeof(IdxDataType);
  553  const size_t data_block_size = idx_data_type_size / 
sizeof(
int);
 
  554 
  555  if (
sizeof(IdxDataType) % 
sizeof(
int)) {
 
  557  }
  558  std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
  561 
  562  
  563  
  564  for (int ss = 0; ss != loop_size; ++ss) {
  565 
  566    DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
  567        hi_miit;
  568    hi_miit = dofs_ptr[ss]->get<0>().end();
  569 
  570    boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
  571        boost::shared_ptr<std::vector<NumeredDofEntity>>(
  572            new std::vector<NumeredDofEntity>());
  573    int nb_dofs_to_add = 0;
  574    miit = dofs_ptr[ss]->get<0>().begin();
  575    for (; miit != hi_miit; ++miit) {
  576      
  577      if (!(miit->get()->getActive()))
  578        continue;
  579      ++nb_dofs_to_add;
  580    }
  581    dofs_array->reserve(nb_dofs_to_add);
  582    if (ss == 0) {
  583      problem_ptr->getRowDofsSequence()->push_back(dofs_array);
  584    } else {
  585      problem_ptr->getColDofsSequence()->push_back(dofs_array);
  586    }
  587 
  588    int &local_idx = *local_nbdof_ptr[ss];
  589    miit = dofs_ptr[ss]->get<0>().begin();
  590    for (; miit != hi_miit; ++miit) {
  591 
  592      
  593      if (!(miit->get()->getActive()))
  594        continue;
  595 
  596      dofs_array->emplace_back(*miit);
  597 
  598      int owner_proc = dofs_array->back().getOwnerProc();
  599      if (owner_proc < 0) {
  601                "data inconsistency");
  602      }
  603 
  605        dofs_array->back().pArt = owner_proc;
  606        dofs_array->back().dofIdx = -1;
  607        dofs_array->back().petscGloablDofIdx = -1;
  608        dofs_array->back().petscLocalDofIdx = -1;
  609      } else {
  610 
  611        
  612        int glob_idx = start_ranges[ss] + local_idx;
  613        dofs_array->back().pArt = owner_proc;
  614        dofs_array->back().dofIdx = glob_idx;
  615        dofs_array->back().petscGloablDofIdx = glob_idx;
  616        dofs_array->back().petscLocalDofIdx = local_idx;
  617        local_idx++;
  618 
  619        unsigned char pstatus = dofs_array->back().getPStatus();
  620        
  621        
  622        if (pstatus > 0) {
  623 
  624          for (int proc = 0;
  625               proc < MAX_SHARING_PROCS &&
  626               -1 != dofs_array->back().getSharingProcsPtr()[proc];
  627               proc++) {
  628            
  629            if (ss == 0) {
  630              ids_data_packed_rows[dofs_array->back()
  631                                       .getSharingProcsPtr()[proc]]
  632                  .emplace_back(dofs_array->back().getGlobalUniqueId(),
  633                                glob_idx);
  634            } else {
  635              ids_data_packed_cols[dofs_array->back()
  636                                       .getSharingProcsPtr()[proc]]
  637                  .emplace_back(dofs_array->back().getGlobalUniqueId(),
  638                                glob_idx);
  639            }
  640            if (!(pstatus & PSTATUS_MULTISHARED)) {
  641              break;
  642            }
  643          }
  644        }
  645      }
  646    }
  647 
  648    auto hint = numered_dofs_ptr[ss]->end();
  649    for (
auto &
v : *dofs_array)
 
  650      hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &
v);
 
  651  }
  652  if (square_matrix) {
  653    local_nbdof_ptr[1] = local_nbdof_ptr[0];
  654  }
  655 
  656  int nsends_rows = 0, nsends_cols = 0;
  657  
  660  lengths_rows.clear();
  661  lengths_cols.clear();
  663    lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
  664    lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
  665    if (!ids_data_packed_rows[proc].empty())
  666      nsends_rows++;
  667    if (!ids_data_packed_cols[proc].empty())
  668      nsends_cols++;
  669  }
  670 
  671  MPI_Status *status;
  673 
  674  
  675  int nrecvs_rows;    
  676  int *onodes_rows;   
  677  int *olengths_rows; 
  678  int **rbuf_row;     
  679 
  680  
  681  MPI_Comm comm;
  683 
  684  {
  685 
  686    
  687 
  688    
  689    CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
 
  690                                       &nrecvs_rows);
  691    
  692 
  693    
  694    
  695    CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
 
  696                                     &lengths_rows[0], &onodes_rows,
  697                                     &olengths_rows);
  698 
  699    
  700    
  701    
  702    
  703    int tag_row;
  704    CHKERR PetscCommGetNewTag(comm, &tag_row);
 
  705 
  706    
  707    
  708    MPI_Request *r_waits_row; 
  709    
  710    
  711    
  712 
  713    CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
 
  714                             olengths_rows, &rbuf_row, &r_waits_row);
  715    CHKERR PetscFree(onodes_rows);
 
  716 
  717    MPI_Request *s_waits_row; 
  718    CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
 
  719 
  720    
  721    for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
 
  722      if (!lengths_rows[proc])
  723        continue; 
  724      CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], 
 
  725                       lengths_rows[proc],               
  726                       MPIU_INT, proc,                   
  727                       tag_row, comm, s_waits_row + kk);
  728      kk++;
  729    }
  730 
  731    if (nrecvs_rows) {
  732      CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
 
  733    }
  734    if (nsends_rows) {
  735      CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
 
  736    }
  737 
  738    CHKERR PetscFree(r_waits_row);
 
  739    CHKERR PetscFree(s_waits_row);
 
  740  }
  741 
  742  
  743  int nrecvs_cols = nrecvs_rows;
  744  int *olengths_cols = olengths_rows;
  745  PetscInt **rbuf_col = rbuf_row;
  746  if (!square_matrix) {
  747 
  748    
  749    CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
 
  750                                       &nrecvs_cols);
  751 
  752    
  753    
  754    int *onodes_cols;
  755    CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
 
  756                                     &lengths_cols[0], &onodes_cols,
  757                                     &olengths_cols);
  758 
  759    
  760    int tag_col;
  761    CHKERR PetscCommGetNewTag(comm, &tag_col);
 
  762 
  763    
  764    
  765    MPI_Request *r_waits_col; 
  766    CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
 
  767                             olengths_cols, &rbuf_col, &r_waits_col);
  768    CHKERR PetscFree(onodes_cols);
 
  769 
  770    MPI_Request *s_waits_col; 
  771    CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
 
  772 
  773    
  774    for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
 
  775      if (!lengths_cols[proc])
  776        continue; 
  777      CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], 
 
  778                       lengths_cols[proc],               
  779                       MPIU_INT, proc,                   
  780                       tag_col, comm, s_waits_col + kk);
  781      kk++;
  782    }
  783 
  784    if (nrecvs_cols) {
  785      CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
 
  786    }
  787    if (nsends_cols) {
  788      CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
 
  789    }
  790 
  791    CHKERR PetscFree(r_waits_col);
 
  792    CHKERR PetscFree(s_waits_col);
 
  793  }
  794 
  795  CHKERR PetscCommDestroy(&comm);
 
  797 
  799  auto hint = dofs_glob_uid_view.begin();
  800  for (auto dof : *m_field.get_dofs())
  801    dofs_glob_uid_view.emplace_hint(hint, dof);
  802 
  803  
  804  for (int ss = 0; ss != loop_size; ++ss) {
  805 
  806    int nrecvs;
  807    int *olengths;
  808    int **data_procs;
  809    if (ss == 0) {
  810      nrecvs = nrecvs_rows;
  811      olengths = olengths_rows;
  812      data_procs = rbuf_row;
  813    } else {
  814      nrecvs = nrecvs_cols;
  815      olengths = olengths_cols;
  816      data_procs = rbuf_col;
  817    }
  818 
  820    for (int kk = 0; kk != nrecvs; ++kk) {
  821      int len = olengths[kk];
  822      int *data_from_proc = data_procs[kk];
  823      for (
int dd = 0; 
dd < len; 
dd += data_block_size) {
 
  824        uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
  825        auto ddit = dofs_glob_uid_view.find(uid);
  826        const auto owner_proc = FieldEntity::getOwnerFromUniqueId(uid);
  827 
  829            PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
  830 
  831#ifndef NDEBUG
  834              << "DOF is shared or multishared between processors. For example "
  835                 "if order of field on given entity is set in inconsistently, "
  836                 "has different value on two processor, error such as this is "
  837                 "triggered";
  838 
  839          MOFEM_LOG(
"SELF", Sev::error) << 
"UId " << uid << 
" is not found";
 
  841              << "Problematic UId owner proc is " << owner_proc;
  842          const auto uid_handle = FieldEntity::getHandleFromUniqueId(uid);
  844              << "Problematic UId entity owning processor handle is "
  845              << uid_handle << " entity type "
  847          const auto uid_bit_number =
  848              FieldEntity::getFieldBitNumberFromUniqueId(uid);
  850              << "Problematic UId field is "
  853 
  855          field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
  856          auto fe_it = field_view.find(FieldEntity::getGlobalUniqueIdCalculate(
  857              owner_proc, uid_bit_number, uid_handle));
  858          if (fe_it == field_view.end()) {
  860                << "Also, no field entity in database for given global UId";
  861          } else {
  862            MOFEM_LOG(
"SELF", Sev::error) << 
"Field entity in databse exist " 
  863                                             "(but have no DOF wih give UId";
  864            MOFEM_LOG(
"SELF", Sev::error) << **fe_it;
 
  865 
  866            
  867            auto error_file_name =
  868                "error_with_missing_entity_" +
  870                ".vtk";
  872                << "Look to file < " << error_file_name
  873                << " > it contains entity with missing DOF.";
  874 
  876            const auto local_fe_ent = (*fe_it)->getEnt();
  877            CHKERR m_field.
get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
 
  878            CHKERR m_field.
get_moab().write_file(error_file_name.c_str(), 
"VTK",
 
  879                                                 "", tmp_msh->get_ptr(), 1);
  880          }
  881#endif
  882 
  884                  "DOF with global UId not found (Compile code in Debug to "
  885                  "learn more about problem");
  886        }
  887 
  888        if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
  889          continue;
  890        }
  891 
  892        auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
  893 
  894        if (dit != numered_dofs_ptr[ss]->end()) {
  895 
  896          int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
  897#ifndef NDEBUG
  898          if (PetscUnlikely(global_idx < 0))
  900                    "received negative dof");
  901#endif
  902          bool success;
  903          success = numered_dofs_ptr[ss]->modify(
  904              dit, NumeredDofEntity_mofem_index_change(global_idx));
  905 
  906#ifndef NDEBUG
  907          if (PetscUnlikely(!success))
  909                    "modification unsuccessful");
  910#endif
  911          success = numered_dofs_ptr[ss]->modify(
  912              dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
  913                                                             global_idx));
  914#ifndef NDEBUG
  915          if (PetscUnlikely(!success))
  917                    "modification unsuccessful");
  918#endif
  919        };
  920 
  921#ifndef NDEBUG
  922        if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
  923 
  924          
  925          
  926          
  927          
  928 
  929          std::ostringstream zz;
  930          zz << **ddit << std::endl;
  932                   "data inconsistency, dofs is not shared, but received "
  933                   "from other proc\n"
  934                   "%s",
  935                   zz.str().c_str());
  936        }
  937#endif
  938      }
  939    }
  940  }
  941 
  942  if (square_matrix) {
  943    (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
  944    (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
  945  }
  946 
  947  CHKERR PetscFree(olengths_rows);
 
  948  CHKERR PetscFree(rbuf_row[0]);
 
  949  CHKERR PetscFree(rbuf_row);
 
  950  if (!square_matrix) {
  951    CHKERR PetscFree(olengths_cols);
 
  952    CHKERR PetscFree(rbuf_col[0]);
 
  953    CHKERR PetscFree(rbuf_col);
 
  954  }
  955 
  956  if (square_matrix) {
  957    if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
  959              "data inconsistency for square_matrix %ld!=%ld",
  960              numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
  961    }
  962    if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
  964              "data inconsistency for square_matrix");
  965    }
  966  }
  967 
  968  CHKERR printPartitionedProblem(problem_ptr, verb);
 
  969  CHKERR debugPartitionedProblem(problem_ptr, verb);
 
  970 
  971  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
  972 
  974}
  975 
  977    const std::string out_name,
  978 
  979    const std::vector<std::string> &fields_row,
  980    const std::vector<std::string> &fields_col,
  981 
  982    const std::string main_problem, const bool square_matrix,
  983 
  984    const map<std::string, boost::shared_ptr<Range>> *entityMapRow,
  985    const map<std::string, boost::shared_ptr<Range>> *entityMapCol,
  986 
  987    int verb) {
  991 
  993 
  994  
  995  using ProblemByName = decltype(problems_ptr->get<Problem_mi_tag>());
  996  auto &problems_by_name =
  997      const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
  998 
  999  
 1000  auto out_problem_it = problems_by_name.find(out_name);
 1001  if (out_problem_it == problems_by_name.end()) {
 1003             "subproblem with name < %s > not defined (top tip check spelling)",
 1004             out_name.c_str());
 1005  }
 1006  
 1007  
 1008  auto main_problem_it = problems_by_name.find(main_problem);
 1009  if (main_problem_it == problems_by_name.end()) {
 1011             "problem of subproblem with name < %s > not defined (top tip "
 1012             "check spelling)",
 1013             main_problem.c_str());
 1014  }
 1015 
 1016  
 1017  boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
 1018      out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
 1019  
 1020  boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
 1021      main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
 1022  
 1023  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
 1024                          &out_problem_it->nbLocDofsCol};
 1025  
 1026  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
 1027 
 1028  
 1029  {
 1030    out_problem_it->nbGhostDofsRow = 0;
 1031    out_problem_it->nbGhostDofsCol = 0;
 1032  }
 1033 
 1034  
 1035  std::vector<std::string> fields[] = {fields_row, fields_col};
 1036  const map<std::string, boost::shared_ptr<Range>> *entityMap[] = {
 1037      entityMapRow, entityMapCol};
 1038 
 1039  
 1040  out_problem_it->subProblemData =
 1041      boost::make_shared<Problem::SubProblemData>();
 1042 
 1043  
 1044  for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
 1045 
 1046    
 1047    (*nb_local_dofs[ss]) = 0;
 1048    (*nb_dofs[ss]) = 0;
 1049    
 1050    out_problem_dofs[ss]->clear();
 1051 
 1052    
 1053    out_problem_it->numeredFiniteElementsPtr->clear();
 1054 
 1055    
 1056    for (auto field : fields[ss]) {
 1057 
 1058      
 1059      
 1060 
 1061      
 1062      boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
 1063          boost::make_shared<std::vector<NumeredDofEntity>>();
 1064      
 1065      if (!ss)
 1066        out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
 1067      else
 1068        out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
 1069 
 1070      
 1072 
 1073      auto add_dit_to_dofs_array = [&](auto &dit) {
 1074        if (dit->get()->getPetscGlobalDofIdx() >= 0)
 1075          dofs_array->emplace_back(
 1076              dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
 1077              dit->get()->getPetscGlobalDofIdx(),
 1078              dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
 1079      };
 1080 
 1081      auto get_dafult_dof_range = [&]() {
 1082        auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
 1083            FieldEntity::getLoBitNumberUId(bit_number));
 1084        auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
 1085            FieldEntity::getHiBitNumberUId(bit_number));
 1086        return std::make_pair(dit, hi_dit);
 1087      };
 1088 
 1089      auto check = [&](auto &field) -> boost::shared_ptr<Range> {
 1090        if (entityMap[ss]) {
 1091          auto mit = entityMap[ss]->find(field);
 1092          if (mit != entityMap[ss]->end()) {
 1093            return mit->second;
 1094          } else {
 1095            return nullptr;
 1096          }
 1097        } else {
 1098          return nullptr;
 1099        }
 1100      };
 1101 
 1102      if (auto r_ptr = check(field)) {
 1103        for (auto p = r_ptr->pair_begin(); p != r_ptr->pair_end(); ++p) {
 1104          const auto lo_ent = p->first;
 1105          const auto hi_ent = p->second;
 1106          auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
 1107              DofEntity::getLoFieldEntityUId(bit_number, lo_ent));
 1108          auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
 1109              DofEntity::getHiFieldEntityUId(bit_number, hi_ent));
 1110          dofs_array->reserve(std::distance(dit, hi_dit));
 1111          for (; dit != hi_dit; dit++) {
 1112            add_dit_to_dofs_array(dit);
 1113          }
 1114        }
 1115      } else {
 1116        auto [dit, hi_dit] = get_dafult_dof_range();
 1117        dofs_array->reserve(std::distance(dit, hi_dit));
 1118        for (; dit != hi_dit; dit++)
 1119          add_dit_to_dofs_array(dit);
 1120      }
 1121 
 1122      
 1123      auto hint = out_problem_dofs[ss]->end();
 1124      for (
auto &
v : *dofs_array)
 
 1125        hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &
v);
 
 1126    }
 1127    
 1128    {
 1129      auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
 1130      auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
 1131      for (; dit != hi_dit; dit++) {
 1132        int idx = -1; 
 1133        if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
 
 1134          idx = (*nb_local_dofs[ss])++;
 1135        }
 1136        bool success = out_problem_dofs[ss]->modify(
 1137            out_problem_dofs[ss]->project<0>(dit),
 1138            NumeredDofEntity_local_idx_change(idx));
 1139        if (!success) {
 1141                  "operation unsuccessful");
 1142        }
 1143      };
 1144    }
 1145    
 1146    {
 1147      auto dit =
 1148          out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
 1149      auto hi_dit =
 1150          out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
 1151              out_problem_dofs[ss]->size());
 1152      const int nb = std::distance(dit, hi_dit);
 1153      
 1154      std::vector<int> main_indices(nb);
 1155      for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
 1156        *it = dit->get()->getPetscGlobalDofIdx();
 1157      }
 1158      
 1159      IS is;
 1160      CHKERR ISCreateGeneral(m_field.
get_comm(), nb, &*main_indices.begin(),
 
 1161                             PETSC_USE_POINTER, &is);
 1162      
 1163      
 1164      AO ao;
 1165      CHKERR AOCreateMappingIS(is, PETSC_NULLPTR, &ao);
 
 1166      if (ss == 0) {
 1167        IS is_dup;
 1168        CHKERR ISDuplicate(is, &is_dup);
 
 1169        out_problem_it->getSubData()->rowIs = SmartPetscObj<IS>(is_dup, false);
 1170        out_problem_it->getSubData()->rowMap = SmartPetscObj<AO>(ao, true);
 1171      } else {
 1172        IS is_dup;
 1173        CHKERR ISDuplicate(is, &is_dup);
 
 1174        out_problem_it->getSubData()->colIs = SmartPetscObj<IS>(is_dup, false);
 1175        out_problem_it->getSubData()->colMap = SmartPetscObj<AO>(ao, true);
 1176      }
 1177      CHKERR AOApplicationToPetscIS(ao, is);
 
 1178      
 1179      CHKERR ISGetSize(is, nb_dofs[ss]);
 
 1181      
 1182      dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
 1183      for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
 1184           dit++, it++) {
 1185        bool success = out_problem_dofs[ss]->modify(
 1186            out_problem_dofs[ss]->project<0>(dit),
 1187            NumeredDofEntity_part_and_all_indices_change(
 1188                dit->get()->getPart(), *it, *it,
 1189                dit->get()->getPetscLocalDofIdx()));
 1190        if (!success) {
 1192                  "operation unsuccessful");
 1193        }
 1194      }
 1195      
 1196      {
 1197        NumeredDofEntityByLocalIdx::iterator dit =
 1198            out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
 1199        NumeredDofEntityByLocalIdx::iterator hi_dit =
 1200            out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
 1201        const int nb = std::distance(dit, hi_dit);
 1202        std::vector<int> main_indices_non_local(nb);
 1203        for (auto it = main_indices_non_local.begin(); dit != hi_dit;
 1204             dit++, it++) {
 1205          *it = dit->get()->getPetscGlobalDofIdx();
 1206        }
 1207        IS is;
 1209                               &*main_indices_non_local.begin(),
 1210                               PETSC_USE_POINTER, &is);
 1211        CHKERR AOApplicationToPetscIS(ao, is);
 
 1213        dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
 1214        for (auto it = main_indices_non_local.begin(); dit != hi_dit;
 1215             dit++, it++) {
 1216          bool success = out_problem_dofs[ss]->modify(
 1217              out_problem_dofs[ss]->project<0>(dit),
 1218              NumeredDofEntity_part_and_all_indices_change(
 1219                  dit->get()->getPart(), dit->get()->getDofIdx(), *it,
 1220                  dit->get()->getPetscLocalDofIdx()));
 1221          if (!success) {
 1223                    "operation unsuccessful");
 1224          }
 1225        }
 1226      }
 1228    }
 1229  }
 1230 
 1231  if (square_matrix) {
 1232    out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
 1233    out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
 1234    out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
 1235    out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
 1236    out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
 1237  }
 1238 
 1239  CHKERR printPartitionedProblem(&*out_problem_it, verb);
 
 1240  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
 
 1241 
 1243}
 1244 
 1246    const std::string out_name, const std::vector<std::string> add_row_problems,
 1247    const std::vector<std::string> add_col_problems, const bool square_matrix,
 1248    int verb) {
 1249  if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
 1251  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
 1253  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
 1258 
 1260  
 1261  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
 1262  ProblemByName &problems_by_name =
 1263      const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
 1264 
 1265  
 1266  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
 1267  if (out_problem_it == problems_by_name.end()) {
 1269             "problem with name < %s > not defined (top tip check spelling)",
 1270             out_name.c_str());
 1271  }
 1272  
 1273  out_problem_it->composedProblemsData =
 1274      boost::make_shared<ComposedProblemsData>();
 1275  boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
 1276      out_problem_it->getComposedProblemsData();
 1277 
 1278  const std::vector<std::string> *add_prb[] = {&add_row_problems,
 1279                                               &add_col_problems};
 1280  std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
 1281                                                 &cmp_prb_data->colProblemsAdd};
 1282  std::vector<SmartPetscObj<IS>> *add_prb_is[] = {&cmp_prb_data->rowIs,
 1283                                                  &cmp_prb_data->colIs};
 1284 
 1285  
 1286  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
 1287                          &out_problem_it->nbLocDofsCol};
 1288  
 1289  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
 1290 
 1291  
 1292  {
 1293    out_problem_it->nbDofsRow = 0;
 1294    out_problem_it->nbDofsCol = 0;
 1295    out_problem_it->nbLocDofsRow = 0;
 1296    out_problem_it->nbLocDofsCol = 0;
 1297    out_problem_it->nbGhostDofsRow = 0;
 1298    out_problem_it->nbGhostDofsCol = 0;
 1299  }
 1300  int nb_dofs_reserve[] = {0, 0};
 1301 
 1302  
 1303  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
 1304    add_prb_ptr[ss]->reserve(add_prb[ss]->size());
 1305    add_prb_is[ss]->reserve(add_prb[ss]->size());
 1306    for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
 1307         vit != add_prb[ss]->end(); vit++) {
 1308      ProblemByName::iterator prb_it = problems_by_name.find(*vit);
 1309      if (prb_it == problems_by_name.end()) {
 1310        SETERRQ(
 1312            "problem with name < %s > not defined (top tip check spelling)",
 1313            vit->c_str());
 1314      }
 1315      add_prb_ptr[ss]->push_back(&*prb_it);
 1316      
 1317      if (ss == 0) {
 1318        
 1319        *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
 1320        *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
 1321        nb_dofs_reserve[ss] +=
 1322            add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
 1323      } else {
 1324        
 1325        *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
 1326        *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
 1327        nb_dofs_reserve[ss] +=
 1328            add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
 1329      }
 1330    }
 1331  }
 1332  
 1333  if (square_matrix) {
 1334    add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
 1335    add_prb_is[1]->reserve(add_prb_ptr[0]->size());
 1336    out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
 1337    *nb_dofs[1] = *nb_dofs[0];
 1338    *nb_local_dofs[1] = *nb_local_dofs[0];
 1339  }
 1340 
 1341  
 1342  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
 1343  
 1344  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
 1345    dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
 1346    dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
 1347    if (!ss)
 1348      out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
 1349    else
 1350      out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
 1351  }
 1352 
 1353  
 1354  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
 1355    NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
 1356        dit,
 1357        hi_dit;
 1358    int shift_glob = 0;
 1359    int shift_loc = 0;
 1360    for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
 1361      PetscInt *dofs_out_idx_ptr;
 1362      int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
 1363      CHKERR PetscMalloc(nb_local_dofs * 
sizeof(
int), &dofs_out_idx_ptr);
 
 1364      if (ss == 0) {
 1365        dit = (*add_prb_ptr[ss])[pp]
 1366                  ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
 1367                  .begin();
 1368        hi_dit = (*add_prb_ptr[ss])[pp]
 1369                     ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
 1370                     .end();
 1371      } else {
 1372        dit = (*add_prb_ptr[ss])[pp]
 1373                  ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
 1374                  .begin();
 1375        hi_dit = (*add_prb_ptr[ss])[pp]
 1376                     ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
 1377                     .end();
 1378      }
 1379      int is_nb = 0;
 1380      for (; dit != hi_dit; dit++) {
 1381        const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
 
 1382        const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
 
 1383        const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
 
 1384        if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
 1385          continue;
 1387        const int part = dit->get()->getPart();
 1388        const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
 1389        const int loc_idx =
 1390            (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
 1391                           : -1;
 1392        dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
 1393                                     glob_idx, loc_idx, part);
 1394        if (part == rank) {
 1395          dofs_out_idx_ptr[is_nb++] = glob_idx;
 1396        }
 1397      }
 1398      if (is_nb > nb_local_dofs) {
 1400                "Data inconsistency");
 1401      }
 1402      IS is;
 1403      CHKERR ISCreateGeneral(m_field.
get_comm(), is_nb, dofs_out_idx_ptr,
 
 1404                             PETSC_OWN_POINTER, &is);
 1405      auto smart_is = SmartPetscObj<IS>(is);
 1406      (*add_prb_is[ss]).push_back(smart_is);
 1407      if (ss == 0) {
 1408        shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
 1409        shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
 1410      } else {
 1411        shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
 1412        shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
 1413      }
 1414      if (square_matrix) {
 1415        (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
 1416        (*add_prb_is[1]).push_back(smart_is);
 1417      }
 1418    }
 1419  }
 1420 
 1421  if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
 1423  }
 1424 
 1425  
 1426  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
 1427    auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
 1428                          : out_problem_it->numeredColDofsPtr->end();
 1429    for (
auto &
v : *dofs_array[ss])
 
 1430      hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
 1431                             hint, dofs_array[ss], &
v)
 
 1432                       : out_problem_it->numeredColDofsPtr->emplace_hint(
 1433                             hint, dofs_array[ss], &
v);
 
 1434  }
 1435 
 1436  
 1437  *nb_dofs[0] = 0;
 1438  *nb_dofs[1] = 0;
 1439  *nb_local_dofs[0] = 0;
 1440  *nb_local_dofs[1] = 0;
 1441  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
 1442 
 1443    boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
 1444    if (ss == 0) {
 1445      dofs_ptr = out_problem_it->numeredRowDofsPtr;
 1446    } else {
 1447      dofs_ptr = out_problem_it->numeredColDofsPtr;
 1448    }
 1449    NumeredDofEntityByUId::iterator dit, hi_dit;
 1450    dit = dofs_ptr->get<Unique_mi_tag>().begin();
 1451    hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
 1452    std::vector<int> idx;
 1453    idx.reserve(std::distance(dit, hi_dit));
 1454    
 1455    for (; dit != hi_dit; dit++) {
 1456      if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
 
 1457        bool success = dofs_ptr->get<Unique_mi_tag>().modify(
 1458            dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
 1459        if (!success) {
 1461                  "modification unsuccessful");
 1462        }
 1463        idx.push_back(dit->get()->getPetscGlobalDofIdx());
 1464      } else {
 1465        if (dit->get()->getPetscLocalDofIdx() != -1) {
 1467                  "local index should be negative");
 1468        }
 1469      }
 1470    }
 1471    if (square_matrix) {
 1472      *nb_local_dofs[1] = *nb_local_dofs[0];
 1473    }
 1474 
 1475    
 1476    IS is;
 1477    CHKERR ISCreateGeneral(m_field.
get_comm(), idx.size(), &*idx.begin(),
 
 1478                           PETSC_USE_POINTER, &is);
 1479    CHKERR ISGetSize(is, nb_dofs[ss]);
 
 1480    if (square_matrix) {
 1481      *nb_dofs[1] = *nb_dofs[0];
 1482    }
 1483 
 1484    AO ao;
 1485    CHKERR AOCreateMappingIS(is, PETSC_NULLPTR, &ao);
 
 1486    for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
 1487      CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
 
 1488 
 1489    
 1490    {
 1491      std::vector<int> idx_new;
 1492      idx_new.reserve(dofs_ptr->size());
 1493      for (NumeredDofEntityByUId::iterator dit =
 1494               dofs_ptr->get<Unique_mi_tag>().begin();
 1495           dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
 1496        idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
 1497      }
 1498      
 1499      IS is_new;
 1501                             &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
 1502      CHKERR AOApplicationToPetscIS(ao, is_new);
 
 1503      
 1504      std::vector<int>::iterator vit = idx_new.begin();
 1505      for (NumeredDofEntityByUId::iterator dit =
 1506               dofs_ptr->get<Unique_mi_tag>().begin();
 1507           dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
 1508        bool success =
 1509            dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
 1510                                      dit->get()->getPart(), *(vit++)));
 1511        if (!success) {
 1513                  "modification unsuccessful");
 1514        }
 1515      }
 1516      CHKERR ISDestroy(&is_new);
 
 1517    }
 1520  }
 1521 
 1522  CHKERR printPartitionedProblem(&*out_problem_it, verb);
 
 1523  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
 
 1524 
 1525  
 1526  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM;
 1527  cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
 1528 
 1530}
 1531 
 1532MoFEMErrorCode ProblemsManager::partitionSimpleProblem(
const std::string name,
 
 1533                                                       int verb) {
 1534 
 1538 
 1539  if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
 1541  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
 1543  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
 1545  if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
 1547  MOFEM_LOG(
"WORLD", Sev::verbose) << 
"Simple partition problem " << name;
 
 1548 
 1549  
 1550  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
 1551  ProblemByName &problems_set =
 1552      const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
 1553  ProblemByName::iterator p_miit = problems_set.find(name);
 1554  if (p_miit == problems_set.end()) {
 1555    SETERRQ(PETSC_COMM_SELF, 1,
 1556             "problem < %s > is not found (top tip: check spelling)",
 1557             name.c_str());
 1558  }
 1560                                    Idx_mi_tag>::type NumeredDofEntitysByIdx;
 1561  NumeredDofEntitysByIdx &dofs_row_by_idx =
 1562      p_miit->numeredRowDofsPtr->get<Idx_mi_tag>();
 1563  NumeredDofEntitysByIdx &dofs_col_by_idx =
 1564      p_miit->numeredColDofsPtr->get<Idx_mi_tag>();
 1566                            Idx_mi_tag>::type::iterator miit_row,
 1567      hi_miit_row;
 1569                            Idx_mi_tag>::type::iterator miit_col,
 1570      hi_miit_col;
 1571  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
 
 1572  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
 
 1573  nb_row_local_dofs = 0;
 1574  nb_row_ghost_dofs = 0;
 1575  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
 
 1576  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
 
 1577  nb_col_local_dofs = 0;
 1578  nb_col_ghost_dofs = 0;
 1579 
 1580  bool square_matrix = false;
 1581  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
 1582    square_matrix = true;
 1583  }
 1584 
 1585  
 1586  PetscLayout layout_row;
 1587  const int *ranges_row;
 1588 
 1589  DofIdx nb_dofs_row = dofs_row_by_idx.size();
 
 1591  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
 
 1592  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
 
 1593  CHKERR PetscLayoutSetUp(layout_row);
 
 1594  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
 
 1595  
 1596  PetscLayout layout_col;
 1597  const int *ranges_col;
 1598  if (!square_matrix) {
 1599    DofIdx nb_dofs_col = dofs_col_by_idx.size();
 
 1601    CHKERR PetscLayoutSetBlockSize(layout_col, 1);
 
 1602    CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
 
 1603    CHKERR PetscLayoutSetUp(layout_col);
 
 1604    CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
 
 1605  }
 1607       part++) {
 1608    miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
 1609    hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
 1610    if (std::distance(miit_row, hi_miit_row) !=
 1611        ranges_row[part + 1] - ranges_row[part]) {
 1612      SETERRQ(
 1613          PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
 1614          "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
 1615          "rstart (%ld != %d - %d = %d) ",
 1616          std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
 1617          ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
 1618    }
 1619    
 1620    for (; miit_row != hi_miit_row; miit_row++) {
 1621      bool success = dofs_row_by_idx.modify(
 1622          miit_row,
 1623          NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
 1624      if (!success)
 1626                "modification unsuccessful");
 1628        success = dofs_row_by_idx.modify(
 1629            miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
 1630        if (!success)
 1632                  "modification unsuccessful");
 1633      }
 1634    }
 1635    if (!square_matrix) {
 1636      miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
 1637      hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
 1638      if (std::distance(miit_col, hi_miit_col) !=
 1639          ranges_col[part + 1] - ranges_col[part]) {
 1640        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
 1641                 "data inconsistency, std::distance(miit_col,hi_miit_col) != "
 1642                 "rend - "
 1643                 "rstart (%ld != %d - %d = %d) ",
 1644                 std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
 1645                 ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
 1646      }
 1647      
 1648      for (; miit_col != hi_miit_col; miit_col++) {
 1649        bool success = dofs_col_by_idx.modify(
 1650            miit_col, NumeredDofEntity_part_and_glob_idx_change(
 1651                          part, (*miit_col)->dofIdx));
 1652        if (!success)
 1654                  "modification unsuccessful");
 1656          success = dofs_col_by_idx.modify(
 1657              miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
 1658          if (!success)
 1660                    "modification unsuccessful");
 1661        }
 1662      }
 1663    }
 1664  }
 1665  CHKERR PetscLayoutDestroy(&layout_row);
 
 1666  if (!square_matrix) {
 1667    CHKERR PetscLayoutDestroy(&layout_col);
 
 1668  }
 1669  if (square_matrix) {
 1670    nb_col_local_dofs = nb_row_local_dofs;
 1671    nb_col_ghost_dofs = nb_row_ghost_dofs;
 1672  }
 1673  CHKERR printPartitionedProblem(&*p_miit, verb);
 
 1674  cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
 1676}
 1677 
 1678MoFEMErrorCode ProblemsManager::partitionProblem(
const std::string name,
 
 1679                                                 int verb) {
 1683 
 1684  MOFEM_LOG(
"WORLD", Sev::noisy) << 
"Partition problem " << name;
 
 1685 
 1686  using NumeredDofEntitysByIdx =
 1687      NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type;
 1688  using ProblemsByName = Problem_multiIndex::index<Problem_mi_tag>::type;
 1689 
 1690  
 1691  auto &problems_set =
 1692      const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
 1693  auto p_miit = problems_set.find(name);
 1694  if (p_miit == problems_set.end()) {
 1696             "problem with name %s not defined (top tip check spelling)",
 1697             name.c_str());
 1698  }
 1699  int nb_dofs_row = p_miit->getNbDofsRow();
 1700 
 1702 
 1703    if (!(cOre.getBuildMoFEM() & (1 << 0)))
 1705    if (!(cOre.getBuildMoFEM() & (1 << 1)))
 1707    if (!(cOre.getBuildMoFEM() & (1 << 2)))
 1709              "entFEAdjacencies not build");
 1710    if (!(cOre.getBuildMoFEM() & (1 << 3)))
 1712 
 1713    Mat Adj;
 1715        ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
 1716 
 1720      MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
 1721 
 1722    
 1723    MatPartitioning part;
 1724    IS is;
 1726    CHKERR MatPartitioningSetAdjacency(part, Adj);
 
 1727    CHKERR MatPartitioningSetFromOptions(part);
 
 1729    CHKERR MatPartitioningApply(part, &is);
 
 1731      ISView(is, PETSC_VIEWER_STDOUT_WORLD);
 1732 
 1733    
 1734    IS is_gather, is_num, is_gather_num;
 1735    CHKERR ISAllGather(is, &is_gather);
 
 1736    CHKERR ISPartitioningToNumbering(is, &is_num);
 
 1737    CHKERR ISAllGather(is_num, &is_gather_num);
 
 1738    const int *part_number, *petsc_idx;
 1739    CHKERR ISGetIndices(is_gather, &part_number);
 
 1740    CHKERR ISGetIndices(is_gather_num, &petsc_idx);
 
 1741    int size_is_num, size_is_gather;
 1742    CHKERR ISGetSize(is_gather, &size_is_gather);
 
 1743    if (size_is_gather != (int)nb_dofs_row)
 1744      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
 1745               "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
 1746 
 1747    CHKERR ISGetSize(is_num, &size_is_num);
 
 1748    if (size_is_num != (int)nb_dofs_row)
 1749      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
 1750               "data inconsistency %d != %d", size_is_num, nb_dofs_row);
 1751 
 1752    bool square_matrix = false;
 1753    if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
 1754      square_matrix = true;
 1755 
 1756    
 1757    
 1758    
 1759    
 1760    
 1761    
 1762    
 1763    
 1764    
 1765    
 1766    
 1767    
 1768    
 1769    
 1770    
 1771    
 1772    
 1773    
 1774    
 1775    
 1776    
 1777    
 1778    
 1779    
 1780 
 1781    auto number_dofs = [&](auto &dofs_idx, auto &counter) {
 1783      for (auto miit_dofs_row = dofs_idx.begin();
 1784           miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
 1785        const int part = part_number[(*miit_dofs_row)->dofIdx];
 1787          const bool success = dofs_idx.modify(
 1788              miit_dofs_row,
 1789              NumeredDofEntity_part_and_indices_change(
 1790                  part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
 1791          if (!success) {
 1793                    "modification unsuccessful");
 1794          }
 1795        } else {
 1796          const bool success = dofs_idx.modify(
 1797              miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
 1798                                 part, petsc_idx[(*miit_dofs_row)->dofIdx]));
 1799          if (!success) {
 1801                    "modification unsuccessful");
 1802          }
 1803        }
 1804      }
 1806    };
 1807 
 1808    
 1809    auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
 1810        p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
 1811    int &nb_row_local_dofs = p_miit->nbLocDofsRow;
 1812    int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
 1813    nb_row_local_dofs = 0;
 1814    nb_row_ghost_dofs = 0;
 1815 
 1816    CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
 
 1817 
 1818    int &nb_col_local_dofs = p_miit->nbLocDofsCol;
 1819    int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
 1820    if (square_matrix) {
 1821      nb_col_local_dofs = nb_row_local_dofs;
 1822      nb_col_ghost_dofs = nb_row_ghost_dofs;
 1823    } else {
 1824      NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
 1825          const_cast<NumeredDofEntitysByIdx &>(
 1826              p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
 1827      nb_col_local_dofs = 0;
 1828      nb_col_ghost_dofs = 0;
 1829      CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
 
 1830    }
 1831 
 1832    CHKERR ISRestoreIndices(is_gather, &part_number);
 
 1833    CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
 
 1834    CHKERR ISDestroy(&is_num);
 
 1835    CHKERR ISDestroy(&is_gather_num);
 
 1836    CHKERR ISDestroy(&is_gather);
 
 1838    CHKERR MatPartitioningDestroy(&part);
 
 1840    CHKERR printPartitionedProblem(&*p_miit, verb);
 
 1841  } else {
 1842 
 1843    auto number_dofs = [&](auto &dof_idx, auto &counter) {
 1845      for (auto miit_dofs_row = dof_idx.begin(); miit_dofs_row != dof_idx.end();
 1846           miit_dofs_row++) {
 1847        const bool success = dof_idx.modify(
 1848            miit_dofs_row,
 1849            NumeredDofEntity_part_and_indices_change(0, counter, counter));
 1850        ++counter;
 1851        if (!success) {
 1853                  "modification unsuccessful");
 1854        }
 1855      }
 1857    };
 1858 
 1859    auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
 1860        p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
 1861    int &nb_row_local_dofs = p_miit->nbLocDofsRow;
 1862    int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
 1863    nb_row_local_dofs = 0;
 1864    nb_row_ghost_dofs = 0;
 1865 
 1866    CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
 
 1867 
 1868    bool square_matrix = false;
 1869    if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
 1870      square_matrix = true;
 1871 
 1872    int &nb_col_local_dofs = p_miit->nbLocDofsCol;
 1873    int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
 1874    if (square_matrix) {
 1875      nb_col_local_dofs = nb_row_local_dofs;
 1876      nb_col_ghost_dofs = nb_row_ghost_dofs;
 1877    } else {
 1878      NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
 1879          const_cast<NumeredDofEntitysByIdx &>(
 1880              p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
 1881      nb_col_local_dofs = 0;
 1882      nb_col_ghost_dofs = 0;
 1883      CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
 
 1884    }
 1885  }
 1886 
 1887  cOre.getBuildMoFEM() |= 1 << 4;
 1889}
 1890 
 1892    const std::string name, const std::string problem_for_rows, bool copy_rows,
 1893    const std::string problem_for_cols, bool copy_cols, int verb) {
 1897 
 1898  if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
 1900 
 1901  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
 1902 
 1903  
 1904  ProblemByName &problems_by_name =
 1905      const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
 1906  ProblemByName::iterator p_miit = problems_by_name.find(name);
 1907  if (p_miit == problems_by_name.end()) {
 1909             "problem with name < %s > not defined (top tip check spelling)",
 1910             name.c_str());
 1911  }
 1914        << p_miit->getName() << " from rows of " << problem_for_rows
 1915        << " and columns of " << problem_for_cols;
 1916 
 1917  
 1918  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
 1919  if (p_miit_row == problems_by_name.end()) {
 1921             "problem with name < %s > not defined (top tip check spelling)",
 1922             problem_for_rows.c_str());
 1923  }
 1924  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
 1925      p_miit_row->numeredRowDofsPtr;
 1926 
 1927  
 1928  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
 1929  if (p_miit_col == problems_by_name.end()) {
 1931             "problem with name < %s > not defined (top tip check spelling)",
 1932             problem_for_cols.c_str());
 1933  }
 1934  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
 1935      p_miit_col->numeredColDofsPtr;
 1936 
 1937  bool copy[] = {copy_rows, copy_cols};
 1938  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
 1939      p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
 1940 
 1941  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
 1942  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
 1943  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
 1944                                                                  dofs_col};
 1945 
 1946  for (int ss = 0; ss < 2; ss++) {
 1947 
 1948    
 1949    *nb_local_dofs[ss] = 0;
 1950    if (!copy[ss]) {
 1951 
 1952      
 1953      std::vector<int> is_local, is_new;
 1954 
 1956          copied_dofs[ss]->get<Unique_mi_tag>();
 1957      for (NumeredDofEntity_multiIndex::iterator dit =
 1958               composed_dofs[ss]->begin();
 1959           dit != composed_dofs[ss]->end(); dit++) {
 1960        NumeredDofEntityByUId::iterator diit =
 1961            dofs_by_uid.find((*dit)->getLocalUniqueId());
 1962        if (diit == dofs_by_uid.end()) {
 1963          SETERRQ(
 1965              "data inconsistency, could not find dof in composite problem");
 1966        }
 1967        int part_number = (*diit)->getPart(); 
 1968        int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
 1969        bool success;
 1970        success = composed_dofs[ss]->modify(
 1971            dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
 1972                                                           petsc_global_dof));
 1973        if (!success) {
 1975                  "modification unsuccessful");
 1976        }
 1977        if ((*dit)->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
 
 1978          success = composed_dofs[ss]->modify(
 1979              dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
 1980          if (!success) {
 1982                    "modification unsuccessful");
 1983          }
 1984          is_local.push_back(petsc_global_dof);
 1985        }
 1986      }
 1987 
 1988      AO ao;
 1989      CHKERR AOCreateMapping(m_field.
get_comm(), is_local.size(), &is_local[0],
 
 1990                             NULL, &ao);
 1991 
 1992      
 1993      is_local.resize(0);
 1994      for (NumeredDofEntity_multiIndex::iterator dit =
 1995               composed_dofs[ss]->begin();
 1996           dit != composed_dofs[ss]->end(); dit++) {
 1997        is_local.push_back((*dit)->getPetscGlobalDofIdx());
 1998      }
 1999      CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
 
 2000      int idx2 = 0;
 2001      for (NumeredDofEntity_multiIndex::iterator dit =
 2002               composed_dofs[ss]->begin();
 2003           dit != composed_dofs[ss]->end(); dit++) {
 2004        int part_number = (*dit)->getPart(); 
 2005        int petsc_global_dof = is_local[idx2++];
 2006        bool success;
 2007        success = composed_dofs[ss]->modify(
 2008            dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
 2009                                                           petsc_global_dof));
 2010        if (!success) {
 2012                  "modification unsuccessful");
 2013        }
 2014      }
 2015 
 2017 
 2018    } else {
 2019 
 2020      for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
 2021           dit != copied_dofs[ss]->end(); dit++) {
 2022        std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
 2023        p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
 2024            new NumeredDofEntity((*dit)->getDofEntityPtr())));
 2025        if (p.second) {
 2026          (*nb_dofs[ss])++;
 2027        }
 2028        int dof_idx = (*dit)->getDofIdx();
 2029        int part_number = (*dit)->getPart(); 
 2030        int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
 2032          const bool success = composed_dofs[ss]->modify(
 2033              p.first, NumeredDofEntity_part_and_all_indices_change(
 2034                           part_number, dof_idx, petsc_global_dof,
 2035                           (*nb_local_dofs[ss])++));
 2036          if (!success) {
 2038                    "modification unsuccessful");
 2039          }
 2040        } else {
 2041          const bool success = composed_dofs[ss]->modify(
 2042              p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
 2043                           part_number, dof_idx, petsc_global_dof));
 2044          if (!success) {
 2046                    "modification unsuccessful");
 2047          }
 2048        }
 2049      }
 2050    }
 2051  }
 2052 
 2053  CHKERR printPartitionedProblem(&*p_miit, verb);
 
 2054  CHKERR debugPartitionedProblem(&*p_miit, verb);
 
 2055 
 2057}
 2058 
 2060ProblemsManager::printPartitionedProblem(const Problem *problem_ptr, int verb) {
 2063 
 2065 
 2067        << problem_ptr->getName() << " Nb. local dof "
 2068        << problem_ptr->getNbLocalDofsRow() << " by "
 2069        << problem_ptr->getNbLocalDofsCol() << " nb global dofs "
 2070        << problem_ptr->getNbDofsRow() << " by " << problem_ptr->getNbDofsCol();
 2071 
 2073  }
 2074 
 2076}
 2077 
 2079ProblemsManager::debugPartitionedProblem(const Problem *problem_ptr, int verb) {
 2082 
 2083  auto save_ent = [](moab::Interface &moab, const std::string name,
 2087    CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
 
 2088    CHKERR moab.add_entities(out_meshset, &ent, 1);
 
 2089    CHKERR moab.write_file(name.c_str(), 
"VTK", 
"", &out_meshset, 1);
 
 2090    CHKERR moab.delete_entities(&out_meshset, 1);
 
 2092  };
 2093 
 2095 
 2096    typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
 2097        NumeredDofEntitysByIdx;
 2098    NumeredDofEntitysByIdx::iterator dit, hi_dit;
 2099    const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
 2100        &(problem_ptr->numeredRowDofsPtr->get<Idx_mi_tag>()),
 2101        &(problem_ptr->numeredColDofsPtr->get<Idx_mi_tag>())};
 2102 
 2103    int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
 2104    int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
 2105                              &problem_ptr->nbLocDofsCol};
 2106 
 2107    for (int ss = 0; ss < 2; ss++) {
 2108 
 2109      dit = numered_dofs_ptr[ss]->begin();
 2110      hi_dit = numered_dofs_ptr[ss]->end();
 2111      for (; dit != hi_dit; dit++) {
 2113          if ((*dit)->getPetscLocalDofIdx() < 0) {
 2114            std::ostringstream zz;
 2117                     "local dof index for %d (0-row, 1-col) not set, i.e. has "
 2118                     "negative value\n %s",
 2119                     ss, zz.str().c_str());
 2120          }
 2121          if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
 2122            std::ostringstream zz;
 2125                     "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
 2126                     zz.str().c_str());
 2127          }
 2128        } else {
 2129          if ((*dit)->getPetscGlobalDofIdx() < 0) {
 2130 
 2134                "debug_part" +
 2136                    "_negative_global_index.vtk",
 2137                ent);
 2138 
 2139            std::ostringstream zz;
 2141               << dit->get()->getBitRefLevel() << " " << **dit;
 2143                     "global dof index for %d (0-row, 1-col) row not set, i.e. "
 2144                     "has negative value\n %s",
 2145                     ss, zz.str().c_str());
 2146          }
 2147          if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
 2148            std::ostringstream zz;
 2150               << *nbdof_ptr[ss] << " " << **dit;
 2152                     "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
 2153                     zz.str().c_str());
 2154          }
 2155        }
 2156      }
 2157    }
 2158  }
 2160}
 2161 
 2162MoFEMErrorCode ProblemsManager::partitionFiniteElements(
const std::string name,
 
 2163                                                        bool part_from_moab,
 2164                                                        int low_proc,
 2165                                                        int hi_proc, int verb) {
 2170 
 2171  if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
 2173 
 2174  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
 2176 
 2177  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
 2179            "adjacencies not build");
 2180 
 2181  if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
 2183 
 2184  if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
 2186            "problem not partitioned");
 2187 
 2188  if (low_proc == -1)
 2190  if (hi_proc == -1)
 2192 
 2193  
 2194  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
 2195  auto &problems =
 2196      const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
 2197  ProblemByName::iterator p_miit = problems.find(name);
 2198  if (p_miit == problems.end()) {
 2200             "problem < %s > not found (top tip: check spelling)",
 2201             name.c_str());
 2202  }
 2203 
 2204  
 2206      *p_miit->numeredFiniteElementsPtr;
 2207 
 2208  
 2209  problem_finite_elements.clear();
 2210 
 2211  
 2212  
 2213  bool do_cols_prob = true;
 2214  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
 2215    do_cols_prob = false;
 2216  }
 2217 
 2218  auto get_good_elems = [&]() {
 2219    auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
 2220    good_elems.reserve(fe_ent_ptr->size());
 2221 
 2222    const auto prb_bit = p_miit->getBitRefLevel();
 2223    const auto prb_mask = p_miit->getBitRefLevelMask();
 2224 
 2225    
 2226    
 2227    for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
 2228 
 2229      
 2230      if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
 2231 
 2232        const auto &fe_bit = (*efit)->getBitRefLevel();
 2233 
 2234        
 2235        if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit).any())
 2236          good_elems.emplace_back(efit);
 2237      }
 2238    }
 2239 
 2240    return good_elems;
 2241  };
 2242 
 2243  auto good_elems = get_good_elems();
 2244 
 2245  auto numbered_good_elems_ptr =
 2246      boost::make_shared<std::vector<NumeredEntFiniteElement>>();
 2247  numbered_good_elems_ptr->reserve(good_elems.size());
 2248  for (auto &efit : good_elems)
 2249    numbered_good_elems_ptr->emplace_back(NumeredEntFiniteElement(*efit));
 2250 
 2251  if (!do_cols_prob) {
 2252    for (auto &fe : *numbered_good_elems_ptr) {
 2253      if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
 2254        fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
 2255      }
 2256    }
 2257  }
 2258 
 2259  if (part_from_moab) {
 2260    for (auto &fe : *numbered_good_elems_ptr) {
 2261      if (fe.getEntType() == MBENTITYSET) {
 2263      } else {
 2264        
 2265        int proc = fe.getPartProc();
 2266        if (proc == -1 && fe.getEntType() == MBVERTEX)
 2267          proc = fe.getOwnerProc();
 2268        fe.part = proc;
 2269      }
 2270    }
 2271  }
 2272 
 2273  for (auto &fe : *numbered_good_elems_ptr) {
 2274 
 2276    CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
 
 2277 
 2278    if (!part_from_moab) {
 2279      if (fe.getEntType() == MBENTITYSET) {
 2281      } else {
 2283        for (auto &dof_ptr : rows_view)
 2284          parts[dof_ptr->pArt]++;
 2285        std::vector<int>::iterator pos =
 2286            max_element(parts.begin(), parts.end());
 2287        const auto max_part = std::distance(parts.begin(), pos);
 2288        fe.part = max_part;
 2289      }
 2290    }
 2291  }
 2292 
 2293  for (auto &fe : *numbered_good_elems_ptr) {
 2294 
 2295    auto check_fields_and_dofs = [&]() {
 2296      if (!part_from_moab) {
 2297        if (fe.getBitFieldIdRow().none() && m_field.
get_comm_size() == 0) {
 
 2299              << "At least one field has to be added to element row to "
 2300                 "determine partition of finite element. Check element " +
 2301                     boost::lexical_cast<std::string>(fe.getName());
 2302        }
 2303      }
 2304 
 2305      return true;
 2306    };
 2307 
 2308    if (check_fields_and_dofs()) {
 2309      
 2310      auto p = problem_finite_elements.insert(
 2311          boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
 2312                                                     &fe));
 2313      if (!p.second)
 2315    }
 2316  }
 2317 
 2319    auto elements_on_rank =
 2320        problem_finite_elements.get<Part_mi_tag>().equal_range(
 2323        << p_miit->getName() << " nb. elems "
 2324        << std::distance(elements_on_rank.first, elements_on_rank.second);
 2326    for (auto &fe : *fe_ptr) {
 2327      auto e_range =
 2328          problem_finite_elements.get<Composite_Name_And_Part_mi_tag>()
 2329              .equal_range(
 2332          << "Element " << fe->getName() << " nb. elems "
 2333          << std::distance(e_range.first, e_range.second);
 2334    }
 2335 
 2337  }
 2338 
 2339  cOre.getBuildMoFEM() |= Core::PARTITION_FE;
 2341}
 2342 
 2343MoFEMErrorCode ProblemsManager::partitionGhostDofs(
const std::string name,
 
 2344                                                   int verb) {
 2348 
 2349  if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
 2351            "partition of problem not build");
 2352  if (!(cOre.getBuildMoFEM() & Core::PARTITION_FE))
 2354            "partitions finite elements not build");
 2355 
 2356  
 2357  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
 2358  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
 2360             name.c_str());
 2361 
 2362  
 2363  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
 2364  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
 2365  nb_row_ghost_dofs = 0;
 2366  nb_col_ghost_dofs = 0;
 2367 
 2368  
 2370 
 2372        ghost_idx_col_view;
 2373 
 2374    
 2375    auto fe_range =
 2376        p_miit->numeredFiniteElementsPtr->get<Part_mi_tag>().equal_range(
 2378 
 2379    
 2380 
 2381    struct Inserter {
 2382      using Vec = std::vector<boost::shared_ptr<NumeredDofEntity>>;
 
 2383      using It = Vec::iterator;
 2384      It operator()(Vec &dofs_view, It &hint,
 2385                    boost::shared_ptr<NumeredDofEntity> &&dof) {
 2386        dofs_view.emplace_back(dof);
 2387        return dofs_view.end();
 2388      }
 2389    };
 2390 
 2391    
 2392    std::vector<boost::shared_ptr<NumeredDofEntity>> fe_vec_view;
 2393    auto hint_r = ghost_idx_row_view.begin();
 2394    for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
 2395 
 2396      fe_vec_view.clear();
 2397      CHKERR EntFiniteElement::getDofView((*fe_ptr)->getRowFieldEnts(),
 
 2398                                          *(p_miit->getNumeredRowDofsPtr()),
 2399                                          fe_vec_view, Inserter());
 2400 
 2401      for (auto &dof_ptr : fe_vec_view) {
 2402        if (dof_ptr->getPart() != (
unsigned int)m_field.
get_comm_rank()) {
 
 2403          hint_r = ghost_idx_row_view.emplace_hint(hint_r, dof_ptr);
 2404        }
 2405      }
 2406    }
 2407 
 2408    
 2409    if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
 2410 
 2411      auto hint_c = ghost_idx_col_view.begin();
 2412      for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
 2413 
 2414        fe_vec_view.clear();
 2415        CHKERR EntFiniteElement::getDofView((*fe_ptr)->getColFieldEnts(),
 
 2416                                            *(p_miit->getNumeredColDofsPtr()),
 2417                                            fe_vec_view, Inserter());
 2418 
 2419        for (auto &dof_ptr : fe_vec_view) {
 2420          if (dof_ptr->getPart() != (
unsigned int)m_field.
get_comm_rank()) {
 
 2421            hint_c = ghost_idx_col_view.emplace_hint(hint_c, dof_ptr);
 2422          }
 2423        }
 2424      }
 2425    }
 2426 
 2427    int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
 2428    int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
 2429 
 2431        &ghost_idx_row_view, &ghost_idx_col_view};
 2433        &p_miit->numeredRowDofsPtr->get<Unique_mi_tag>(),
 2434        &p_miit->numeredColDofsPtr->get<Unique_mi_tag>()};
 2435 
 2436    int loop_size = 2;
 2437    if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
 2438      loop_size = 1;
 2439    }
 2440 
 2441    
 2442    for (int ss = 0; ss != loop_size; ++ss) {
 2443      for (auto &gid : *ghost_idx_view[ss]) {
 2444        NumeredDofEntityByUId::iterator dof =
 2445            dof_by_uid_no_const[ss]->find(gid->getLocalUniqueId());
 2446        if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
 2448                  "inconsistent data, ghost dof already set");
 2449        bool success = dof_by_uid_no_const[ss]->modify(
 2450            dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
 2451        if (PetscUnlikely(!success))
 2453                  "modification unsuccessful");
 2454        (*nb_ghost_dofs[ss])++;
 2455      }
 2456    }
 2457    if (loop_size == 1) {
 2458      (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
 2459    }
 2460  }
 2461 
 2464        << " FEs ghost dofs on problem " << p_miit->getName()
 2465        << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
 2466        << p_miit->getNbGhostDofsCol() << " Nb. local dof "
 2467        << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
 2468 
 2470  }
 2471 
 2472  cOre.getBuildMoFEM() |= Core::PARTITION_GHOST_DOFS;
 2474}
 2475 
 2477ProblemsManager::partitionGhostDofsOnDistributedMesh(const std::string name,
 2478                                                     int verb) {
 2482 
 2483  if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
 2485            "partition of problem not build");
 2486  if (!(cOre.getBuildMoFEM() & Core::PARTITION_FE))
 2488            "partitions finite elements not build");
 2489 
 2490  
 2491  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
 2492  
 2493  ProblemsByName &problems_set =
 2494      const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
 2495  ProblemsByName::iterator p_miit = problems_set.find(name);
 2496 
 2497  
 2498  
 2499  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
 
 2500                              &(p_miit->nbGhostDofsCol)};
 2501  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
 
 2502  for (int ss = 0; ss != 2; ++ss) {
 2503    (*nb_ghost_dofs[ss]) = 0;
 2504  }
 2505 
 2506  
 2508    
 2509    int loop_size = 2;
 2510    if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
 2511      loop_size = 1;
 2512    }
 2513 
 2514    typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
 2515    NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
 2516                                           p_miit->numeredColDofsPtr};
 2517 
 2518    
 2519    for (int ss = 0; ss != loop_size; ++ss) {
 2520 
 2521      
 2522      auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
 
 2523 
 2524      std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
 2525      ghost_idx_view.reserve(std::distance(
r.first, 
r.second));
 
 2526      for (; 
r.first != 
r.second; ++
r.first)
 
 2527        ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(
r.first));
 
 2528 
 2529      auto cmp = [](
auto a, 
auto b) {
 
 2530        return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
 2531      };
 2532      sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
 2533 
 2534      
 2535      for (auto gid_it : ghost_idx_view) {
 2536        bool success = numered_dofs[ss]->modify(
 2537            gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
 2538        if (!success)
 2540                  "modification unsuccessful");
 2541        ++(*nb_ghost_dofs[ss]);
 2542      }
 2543    }
 2544    if (loop_size == 1) {
 2545      (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
 2546    }
 2547  }
 2548 
 2551        << " FEs ghost dofs on problem " << p_miit->getName()
 2552        << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
 2553        << p_miit->getNbGhostDofsCol() << " Nb. local dof "
 2554        << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
 2555 
 2557  }
 2558 
 2559  cOre.getBuildMoFEM() |= Core::PARTITION_GHOST_DOFS;
 2561}
 2562 
 2563MoFEMErrorCode ProblemsManager::getFEMeshset(
const std::string prb_name,
 
 2564                                             const std::string &fe_name,
 2567  const Problem *problem_ptr;
 2570 
 2574 
 2575  auto fe_miit = fes_ptr->get<FiniteElement_name_mi_tag>().find(fe_name);
 2576  if (fe_miit != fes_ptr->get<FiniteElement_name_mi_tag>().end()) {
 2577    auto fit =
 2578        problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
 2579            EntFiniteElement::getLocalUniqueIdCalculate(
 2580                0, (*fe_miit)->getFEUId()));
 2581    auto hi_fe_it =
 2582        problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
 2583            EntFiniteElement::getLocalUniqueIdCalculate(
 2584                get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
 2585    std::vector<EntityHandle> fe_vec;
 2586    fe_vec.reserve(std::distance(fit, hi_fe_it));
 2587    for (; fit != hi_fe_it; fit++)
 2588      fe_vec.push_back(fit->get()->getEnt());
 2589    CHKERR m_field.
get_moab().add_entities(*meshset, &*fe_vec.begin(),
 
 2590                                           fe_vec.size());
 2591  }
 2592 
 2594}
 2595 
 2597ProblemsManager::getProblemElementsLayout(const std::string name,
 2598                                          const std::string &fe_name,
 2599                                          PetscLayout *layout) const {
 2601  const Problem *problem_ptr;
 2605                                                       fe_name, layout);
 2607}
 2608 
 2610 
 2611    std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view,
 2612 
 2614    const Range ents, 
int bridge_dim, 
const int lo_coeff, 
const int hi_coeff,
 
 2615    const int lo_order, 
const int hi_order, 
int verb, 
const bool debug 
 2616 
 2617) const {
 2620 
 2621  switch (rc) {
 2624    break;
 2625  default:
 2627    break;
 2628  }
 2629 
 2630  const Problem *prb_ptr;
 2632 
 2633  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
 2635  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
 2636    numered_dofs[1] = prb_ptr->numeredColDofsPtr;
 2637 
 2640      ents, bridge_dim, false, bridge_ents, moab::Interface::UNION);
 2641 
 2644  const auto nb_coeffs =
 2648 
 2649  using NumeredDofEntity_it_view_multiIndex = multi_index_container<
 2650 
 2651      NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
 2652 
 2653      >;
 2654  NumeredDofEntity_it_view_multiIndex dofs_it_view;
 2655 
 2656  for (auto pit = bridge_ents.const_pair_begin();
 2657       pit != bridge_ents.const_pair_end(); ++pit) {
 2658    auto lo = numered_dofs[rc]->get<Unique_mi_tag>().lower_bound(
 2659        DofEntity::getLoFieldEntityUId(bit_number, pit->first));
 2660    auto hi = numered_dofs[rc]->get<Unique_mi_tag>().upper_bound(
 2661        DofEntity::getHiFieldEntityUId(bit_number, pit->second));
 2662 
 2663    auto bride_type = moab.type_from_handle(pit->first);
 2664 
 2665    for (; lo != hi; ++lo) {
 2666      if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
 2667          (*lo)->getDofCoeffIdx() <= hi_coeff &&
 2668          (*lo)->getDofOrder() >= lo_order &&
 2669          (*lo)->getDofOrder() <= hi_order) {
 2670        auto ent_dof_index = (*lo)->getEntDofIdx();
 2671        auto side_it = side_dof_map.at(bride_type)
 2672                           .get<EntDofIdx_mi_tag>()
 2673                           .find(std::floor(static_cast<double>(ent_dof_index) /
 2674                                            nb_coeffs));
 2675        if (side_it !=
 2676            side_dof_map.at(bride_type).get<EntDofIdx_mi_tag>().end()) {
 2677          auto bridge_ent = (*lo)->getEnt();
 2678          auto type = side_it->type;
 
 2679          auto side = side_it->side;
 2680          auto dim = CN::Dimension(type);
 2682          CHKERR moab.side_element(bridge_ent, dim, side, side_ent);
 
 2683          if (ents.find(side_ent) != ents.end()) {
 2684            dofs_it_view.emplace_back(numered_dofs[rc]->project<0>(lo));
 2685          }
 2686        } else {
 2687#ifndef NDEBUG
 2691              << "side not found for entity " << CN::EntityTypeName(bride_type);
 2693              << "ent_dof_index / nb_coeffs "
 2694              << std::floor(static_cast<double>(ent_dof_index) / nb_coeffs);
 2696              << "side_dof_map.size() " << side_dof_map.size();
 2697#endif 
 2699                  "side not found - you will get more information in debug");
 2700        }
 2701      }
 2702    }
 2703  }
 2704 
 2706    for (auto &dof : dofs_it_view)
 2709  }
 2710 
 2711  vec_dof_view.reserve(vec_dof_view.size() + dofs_it_view.size());
 2712  for (auto dit : dofs_it_view)
 2713    vec_dof_view.push_back(*dit);
 2714 
 2716}
 2717 
 2719    const std::string problem_name, 
RowColData rc,
 
 2720    std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view, int verb,
 2722 
 2725 
 2726  switch (rc) {
 2729    break;
 2730  default:
 2732  }
 2733 
 2734  const Problem *prb_ptr;
 2736 
 2737  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
 2739  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
 2740    numered_dofs[1] = prb_ptr->numeredColDofsPtr;
 2741 
 2742  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
 2743  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
 2744  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
 2745 
 2746  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
 2747  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
 2748  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
 2749  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
 2750  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
 2751 
 2752  if (numered_dofs[rc]) {
 2753 
 2756                  "Number of DOFs in multi-index %d and to delete %d\n",
 2757                  numered_dofs[rc]->size(), vec_dof_view.size());
 2758 
 2759    
 2760    for (auto weak_dit : vec_dof_view)
 2761      if (auto dit = weak_dit.lock()) {
 2762        numered_dofs[rc]->erase(dit->getLocalUniqueId());
 2763      }
 2764 
 2767                  "Number of DOFs in multi-index after delete %d\n",
 2768                  numered_dofs[rc]->size());
 2769 
 2770    
 2771    int nb_local_dofs = 0;
 2772    int nb_ghost_dofs = 0;
 2773    for (auto dit = numered_dofs[rc]->get<PetscLocalIdx_mi_tag>().begin();
 2774         dit != numered_dofs[rc]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
 2775      if ((*dit)->getPetscLocalDofIdx() >= 0 &&
 2776          (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[rc]))
 2777        ++nb_local_dofs;
 2778      else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[rc]))
 2779        ++nb_ghost_dofs;
 2780      else
 2782                 "Impossible case. You could run problem on no distributed "
 2783                 "mesh. That is not implemented. Dof local index is %d",
 2784                 (*dit)->getPetscLocalDofIdx());
 2785    }
 2786 
 2787    
 2788    auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
 2789      const int nb_dofs = numered_dofs[rc]->size();
 2790      indices.clear();
 2791      indices.reserve(nb_dofs);
 2792      for (auto dit = numered_dofs[rc]->get<decltype(tag)>().begin();
 2793           dit != numered_dofs[rc]->get<decltype(tag)>().end(); ++dit) {
 2794        bool add = true;
 2795        if (only_local) {
 2796          if ((*dit)->getPetscLocalDofIdx() < 0 ||
 2797              (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[rc])) {
 2798            add = false;
 2799          }
 2800        }
 2801        if (add)
 2802          indices.push_back(decltype(tag)::get_index(dit));
 2803      }
 2804    };
 2805 
 2806    auto get_indices_by_uid = [&](auto tag, auto &indices) {
 2807      const int nb_dofs = numered_dofs[rc]->size();
 2808      indices.clear();
 2809      indices.reserve(nb_dofs);
 2810      for (auto dit = numered_dofs[rc]->begin(); dit != numered_dofs[rc]->end();
 2811           ++dit)
 2812        indices.push_back(decltype(tag)::get_index(dit));
 2813    };
 2814 
 2815    auto get_sub_ao = [&](auto sub_data) {
 2816      if (rc == 0) {
 2817        return sub_data->getSmartRowMap();
 2818      } else {
 2819        return sub_data->getSmartColMap();
 2820      }
 2821    };
 2822 
 2823    auto set_sub_is_and_ao = [&rc, &prb_ptr](auto sub_data, auto is, auto ao) {
 2825        sub_data->rowIs = is;
 2826        sub_data->rowMap = ao;
 2827        sub_data->colIs = is; 
 2828        sub_data->colMap = ao;
 2829      } else {
 2830        sub_data->colIs = is;
 2831        sub_data->colMap = ao;
 2832      }
 2833    };
 2834 
 2835    auto apply_symmetry = [&rc, &prb_ptr](auto sub_data) {
 2837        if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
 2838          sub_data->colIs = sub_data->getSmartRowIs();
 2839          sub_data->colMap = sub_data->getSmartRowMap();
 2840        }
 2841      }
 2842    };
 2843 
 2844    auto concatenate_dofs = [&](auto tag, auto &indices,
 2845                                const auto local_only) {
 2847      get_indices_by_tag(tag, indices, local_only);
 2848 
 2849      SmartPetscObj<AO> ao;
 2850      
 2851      
 2852      if (local_only)
 2854                             &*indices.begin(), PETSC_NULLPTR);
 2855      else
 2856        ao = 
createAOMapping(PETSC_COMM_SELF, indices.size(), &*indices.begin(),
 
 2857                             PETSC_NULLPTR);
 2858 
 2859      
 2860      if (local_only) {
 2861        if (auto sub_data = prb_ptr->getSubData()) {
 2862          
 2864                                        &*indices.begin(), PETSC_COPY_VALUES);
 2865          
 2866          
 2867          auto sub_ao = get_sub_ao(sub_data);
 2868          CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
 
 2870          
 2871          set_sub_is_and_ao(sub_data, sub_is, sub_ao);
 2872          apply_symmetry(sub_data);
 2873        } else {
 2874          
 2875          prb_ptr->getSubData() = boost::make_shared<Problem::SubProblemData>();
 2877                                        &*indices.begin(), PETSC_COPY_VALUES);
 2878          
 2879          set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
 2880          apply_symmetry(prb_ptr->getSubData());
 2881        }
 2882      }
 2883 
 2884      get_indices_by_uid(tag, indices);
 2885      CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
 
 2886 
 2888    };
 2889 
 2890    
 2891    auto set_concatenated_indices = [&]() {
 2892      std::vector<int> global_indices;
 2893      std::vector<int> local_indices;
 2895      CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, 
true);
 
 2896      CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, 
false);
 
 2897      auto gi = global_indices.begin();
 2898      auto li = local_indices.begin();
 2899      for (auto dit = numered_dofs[rc]->begin(); dit != numered_dofs[rc]->end();
 2900           ++dit) {
 2901        auto mod = NumeredDofEntity_part_and_all_indices_change(
 2902            (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
 2903        bool success = numered_dofs[rc]->modify(dit, mod);
 2904        if (!success)
 2906                  "can not set negative indices");
 2907        ++gi;
 2908        ++li;
 2909      }
 2911    };
 2912    CHKERR set_concatenated_indices();
 
 2913 
 2914    MPI_Allreduce(&nb_local_dofs, nbdof_ptr[rc], 1, MPI_INT, MPI_SUM,
 2916    *(local_nbdof_ptr[rc]) = nb_local_dofs;
 2917    *(ghost_nbdof_ptr[rc]) = nb_ghost_dofs;
 2918 
 2920      for (auto dof : (*numered_dofs[rc])) {
 2921        if (dof->getPetscGlobalDofIdx() < 0) {
 2923                  "Negative global idx");
 2924        }
 2925        if (dof->getPetscLocalDofIdx() < 0) {
 2927                  "Negative local idx");
 2928        }
 2929      }
 2930 
 2931  } else {
 2932 
 2933    *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
 2934    *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
 2935    *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
 2936  }
 2937 
 2940        "WORLD", Sev::inform,
 2941        "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
 2942        prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
 2943        prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
 2945                "Removed DOFs from problem %s dofs [ %d / %d  "
 2946                "(before %d / %d) local, %d / %d (before %d / %d)]",
 2947                prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
 2948                prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
 2949                nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
 2950                prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
 2951                nb_init_ghost_col_dofs);
 2953  }
 2954 
 2956}
 2957 
 2959    const std::string problem_name, 
const std::string 
field_name,
 
 2960    const Range ents, 
const int lo_coeff, 
const int hi_coeff,
 
 2961    const int lo_order, 
const int hi_order, 
int verb, 
const bool debug) {
 
 2962 
 2965 
 2966  const Problem *prb_ptr;
 2968 
 2969  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
 2971  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
 2972    numered_dofs[1] = prb_ptr->numeredColDofsPtr;
 2973 
 2974  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
 2975  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
 2976  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
 2977 
 2978  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
 2979  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
 2980  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
 2981  
 2982  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
 2983  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
 2984 
 2985  for (int s = 0; s != 2; ++s)
 2986    if (numered_dofs[s]) {
 2987 
 2988      using NumeredDofEntity_it_view_multiIndex = multi_index_container<
 2989 
 2990          NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
 2991 
 2992          >;
 2993 
 2995      NumeredDofEntity_it_view_multiIndex dofs_it_view;
 2996 
 2997      
 2998      for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
 2999           ++pit) {
 3000        auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
 3001            DofEntity::getLoFieldEntityUId(bit_number, pit->first));
 3002        auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
 3003            DofEntity::getHiFieldEntityUId(bit_number, pit->second));
 3004 
 3005        for (; lo != hi; ++lo)
 3006          if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
 3007              (*lo)->getDofCoeffIdx() <= hi_coeff &&
 3008              (*lo)->getDofOrder() >= lo_order &&
 3009              (*lo)->getDofOrder() <= hi_order)
 3010            dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
 3011      }
 3012 
 3014        for (auto &dof : dofs_it_view)
 3017      }
 3018 
 3019      
 3020      std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
 3021      dofs_weak_view.reserve(dofs_it_view.size());
 3022      for (auto dit : dofs_it_view)
 3023        dofs_weak_view.push_back(*dit);
 3024 
 3027                    "Number of DOFs in multi-index %d and to delete %d\n",
 3028                    numered_dofs[s]->size(), dofs_it_view.size());
 3029 
 3030      
 3031      for (auto weak_dit : dofs_weak_view)
 3032        if (auto dit = weak_dit.lock()) {
 3033          numered_dofs[s]->erase(dit->getLocalUniqueId());
 3034        }
 3035 
 3038                    "Number of DOFs in multi-index after delete %d\n",
 3039                    numered_dofs[s]->size());
 3040 
 3041      
 3042      int nb_local_dofs = 0;
 3043      int nb_ghost_dofs = 0;
 3044      for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
 3045           dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
 3046        if ((*dit)->getPetscLocalDofIdx() >= 0 &&
 3047            (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
 3048          ++nb_local_dofs;
 3049        else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
 3050          ++nb_ghost_dofs;
 3051        else
 3053                   "Impossible case. You could run problem on no distributed "
 3054                   "mesh. That is not implemented. Dof local index is %d",
 3055                   (*dit)->getPetscLocalDofIdx());
 3056      }
 3057 
 3058      
 3059      auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
 3060        const int nb_dofs = numered_dofs[s]->size();
 3061        indices.clear();
 3062        indices.reserve(nb_dofs);
 3063        for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
 3064             dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
 3065          bool add = true;
 3066          if (only_local) {
 3067            if ((*dit)->getPetscLocalDofIdx() < 0 ||
 3068                (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
 3069              add = false;
 3070            }
 3071          }
 3072          if (add)
 3073            indices.push_back(decltype(tag)::get_index(dit));
 3074        }
 3075      };
 3076 
 3077      auto get_indices_by_uid = [&](auto tag, auto &indices) {
 3078        const int nb_dofs = numered_dofs[s]->size();
 3079        indices.clear();
 3080        indices.reserve(nb_dofs);
 3081        for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
 3082             ++dit)
 3083          indices.push_back(decltype(tag)::get_index(dit));
 3084      };
 3085 
 3086      auto get_sub_ao = [&](auto sub_data) {
 3087        if (s == 0) {
 3088          return sub_data->getSmartRowMap();
 3089        } else {
 3090          return sub_data->getSmartColMap();
 3091        }
 3092      };
 3093 
 3094      auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
 3095        if (s == 0) {
 3096          sub_data->rowIs = is;
 3097          sub_data->rowMap = ao;
 3098          sub_data->colIs = is; 
 3099          sub_data->colMap = ao;
 3100        } else {
 3101          sub_data->colIs = is;
 3102          sub_data->colMap = ao;
 3103        }
 3104      };
 3105 
 3106      auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
 3107        if (s == 0) {
 3108          if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
 3109            sub_data->colIs = sub_data->getSmartRowIs();
 3110            sub_data->colMap = sub_data->getSmartRowMap();
 3111          }
 3112        }
 3113      };
 3114 
 3115      auto concatenate_dofs = [&](auto tag, auto &indices,
 3116                                  const auto local_only) {
 3118        get_indices_by_tag(tag, indices, local_only);
 3119 
 3120        SmartPetscObj<AO> ao;
 3121        
 3122        
 3123        if (local_only)
 3125                               &*indices.begin(), PETSC_NULLPTR);
 3126        else
 3128                               &*indices.begin(), PETSC_NULLPTR);
 3129 
 3130        
 3131        if (local_only) {
 3132          if (auto sub_data = prb_ptr->getSubData()) {
 3133            
 3135                                          &*indices.begin(), PETSC_COPY_VALUES);
 3136            
 3137            
 3138            auto sub_ao = get_sub_ao(sub_data);
 3139            CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
 
 3141            
 3142            set_sub_is_and_ao(sub_data, sub_is, sub_ao);
 3143            apply_symmetry(sub_data);
 3144          } else {
 3145            
 3146            prb_ptr->getSubData() =
 3147                boost::make_shared<Problem::SubProblemData>();
 3149                                          &*indices.begin(), PETSC_COPY_VALUES);
 3150            
 3151            set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
 3152            apply_symmetry(prb_ptr->getSubData());
 3153          }
 3154        }
 3155 
 3156        get_indices_by_uid(tag, indices);
 3157        CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
 
 3158 
 3160      };
 3161 
 3162      
 3163      auto set_concatenated_indices = [&]() {
 3164        std::vector<int> global_indices;
 3165        std::vector<int> local_indices;
 3167        CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, 
true);
 
 3168        CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, 
false);
 
 3169        auto gi = global_indices.begin();
 3170        auto li = local_indices.begin();
 3171        for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
 3172             ++dit) {
 3173          auto mod = NumeredDofEntity_part_and_all_indices_change(
 3174              (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
 3175          bool success = numered_dofs[s]->modify(dit, mod);
 3176          if (!success)
 3178                    "can not set negative indices");
 3179          ++gi;
 3180          ++li;
 3181        }
 3183      };
 3184      CHKERR set_concatenated_indices();
 
 3185 
 3186      MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
 3188      *(local_nbdof_ptr[s]) = nb_local_dofs;
 3189      *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
 3190 
 3192        for (auto dof : (*numered_dofs[s])) {
 3193          if (dof->getPetscGlobalDofIdx() < 0) {
 3195                    "Negative global idx");
 3196          }
 3197          if (dof->getPetscLocalDofIdx() < 0) {
 3199                    "Negative local idx");
 3200          }
 3201        }
 3202 
 3203    } else {
 3204 
 3205      *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
 3206      *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
 3207      *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
 3208    }
 3209 
 3212        "WORLD", Sev::inform,
 3213        "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
 3214        prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
 3215        prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
 3217                "Removed DOFs from problem %s dofs [ %d / %d  "
 3218                "(before %d / %d) local, %d / %d (before %d / %d)]",
 3219                prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
 3220                prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
 3221                nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
 3222                prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
 3223                nb_init_ghost_col_dofs);
 3225  }
 3226 
 3228}
 3229 
 3230MoFEMErrorCode ProblemsManager::removeDofsOnEntitiesNotDistributed(
 
 3231    const std::string problem_name, 
const std::string 
field_name,
 
 3232    const Range ents, 
const int lo_coeff, 
const int hi_coeff,
 
 3233    const int lo_order, 
const int hi_order, 
int verb, 
const bool debug) {
 
 3234 
 3237 
 3238  const Problem *prb_ptr;
 3240 
 3241  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
 3243  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
 3244    numered_dofs[1] = prb_ptr->numeredColDofsPtr;
 3245 
 3246  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
 3247  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
 3248  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
 3249 
 3250  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
 3251  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
 3252  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
 3253  
 3254  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
 3255  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
 3256 
 3257  const std::array<int, 2> nb_init_dofs = {nb_init_row_dofs, nb_init_col_dofs};
 3258 
 3259  for (int s = 0; s != 2; ++s)
 3260    if (numered_dofs[s]) {
 3261 
 3262      typedef multi_index_container<
 3263 
 3264          NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
 3265 
 3266          >
 3267          NumeredDofEntity_it_view_multiIndex;
 3268 
 3270      NumeredDofEntity_it_view_multiIndex dofs_it_view;
 3271 
 3272      
 3273      for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
 3274           ++pit) {
 3275        auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
 3276            DofEntity::getLoFieldEntityUId(bit_number, pit->first));
 3277        auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
 3278            DofEntity::getHiFieldEntityUId(bit_number, pit->second));
 3279 
 3280        for (; lo != hi; ++lo)
 3281          if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
 3282              (*lo)->getDofCoeffIdx() <= hi_coeff &&
 3283              (*lo)->getDofOrder() >= lo_order &&
 3284              (*lo)->getDofOrder() <= hi_order)
 3285            dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
 3286      }
 3287 
 3289        for (auto &dof : dofs_it_view)
 3291      }
 3292 
 3293      
 3294      auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
 3295      for (auto dit : dofs_it_view) {
 3296        bool success = numered_dofs[s]->modify(dit, mod);
 3297        if (!success)
 3299                  "can not set negative indices");
 3300      }
 3301 
 3302      
 3303      std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
 3304      dosf_weak_view.reserve(dofs_it_view.size());
 3305      for (auto dit : dofs_it_view)
 3306        dosf_weak_view.push_back(*dit);
 3307 
 3310                    "Number of DOFs in multi-index %d and to delete %d\n",
 3311                    numered_dofs[s]->size(), dofs_it_view.size());
 3312 
 3313      
 3314      for (auto weak_dit : dosf_weak_view)
 3315        if (auto dit = weak_dit.lock()) {
 3316          numered_dofs[s]->erase(dit->getLocalUniqueId());
 3317        }
 3318 
 3321                    "Number of DOFs in multi-index after delete %d\n",
 3322                    numered_dofs[s]->size());
 3323 
 3324      
 3325      int nb_global_dof = 0;
 3326      int nb_local_dofs = 0;
 3327      int nb_ghost_dofs = 0;
 3328 
 3329      for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
 3330           ++dit) {
 3331 
 3332        if ((*dit)->getDofIdx() >= 0) {
 3333 
 3334          if ((*dit)->getPetscLocalDofIdx() >= 0 &&
 3335              (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
 3336            ++nb_local_dofs;
 3337          else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
 3338            ++nb_ghost_dofs;
 3339 
 3340          ++nb_global_dof;
 3341        }
 3342      }
 3343 
 3345        MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
 3347        if (*(nbdof_ptr[s]) != nb_global_dof)
 3349                   "Number of local DOFs do not add up %d != %d",
 3350                   *(nbdof_ptr[s]), nb_global_dof);
 3351      }
 3352 
 3353      *(nbdof_ptr[s]) = nb_global_dof;
 3354      *(local_nbdof_ptr[s]) = nb_local_dofs;
 3355      *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
 3356 
 3357      
 3358      auto get_indices_by_tag = [&](auto tag) {
 3359        std::vector<int> indices;
 3360        indices.resize(nb_init_dofs[s], -1);
 3361        for (auto dit = numered_dofs[s]->get<Idx_mi_tag>().lower_bound(0);
 3362             dit != numered_dofs[s]->get<Idx_mi_tag>().end(); ++dit) {
 3363          indices[(*dit)->getDofIdx()] = decltype(tag)::get_index(dit);
 3364        }
 3365        return indices;
 3366      };
 3367 
 3368      auto renumber = [&](auto tag, auto &indices) {
 3370        int idx = 0;
 3371        for (auto dit = numered_dofs[s]->get<decltype(tag)>().lower_bound(0);
 3372             dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
 3373          indices[(*dit)->getDofIdx()] = idx++;
 3374        }
 3376      };
 3377 
 3378      auto get_sub_ao = [&](auto sub_data) {
 3379        if (s == 0) {
 3380          return sub_data->getSmartRowMap();
 3381        } else {
 3382          return sub_data->getSmartColMap();
 3383        }
 3384      };
 3385 
 3386      auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
 3387        if (s == 0) {
 3388          sub_data->rowIs = is;
 3389          sub_data->rowMap = ao;
 3390          sub_data->colIs = is;
 3391          sub_data->colMap = ao;
 3392        } else {
 3393          sub_data->colIs = is;
 3394          sub_data->colMap = ao;
 3395        }
 3396      };
 3397 
 3398      auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
 3399        if (s == 0) {
 3400          if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
 3401            sub_data->colIs = sub_data->getSmartRowIs();
 3402            sub_data->colMap = sub_data->getSmartRowMap();
 3403          }
 3404        }
 3405      };
 3406 
 3407      auto set_sub_data = [&](auto &indices) {
 3409        if (auto sub_data = prb_ptr->getSubData()) {
 3410          
 3412                                        &*indices.begin(), PETSC_COPY_VALUES);
 3413          
 3414          
 3415          auto sub_ao = get_sub_ao(sub_data);
 3416          CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
 
 3418          
 3419          set_sub_is_and_ao(sub_data, sub_is, sub_ao);
 3420          apply_symmetry(sub_data);
 3421        } else {
 3422          prb_ptr->getSubData() = boost::make_shared<Problem::SubProblemData>();
 3424                                        &*indices.begin(), PETSC_COPY_VALUES);
 3426          
 3427          set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, sub_ao);
 3428          apply_symmetry(prb_ptr->getSubData());
 3429        }
 3431      };
 3432 
 3433      auto global_indices = get_indices_by_tag(PetscGlobalIdx_mi_tag());
 3434      auto local_indices = get_indices_by_tag(PetscLocalIdx_mi_tag());
 3435      CHKERR set_sub_data(global_indices);
 
 3436      CHKERR renumber(PetscGlobalIdx_mi_tag(), global_indices);
 
 3437      CHKERR renumber(PetscLocalIdx_mi_tag(), local_indices);
 
 3438 
 3440      for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
 3441           ++dit) {
 3442        auto idx = (*dit)->getDofIdx();
 3443        if (idx >= 0) {
 3444          auto mod = NumeredDofEntity_part_and_all_indices_change(
 3445              (*dit)->getPart(), 
i++, global_indices[idx], local_indices[idx]);
 
 3446          bool success = numered_dofs[s]->modify(dit, mod);
 3447          if (!success)
 3449                    "can not set negative indices");
 3450        } else {
 3451          auto mod = NumeredDofEntity_part_and_all_indices_change(
 3452              (*dit)->getPart(), -1, -1, -1);
 3453          bool success = numered_dofs[s]->modify(dit, mod);
 3454          if (!success)
 3456                    "can not set negative indices");
 3457        }
 3458      };
 3459 
 3461        for (auto dof : (*numered_dofs[s])) {
 3462          if (dof->getDofIdx() >= 0 && dof->getPetscGlobalDofIdx() < 0) {
 3464                    "Negative global idx");
 3465          }
 3466        }
 3467      }
 3468 
 3469    } else {
 3470 
 3471      *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
 3472      *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
 3473      *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
 3474    }
 3475 
 3478 
 3481        "WORLD", Sev::inform,
 3482        "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
 3483        prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
 3484        prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
 3486                "Removed DOFs from problem %s dofs [ %d / %d  "
 3487                "(before %d / %d) local, %d / %d (before %d / %d)]",
 3488                prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
 3489                prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
 3490                nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
 3491                prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
 3492                nb_init_ghost_col_dofs);
 3494  }
 3496}
 3497 
 3499    const std::string problem_name, 
const std::string 
field_name,
 
 3500    const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
 3501    Range *ents_ptr, 
const int lo_coeff, 
const int hi_coeff, 
const int lo_order,
 
 3502    const int hi_order, 
int verb, 
const bool debug) {
 
 3505 
 3506  auto bit_manager = m_field.
getInterface<BitRefManager>();
 
 3507 
 3509  if (ents_ptr) {
 3510    ents = *ents_ptr;
 3511    CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
 
 3512                                                 ents, verb);
 3513  } else {
 3514    CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
 
 3515                                              verb);
 3516  }
 3517 
 3519                              hi_coeff, lo_order, hi_order, verb, 
debug);
 
 3520 
 3522}
 3523 
 3524MoFEMErrorCode ProblemsManager::removeDofsOnEntitiesNotDistributed(
 
 3525    const std::string problem_name, 
const std::string 
field_name,
 
 3526    const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
 3527    Range *ents_ptr, 
const int lo_coeff, 
const int hi_coeff, 
const int lo_order,
 
 3528    const int hi_order, 
int verb, 
const bool debug) {
 
 3531 
 3532  auto bit_manager = m_field.
getInterface<BitRefManager>();
 
 3533 
 3535  if (ents_ptr) {
 3536    ents = *ents_ptr;
 3537    CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
 
 3538                                                 ents, verb);
 3539  } else {
 3540    CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
 
 3541                                              verb);
 3542  }
 3543 
 3545                                            lo_coeff, hi_coeff, lo_order,
 3546                                            hi_order, verb, 
debug);
 
 3547 
 3549}
 3550 
 3552ProblemsManager::markDofs(
const std::string problem_name, 
RowColData rc,
 
 3553                          const enum MarkOP op, 
const Range ents,
 
 3554                          std::vector<unsigned char> &
marker)
 const {
 
 3555 
 3557  const Problem *problem_ptr;
 3560  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
 3561  switch (rc) {
 3564    break;
 3566    dofs = problem_ptr->getNumeredColDofsPtr();
 3567  default:
 3569  }
 3570  marker.resize(dofs->size(), 0);
 
 3571  std::vector<unsigned char> marker_tmp;
 3572 
 3573  switch (op) {
 3574  case MarkOP::OR:
 3575    for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
 3576      auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
 3577      auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
 3578      for (; lo != hi; ++lo) {
 3579        if ((*lo)->getPetscLocalDofIdx() < 0)
 3580          continue;
 3581        marker[(*lo)->getPetscLocalDofIdx()] |= 1;
 
 3582      }
 3583    }
 3584    break;
 3585  case MarkOP::AND:
 3586    marker_tmp.resize(dofs->size(), 0);
 3587    for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
 3588      auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
 3589      auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
 3590      for (; lo != hi; ++lo) {
 3591        if ((*lo)->getPetscLocalDofIdx() < 0)
 3592          continue;
 3593        marker_tmp[(*lo)->getPetscLocalDofIdx()] = 1;
 3594      }
 3595    }
 3596    for (
int i = 0; 
i != 
marker.size(); ++
i) {
 
 3598    }
 3599    break;
 3600  }
 3601 
 3603}
 3604 
 3607    const int lo, const int hi, const enum ProblemsManager::MarkOP op,
 3608    const unsigned char c, std::vector<unsigned char> &
marker)
 const {
 
 3609 
 3611  const Problem *problem_ptr;
 3614  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
 3615  switch (rc) {
 3618    break;
 3620    dofs = problem_ptr->getNumeredColDofsPtr();
 3621  default:
 3623  }
 3624  marker.resize(dofs->size(), 0);
 
 3625 
 3626  auto dof_lo = dofs->get<Unique_mi_tag>().lower_bound(
 3627      FieldEntity::getLoBitNumberUId(m_field.get_field_bit_number(
field_name)));
 
 3628  auto dof_hi = dofs->get<Unique_mi_tag>().upper_bound(
 3629      FieldEntity::getHiBitNumberUId(m_field.get_field_bit_number(
field_name)));
 
 3630 
 3631  switch (op) {
 3632  case MarkOP::OR:
 3633    for (; dof_lo != dof_hi; ++dof_lo) {
 3634      if ((*dof_lo)->getPetscLocalDofIdx() < 0)
 3635        continue;
 3636      if ((*dof_lo)->getDofCoeffIdx() >= lo &&
 3637          (*dof_lo)->getDofCoeffIdx() <= hi)
 3638        marker[(*dof_lo)->getPetscLocalDofIdx()] |= 
c;
 
 3639    }
 3640    break;
 3641  case MarkOP::AND:
 3642    for (; dof_lo != dof_hi; ++dof_lo) {
 3643      if ((*dof_lo)->getPetscLocalDofIdx() < 0)
 3644        continue;
 3645      if ((*dof_lo)->getDofCoeffIdx() >= lo &&
 3646          (*dof_lo)->getDofCoeffIdx() <= hi)
 3647        marker[(*dof_lo)->getPetscLocalDofIdx()] &= 
c;
 
 3648    }
 3649    break;
 3650  }
 3651 
 3653}
 3654 
 3656 
 3657    const std::string problem_name, 
RowColData rc,
 
 3658    std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view,
 3659    const enum MarkOP op, std::vector<unsigned char> &
marker 
 3660 
 3661) const {
 3664 
 3665  const Problem *problem_ptr;
 3667  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
 3668  switch (rc) {
 3671    break;
 3673    dofs = problem_ptr->getNumeredColDofsPtr();
 3674    break;
 3675  default:
 3677  }
 3678  marker.resize(dofs->size(), 0);
 
 3679 
 3680  for (auto &dof : vec_dof_view) {
 3681    if (auto dof_ptr = dof.lock()) {
 3682      if (dof_ptr->getPetscLocalDofIdx() < 0)
 3683        continue;
 3684      if (op == MarkOP::OR)
 3685        marker[dof_ptr->getPetscLocalDofIdx()] |= 1;
 
 3686      else
 3687        marker[dof_ptr->getPetscLocalDofIdx()] &= 1;
 
 3688    }
 3689  }
 3690 
 3692}
 3693 
 3695ProblemsManager::addFieldToEmptyFieldBlocks(const std::string problem_name,
 3696                                            const std::string row_field,
 3697                                            const std::string col_field) const {
 3698 
 3701 
 3702  const auto problem_ptr = m_field.
get_problem(problem_name);
 
 3703  auto get_field_id = [&](
const std::string 
field_name) {
 
 3705  };
 3706  const auto row_id = get_field_id(row_field);
 3707  const auto col_id = get_field_id(col_field);
 3708 
 3709  problem_ptr->addFieldToEmptyFieldBlocks(
BlockFieldPair(row_id, col_id));
 
 3710 
 3712}
 3713 
 3714} 
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define ProblemManagerFunctionBegin
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FieldEntity, UId, &FieldEntity::getGlobalUniqueId > > > > FieldEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getLocalUniqueId > >, ordered_non_unique< const_mem_fun< DofEntity, char, &DofEntity::getActive > > > > DofEntity_multiIndex_active_view
multi-index view on DofEntity activity
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getGlobalUniqueId > > > > DofEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
auto marker
set bit to marker
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
auto createISGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
Creates a data structure for an index set containing a list of integers.
auto createAOMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[])
Creates an application mapping using two integer arrays.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
MoFEM::LogManager::SeverityLevel Sev
Problem::BlockFieldPair BlockFieldPair
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
DeprecatedCoreInterface Interface
constexpr auto field_name
FTensor::Index< 'm', 3 > m
virtual int get_comm_size() const =0
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
virtual std::string get_field_name(const BitFieldId id) const =0
get field name from id
virtual int get_comm_rank() const =0
Deprecated interface functions.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients per DOF.
std::map< int, BaseFunction::DofsSideMap > & getDofSideMap() const
Get DOF side mapping for broken spaces.
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
BitProblemId getId() const
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.