v0.10.0
Files | Functions
ProblemsManager

Adding and managing problems. More...

Collaboration diagram for ProblemsManager:

Files

file  ProblemsManager.cpp
 Managing complexities for problem.
 
file  ProblemsManager.hpp
 Interface managing problems.
 

Functions

MoFEMErrorCode MoFEM::ProblemsManager::partitionMesh (const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
 Set partition tag to each finite element in the problem. More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (Problem *problem_ptr, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective) More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (Problem *problem_ptr, const bool square_matrix=true, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective) More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildSubProblem (const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, std::pair< EntityType, EntityType >> *entityMapRow=nullptr, const map< std::string, std::pair< EntityType, EntityType >> *entityMapCol=nullptr, int verb=VERBOSE)
 build sub problem More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildCompsedProblem (const std::string out_name, const std::vector< std::string > add_row_problems, const std::vector< std::string > add_col_problems, const bool square_matrix=true, int verb=1)
 build composite problem More...
 
MoFEMErrorCode MoFEM::ProblemsManager::inheritPartition (const std::string name, const std::string problem_for_rows, bool copy_rows, const std::string problem_for_cols, bool copy_cols, int verb=VERBOSE)
 build indexing and partition problem inheriting indexing and partitioning from two other problems More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionSimpleProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs (collective) More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionFiniteElements (const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
 partition finite elements More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofs (const std::string name, int verb=VERBOSE)
 determine ghost nodes More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh (const std::string name, int verb=VERBOSE)
 determine ghost nodes on distributed meshes More...
 
MoFEMErrorCode MoFEM::ProblemsManager::getFEMeshset (const std::string prb_name, const std::string fe_name, EntityHandle *meshset) const
 create add entities of finite element in the problem More...
 
MoFEMErrorCode MoFEM::ProblemsManager::getProblemElementsLayout (const std::string name, const std::string fe_name, PetscLayout *layout) const
 Get layout of elements in the problem. More...
 
MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities (const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, int verb=VERBOSE, const bool debug=false)
 Remove DOFs from problem. More...
 

Detailed Description

Adding and managing problems.

Function Documentation

◆ buildCompsedProblem()

MoFEMErrorCode MoFEM::ProblemsManager::buildCompsedProblem ( const std::string  out_name,
const std::vector< std::string >  add_row_problems,
const std::vector< std::string >  add_col_problems,
const bool  square_matrix = true,
int  verb = 1 
)

build composite problem

Parameters
out_namename of build problem
add_row_problemsvector of add row problems
add_col_problemsvector of add col problems
square_matrixtrue if structurally squared matrix
verbverbosity level
Returns
error code

Definition at line 1521 of file ProblemsManager.cpp.

1524  {
1526  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1527  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1528  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1529  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1530  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1531  MoFEM::Interface &m_field = cOre;
1532  auto problems_ptr = m_field.get_problems();
1534 
1535  CHKERR m_field.clear_problem(out_name);
1536  // get reference to all problems
1538  ProblemByName &problems_by_name =
1539  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1540 
1541  // Get iterators to out problem, i.e. build problem
1542  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1543  if (out_problem_it == problems_by_name.end()) {
1544  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1545  "problem with name < %s > not defined (top tip check spelling)",
1546  out_name.c_str());
1547  }
1548  // Make data structure for composed-problem data
1549  out_problem_it->composedProblemsData =
1550  boost::make_shared<ComposedProblemsData>();
1551  boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1552  out_problem_it->getComposedProblemsData();
1553 
1554  const std::vector<std::string> *add_prb[] = {&add_row_problems,
1555  &add_col_problems};
1556  std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1557  &cmp_prb_data->colProblemsAdd};
1558  std::vector<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1559 
1560  // Get local indices counter
1561  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1562  &out_problem_it->nbLocDofsCol};
1563  // Get global indices counter
1564  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1565 
1566  // Set number of ghost nodes to zero
1567  {
1568  out_problem_it->nbDofsRow = 0;
1569  out_problem_it->nbDofsCol = 0;
1570  out_problem_it->nbLocDofsRow = 0;
1571  out_problem_it->nbLocDofsCol = 0;
1572  out_problem_it->nbGhostDofsRow = 0;
1573  out_problem_it->nbGhostDofsCol = 0;
1574  }
1575  int nb_dofs_reserve[] = {0, 0};
1576 
1577  // Loop over rows and columns in the main problem and sub-problems
1578  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1579  // cerr << "SS " << ss << endl;
1580  // cerr << add_prb[ss]->size() << endl;
1581  // cerr << add_prb_ptr[ss]->size() << endl;
1582  add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1583  add_prb_is[ss]->reserve(add_prb[ss]->size());
1584  for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1585  vit != add_prb[ss]->end(); vit++) {
1586  ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1587  if (prb_it == problems_by_name.end()) {
1588  SETERRQ1(
1589  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1590  "problem with name < %s > not defined (top tip check spelling)",
1591  vit->c_str());
1592  }
1593  add_prb_ptr[ss]->push_back(&*prb_it);
1594  // set number of dofs on rows and columns
1595  if (ss == 0) {
1596  // row
1597  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1598  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1599  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1600  } else {
1601  // column
1602  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1603  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1604  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1605  }
1606  }
1607  }
1608  // if squre problem, rows and columns are the same
1609  if (square_matrix) {
1610  add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1611  add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1612  out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1613  *nb_dofs[1] = *nb_dofs[0];
1614  *nb_local_dofs[1] = *nb_local_dofs[0];
1615  }
1616 
1617  // reserve memory for dofs
1618  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1619  // Reserve memory
1620  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1621  dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1622  dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1623  if (!ss)
1624  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1625  else
1626  out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1627  }
1628 
1629  // Push back DOFs
1630  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1631  NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1632  dit,
1633  hi_dit;
1634  int shift_glob = 0;
1635  int shift_loc = 0;
1636  for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1637  PetscInt *dofs_out_idx_ptr;
1638  int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1639  CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1640  if (ss == 0) {
1641  dit = (*add_prb_ptr[ss])[pp]
1642  ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1643  .begin();
1644  hi_dit = (*add_prb_ptr[ss])[pp]
1645  ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1646  .end();
1647  } else {
1648  dit = (*add_prb_ptr[ss])[pp]
1649  ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1650  .begin();
1651  hi_dit = (*add_prb_ptr[ss])[pp]
1652  ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1653  .end();
1654  }
1655  int is_nb = 0;
1656  for (; dit != hi_dit; dit++) {
1657  BitRefLevel prb_bit = out_problem_it->getBitRefLevel();
1658  BitRefLevel prb_mask = out_problem_it->getMaskBitRefLevel();
1659  BitRefLevel dof_bit = dit->get()->getBitRefLevel();
1660  if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1661  continue;
1662  const int rank = m_field.get_comm_rank();
1663  const int part = dit->get()->getPart();
1664  const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1665  const int loc_idx =
1666  (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1667  : -1;
1668  dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1669  glob_idx, loc_idx, part);
1670  if (part == rank) {
1671  dofs_out_idx_ptr[is_nb++] = glob_idx;
1672  }
1673  }
1674  if (is_nb > nb_local_dofs) {
1675  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1676  "Data inconsistency");
1677  }
1678  IS is;
1679  CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1680  PETSC_OWN_POINTER, &is);
1681  (*add_prb_is[ss]).push_back(is);
1682  if (ss == 0) {
1683  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1684  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1685  } else {
1686  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1687  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1688  }
1689  if (square_matrix) {
1690  (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1691  (*add_prb_is[1]).push_back(is);
1692  CHKERR PetscObjectReference((PetscObject)is);
1693  }
1694  }
1695  }
1696 
1697  if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1698  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1699  }
1700 
1701  // Insert DOFs to problem multi-index
1702  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1703  auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1704  : out_problem_it->numeredColDofsPtr->end();
1705  for (auto &v : *dofs_array[ss])
1706  hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1707  hint, dofs_array[ss], &v)
1708  : out_problem_it->numeredColDofsPtr->emplace_hint(
1709  hint, dofs_array[ss], &v);
1710  }
1711 
1712  // Compress DOFs
1713  *nb_dofs[0] = 0;
1714  *nb_dofs[1] = 0;
1715  *nb_local_dofs[0] = 0;
1716  *nb_local_dofs[1] = 0;
1717  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1718 
1719  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1720  if (ss == 0) {
1721  dofs_ptr = out_problem_it->numeredRowDofsPtr;
1722  } else {
1723  dofs_ptr = out_problem_it->numeredColDofsPtr;
1724  }
1725  NumeredDofEntityByUId::iterator dit, hi_dit;
1726  dit = dofs_ptr->get<Unique_mi_tag>().begin();
1727  hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1728  std::vector<int> idx;
1729  idx.reserve(std::distance(dit, hi_dit));
1730  // set dofs in order entity and dof number on entity
1731  for (; dit != hi_dit; dit++) {
1732  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1733  bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1734  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1735  if (!success) {
1736  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1737  "modification unsuccessful");
1738  }
1739  idx.push_back(dit->get()->getPetscGlobalDofIdx());
1740  } else {
1741  if (dit->get()->getPetscLocalDofIdx() != -1) {
1742  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1743  "local index should be negative");
1744  }
1745  }
1746  }
1747  if (square_matrix) {
1748  *nb_local_dofs[1] = *nb_local_dofs[0];
1749  }
1750 
1751  // set new dofs mapping
1752  IS is;
1753  CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1754  PETSC_USE_POINTER, &is);
1755  CHKERR ISGetSize(is, nb_dofs[ss]);
1756  if (square_matrix) {
1757  *nb_dofs[1] = *nb_dofs[0];
1758  }
1759 
1760  AO ao;
1761  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1762  for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1763  CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1764 
1765  // Set DOFs numeration
1766  {
1767  std::vector<int> idx_new;
1768  idx_new.reserve(dofs_ptr->size());
1769  for (NumeredDofEntityByUId::iterator dit =
1770  dofs_ptr->get<Unique_mi_tag>().begin();
1771  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1772  idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1773  }
1774  // set new global dofs numeration
1775  IS is_new;
1776  CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1777  &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1778  CHKERR AOApplicationToPetscIS(ao, is_new);
1779  // set global indices to multi-index
1780  std::vector<int>::iterator vit = idx_new.begin();
1781  for (NumeredDofEntityByUId::iterator dit =
1782  dofs_ptr->get<Unique_mi_tag>().begin();
1783  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1784  bool success =
1785  dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1786  dit->get()->getPart(), *(vit++)));
1787  if (!success) {
1788  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1789  "modification unsuccessful");
1790  }
1791  }
1792  CHKERR ISDestroy(&is_new);
1793  }
1794  CHKERR ISDestroy(&is);
1795  CHKERR AODestroy(&ao);
1796  }
1797 
1798  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1799  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1800 
1801  // Inidcate that porble has been build
1804 
1806 }
#define ProblemManagerFunctionBegin
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:126
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221
Deprecated interface functions.
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)

