v0.10.0
Files | Classes | Functions
ProblemsManager

Adding and managing problems. More...

Collaboration diagram for ProblemsManager:

Files

file  ProblemsManager.cpp
 Managing complexities for problem.
 
file  ProblemsManager.hpp
 Interface managing problemsManaging problems, build and partitioning.
 

Classes

struct  MoFEM::ProblemsManager
 Problem manager is used to build and partition problems. More...
 

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 problemThis 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 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)Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities. 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)Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities. 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 elementsFunction which partition finite elements based on dofs partitioning.
In addition it sets information about local row and cols dofs at given element on partition. 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 problemIn 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 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 problemRemove 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. 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 1523 of file ProblemsManager.cpp.

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

◆ 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
Examples
build_problems.cpp, continuity_check_on_contact_prism_side_ele.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, and quad_polynomial_approximation.cpp.

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 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
#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 * get_problem(const std::string &problem_name) const =0
Get the problem object.
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221

◆ 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 ProblemManagerFunctionBegin
Deprecated interface functions.
PetscLogEvent MOFEM_EVENT_ProblemsManager
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
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
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:305
virtual MPI_Comm & get_comm() const =0
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340

◆ 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.

Examples
remove_entities_from_problem.cpp.

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 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
#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
MoFEMErrorCode 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,...
virtual const Problem * get_problem(const std::string &problem_name) const =0
Get the problem object.
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221

◆ 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 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
PetscLogEvent MOFEM_EVENT_ProblemsManager
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual int get_comm_size() const =0
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getGlobalUniqueId > > > > DofEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< 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
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
virtual int get_comm_rank() const =0
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
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:305
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
MoFEMErrorCode getOptions()
virtual MPI_Comm & get_comm() const =0
uint128_t UId
Unique Id.
Definition: Types.hpp:42
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
static UId getHiBitNumberUId(const FieldBitNumber bit_number)

◆ 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  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->rowIs));
1442  // CHKERR ISSort(out_problem_it->getSubData()->rowIs);
1443  out_problem_it->getSubData()->rowMap = ao;
1444  CHKERR PetscObjectReference((PetscObject)ao);
1445  } else {
1446  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->colIs));
1447  // CHKERR ISSort(out_problem_it->getSubData()->colIs);
1448  out_problem_it->getSubData()->colMap = ao;
1449  CHKERR PetscObjectReference((PetscObject)ao);
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  CHKERR PetscObjectReference(
1514  (PetscObject)out_problem_it->getSubData()->rowMap);
1515  }
1516 
1517  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1518  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1519 
1521 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
virtual int get_comm_rank() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
#define CHKERR
Inline error check.
Definition: definitions.h:604
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
static UId getHiBitNumberUId(const FieldBitNumber 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 2821 of file ProblemsManager.cpp.

2823  {
2824  MoFEM::Interface &m_field = cOre;
2825  const Problem *problem_ptr;
2827 
2828  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2829  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2830  auto fit =
2831  problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2832  .lower_bound(fe_name);
2833  auto hi_fe_it =
2834  problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2835  .upper_bound(fe_name);
2836  std::vector<EntityHandle> fe_vec;
2837  fe_vec.reserve(std::distance(fit, hi_fe_it));
2838  for (; fit != hi_fe_it; fit++)
2839  fe_vec.push_back(fit->get()->getEnt());
2840  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2841  fe_vec.size());
2843 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual const Problem * get_problem(const std::string &problem_name) const =0
Get the problem object.

◆ getProblemElementsLayout()

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

Get layout of elements in the problemIn 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 2846 of file ProblemsManager.cpp.

2848  {
2849  MoFEM::Interface &m_field = cOre;
2850  const Problem *problem_ptr;
2852  CHKERR m_field.get_problem(name, &problem_ptr);
2853  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2854  fe_name, layout);
2856 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
#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 * get_problem(const std::string &problem_name) const =0
Get the problem object.

◆ 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.

Examples
build_problems.cpp.

Definition at line 2168 of file ProblemsManager.cpp.

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

◆ 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 elementsFunction 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
Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, and remove_entities_from_problem.cpp.

Definition at line 2418 of file ProblemsManager.cpp.

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

◆ 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.

Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, and remove_entities_from_problem.cpp.

Definition at line 2601 of file ProblemsManager.cpp.

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

◆ 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 2735 of file ProblemsManager.cpp.

2736  {
2737  MoFEM::Interface &m_field = cOre;
2738  auto problems_ptr = m_field.get_problems();
2740 
2742  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2743  "partition of problem not build");
2745  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2746  "partitions finite elements not build");
2747 
2748  // get problem pointer
2749  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2750  // find p_miit
2751  ProblemsByName &problems_set =
2752  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2753  ProblemsByName::iterator p_miit = problems_set.find(name);
2754 
2755  // get reference to number of ghost dofs
2756  // get number of local dofs
2757  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2758  &(p_miit->nbGhostDofsCol)};
2759  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2760  for (int ss = 0; ss != 2; ++ss) {
2761  (*nb_ghost_dofs[ss]) = 0;
2762  }
2763 
2764  // do work if more than one processor
2765  if (m_field.get_comm_size() > 1) {
2766  // determine if rows on columns are different from dofs on rows
2767  int loop_size = 2;
2768  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2769  loop_size = 1;
2770  }
2771 
2772  typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2773  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2774  p_miit->numeredColDofsPtr};
2775 
2776  // interate over dofs on rows and dofs on columns
2777  for (int ss = 0; ss != loop_size; ++ss) {
2778 
2779  // create dofs view by uid
2780  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2781 
2782  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2783  ghost_idx_view.reserve(std::distance(r.first, r.second));
2784  for (; r.first != r.second; ++r.first)
2785  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2786 
2787  auto cmp = [](auto a, auto b) {
2788  return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2789  };
2790  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2791 
2792  // intare over dofs which have negative local index
2793  for (auto gid_it : ghost_idx_view) {
2794  bool success = numered_dofs[ss]->modify(
2795  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2796  if (!success)
2797  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2798  "modification unsuccessful");
2799  ++(*nb_ghost_dofs[ss]);
2800  }
2801  }
2802  if (loop_size == 1) {
2803  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2804  }
2805  }
2806 
2807  if (verb > QUIET) {
2808  MOFEM_LOG("SYNC", Sev::inform)
2809  << " FEs ghost dofs on problem " << p_miit->getName()
2810  << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2811  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2812  << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2813 
2814  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2815  }
2816 
2819 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
virtual int get_comm_size() const =0
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
const double r
rate factor
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
int DofIdx
Index of DOF.
Definition: Types.hpp:29
virtual MPI_Comm & get_comm() const =0
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340

◆ 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 problemThis 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 }
#define ProblemManagerFunctionBegin
Deprecated interface functions.
FTensor::Index< 'i', 3 > i
virtual moab::Interface & get_moab()=0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
virtual int get_comm_size() const =0
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
const int dim
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
static const bool debug
FTensor::Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual MPI_Comm & get_comm() const =0
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:221
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340

◆ partitionProblem()

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

partition problem dofs (collective)

Parameters
nameproblem name
Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, and mesh_insert_interface_atom.cpp.

Definition at line 1956 of file ProblemsManager.cpp.

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

◆ partitionSimpleProblem()

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

partition problem dofs

Parameters
nameproblem name
Examples
forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, and quad_polynomial_approximation.cpp.

Definition at line 1810 of file ProblemsManager.cpp.

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

◆ 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 problemRemove 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
Examples
EshelbianPlasticity.cpp, and remove_entities_from_problem.cpp.

Definition at line 2858 of file ProblemsManager.cpp.

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