◆ buildProblem() [1/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblem ( const std::string  name,
const bool  square_matrix,
int  verb = VERBOSE 
)

build problem data structures

Note
If square_matrix is set to true, that indicate that problem is structurally symmetric, i.e. rows and columns have the same dofs and are indexed in the same way.
Parameters
nameproblem name
square_matrixstructurally symmetric problem
verbverbosity level
Returns
error code

Definition at line 428 of file ProblemsManager.cpp.

430  {
431 
432  MoFEM::Interface &m_field = cOre;
434  if (!(cOre.getBuildMoFEM() & (1 << 0)))
435  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
436  if (!(cOre.getBuildMoFEM() & (1 << 1)))
437  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
438  if (!(cOre.getBuildMoFEM() & (1 << 2)))
439  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
440  const Problem *problem_ptr;
441  CHKERR m_field.get_problem(name, &problem_ptr);
442  CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
443  cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
444  // function knows what he is doing
446 }
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures

◆ buildProblem() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblem ( Problem problem_ptr,
const bool  square_matrix,
int  verb = VERBOSE 
)

build problem data structures

Note
If square_matrix is set to true, that indicate that problem is structurally symmetric, i.e. rows and columns have the same dofs and are indexed in the same way.
Parameters
problempointer
square_matrixstructurally symmetric problem
verbverbosity level
Returns
error code

Definition at line 448 of file ProblemsManager.cpp.

450  {
451  MoFEM::Interface &m_field = cOre;
452  auto dofs_field_ptr = m_field.get_dofs();
453  auto ents_field_ptr = m_field.get_field_ents();
454  auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
456  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
457 
458  // Note: Only allowed changes on problem_ptr structure which not influence
459  // multi-index.
460 
461  if (problem_ptr->getBitRefLevel().none()) {
462  SETERRQ1(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
463  problem_ptr->getName().c_str());
464  }
465  CHKERR m_field.clear_problem(problem_ptr->getName());
466 
467  // zero finite elements
468  problem_ptr->numeredFiniteElementsPtr->clear();
469 
470  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
471  auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
473  for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
474 
475  const auto uid = (*it)->getLocalUniqueId();
476 
477  auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
478  for (auto lo = r.first; lo != r.second; ++lo) {
479 
480  if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
481  std::array<bool, 2> row_col = {false, false};
482 
483  const BitRefLevel prb_bit = problem_ptr->getBitRefLevel();
484  const BitRefLevel prb_mask = problem_ptr->getMaskBitRefLevel();
485  const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
486 
487  // if entity is not problem refinement level
488  if ((fe_bit & prb_mask) != fe_bit)
489  continue;
490  if ((fe_bit & prb_bit) != prb_bit)
491  continue;
492 
493  auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
494  auto dit = nb_dofs->lower_bound(uid);
495  decltype(dit) hi_dit;
496  if (dit != nb_dofs->end()) {
497  hi_dit = nb_dofs->upper_bound(
498  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
499  view.insert(dit, hi_dit);
500  row_col[rc] = true;
501  }
502  };
503 
504  if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
505  add_to_view(dofs_field_ptr, dofs_rows, ROW);
506 
507  if (!square_matrix)
508  if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
509  add_to_view(dofs_field_ptr, dofs_cols, COL);
510 
511  if (square_matrix && row_col[ROW])
512  break;
513  else if (row_col[ROW] && row_col[COL])
514  break;
515  }
516  }
517  }
519  };
520 
521  CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
522 
523  // Add row dofs to problem
524  {
525  // zero rows
526  problem_ptr->nbDofsRow = 0;
527  problem_ptr->nbLocDofsRow = 0;
528  problem_ptr->nbGhostDofsRow = 0;
529 
530  // add dofs for rows
531  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
532  hi_miit;
533  hi_miit = dofs_rows.get<0>().end();
534 
535  int count_dofs = dofs_rows.get<1>().count(true);
536  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
537  boost::shared_ptr<std::vector<NumeredDofEntity>>(
538  new std::vector<NumeredDofEntity>());
539  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
540  dofs_array->reserve(count_dofs);
541  miit = dofs_rows.get<0>().begin();
542  for (; miit != hi_miit; miit++) {
543  if ((*miit)->getActive()) {
544  dofs_array->emplace_back(*miit);
545  dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
546  }
547  }
548  auto hint = problem_ptr->numeredRowDofsPtr->end();
549  for (auto &v : *dofs_array) {
550  hint = problem_ptr->numeredRowDofsPtr->emplace_hint(hint, dofs_array, &v);
551  }
552  }
553 
554  // Add col dofs to problem
555  if (!square_matrix) {
556  // zero cols
557  problem_ptr->nbDofsCol = 0;
558  problem_ptr->nbLocDofsCol = 0;
559  problem_ptr->nbGhostDofsCol = 0;
560 
561  // add dofs for cols
562  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
563  hi_miit;
564  hi_miit = dofs_cols.get<0>().end();
565 
566  int count_dofs = 0;
567  miit = dofs_cols.get<0>().begin();
568  for (; miit != hi_miit; miit++) {
569  if (!(*miit)->getActive()) {
570  continue;
571  }
572  count_dofs++;
573  }
574 
575  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
576  boost::shared_ptr<std::vector<NumeredDofEntity>>(
577  new std::vector<NumeredDofEntity>());
578  problem_ptr->getColDofsSequence()->push_back(dofs_array);
579  dofs_array->reserve(count_dofs);
580  miit = dofs_cols.get<0>().begin();
581  for (; miit != hi_miit; miit++) {
582  if (!(*miit)->getActive()) {
583  continue;
584  }
585  dofs_array->emplace_back(*miit);
586  dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
587  }
588  auto hint = problem_ptr->numeredColDofsPtr->end();
589  for (auto &v : *dofs_array) {
590  hint = problem_ptr->numeredColDofsPtr->emplace_hint(hint, dofs_array, &v);
591  }
592  } else {
593  problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
594  problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
595  problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
596  }
597 
598  // job done, some debugging and postprocessing
599  if (verb >= VERBOSE) {
600  MOFEM_LOG("SYNC", Sev::verbose)
601  << problem_ptr->getName() << " Nb. local dofs "
602  << problem_ptr->numeredRowDofsPtr->size() << " by "
603  << problem_ptr->numeredColDofsPtr->size();
604  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
605  }
606 
607  if (verb >= NOISY) {
608  MOFEM_LOG("SYNC", Sev::noisy)
609  << "FEs row dofs " << *problem_ptr << " Nb. row dof "
610  << problem_ptr->getNbDofsRow();
611  for (auto &miit : *problem_ptr->numeredRowDofsPtr)
612  MOFEM_LOG("SYNC", Sev::noisy) << *miit;
613 
614  MOFEM_LOG("SYNC", Sev::noisy)
615  << "FEs col dofs " << *problem_ptr << " Nb. col dof "
616  << problem_ptr->getNbDofsCol();
617  for (auto &miit : *problem_ptr->numeredColDofsPtr)
618  MOFEM_LOG("SYNC", Sev::noisy) << *miit;
619  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
620  }
621 
622  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
623  // uses this function knows
624  // what he is doing
625 
626  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
627 
629 }
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:349
@ NOISY
Definition: definitions.h:280
@ VERBOSE
Definition: definitions.h:278
@ COL
Definition: definitions.h:192
@ ROW
Definition: definitions.h:192
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:305
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
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
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 DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
uint128_t UId
Unique Id.
Definition: Types.hpp:42
const double r
rate factor
PetscLogEvent MOFEM_EVENT_ProblemsManager

◆ buildProblemOnDistributedMesh() [1/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh ( const std::string  name,
const bool  square_matrix,
int  verb = VERBOSE 
)

build problem data structures, assuming that mesh is distributed (collective)

Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities.

Collective - need to be run on all processors in communicator, i.e. each function has to call this function.

Definition at line 631 of file ProblemsManager.cpp.

632  {
633  MoFEM::Interface &m_field = cOre;
635 
636  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
637  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
638  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
639  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
640  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
641  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
642 
643  const Problem *problem_ptr;
644  CHKERR m_field.get_problem(name, &problem_ptr);
645  CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
646  square_matrix, verb);
647 
650 
652 }
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)

◆ buildProblemOnDistributedMesh() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh ( Problem problem_ptr,
const bool  square_matrix = true,
int  verb = VERBOSE 
)

build problem data structures, assuming that mesh is distributed (collective)

Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities.

Collective - need to be run on all processors in communicator, i.e. each function has to call this function.

Definition at line 654 of file ProblemsManager.cpp.

655  {
656  MoFEM::Interface &m_field = cOre;
657  auto fields_ptr = m_field.get_fields();
658  auto fe_ptr = m_field.get_finite_elements();
659  auto fe_ent_ptr = m_field.get_ents_finite_elements();
660  auto dofs_field_ptr = m_field.get_dofs();
662  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
663 
664  // clear data structures
665  CHKERR m_field.clear_problem(problem_ptr->getName());
666 
667  CHKERR getOptions();
668 
669  if (problem_ptr->getBitRefLevel().none())
670  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
671  "problem <%s> refinement level not set",
672  problem_ptr->getName().c_str());
673 
674  int loop_size = 2;
675  if (square_matrix) {
676  loop_size = 1;
677  problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
678  } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
679  problem_ptr->numeredColDofsPtr =
680  boost::shared_ptr<NumeredDofEntity_multiIndex>(
682  }
683 
684  const BitRefLevel prb_bit = problem_ptr->getBitRefLevel();
685  const BitRefLevel prb_mask = problem_ptr->getMaskBitRefLevel();
686 
687  // // get rows and cols dofs view based on data on elements
688  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
689 
690  // Add DOFs to problem by visiting all elements and adding DOFs from
691  // elements to the problem
692  if (buildProblemFromFields == PETSC_FALSE) {
693 
694  auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
695  auto ents_field_ptr = m_field.get_field_ents();
696  auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
698  for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
699  ++it) {
700 
701  const auto uid = (*it)->getLocalUniqueId();
702 
703  auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
704  for (auto lo = r.first; lo != r.second; ++lo) {
705 
706  if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
707  std::array<bool, 2> row_col = {false, false};
708 
709  const BitRefLevel prb_bit = problem_ptr->getBitRefLevel();
710  const BitRefLevel prb_mask = problem_ptr->getMaskBitRefLevel();
711  const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
712 
713  // if entity is not problem refinement level
714  if ((fe_bit & prb_mask) != fe_bit)
715  continue;
716  if ((fe_bit & prb_bit) != prb_bit)
717  continue;
718 
719  auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
720  auto dit = nb_dofs->lower_bound(uid);
721  decltype(dit) hi_dit;
722  if (dit != nb_dofs->end()) {
723  hi_dit = nb_dofs->upper_bound(
724  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
725  view.insert(dit, hi_dit);
726  row_col[rc] = true;
727  }
728  };
729 
730  if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
731  add_to_view(dofs_field_ptr, dofs_rows, ROW);
732 
733  if (!square_matrix)
734  if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
735  add_to_view(dofs_field_ptr, dofs_cols, COL);
736 
737  if (square_matrix && row_col[ROW])
738  break;
739  else if (row_col[ROW] && row_col[COL])
740  break;
741  }
742  }
743  }
745  };
746 
747  CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
748  }
749 
750  // Add DOFS to the problem by searching all the filedes, and adding to
751  // problem owned or shared DOFs
752  if (buildProblemFromFields == PETSC_TRUE) {
753  // Get fields IDs on elements
754  BitFieldId fields_ids_row, fields_ids_col;
755  for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
756  fit != fe_ptr->end(); fit++) {
757  if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
758  fields_ids_row |= fit->get()->getBitFieldIdRow();
759  fields_ids_col |= fit->get()->getBitFieldIdCol();
760  }
761  }
762  // Get fields DOFs
763  for (Field_multiIndex::iterator fit = fields_ptr->begin();
764  fit != fields_ptr->end(); fit++) {
765  if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
766 
767  auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
768  FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
769  auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
770  FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
771 
772  for (; dit != hi_dit; dit++) {
773 
774  const int owner_proc = dit->get()->getOwnerProc();
775  if (owner_proc != m_field.get_comm_rank()) {
776  const unsigned char pstatus = dit->get()->getPStatus();
777  if (pstatus == 0) {
778  continue;
779  }
780  }
781  const BitRefLevel dof_bit = (*dit)->getBitRefLevel();
782  // if entity is not problem refinement level
783  if ((dof_bit & prb_mask) != dof_bit)
784  continue;
785  if ((dof_bit & prb_bit) != prb_bit)
786  continue;
787  if ((fit->get()->getId() & fields_ids_row).any()) {
788  dofs_rows.insert(*dit);
789  }
790  if (!square_matrix) {
791  if ((fit->get()->getId() & fields_ids_col).any()) {
792  dofs_cols.insert(*dit);
793  }
794  }
795  }
796  }
797  }
798  }
799 
801  // Get fields IDs on elements
802  BitFieldId fields_ids_row, fields_ids_col;
803  BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
804  for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
805  fit != fe_ptr->end(); fit++) {
806  if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
807  fields_ids_row |= fit->get()->getBitFieldIdRow();
808  fields_ids_col |= fit->get()->getBitFieldIdCol();
809  }
810  }
811 
812  DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
813  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
814  DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
815  hi_miit;
816  miit = dofs_ptr[ss]->get<1>().lower_bound(1);
817  hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
818  Range ents_to_synchronise;
819  for (; miit != hi_miit; ++miit) {
820  if (miit->get()->getEntDofIdx() != 0)
821  continue;
822  ents_to_synchronise.insert(miit->get()->getEnt());
823  }
824  Range tmp_ents = ents_to_synchronise;
825  CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
826  ents_to_synchronise, verb);
827  ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
828  for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
829  if ((fit->get()->getId() & *fields_ids[ss]).any()) {
830  const auto bit_number = (*fit)->getBitNumber();
831  for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
832  pit != ents_to_synchronise.pair_end(); ++pit) {
833  const auto f = pit->first;
834  const auto s = pit->second;
835  const auto lo_uid =
837  const auto hi_uid =
839 
840  auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
841  auto hi_dit =
842  dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
843 
844  dofs_ptr[ss]->insert(dit, hi_dit);
845  }
846  }
847  }
848  }
849  }
850 
851  // add dofs for rows and cols and set ownership
852  DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
853  boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
854  problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
855  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
856  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
857  &problem_ptr->nbLocDofsCol};
858  int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
859  &problem_ptr->nbGhostDofsCol};
860  for (int ss = 0; ss < 2; ss++) {
861  *(nbdof_ptr[ss]) = 0;
862  *(local_nbdof_ptr[ss]) = 0;
863  *(ghost_nbdof_ptr[ss]) = 0;
864  }
865 
866  // Loop over dofs on rows and columns and add to multi-indices in dofs
867  // problem structure, set partition for each dof
868  int nb_local_dofs[] = {0, 0};
869  for (int ss = 0; ss < loop_size; ss++) {
870  DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
871  hi_miit;
872  miit = dofs_ptr[ss]->get<1>().lower_bound(1);
873  hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
874  for (; miit != hi_miit; miit++) {
875  int owner_proc = (*miit)->getOwnerProc();
876  if (owner_proc == m_field.get_comm_rank()) {
877  nb_local_dofs[ss]++;
878  }
879  }
880  }
881 
882  // get layout
883  int start_ranges[2], end_ranges[2];
884  for (int ss = 0; ss != loop_size; ss++) {
885  PetscLayout layout;
886  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
887  CHKERR PetscLayoutSetBlockSize(layout, 1);
888  CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
889  CHKERR PetscLayoutSetUp(layout);
890  CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
891  CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
892  CHKERR PetscLayoutDestroy(&layout);
893  }
894  if (square_matrix) {
895  nbdof_ptr[1] = nbdof_ptr[0];
896  nb_local_dofs[1] = nb_local_dofs[0];
897  start_ranges[1] = start_ranges[0];
898  end_ranges[1] = end_ranges[0];
899  }
900 
901  // if(sizeof(UId) != SIZEOFUID) {
902  // SETERRQ2(
903  // PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,
904  // "check size of UId, size of UId is %u != %u",
905  // sizeof(UId),SIZEOFUID
906  // );
907  // }
908 
909  // set local and global indices on own dofs
910  const size_t idx_data_type_size = sizeof(IdxDataType);
911  const size_t data_block_size = idx_data_type_size / sizeof(int);
912 
913  if (sizeof(IdxDataType) % sizeof(int)) {
914  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
915  }
916  std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
917  m_field.get_comm_size()),
918  ids_data_packed_cols(m_field.get_comm_size());
919 
920  // Loop over dofs on this processor and prepare those dofs to send on
921  // another proc
922  for (int ss = 0; ss != loop_size; ++ss) {
923 
924  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
925  hi_miit;
926  hi_miit = dofs_ptr[ss]->get<0>().end();
927 
928  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
929  boost::shared_ptr<std::vector<NumeredDofEntity>>(
930  new std::vector<NumeredDofEntity>());
931  int nb_dofs_to_add = 0;
932  miit = dofs_ptr[ss]->get<0>().begin();
933  for (; miit != hi_miit; ++miit) {
934  // Only set global idx for dofs on this processor part
935  if (!(miit->get()->getActive()))
936  continue;
937  ++nb_dofs_to_add;
938  }
939  dofs_array->reserve(nb_dofs_to_add);
940  if (ss == 0) {
941  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
942  } else {
943  problem_ptr->getColDofsSequence()->push_back(dofs_array);
944  }
945 
946  int &local_idx = *local_nbdof_ptr[ss];
947  miit = dofs_ptr[ss]->get<0>().begin();
948  for (; miit != hi_miit; ++miit) {
949 
950  // Only set global idx for dofs on this processor part
951  if (!(miit->get()->getActive()))
952  continue;
953 
954  dofs_array->emplace_back(*miit);
955 
956  int owner_proc = dofs_array->back().getOwnerProc();
957  if (owner_proc < 0) {
958  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
959  "data inconsistency");
960  }
961 
962  if (owner_proc != m_field.get_comm_rank()) {
963  dofs_array->back().pArt = owner_proc;
964  dofs_array->back().dofIdx = -1;
965  dofs_array->back().petscGloablDofIdx = -1;
966  dofs_array->back().petscLocalDofIdx = -1;
967  } else {
968 
969  // Set part and indexes
970  int glob_idx = start_ranges[ss] + local_idx;
971  dofs_array->back().pArt = owner_proc;
972  dofs_array->back().dofIdx = glob_idx;
973  dofs_array->back().petscGloablDofIdx = glob_idx;
974  dofs_array->back().petscLocalDofIdx = local_idx;
975  local_idx++;
976 
977  unsigned char pstatus = dofs_array->back().getPStatus();
978  // check id dof is shared, if that is a case global idx added to data
979  // structure and is sended to other processors
980  if (pstatus > 0) {
981 
982  for (int proc = 0;
983  proc < MAX_SHARING_PROCS &&
984  -1 != dofs_array->back().getSharingProcsPtr()[proc];
985  proc++) {
986  // make it different for rows and columns when needed
987  if (ss == 0) {
988  ids_data_packed_rows[dofs_array->back()
989  .getSharingProcsPtr()[proc]]
990  .emplace_back(dofs_array->back().getGlobalUniqueId(),
991  glob_idx);
992  } else {
993  ids_data_packed_cols[dofs_array->back()
994  .getSharingProcsPtr()[proc]]
995  .emplace_back(dofs_array->back().getGlobalUniqueId(),
996  glob_idx);
997  }
998  if (!(pstatus & PSTATUS_MULTISHARED)) {
999  break;
1000  }
1001  }
1002  }
1003  }
1004  }
1005 
1006  auto hint = numered_dofs_ptr[ss]->end();
1007  for (auto &v : *dofs_array)
1008  hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
1009  }
1010  if (square_matrix) {
1011  local_nbdof_ptr[1] = local_nbdof_ptr[0];
1012  }
1013 
1014  int nsends_rows = 0, nsends_cols = 0;
1015  // Non zero lengths[i] represent a message to i of length lengths[i].
1016  std::vector<int> lengths_rows(m_field.get_comm_size()),
1017  lengths_cols(m_field.get_comm_size());
1018  lengths_rows.clear();
1019  lengths_cols.clear();
1020  for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
1021  lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
1022  lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
1023  if (!ids_data_packed_rows[proc].empty())
1024  nsends_rows++;
1025  if (!ids_data_packed_cols[proc].empty())
1026  nsends_cols++;
1027  }
1028 
1029  MPI_Status *status;
1030  CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
1031 
1032  // Do rows
1033  int nrecvs_rows; // number of messages received
1034  int *onodes_rows; // list of node-ids from which messages are expected
1035  int *olengths_rows; // corresponding message lengths
1036  int **rbuf_row; // must bee freed by user
1037 
1038  {
1039 
1040  // make sure it is a PETSc comm
1041  CHKERR PetscCommDuplicate(m_field.get_comm(), &m_field.get_comm(), NULL);
1042 
1043  // rows
1044 
1045  // Computes the number of messages a node expects to receive
1046  CHKERR PetscGatherNumberOfMessages(m_field.get_comm(), NULL,
1047  &lengths_rows[0], &nrecvs_rows);
1048  // std::cerr << nrecvs_rows << std::endl;
1049 
1050  // Computes info about messages that a MPI-node will receive, including
1051  // (from-id,length) pairs for each message.
1052  CHKERR PetscGatherMessageLengths(m_field.get_comm(), nsends_rows,
1053  nrecvs_rows, &lengths_rows[0],
1054  &onodes_rows, &olengths_rows);
1055 
1056  // Gets a unique new tag from a PETSc communicator. All processors that
1057  // share the communicator MUST call this routine EXACTLY the same number
1058  // of times. This tag should only be used with the current objects
1059  // communicator; do NOT use it with any other MPI communicator.
1060  int tag_row;
1061  CHKERR PetscCommGetNewTag(m_field.get_comm(), &tag_row);
1062 
1063  // Allocate a buffer sufficient to hold messages of size specified in
1064  // olengths. And post Irecvs on these buffers using node info from onodes
1065  MPI_Request *r_waits_row; // must bee freed by user
1066  // rbuf has a pointers to messeges. It has size of of nrecvs (number of
1067  // messages) +1. In the first index a block is allocated,
1068  // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
1069 
1070  CHKERR PetscPostIrecvInt(m_field.get_comm(), tag_row, nrecvs_rows,
1071  onodes_rows, olengths_rows, &rbuf_row,
1072  &r_waits_row);
1073  CHKERR PetscFree(onodes_rows);
1074 
1075  MPI_Request *s_waits_row; // status of sens messages
1076  CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
1077 
1078  // Send messeges
1079  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
1080  if (!lengths_rows[proc])
1081  continue; // no message to send to this proc
1082  CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
1083  lengths_rows[proc], // message length
1084  MPIU_INT, proc, // to proc
1085  tag_row, m_field.get_comm(), s_waits_row + kk);
1086  kk++;
1087  }
1088 
1089  if (nrecvs_rows) {
1090  CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
1091  }
1092  if (nsends_rows) {
1093  CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
1094  }
1095 
1096  CHKERR PetscFree(r_waits_row);
1097  CHKERR PetscFree(s_waits_row);
1098  }
1099 
1100  // cols
1101  int nrecvs_cols = nrecvs_rows;
1102  int *olengths_cols = olengths_rows;
1103  PetscInt **rbuf_col = rbuf_row;
1104  if (!square_matrix) {
1105 
1106  // Computes the number of messages a node expects to receive
1107  CHKERR PetscGatherNumberOfMessages(m_field.get_comm(), NULL,
1108  &lengths_cols[0], &nrecvs_cols);
1109 
1110  // Computes info about messages that a MPI-node will receive, including
1111  // (from-id,length) pairs for each message.
1112  int *onodes_cols;
1113  CHKERR PetscGatherMessageLengths(m_field.get_comm(), nsends_cols,
1114  nrecvs_cols, &lengths_cols[0],
1115  &onodes_cols, &olengths_cols);
1116 
1117  // Gets a unique new tag from a PETSc communicator.
1118  int tag_col;
1119  CHKERR PetscCommGetNewTag(m_field.get_comm(), &tag_col);
1120 
1121  // Allocate a buffer sufficient to hold messages of size specified in
1122  // olengths. And post Irecvs on these buffers using node info from onodes
1123  MPI_Request *r_waits_col; // must bee freed by user
1124  CHKERR PetscPostIrecvInt(m_field.get_comm(), tag_col, nrecvs_cols,
1125  onodes_cols, olengths_cols, &rbuf_col,
1126  &r_waits_col);
1127  CHKERR PetscFree(onodes_cols);
1128 
1129  MPI_Request *s_waits_col; // status of sens messages
1130  CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
1131 
1132  // Send messages
1133  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
1134  if (!lengths_cols[proc])
1135  continue; // no message to send to this proc
1136  CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
1137  lengths_cols[proc], // message length
1138  MPIU_INT, proc, // to proc
1139  tag_col, m_field.get_comm(), s_waits_col + kk);
1140  kk++;
1141  }
1142 
1143  if (nrecvs_cols) {
1144  CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
1145  }
1146  if (nsends_cols) {
1147  CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
1148  }
1149 
1150  CHKERR PetscFree(r_waits_col);
1151  CHKERR PetscFree(s_waits_col);
1152  }
1153 
1154  CHKERR PetscFree(status);
1155 
1156  DofEntity_multiIndex_global_uid_view dofs_glob_uid_view;
1157  auto hint = dofs_glob_uid_view.begin();
1158  for (auto dof : *m_field.get_dofs())
1159  dofs_glob_uid_view.emplace_hint(hint, dof);
1160 
1161  // set values received from other processors
1162  for (int ss = 0; ss != loop_size; ++ss) {
1163 
1164  int nrecvs;
1165  int *olengths;
1166  int **data_procs;
1167  if (ss == 0) {
1168  nrecvs = nrecvs_rows;
1169  olengths = olengths_rows;
1170  data_procs = rbuf_row;
1171  } else {
1172  nrecvs = nrecvs_cols;
1173  olengths = olengths_cols;
1174  data_procs = rbuf_col;
1175  }
1176 
1177  UId uid;
1178  for (int kk = 0; kk != nrecvs; ++kk) {
1179  int len = olengths[kk];
1180  int *data_from_proc = data_procs[kk];
1181  for (int dd = 0; dd < len; dd += data_block_size) {
1182  uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
1183  auto ddit = dofs_glob_uid_view.find(uid);
1184  if (ddit == dofs_glob_uid_view.end()) {
1185  std::ostringstream zz;
1186  zz << uid << std::endl;
1187  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1188  "no such dof %s in mofem database", zz.str().c_str());
1189  }
1190 
1191  auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
1192 
1193  if (dit != numered_dofs_ptr[ss]->end()) {
1194 
1195  int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
1196  if (global_idx < 0)
1197  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1198  "received negative dof");
1199  bool success;
1200  success = numered_dofs_ptr[ss]->modify(
1201  dit, NumeredDofEntity_mofem_index_change(global_idx));
1202  if (!success)
1203  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1204  "modification unsuccessful");
1205  success = numered_dofs_ptr[ss]->modify(
1206  dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
1207  global_idx));
1208  if (!success)
1209  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1210  "modification unsuccessful");
1211 
1212  } else if (ddit->get()->getPStatus() == 0) {
1213 
1214  // Dof is shared on this processor, however there is no element
1215  // which have this dof. If DOF is not shared and received from other
1216  // processor, but not marked as a shared on other that means that is
1217  // data inconstancy and error should be thorwed.
1218 
1219  std::ostringstream zz;
1220  zz << **ddit << std::endl;
1221  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1222  "data inconsistency, dofs is not shared, but received "
1223  "from other proc\n"
1224  "%s",
1225  zz.str().c_str());
1226  }
1227  }
1228  }
1229  }
1230 
1231  if (square_matrix) {
1232  (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
1233  (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
1234  }
1235 
1236  CHKERR PetscFree(olengths_rows);
1237  CHKERR PetscFree(rbuf_row[0]);
1238  CHKERR PetscFree(rbuf_row);
1239  if (!square_matrix) {
1240  CHKERR PetscFree(olengths_cols);
1241  CHKERR PetscFree(rbuf_col[0]);
1242  CHKERR PetscFree(rbuf_col);
1243  }
1244 
1245  if (square_matrix) {
1246  if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
1247  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1248  "data inconsistency for square_matrix %d!=%d",
1249  numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
1250  }
1251  if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
1252  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1253  "data inconsistency for square_matrix");
1254  }
1255  }
1256 
1257  CHKERR printPartitionedProblem(problem_ptr, verb);
1258  CHKERR debugPartitionedProblem(problem_ptr, verb);
1259 
1260  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
1261 
1263 }
@ MOFEM_NOT_FOUND
Definition: definitions.h:125
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
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
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)
Definition: ddTensor0.hpp:33
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
virtual int get_comm_size() const =0
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
MoFEMErrorCode getOptions()
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.

◆ buildSubProblem()

MoFEMErrorCode MoFEM::ProblemsManager::buildSubProblem ( const std::string  out_name,
const std::vector< std::string > &  fields_row,
const std::vector< std::string > &  fields_col,
const std::string  main_problem,
const bool  square_matrix = true,
const map< std::string, std::pair< EntityType, EntityType >> *  entityMapRow = nullptr,
const map< std::string, std::pair< EntityType, EntityType >> *  entityMapCol = nullptr,
int  verb = VERBOSE 
)

build sub problem

Parameters
out_nameproblem
fields_rowvector of fields composing problem
fields_colvector of fields composing problem
main_problemmain problem
Returns
error code

Definition at line 1265 of file ProblemsManager.cpp.

1271  {
1272  MoFEM::Interface &m_field = cOre;
1273  auto problems_ptr = m_field.get_problems();
1275 
1276  CHKERR m_field.clear_problem(out_name);
1277 
1278  // get reference to all problems
1280  ProblemByName &problems_by_name =
1281  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1282 
1283  // get iterators to out problem, i.e. build problem
1284  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1285  if (out_problem_it == problems_by_name.end()) {
1286  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1287  "subproblem with name < %s > not defined (top tip check spelling)",
1288  out_name.c_str());
1289  }
1290  // get iterator to main problem, i.e. out problem is subproblem of main
1291  // problem
1292  ProblemByName::iterator main_problem_it = problems_by_name.find(main_problem);
1293  if (main_problem_it == problems_by_name.end()) {
1294  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1295  "problem of subproblem with name < %s > not defined (top tip "
1296  "check spelling)",
1297  main_problem.c_str());
1298  }
1299 
1300  // get dofs for row & columns for out problem,
1301  boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1302  out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1303  // get dofs for row & columns for main problem
1304  boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1305  main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1306  // get local indices counter
1307  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1308  &out_problem_it->nbLocDofsCol};
1309  // get global indices counter
1310  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1311 
1312  // set number of ghost nodes to zero
1313  {
1314  out_problem_it->nbGhostDofsRow = 0;
1315  out_problem_it->nbGhostDofsCol = 0;
1316  }
1317 
1318  // put rows & columns field names in array
1319  std::vector<std::string> fields[] = {fields_row, fields_col};
1320  const map<std::string, std::pair<EntityType, EntityType>> *entityMap[] = {
1321  entityMapRow, entityMapCol};
1322 
1323  // make data structure fos sub-problem data
1324  out_problem_it->subProblemData =
1325  boost::make_shared<Problem::SubProblemData>();
1326 
1327  // Loop over rows and columns
1328  for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1329 
1330  // reset dofs and columns counters
1331  (*nb_local_dofs[ss]) = 0;
1332  (*nb_dofs[ss]) = 0;
1333  // clear arrays
1334  out_problem_dofs[ss]->clear();
1335 
1336  // If DOFs are cleared clear finite elements too.
1337  out_problem_it->numeredFiniteElementsPtr->clear();
1338 
1339  // get dofs by field name and insert them in out problem multi-indices
1340  for (auto field : fields[ss]) {
1341 
1342  // Following reserve memory in sequences, only two allocations are here,
1343  // once for array of objects, next for array of shared pointers
1344 
1345  // aliased sequence of pointer is killed with element
1346  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1347  boost::make_shared<std::vector<NumeredDofEntity>>();
1348  // reserve memory for field dofs
1349  if (!ss)
1350  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1351  else
1352  out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1353 
1354  // create elements objects
1355  auto bit_number = m_field.get_field_bit_number(field);
1356  auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1357  FieldEntity::getLoBitNumberUId(bit_number));
1358  auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1359  FieldEntity::getHiBitNumberUId(bit_number));
1360 
1361  auto add_dit_to_dofs_array = [&](auto &dit) {
1362  dofs_array->emplace_back(
1363  dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1364  dit->get()->getPetscGlobalDofIdx(),
1365  dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1366  };
1367 
1368  if (entityMap[ss]) {
1369  auto mit = entityMap[ss]->find(field);
1370  if (mit != entityMap[ss]->end()) {
1371  EntityType lo_type = mit->second.first;
1372  EntityType hi_type = mit->second.second;
1373  int count = 0;
1374  for (auto diit = dit; diit != hi_dit; ++diit) {
1375  EntityType ent_type = (*diit)->getEntType();
1376  if (ent_type >= lo_type && ent_type <= hi_type)
1377  ++count;
1378  }
1379  dofs_array->reserve(count);
1380  for (; dit != hi_dit; ++dit) {
1381  EntityType ent_type = (*dit)->getEntType();
1382  if (ent_type >= lo_type && ent_type <= hi_type)
1383  add_dit_to_dofs_array(dit);
1384  }
1385  } else {
1386  dofs_array->reserve(std::distance(dit, hi_dit));
1387  for (; dit != hi_dit; dit++)
1388  add_dit_to_dofs_array(dit);
1389  }
1390  } else {
1391  dofs_array->reserve(std::distance(dit, hi_dit));
1392  for (; dit != hi_dit; dit++)
1393  add_dit_to_dofs_array(dit);
1394  }
1395 
1396  // fill multi-index
1397  auto hint = out_problem_dofs[ss]->end();
1398  for (auto &v : *dofs_array)
1399  hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1400  }
1401  // Set local indexes
1402  {
1403  auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1404  auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1405  for (; dit != hi_dit; dit++) {
1406  int idx = -1; // if dof is not part of partition, set local index to -1
1407  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1408  idx = (*nb_local_dofs[ss])++;
1409  }
1410  bool success = out_problem_dofs[ss]->modify(
1411  out_problem_dofs[ss]->project<0>(dit),
1412  NumeredDofEntity_local_idx_change(idx));
1413  if (!success) {
1414  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1415  "operation unsuccessful");
1416  }
1417  };
1418  }
1419  // Set global indexes, compress global indices
1420  {
1421  auto dit =
1422  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1423  auto hi_dit =
1424  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1425  out_problem_dofs[ss]->size());
1426  const int nb = std::distance(dit, hi_dit);
1427  // get main problem global indices
1428  std::vector<int> main_indices(nb);
1429  for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1430  *it = dit->get()->getPetscGlobalDofIdx();
1431  }
1432  // create is with global dofs
1433  IS is;
1434  CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1435  PETSC_USE_POINTER, &is);
1436  // create map form main problem global indices to out problem global
1437  // indices
1438  AO ao;
1439  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1440  if (ss == 0) {
1441  IS is_dup;
1442  CHKERR ISDuplicate(is, &is_dup);
1443  out_problem_it->getSubData()->rowIs = SmartPetscObj<IS>(is_dup, false);
1444  out_problem_it->getSubData()->rowMap = SmartPetscObj<AO>(ao, true);
1445  } else {
1446  IS is_dup;
1447  CHKERR ISDuplicate(is, &is_dup);
1448  out_problem_it->getSubData()->colIs = SmartPetscObj<IS>(is_dup,false);
1449  out_problem_it->getSubData()->colMap = SmartPetscObj<AO>(ao, true);
1450  }
1451  CHKERR AOApplicationToPetscIS(ao, is);
1452  // set global number of DOFs
1453  CHKERR ISGetSize(is, nb_dofs[ss]);
1454  CHKERR ISDestroy(&is);
1455  // set out problem global indices after applying map
1456  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1457  for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1458  dit++, it++) {
1459  bool success = out_problem_dofs[ss]->modify(
1460  out_problem_dofs[ss]->project<0>(dit),
1461  NumeredDofEntity_part_and_all_indices_change(
1462  dit->get()->getPart(), *it, *it,
1463  dit->get()->getPetscLocalDofIdx()));
1464  if (!success) {
1465  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1466  "operation unsuccessful");
1467  }
1468  }
1469  // set global indices to nodes not on this part
1470  {
1471  NumeredDofEntityByLocalIdx::iterator dit =
1472  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1473  NumeredDofEntityByLocalIdx::iterator hi_dit =
1474  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1475  const int nb = std::distance(dit, hi_dit);
1476  std::vector<int> main_indices_non_local(nb);
1477  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1478  dit++, it++) {
1479  *it = dit->get()->getPetscGlobalDofIdx();
1480  }
1481  IS is;
1482  CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1483  &*main_indices_non_local.begin(),
1484  PETSC_USE_POINTER, &is);
1485  CHKERR AOApplicationToPetscIS(ao, is);
1486  CHKERR ISDestroy(&is);
1487  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1488  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1489  dit++, it++) {
1490  bool success = out_problem_dofs[ss]->modify(
1491  out_problem_dofs[ss]->project<0>(dit),
1492  NumeredDofEntity_part_and_all_indices_change(
1493  dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1494  dit->get()->getPetscLocalDofIdx()));
1495  if (!success) {
1496  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1497  "operation unsuccessful");
1498  }
1499  }
1500  }
1501  CHKERR AODestroy(&ao);
1502  }
1503  }
1504 
1505  if (square_matrix) {
1506  out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1507  out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1508  out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1509  out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1510  out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1511  CHKERR PetscObjectReference(
1512  (PetscObject)out_problem_it->getSubData()->rowIs);
1513  }
1514 
1515  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1516  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1517 
1519 }
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number

◆ getFEMeshset()

MoFEMErrorCode MoFEM::ProblemsManager::getFEMeshset ( const std::string  prb_name,
const std::string  fe_name,
EntityHandle *  meshset 
) const

create add entities of finite element in the problem

\create add entities of finite element in the problem

Note
Meshset entity has to be crated
Meshset entity has to be crated
Parameters
prb_namename of the problem
fe_namename of the entity
meshsetpointer meshset handle
Returns
MoFEMErrorCode

Definition at line 2819 of file ProblemsManager.cpp.

2821  {
2822  MoFEM::Interface &m_field = cOre;
2823  const Problem *problem_ptr;
2825 
2826  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2827  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2828  auto fit =
2829  problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2830  .lower_bound(fe_name);
2831  auto hi_fe_it =
2832  problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2833  .upper_bound(fe_name);
2834  std::vector<EntityHandle> fe_vec;
2835  fe_vec.reserve(std::distance(fit, hi_fe_it));
2836  for (; fit != hi_fe_it; fit++)
2837  fe_vec.push_back(fit->get()->getEnt());
2838  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2839  fe_vec.size());
2841 }
virtual moab::Interface & get_moab()=0
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements

◆ getProblemElementsLayout()

MoFEMErrorCode MoFEM::ProblemsManager::getProblemElementsLayout ( const std::string  name,
const std::string  fe_name,
PetscLayout *  layout 
) const

Get layout of elements in the problem.

In layout is stored information how many elements is on each processor, for more information look int petsc documentation http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscLayoutCreate.html#PetscLayoutCreate

Parameters
nameproblem name
fe_namefinite elements
layoutlayout
verbverbosity level
Returns
error code

Definition at line 2844 of file ProblemsManager.cpp.

2846  {
2847  MoFEM::Interface &m_field = cOre;
2848  const Problem *problem_ptr;
2850  CHKERR m_field.get_problem(name, &problem_ptr);
2851  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2852  fe_name, layout);
2854 }

◆ inheritPartition()

MoFEMErrorCode MoFEM::ProblemsManager::inheritPartition ( const std::string  name,
const std::string  problem_for_rows,
bool  copy_rows,
const std::string  problem_for_cols,
bool  copy_cols,
int  verb = VERBOSE 
)

build indexing and partition problem inheriting indexing and partitioning from two other problems

Parameters
nameproblem name
problem_for_rowsproblem used to index rows
copy_rowsjust copy rows dofs
problem_for_colsproblem used to index cols
copy_colsjust copy cols dofs

If copy_rows/copy_cols is set to false only partition is copied between problems.

Definition at line 2166 of file ProblemsManager.cpp.

2168  {
2169  MoFEM::Interface &m_field = cOre;
2170  auto problems_ptr = m_field.get_problems();
2172 
2174  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
2175 
2177 
2178  // find p_miit
2179  ProblemByName &problems_by_name =
2180  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2181  ProblemByName::iterator p_miit = problems_by_name.find(name);
2182  if (p_miit == problems_by_name.end()) {
2183  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2184  "problem with name < %s > not defined (top tip check spelling)",
2185  name.c_str());
2186  }
2187  if (verb > QUIET)
2188  MOFEM_LOG("WORLD", Sev::inform)
2189  << p_miit->getName() << " from rows of " << problem_for_rows
2190  << " and columns of " << problem_for_cols;
2191 
2192  // find p_miit_row
2193  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
2194  if (p_miit_row == problems_by_name.end()) {
2195  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2196  "problem with name < %s > not defined (top tip check spelling)",
2197  problem_for_rows.c_str());
2198  }
2199  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
2200  p_miit_row->numeredRowDofsPtr;
2201 
2202  // find p_mit_col
2203  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
2204  if (p_miit_col == problems_by_name.end()) {
2205  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2206  "problem with name < %s > not defined (top tip check spelling)",
2207  problem_for_cols.c_str());
2208  }
2209  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
2210  p_miit_col->numeredColDofsPtr;
2211 
2212  bool copy[] = {copy_rows, copy_cols};
2213  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
2214  p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
2215 
2216  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
2217  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
2218  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
2219  dofs_col};
2220 
2221  for (int ss = 0; ss < 2; ss++) {
2222 
2223  // build indices
2224  *nb_local_dofs[ss] = 0;
2225  if (!copy[ss]) {
2226 
2227  // only copy indices which are belong to some elements if this problem
2228  std::vector<int> is_local, is_new;
2229 
2230  NumeredDofEntityByUId &dofs_by_uid =
2231  copied_dofs[ss]->get<Unique_mi_tag>();
2232  for (NumeredDofEntity_multiIndex::iterator dit =
2233  composed_dofs[ss]->begin();
2234  dit != composed_dofs[ss]->end(); dit++) {
2235  NumeredDofEntityByUId::iterator diit =
2236  dofs_by_uid.find((*dit)->getLocalUniqueId());
2237  if (diit == dofs_by_uid.end()) {
2238  SETERRQ(
2240  "data inconsistency, could not find dof in composite problem");
2241  }
2242  int part_number = (*diit)->getPart(); // get part number
2243  int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
2244  bool success;
2245  success = composed_dofs[ss]->modify(
2246  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2247  petsc_global_dof));
2248  if (!success) {
2249  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2250  "modification unsuccessful");
2251  }
2252  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2253  success = composed_dofs[ss]->modify(
2254  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
2255  if (!success) {
2256  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2257  "modification unsuccessful");
2258  }
2259  is_local.push_back(petsc_global_dof);
2260  }
2261  }
2262 
2263  AO ao;
2264  CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
2265  NULL, &ao);
2266 
2267  // apply local to global mapping
2268  is_local.resize(0);
2269  for (NumeredDofEntity_multiIndex::iterator dit =
2270  composed_dofs[ss]->begin();
2271  dit != composed_dofs[ss]->end(); dit++) {
2272  is_local.push_back((*dit)->getPetscGlobalDofIdx());
2273  }
2274  CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2275  int idx2 = 0;
2276  for (NumeredDofEntity_multiIndex::iterator dit =
2277  composed_dofs[ss]->begin();
2278  dit != composed_dofs[ss]->end(); dit++) {
2279  int part_number = (*dit)->getPart(); // get part number
2280  int petsc_global_dof = is_local[idx2++];
2281  bool success;
2282  success = composed_dofs[ss]->modify(
2283  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2284  petsc_global_dof));
2285  if (!success) {
2286  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2287  "modification unsuccessful");
2288  }
2289  }
2290 
2291  CHKERR AODestroy(&ao);
2292 
2293  } else {
2294 
2295  for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2296  dit != copied_dofs[ss]->end(); dit++) {
2297  std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2298  p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2299  new NumeredDofEntity((*dit)->getDofEntityPtr())));
2300  if (p.second) {
2301  (*nb_dofs[ss])++;
2302  }
2303  int dof_idx = (*dit)->getDofIdx();
2304  int part_number = (*dit)->getPart(); // get part number
2305  int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2306  if (part_number == (unsigned int)m_field.get_comm_rank()) {
2307  const bool success = composed_dofs[ss]->modify(
2308  p.first, NumeredDofEntity_part_and_all_indices_change(
2309  part_number, dof_idx, petsc_global_dof,
2310  (*nb_local_dofs[ss])++));
2311  if (!success) {
2312  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2313  "modification unsuccessful");
2314  }
2315  } else {
2316  const bool success = composed_dofs[ss]->modify(
2317  p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
2318  part_number, dof_idx, petsc_global_dof));
2319  if (!success) {
2320  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2321  "modification unsuccessful");
2322  }
2323  }
2324  }
2325  }
2326  }
2327 
2328  CHKERR printPartitionedProblem(&*p_miit, verb);
2329  CHKERR debugPartitionedProblem(&*p_miit, verb);
2330 
2332 }
static Index< 'p', 3 > p
@ QUIET
Definition: definitions.h:277
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.

◆ partitionFiniteElements()

MoFEMErrorCode MoFEM::ProblemsManager::partitionFiniteElements ( const std::string  name,
bool  part_from_moab = false,
int  low_proc = -1,
int  hi_proc = -1,
int  verb = VERBOSE 
)

partition finite elements

Function which partition finite elements based on dofs partitioning.
In addition it sets information about local row and cols dofs at given element on partition.

Parameters
nameproblem name

Definition at line 2416 of file ProblemsManager.cpp.

2419  {
2420  MoFEM::Interface &m_field = cOre;
2421  auto problems_ptr = m_field.get_problems();
2422  auto fe_ent_ptr = m_field.get_ents_finite_elements();
2424 
2426  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2427 
2428  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2429  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2430 
2431  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2432  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2433  "adjacencies not build");
2434 
2436  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2437 
2439  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2440  "problem not partitioned");
2441 
2442  if (low_proc == -1)
2443  low_proc = m_field.get_comm_rank();
2444  if (hi_proc == -1)
2445  hi_proc = m_field.get_comm_rank();
2446 
2447  // Find pointer to problem of given name
2449  auto &problems =
2450  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2451  ProblemByName::iterator p_miit = problems.find(name);
2452  if (p_miit == problems.end()) {
2453  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2454  "problem < %s > not found (top tip: check spelling)",
2455  name.c_str());
2456  }
2457 
2458  // Get reference on finite elements multi-index on the problem
2459  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2460  *p_miit->numeredFiniteElementsPtr;
2461 
2462  // Clear all elements and data, build it again
2463  problem_finite_elements.clear();
2464 
2465  // Check if dofs and columns are the same, i.e. structurally symmetric
2466  // problem
2467  bool do_cols_prob = true;
2468  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2469  do_cols_prob = false;
2470  }
2471 
2472  auto get_good_elems = [&]() {
2473  auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
2474  good_elems.reserve(fe_ent_ptr->size());
2475 
2476  const auto prb_bit = p_miit->getBitRefLevel();
2477  const auto prb_mask = p_miit->getMaskBitRefLevel();
2478 
2479  // Loop over all elements in database and if right element is there add it
2480  // to problem finite element multi-index
2481  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2482 
2483  // if element is not part of problem
2484  if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2485 
2486  const auto fe_bit = (*efit)->getBitRefLevel();
2487 
2488  // if entity is not problem refinement level
2489  if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit) == prb_bit)
2490  good_elems.emplace_back(efit);
2491  }
2492  }
2493 
2494  return good_elems;
2495  };
2496 
2497  auto good_elems = get_good_elems();
2498 
2499  auto numbered_good_elems_ptr =
2500  boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2501  numbered_good_elems_ptr->reserve(good_elems.size());
2502  for (auto &efit : good_elems)
2503  numbered_good_elems_ptr->emplace_back(NumeredEntFiniteElement(*efit));
2504 
2505  if (!do_cols_prob) {
2506  for (auto &fe : *numbered_good_elems_ptr) {
2507  if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2508  fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2509  }
2510  }
2511  }
2512 
2513  if (part_from_moab) {
2514  for (auto &fe : *numbered_good_elems_ptr) {
2515  // if partition is taken from moab partition
2516  int proc = fe.getPartProc();
2517  if (proc == -1 && fe.getEntType() == MBVERTEX)
2518  proc = fe.getOwnerProc();
2519  fe.part = proc;
2520  }
2521  }
2522 
2523  for (auto &fe : *numbered_good_elems_ptr) {
2524 
2526  CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2527 
2528  if (!part_from_moab) {
2529  std::vector<int> parts(m_field.get_comm_size(), 0);
2530  for (auto &dof_ptr : rows_view)
2531  parts[dof_ptr->pArt]++;
2532  std::vector<int>::iterator pos = max_element(parts.begin(), parts.end());
2533  const auto max_part = std::distance(parts.begin(), pos);
2534  fe.part = max_part;
2535  }
2536  }
2537 
2538  for (auto &fe : *numbered_good_elems_ptr) {
2539 
2540  auto check_fields_and_dofs = [&]() {
2541  if (!part_from_moab) {
2542 
2543  if (fe.getBitFieldIdRow().none() && m_field.get_comm_size() == 0) {
2544  MOFEM_LOG("WORLD", Sev::warning)
2545  << "At least one field has to be added to element row to "
2546  "determine partition of finite element. Check element " +
2547  boost::lexical_cast<std::string>(fe.getName());
2548  }
2549  }
2550 
2551  // Adding elements if row or column has DOFs, or there is no field set
2552  // to rows and columns. The second case would be used by elements
2553  // performing tasks which do not assemble matrices or vectors, but
2554  // evaluate fields or modify base functions.
2555 
2556  return (!fe.sPtr->getRowFieldEnts().empty() ||
2557  !fe.sPtr->getColFieldEnts().empty())
2558 
2559  ||
2560 
2561  (fe.getBitFieldIdRow().none() || fe.getBitFieldIdCol().none());
2562  };
2563 
2564  if (check_fields_and_dofs()) {
2565  // Add element to the problem
2566  auto p = problem_finite_elements.insert(
2567  boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
2568  &fe));
2569  if (!p.second)
2570  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2571  }
2572  }
2573 
2574  if (verb >= VERBOSE) {
2575  auto elements_on_rank =
2576  problem_finite_elements.get<Part_mi_tag>().equal_range(
2577  m_field.get_comm_rank());
2578  MOFEM_LOG("SYNC", Sev::verbose)
2579  << p_miit->getName() << " nb. elems "
2580  << std::distance(elements_on_rank.first, elements_on_rank.second);
2581  auto fe_ptr = m_field.get_finite_elements();
2582  for (auto &fe : *fe_ptr) {
2583  auto e_range =
2584  problem_finite_elements.get<Composite_Name_And_Part_mi_tag>()
2585  .equal_range(
2586  boost::make_tuple(fe->getName(), m_field.get_comm_rank()));
2587  MOFEM_LOG("SYNC", Sev::noisy)
2588  << "Element " << fe->getName() << " nb. elems "
2589  << std::distance(e_range.first, e_range.second);
2590  }
2591 
2592  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
2593  }
2594 
2597 }
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< FiniteElement_name_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, 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< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered

◆ partitionGhostDofs()

MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofs ( const std::string  name,
int  verb = VERBOSE 
)

determine ghost nodes

Parameters
nameproblem name

DOFs are ghost dofs if are used by elements on given partitition, but not owned by that partitition.

Definition at line 2599 of file ProblemsManager.cpp.

2600  {
2601  MoFEM::Interface &m_field = cOre;
2602  auto problems_ptr = m_field.get_problems();
2604 
2606  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2607  "partition of problem not build");
2609  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2610  "partitions finite elements not build");
2611 
2612  // get problem pointer
2613  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2614  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2615  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2616  name.c_str());
2617 
2618  // get reference to number of ghost dofs
2619  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2620  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2621  nb_row_ghost_dofs = 0;
2622  nb_col_ghost_dofs = 0;
2623 
2624  // do work if more than one processor
2625  if (m_field.get_comm_size() > 1) {
2626 
2628  ghost_idx_col_view;
2629 
2630  // get elements on this partition
2631  auto fe_range =
2632  p_miit->numeredFiniteElementsPtr->get<Part_mi_tag>().equal_range(
2633  m_field.get_comm_rank());
2634 
2635  // get dofs on elements which are not part of this partition
2636 
2637  struct Inserter {
2638  using Vec = std::vector<boost::shared_ptr<NumeredDofEntity>>;
2639  using It = Vec::iterator;
2640  It operator()(Vec &dofs_view, It &hint,
2641  boost::shared_ptr<NumeredDofEntity> &&dof) {
2642  dofs_view.emplace_back(dof);
2643  return dofs_view.end();
2644  }
2645  };
2646 
2647  // rows
2648  std::vector<boost::shared_ptr<NumeredDofEntity>> fe_vec_view;
2649  auto hint_r = ghost_idx_row_view.begin();
2650  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2651 
2652  fe_vec_view.clear();
2653  CHKERR EntFiniteElement::getDofView((*fe_ptr)->getRowFieldEnts(),
2654  *(p_miit->getNumeredRowDofsPtr()),
2655  fe_vec_view, Inserter());
2656 
2657  for (auto &dof_ptr : fe_vec_view) {
2658  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2659  hint_r = ghost_idx_row_view.emplace_hint(hint_r, dof_ptr);
2660  }
2661  }
2662  }
2663 
2664  // columns
2665  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2666 
2667  auto hint_c = ghost_idx_col_view.begin();
2668  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2669 
2670  fe_vec_view.clear();
2672  (*fe_ptr)->getColFieldEnts(), *(p_miit->getNumeredColDofsPtr()),
2673  fe_vec_view, Inserter());
2674 
2675  for (auto &dof_ptr : fe_vec_view) {
2676  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2677  hint_c = ghost_idx_col_view.emplace_hint(hint_c, dof_ptr);
2678  }
2679  }
2680  }
2681  }
2682 
2683  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2684  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2685 
2686  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2687  &ghost_idx_row_view, &ghost_idx_col_view};
2688  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2689  &p_miit->numeredRowDofsPtr->get<Unique_mi_tag>(),
2690  &p_miit->numeredColDofsPtr->get<Unique_mi_tag>()};
2691 
2692  int loop_size = 2;
2693  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2694  loop_size = 1;
2695  }
2696 
2697  // set local ghost dofs indices
2698  for (int ss = 0; ss != loop_size; ++ss) {
2699  for (auto &gid : *ghost_idx_view[ss]) {
2700  NumeredDofEntityByUId::iterator dof =
2701  dof_by_uid_no_const[ss]->find(gid->getLocalUniqueId());
2702  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2703  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2704  "inconsistent data, ghost dof already set");
2705  bool success = dof_by_uid_no_const[ss]->modify(
2706  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2707  if (PetscUnlikely(!success))
2708  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2709  "modification unsuccessful");
2710  (*nb_ghost_dofs[ss])++;
2711  }
2712  }
2713  if (loop_size == 1) {
2714  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2715  }
2716  }
2717 
2718  if (verb > QUIET) {
2719  MOFEM_LOG("SYNC", Sev::inform)
2720  << " FEs ghost dofs on problem " << p_miit->getName()
2721  << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2722  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2723  << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2724 
2725  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2726  }
2727 
2730 }
@ MOFEM_INVALID_DATA
Definition: definitions.h:128
const FTensor::Tensor2< T, Dim, Dim > Vec
int DofIdx
Index of DOF.
Definition: Types.hpp:29
static MoFEMErrorCode getDofView(const FE_ENTS &fe_ents_view, const MOFEM_DOFS &mofem_dofs, MOFEM_DOFS_VIEW &dofs_view, INSERTER &&inserter)

◆ partitionGhostDofsOnDistributedMesh()

MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh ( const std::string  name,
int  verb = VERBOSE 
)

determine ghost nodes on distributed meshes

Parameters
nameproblem name

It is very similar for partitionGhostDofs, however this explits that mesh is distributed.

DOFs are ghosted if are shered but not owned by partition.

Definition at line 2733 of file ProblemsManager.cpp.

2734  {
2735  MoFEM::Interface &m_field = cOre;
2736  auto problems_ptr = m_field.get_problems();
2738 
2740  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2741  "partition of problem not build");
2743  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2744  "partitions finite elements not build");
2745 
2746  // get problem pointer
2747  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2748  // find p_miit
2749  ProblemsByName &problems_set =
2750  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2751  ProblemsByName::iterator p_miit = problems_set.find(name);
2752 
2753  // get reference to number of ghost dofs
2754  // get number of local dofs
2755  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2756  &(p_miit->nbGhostDofsCol)};
2757  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2758  for (int ss = 0; ss != 2; ++ss) {
2759  (*nb_ghost_dofs[ss]) = 0;
2760  }
2761 
2762  // do work if more than one processor
2763  if (m_field.get_comm_size() > 1) {
2764  // determine if rows on columns are different from dofs on rows
2765  int loop_size = 2;
2766  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2767  loop_size = 1;
2768  }
2769 
2770  typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2771  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2772  p_miit->numeredColDofsPtr};
2773 
2774  // interate over dofs on rows and dofs on columns
2775  for (int ss = 0; ss != loop_size; ++ss) {
2776 
2777  // create dofs view by uid
2778  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2779 
2780  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2781  ghost_idx_view.reserve(std::distance(r.first, r.second));
2782  for (; r.first != r.second; ++r.first)
2783  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2784 
2785  auto cmp = [](auto a, auto b) {
2786  return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2787  };
2788  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2789 
2790  // intare over dofs which have negative local index
2791  for (auto gid_it : ghost_idx_view) {
2792  bool success = numered_dofs[ss]->modify(
2793  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2794  if (!success)
2795  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2796  "modification unsuccessful");
2797  ++(*nb_ghost_dofs[ss]);
2798  }
2799  }
2800  if (loop_size == 1) {
2801  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2802  }
2803  }
2804 
2805  if (verb > QUIET) {
2806  MOFEM_LOG("SYNC", Sev::inform)
2807  << " FEs ghost dofs on problem " << p_miit->getName()
2808  << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2809  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2810  << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2811 
2812  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2813  }
2814 
2817 }

◆ partitionMesh()

MoFEMErrorCode MoFEM::ProblemsManager::partitionMesh ( const Range &  ents,
const int  dim,
const int  adj_dim,
const int  n_parts,
Tag *  th_vertex_weights = nullptr,
Tag *  th_edge_weights = nullptr,
Tag *  th_part_weights = nullptr,
int  verb = VERBOSE,
const bool  debug = false 
)

Set partition tag to each finite element in the problem.

This will use one of the mesh partitioning programs available from PETSc See http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPartitioningType.html

Parameters
entsEntities to partition
dimDimension of entities to partition
adj_dimAdjacency dimension
n_partsNumber of partitions
verbVerbosity level
Returns
Error code

Definition at line 104 of file ProblemsManager.cpp.

107  {
108  MoFEM::Interface &m_field = cOre;
110 
111  // get layout
112  int rstart, rend, nb_elems;
113  {
114  PetscLayout layout;
115  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
116  CHKERR PetscLayoutSetBlockSize(layout, 1);
117  CHKERR PetscLayoutSetSize(layout, ents.size());
118  CHKERR PetscLayoutSetUp(layout);
119  CHKERR PetscLayoutGetSize(layout, &nb_elems);
120  CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
121  CHKERR PetscLayoutDestroy(&layout);
122  if (verb >= VERBOSE) {
123  MOFEM_LOG("SYNC", Sev::inform)
124  << "Finite elements in problem: row lower " << rstart << " row upper "
125  << rend << " nb. elems " << nb_elems << " ( " << ents.size() << " )";
127  }
128  }
129 
130  std::vector<EntityHandle> weight_ents;
131  weight_ents.reserve(rend - rstart + 1);
132 
133  struct AdjBridge {
134  EntityHandle ent;
135  std::vector<int> adj;
136  AdjBridge(const EntityHandle ent, std::vector<int> &adj)
137  : ent(ent), adj(adj) {}
138  };
139 
140  typedef multi_index_container<
141  AdjBridge,
142  indexed_by<
143 
144  hashed_unique<member<AdjBridge, EntityHandle, &AdjBridge::ent>>
145 
146  >>
147  AdjBridgeMap;
148 
149  Range all_dim_ents;
150  CHKERR m_field.get_moab().get_adjacencies(ents, adj_dim, true, all_dim_ents,
151  moab::Interface::UNION);
152 
153  AdjBridgeMap adj_bridge_map;
154  auto hint = adj_bridge_map.begin();
155  std::vector<int> adj;
156  for (auto ent : all_dim_ents) {
157  Range adj_ents;
158  CHKERR m_field.get_moab().get_adjacencies(&ent, 1, dim, false, adj_ents);
159  adj_ents = intersect(adj_ents, ents);
160  adj.clear();
161  adj.reserve(adj_ents.size());
162  for (auto a : adj_ents)
163  adj.emplace_back(ents.index(a));
164  hint = adj_bridge_map.emplace_hint(hint, ent, adj);
165  }
166 
167  int *_i;
168  int *_j;
169  {
170  const int nb_loc_elements = rend - rstart;
171  std::vector<int> i(nb_loc_elements + 1, 0), j;
172  {
173  std::vector<int> row_adj;
174  Range::iterator fe_it;
175  int ii, jj;
176  size_t max_row_size;
177  for (
178 
179  fe_it = ents.begin(), ii = 0, jj = 0, max_row_size = 0;
180 
181  fe_it != ents.end(); ++fe_it, ++ii) {
182 
183  if (ii < rstart)
184  continue;
185  if (ii >= rend)
186  break;
187 
188  if (m_field.get_moab().type_from_handle(*fe_it) == MBENTITYSET) {
189  SETERRQ(
190  PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
191  "not yet implemented, don't know what to do for meshset element");
192  } else {
193 
194  Range dim_ents;
195  CHKERR m_field.get_moab().get_adjacencies(&*fe_it, 1, adj_dim, false,
196  dim_ents);
197  dim_ents = intersect(dim_ents, all_dim_ents);
198 
199  row_adj.clear();
200  for (auto e : dim_ents) {
201  auto adj_it = adj_bridge_map.find(e);
202  if (adj_it != adj_bridge_map.end()) {
203 
204  for (const auto idx : adj_it->adj)
205  row_adj.push_back(idx);
206 
207  } else
208  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
209  "Entity not found");
210  }
211 
212  std::sort(row_adj.begin(), row_adj.end());
213  auto end = std::unique(row_adj.begin(), row_adj.end());
214 
215  size_t row_size = std::distance(row_adj.begin(), end);
216  max_row_size = std::max(max_row_size, row_size);
217  if (j.capacity() < (nb_loc_elements - jj) * max_row_size)
218  j.reserve(nb_loc_elements * max_row_size);
219 
220  i[jj] = j.size();
221  auto diag = ents.index(*fe_it);
222  for (auto it = row_adj.begin(); it != end; ++it)
223  if (*it != diag)
224  j.push_back(*it);
225  }
226 
227  ++jj;
228 
229  if (th_vertex_weights != NULL)
230  weight_ents.push_back(*fe_it);
231  }
232 
233  i[jj] = j.size();
234  }
235 
236  CHKERR PetscMalloc(i.size() * sizeof(int), &_i);
237  CHKERR PetscMalloc(j.size() * sizeof(int), &_j);
238  copy(i.begin(), i.end(), _i);
239  copy(j.begin(), j.end(), _j);
240  }
241 
242  // get weights
243  int *vertex_weights = NULL;
244  if (th_vertex_weights != NULL) {
245  CHKERR PetscMalloc(weight_ents.size() * sizeof(int), &vertex_weights);
246  CHKERR m_field.get_moab().tag_get_data(*th_vertex_weights,
247  &*weight_ents.begin(),
248  weight_ents.size(), vertex_weights);
249  }
250 
251  {
252  Mat Adj;
253  // Adjacency matrix used to partition problems, f.e. METIS
254  CHKERR MatCreateMPIAdj(m_field.get_comm(), rend - rstart, nb_elems, _i, _j,
255  PETSC_NULL, &Adj);
256  CHKERR MatSetOption(Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
257 
258  if (debug) {
259  Mat A;
260  CHKERR MatConvert(Adj, MATMPIAIJ, MAT_INITIAL_MATRIX, &A);
261  CHKERR MatView(A, PETSC_VIEWER_DRAW_WORLD);
262  std::string wait;
263  std::cin >> wait;
264  CHKERR MatDestroy(&A);
265  }
266 
267  // run pets to do partitioning
268  MOFEM_LOG("WORLD", Sev::verbose) << "Start";
269 
270  MatPartitioning part;
271  IS is;
272  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
273 
274  CHKERR MatPartitioningSetAdjacency(part, Adj);
275  CHKERR MatPartitioningSetFromOptions(part);
276  CHKERR MatPartitioningSetNParts(part, n_parts);
277  if (th_vertex_weights != NULL) {
278  CHKERR MatPartitioningSetVertexWeights(part, vertex_weights);
279  }
280  PetscBool same;
281  PetscObjectTypeCompare((PetscObject)part, MATPARTITIONINGPARMETIS, &same);
282  if (same) {
283 #ifdef PARMETIS
284  CHKERR MatPartitioningApply_Parmetis_MoFEM(part, &is);
285 #endif
286  } else {
287  CHKERR MatPartitioningApply(part, &is);
288  }
289 
290  MOFEM_LOG("WORLD", Sev::verbose) << "End";
291 
292  // gather
293  IS is_gather, is_num, is_gather_num;
294  CHKERR ISAllGather(is, &is_gather);
295  CHKERR ISPartitioningToNumbering(is, &is_num);
296  CHKERR ISAllGather(is_num, &is_gather_num);
297 
298  const int *part_number, *gids;
299  CHKERR ISGetIndices(is_gather, &part_number);
300  CHKERR ISGetIndices(is_gather_num, &gids);
301 
302  // set partition tag and gid tag to entities
303  ParallelComm *pcomm = ParallelComm::get_pcomm(
304  &m_field.get_moab(), m_field.get_basic_entity_data_ptr()->pcommID);
305  Tag part_tag = pcomm->part_tag();
306  CHKERR m_field.get_moab().tag_set_data(part_tag, ents, part_number);
307  auto get_gid_tag = [&]() {
308  Tag gid_tag;
309  rval = m_field.get_moab().tag_get_handle(GLOBAL_ID_TAG_NAME, gid_tag);
310  if (rval != MB_SUCCESS) {
311  const int zero = 0;
312  CHKERR m_field.get_moab().tag_get_handle(
313  GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag,
314  MB_TAG_DENSE | MB_TAG_CREAT, &zero);
315  }
316  return gid_tag;
317  };
318  Tag gid_tag = get_gid_tag();
319 
320  std::map<int, Range> parts_ents;
321  {
322  // get entities on each part
323  Range::iterator eit = ents.begin();
324  for (int ii = 0; eit != ents.end(); eit++, ii++) {
325  parts_ents[part_number[ii]].insert(*eit);
326  }
327  Range tagged_sets;
328  CHKERR m_field.get_moab().get_entities_by_type_and_tag(
329  0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
330  moab::Interface::UNION);
331  if (!tagged_sets.empty())
332  CHKERR m_field.get_moab().tag_delete_data(part_tag, tagged_sets);
333 
334  if (n_parts > (int)tagged_sets.size()) {
335  // too few partition sets - create missing ones
336  int num_new = n_parts - tagged_sets.size();
337  for (int i = 0; i < num_new; i++) {
338  EntityHandle new_set;
339  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, new_set);
340  tagged_sets.insert(new_set);
341  }
342  } else if (n_parts < (int)tagged_sets.size()) {
343  // too many partition sets - delete extras
344  int num_del = tagged_sets.size() - n_parts;
345  for (int i = 0; i < num_del; i++) {
346  EntityHandle old_set = tagged_sets.pop_back();
347  CHKERR m_field.get_moab().delete_entities(&old_set, 1);
348  }
349  }
350  // write a tag to those sets denoting they're partition sets, with a value
351  // of the proc number
352  std::vector<int> dum_ids(n_parts);
353  for (int i = 0; i < n_parts; i++)
354  dum_ids[i] = i;
355  CHKERR m_field.get_moab().tag_set_data(part_tag, tagged_sets,
356  &*dum_ids.begin());
357  CHKERR m_field.get_moab().clear_meshset(tagged_sets);
358 
359  // get lower dimension entities on each part
360  for (int pp = 0; pp != n_parts; pp++) {
361  Range dim_ents = parts_ents[pp].subset_by_dimension(dim);
362  for (int dd = dim - 1; dd != -1; dd--) {
363  Range adj_ents;
364  if (dim > 0) {
365  CHKERR m_field.get_moab().get_adjacencies(
366  dim_ents, dd, true, adj_ents, moab::Interface::UNION);
367  } else {
368  CHKERR m_field.get_moab().get_connectivity(dim_ents, adj_ents,
369  true);
370  }
371  parts_ents[pp].merge(adj_ents);
372  }
373  }
374  for (int pp = 1; pp != n_parts; pp++) {
375  for (int ppp = 0; ppp != pp; ppp++) {
376  parts_ents[pp] = subtract(parts_ents[pp], parts_ents[ppp]);
377  }
378  }
379  if (debug) {
380  for (int rr = 0; rr != m_field.get_comm_size(); rr++) {
381  ostringstream ss;
382  ss << "out_part_" << rr << ".vtk";
383  EntityHandle meshset;
384  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
385  CHKERR m_field.get_moab().add_entities(meshset, parts_ents[rr]);
386  CHKERR m_field.get_moab().write_file(ss.str().c_str(), "VTK", "",
387  &meshset, 1);
388  CHKERR m_field.get_moab().delete_entities(&meshset, 1);
389  }
390  }
391  for (int pp = 0; pp != n_parts; pp++) {
392  CHKERR m_field.get_moab().add_entities(tagged_sets[pp], parts_ents[pp]);
393  }
394 
395  // set gid to lower dimension entities
396  for (int dd = 0; dd <= dim; dd++) {
397  int gid = 1; // moab indexing from 1
398  for (int pp = 0; pp != n_parts; pp++) {
399  Range dim_ents = parts_ents[pp].subset_by_dimension(dd);
400  // std::endl;
401  for (Range::iterator eit = dim_ents.begin(); eit != dim_ents.end();
402  eit++) {
403  if (dd > 0) {
404  CHKERR m_field.get_moab().tag_set_data(part_tag, &*eit, 1, &pp);
405  }
406  CHKERR m_field.get_moab().tag_set_data(gid_tag, &*eit, 1, &gid);
407  gid++;
408  }
409  }
410  }
411  }
412 
413  CHKERR ISRestoreIndices(is_gather, &part_number);
414  CHKERR ISRestoreIndices(is_gather_num, &gids);
415  CHKERR ISDestroy(&is_num);
416  CHKERR ISDestroy(&is_gather_num);
417  CHKERR ISDestroy(&is_gather);
418  CHKERR ISDestroy(&is);
419  CHKERR MatPartitioningDestroy(&part);
420  CHKERR MatDestroy(&Adj);
421  }
422 
424 
426 }
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
static const bool debug
const int dim
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.

◆ partitionProblem()

MoFEMErrorCode MoFEM::ProblemsManager::partitionProblem ( const std::string  name,
int  verb = VERBOSE 
)

partition problem dofs (collective)

Parameters
nameproblem name

Definition at line 1954 of file ProblemsManager.cpp.

1955  {
1956  MoFEM::Interface &m_field = cOre;
1957  auto problems_ptr = m_field.get_problems();
1959 
1960  MOFEM_LOG("WORLD", Sev::noisy) << "Partition problem " << name;
1961 
1962  using NumeredDofEntitysByIdx =
1964  using ProblemsByName = Problem_multiIndex::index<Problem_mi_tag>::type;
1965 
1966  // Find problem pointer by name
1967  auto &problems_set =
1968  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1969  auto p_miit = problems_set.find(name);
1970  if (p_miit == problems_set.end()) {
1971  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1972  "problem with name %s not defined (top tip check spelling)",
1973  name.c_str());
1974  }
1975  int nb_dofs_row = p_miit->getNbDofsRow();
1976 
1977  if (m_field.get_comm_size() != 1) {
1978 
1979  if (!(cOre.getBuildMoFEM() & (1 << 0)))
1980  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1981  if (!(cOre.getBuildMoFEM() & (1 << 1)))
1982  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1983  if (!(cOre.getBuildMoFEM() & (1 << 2)))
1984  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1985  "entFEAdjacencies not build");
1986  if (!(cOre.getBuildMoFEM() & (1 << 3)))
1987  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1988 
1989  Mat Adj;
1990  CHKERR m_field.getInterface<MatrixManager>()
1991  ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1992 
1993  int m, n;
1994  CHKERR MatGetSize(Adj, &m, &n);
1995  if (verb > VERY_VERBOSE)
1996  MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1997 
1998  // partitioning
1999  MatPartitioning part;
2000  IS is;
2001  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
2002  CHKERR MatPartitioningSetAdjacency(part, Adj);
2003  CHKERR MatPartitioningSetFromOptions(part);
2004  CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
2005  CHKERR MatPartitioningApply(part, &is);
2006  if (verb > VERY_VERBOSE)
2007  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
2008 
2009  // gather
2010  IS is_gather, is_num, is_gather_num;
2011  CHKERR ISAllGather(is, &is_gather);
2012  CHKERR ISPartitioningToNumbering(is, &is_num);
2013  CHKERR ISAllGather(is_num, &is_gather_num);
2014  const int *part_number, *petsc_idx;
2015  CHKERR ISGetIndices(is_gather, &part_number);
2016  CHKERR ISGetIndices(is_gather_num, &petsc_idx);
2017  int size_is_num, size_is_gather;
2018  CHKERR ISGetSize(is_gather, &size_is_gather);
2019  if (size_is_gather != (int)nb_dofs_row)
2020  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
2021  "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
2022 
2023  CHKERR ISGetSize(is_num, &size_is_num);
2024  if (size_is_num != (int)nb_dofs_row)
2025  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
2026  "data inconsistency %d != %d", size_is_num, nb_dofs_row);
2027 
2028  bool square_matrix = false;
2029  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
2030  square_matrix = true;
2031 
2032  // if (!square_matrix) {
2033  // // FIXME: This is for back compatibility, if deprecate interface function
2034  // // build interfaces is removed, this part of the code will be obsolete
2035  // auto mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().begin();
2036  // auto hi_mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().end();
2037  // auto mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().begin();
2038  // auto hi_mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().end();
2039  // for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
2040  // if (mit_col == hi_mit_col) {
2041  // SETERRQ(
2042  // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2043  // "check finite element definition, nb. of rows is not equal to "
2044  // "number for columns");
2045  // }
2046  // if (mit_row->get()->getLocalUniqueId() !=
2047  // mit_col->get()->getLocalUniqueId()) {
2048  // SETERRQ(
2049  // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2050  // "check finite element definition, nb. of rows is not equal to "
2051  // "number for columns");
2052  // }
2053  // }
2054  // }
2055 
2056  auto number_dofs = [&](auto &dofs_idx, auto &counter) {
2058  for (auto miit_dofs_row = dofs_idx.begin();
2059  miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
2060  const int part = part_number[(*miit_dofs_row)->dofIdx];
2061  if (part == (unsigned int)m_field.get_comm_rank()) {
2062  const bool success = dofs_idx.modify(
2063  miit_dofs_row,
2064  NumeredDofEntity_part_and_indices_change(
2065  part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
2066  if (!success) {
2067  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2068  "modification unsuccessful");
2069  }
2070  } else {
2071  const bool success = dofs_idx.modify(
2072  miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
2073  part, petsc_idx[(*miit_dofs_row)->dofIdx]));
2074  if (!success) {
2075  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2076  "modification unsuccessful");
2077  }
2078  }
2079  }
2081  };
2082 
2083  // Set petsc global indices
2084  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
2085  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
2086  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2087  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2088  nb_row_local_dofs = 0;
2089  nb_row_ghost_dofs = 0;
2090 
2091  CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
2092 
2093  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2094  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2095  if (square_matrix) {
2096  nb_col_local_dofs = nb_row_local_dofs;
2097  nb_col_ghost_dofs = nb_row_ghost_dofs;
2098  } else {
2099  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2100  const_cast<NumeredDofEntitysByIdx &>(
2101  p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
2102  nb_col_local_dofs = 0;
2103  nb_col_ghost_dofs = 0;
2104  CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
2105  }
2106 
2107  CHKERR ISRestoreIndices(is_gather, &part_number);
2108  CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
2109  CHKERR ISDestroy(&is_num);
2110  CHKERR ISDestroy(&is_gather_num);
2111  CHKERR ISDestroy(&is_gather);
2112  CHKERR ISDestroy(&is);
2113  CHKERR MatPartitioningDestroy(&part);
2114  CHKERR MatDestroy(&Adj);
2115  CHKERR printPartitionedProblem(&*p_miit, verb);
2116  } else {
2117 
2118  auto number_dofs = [&](auto &dof_idx, auto &counter) {
2120  for (auto miit_dofs_row = dof_idx.begin();
2121  miit_dofs_row != dof_idx.end(); miit_dofs_row++) {
2122  const bool success = dof_idx.modify(
2123  miit_dofs_row,
2124  NumeredDofEntity_part_and_indices_change(0, counter, counter));
2125  ++counter;
2126  if (!success) {
2127  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2128  "modification unsuccessful");
2129  }
2130  }
2132  };
2133 
2134  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
2135  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
2136  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2137  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2138  nb_row_local_dofs = 0;
2139  nb_row_ghost_dofs = 0;
2140 
2141  CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
2142 
2143  bool square_matrix = false;
2144  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
2145  square_matrix = true;
2146 
2147  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2148  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2149  if (square_matrix) {
2150  nb_col_local_dofs = nb_row_local_dofs;
2151  nb_col_ghost_dofs = nb_row_ghost_dofs;
2152  } else {
2153  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2154  const_cast<NumeredDofEntitysByIdx &>(
2155  p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
2156  nb_col_local_dofs = 0;
2157  nb_col_ghost_dofs = 0;
2158  CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
2159  }
2160  }
2161 
2162  cOre.getBuildMoFEM() |= 1 << 4;
2164 }
static Index< 'm', 3 > m
static Index< 'n', 3 > n
@ VERY_VERBOSE
Definition: definitions.h:279
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

◆ partitionSimpleProblem()

MoFEMErrorCode MoFEM::ProblemsManager::partitionSimpleProblem ( const std::string  name,
int  verb = VERBOSE 
)

partition problem dofs

Parameters
nameproblem name

Definition at line 1808 of file ProblemsManager.cpp.

1809  {
1810 
1811  MoFEM::Interface &m_field = cOre;
1812  auto problems_ptr = m_field.get_problems();
1814 
1816  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1817  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1818  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1819  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1820  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1822  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1823  MOFEM_LOG("WORLD", Sev::verbose) << "Simple partition problem " << name;
1824 
1825  // find p_miit
1827  ProblemByName &problems_set =
1828  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1829  ProblemByName::iterator p_miit = problems_set.find(name);
1830  if (p_miit == problems_set.end()) {
1831  SETERRQ1(PETSC_COMM_SELF, 1,
1832  "problem < %s > is not found (top tip: check spelling)",
1833  name.c_str());
1834  }
1835  typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1836  Idx_mi_tag>::type NumeredDofEntitysByIdx;
1837  NumeredDofEntitysByIdx &dofs_row_by_idx =
1838  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>();
1839  NumeredDofEntitysByIdx &dofs_col_by_idx =
1840  p_miit->numeredColDofsPtr->get<Idx_mi_tag>();
1841  boost::multi_index::index<NumeredDofEntity_multiIndex,
1842  Idx_mi_tag>::type::iterator miit_row,
1843  hi_miit_row;
1844  boost::multi_index::index<NumeredDofEntity_multiIndex,
1845  Idx_mi_tag>::type::iterator miit_col,
1846  hi_miit_col;
1847  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1848  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1849  nb_row_local_dofs = 0;
1850  nb_row_ghost_dofs = 0;
1851  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1852  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1853  nb_col_local_dofs = 0;
1854  nb_col_ghost_dofs = 0;
1855 
1856  bool square_matrix = false;
1857  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1858  square_matrix = true;
1859  }
1860 
1861  // get row range of local indices
1862  PetscLayout layout_row;
1863  const int *ranges_row;
1864 
1865  DofIdx nb_dofs_row = dofs_row_by_idx.size();
1866  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1867  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1868  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1869  CHKERR PetscLayoutSetUp(layout_row);
1870  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1871  // get col range of local indices
1872  PetscLayout layout_col;
1873  const int *ranges_col;
1874  if (!square_matrix) {
1875  DofIdx nb_dofs_col = dofs_col_by_idx.size();
1876  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1877  CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1878  CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1879  CHKERR PetscLayoutSetUp(layout_col);
1880  CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1881  }
1882  for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1883  part++) {
1884  miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1885  hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1886  if (std::distance(miit_row, hi_miit_row) !=
1887  ranges_row[part + 1] - ranges_row[part]) {
1888  SETERRQ4(
1889  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1890  "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1891  "rstart (%d != %d - %d = %d) ",
1892  std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1893  ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1894  }
1895  // loop rows
1896  for (; miit_row != hi_miit_row; miit_row++) {
1897  bool success = dofs_row_by_idx.modify(
1898  miit_row,
1899  NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1900  if (!success)
1901  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1902  "modification unsuccessful");
1903  if (part == (unsigned int)m_field.get_comm_rank()) {
1904  success = dofs_row_by_idx.modify(
1905  miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1906  if (!success)
1907  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1908  "modification unsuccessful");
1909  }
1910  }
1911  if (!square_matrix) {
1912  miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1913  hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1914  if (std::distance(miit_col, hi_miit_col) !=
1915  ranges_col[part + 1] - ranges_col[part]) {
1916  SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1917  "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1918  "rend - "
1919  "rstart (%d != %d - %d = %d) ",
1920  std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1921  ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1922  }
1923  // loop cols
1924  for (; miit_col != hi_miit_col; miit_col++) {
1925  bool success = dofs_col_by_idx.modify(
1926  miit_col, NumeredDofEntity_part_and_glob_idx_change(
1927  part, (*miit_col)->dofIdx));
1928  if (!success)
1929  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1930  "modification unsuccessful");
1931  if (part == (unsigned int)m_field.get_comm_rank()) {
1932  success = dofs_col_by_idx.modify(
1933  miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1934  if (!success)
1935  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1936  "modification unsuccessful");
1937  }
1938  }
1939  }
1940  }
1941  CHKERR PetscLayoutDestroy(&layout_row);
1942  if (!square_matrix) {
1943  CHKERR PetscLayoutDestroy(&layout_col);
1944  }
1945  if (square_matrix) {
1946  nb_col_local_dofs = nb_row_local_dofs;
1947  nb_col_ghost_dofs = nb_row_ghost_dofs;
1948  }
1949  CHKERR printPartitionedProblem(&*p_miit, verb);
1952 }

◆ removeDofsOnEntities()

MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities ( const std::string  problem_name,
const std::string  field_name,
const Range  ents,
const int  lo_coeff = 0,
const int  hi_coeff = MAX_DOFS_ON_ENTITY,
int  verb = VERBOSE,
const bool  debug = false 
)

Remove DOFs from problem.

Remove DOFs from problem which are on entities on the given range and given field name. On the finite element level, DOFs can be still accessed however local PETSc indices and global PETSc indices are marked with the index -1.

Note
If the index is marked -1 it is not assembled and dropped by VecSetValues and MatSetValues.
Todo:
Not yet implemented update for AO maps and IS ranges if removed entities in composite problem or sub-problem
Parameters
problem_namename of the problem
field_namename of the field
entsentities on which DOFs are removed
lo_coefflow dof coefficient (rank)
hi_coeffhigh dof coefficient (rank)
verbverbosity level
debugto debug and seek for inconsistencies set to true
Returns
MoFEMErrorCode

Definition at line 2856 of file ProblemsManager.cpp.

2859  {
2860 
2861  MoFEM::Interface &m_field = cOre;
2863 
2864  const Problem *prb_ptr;
2865  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2866 
2867  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {prb_ptr->numeredRowDofsPtr,
2868  nullptr};
2869  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2870  numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2871 
2872  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2873  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2874  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2875 
2876  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2877  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2878  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2879  const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2880  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2881  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2882 
2883  for (int s = 0; s != 2; ++s)
2884  if (numered_dofs[s]) {
2885 
2886  typedef multi_index_container<
2887 
2888  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2889 
2890  >
2891  NumeredDofEntity_it_view_multiIndex;
2892 
2893  const auto bit_number = m_field.get_field_bit_number(field_name);
2894  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2895 
2896  // Set -1 to global and local dofs indices
2897  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2898  ++pit) {
2899  auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2900  DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2901  auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2902  DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2903 
2904  for (; lo != hi; ++lo)
2905  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2906  (*lo)->getDofCoeffIdx() <= hi_coeff)
2907  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2908  }
2909 
2910  if (verb > QUIET) {
2911  for (auto &dof : dofs_it_view)
2912  MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2913 
2914  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
2915  }
2916 
2917  // set negative index
2918  auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2919  for (auto dit : dofs_it_view) {
2920  bool success = numered_dofs[s]->modify(dit, mod);
2921  if (!success)
2922  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2923  "can not set negative indices");
2924  }
2925 
2926  // create weak view
2927  std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2928  dosf_weak_view.reserve(dofs_it_view.size());
2929  for (auto dit : dofs_it_view)
2930  dosf_weak_view.push_back(*dit);
2931 
2932  if (verb >= NOISY)
2933  MOFEM_LOG_C("SYNC", Sev::noisy,
2934  "Number of DOFs in multi-index %d and to delete %d\n",
2935  numered_dofs[s]->size(), dofs_it_view.size());
2936 
2937  // erase dofs from problem
2938  for (auto weak_dit : dosf_weak_view)
2939  if (auto dit = weak_dit.lock()) {
2940  numered_dofs[s]->erase(dit->getLocalUniqueId());
2941  }
2942 
2943  if (verb >= NOISY)
2944  MOFEM_LOG_C("SYNC", Sev::noisy,
2945  "Number of DOFs in multi-index after delete %d\n",
2946  numered_dofs[s]->size());
2947 
2948  // get current number of ghost dofs
2949  int nb_local_dofs = 0;
2950  int nb_ghost_dofs = 0;
2951  for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2952  dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2953  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2954  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2955  ++nb_local_dofs;
2956  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2957  ++nb_ghost_dofs;
2958  else
2959  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2960  "Imposible case. You could run problem on no distributed "
2961  "mesh. That is not implemented.");
2962  }
2963 
2964  // get indices
2965  auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2966  const int nb_dofs = numered_dofs[s]->size();
2967  indices.clear();
2968  indices.reserve(nb_dofs);
2969  for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2970  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2971  bool add = true;
2972  if (only_local) {
2973  if ((*dit)->getPetscLocalDofIdx() < 0 ||
2974  (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2975  add = false;
2976  }
2977  }
2978  if (add)
2979  indices.push_back(decltype(tag)::get_index(dit));
2980  }
2981  };
2982 
2983  auto get_indices_by_uid = [&](auto tag, auto &indices) {
2984  const int nb_dofs = numered_dofs[s]->size();
2985  indices.clear();
2986  indices.reserve(nb_dofs);
2987  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2988  ++dit)
2989  indices.push_back(decltype(tag)::get_index(dit));
2990  };
2991 
2992  auto concatenate_dofs = [&](auto tag, auto &indices,
2993  const auto local_only) {
2995  get_indices_by_tag(tag, indices, local_only);
2996 
2997  AO ao;
2998  if (local_only)
2999  CHKERR AOCreateMapping(m_field.get_comm(), indices.size(),
3000  &*indices.begin(), PETSC_NULL, &ao);
3001  else
3002  CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
3003  &*indices.begin(), PETSC_NULL, &ao);
3004 
3005  get_indices_by_uid(tag, indices);
3006  CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3007  CHKERR AODestroy(&ao);
3009  };
3010 
3011  // set indices index
3012  auto set_concatinated_indices = [&]() {
3013  std::vector<int> global_indices;
3014  std::vector<int> local_indices;
3016  CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3017  CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3018  auto gi = global_indices.begin();
3019  auto li = local_indices.begin();
3020  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3021  ++dit) {
3022  auto mod = NumeredDofEntity_part_and_all_indices_change(
3023  (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3024  bool success = numered_dofs[s]->modify(dit, mod);
3025  if (!success)
3026  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3027  "can not set negative indices");
3028  ++gi;
3029  ++li;
3030  }
3032  };
3033  CHKERR set_concatinated_indices();
3034 
3035  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3036  m_field.get_comm());
3037  *(local_nbdof_ptr[s]) = nb_local_dofs;
3038  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3039 
3040  if (debug)
3041  for (auto dof : (*numered_dofs[s])) {
3042  if (dof->getPetscGlobalDofIdx() < 0) {
3043  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3044  "Negative global idx");
3045  }
3046  if (dof->getPetscLocalDofIdx() < 0) {
3047  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3048  "Negative local idx");
3049  }
3050  }
3051 
3052  } else {
3053 
3054  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3055  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3056  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3057  }
3058 
3059  if (verb > QUIET) {
3060  MOFEM_LOG_C("SYNC", Sev::inform,
3061  "removed ents on rank %d from problem %s dofs [ %d / %d "
3062  "(before %d / "
3063  "%d) local, %d / %d (before %d / %d) "
3064  "ghost, %d / %d (before %d / %d) global]",
3065  m_field.get_comm_rank(), prb_ptr->getName().c_str(),
3066  prb_ptr->getNbLocalDofsRow(), prb_ptr->getNbLocalDofsCol(),
3067  nb_init_row_dofs, nb_init_col_dofs,
3068  prb_ptr->getNbGhostDofsRow(), prb_ptr->getNbGhostDofsCol(),
3069  nb_init_ghost_row_dofs, nb_init_ghost_col_dofs,
3070  prb_ptr->getNbDofsRow(), prb_ptr->getNbDofsCol(),
3071  nb_init_loc_row_dofs, nb_init_loc_col_dofs);
3072  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
3073  }
3074 
3076 }
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:315
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem