v0.9.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 1513 of file ProblemsManager.cpp.

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

460  {
461 
462  MoFEM::Interface &m_field = cOre;
464  if (!(cOre.getBuildMoFEM() & (1 << 0)))
465  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
466  if (!(cOre.getBuildMoFEM() & (1 << 1)))
467  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
468  if (!(cOre.getBuildMoFEM() & (1 << 2)))
469  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
470  const Problem *problem_ptr;
471  CHKERR m_field.get_problem(name, &problem_ptr);
472  CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
473  cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
474  // function knows what he is doing
476 }
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:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

480  {
481  MoFEM::Interface &m_field = cOre;
482  const EntFiniteElement_multiIndex *fe_ent_ptr;
483  const DofEntity_multiIndex *dofs_field_ptr;
485  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
486 
487  // Note: Only allowed changes on problem_ptr structure which not influence
488  // multi-index.
489 
490  if (problem_ptr->getBitRefLevel().none()) {
491  SETERRQ1(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
492  problem_ptr->getName().c_str());
493  }
494  CHKERR m_field.clear_problem(problem_ptr->getName());
495  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
496  CHKERR m_field.get_dofs(&dofs_field_ptr);
497 
498  // zero finite elements
499  problem_ptr->numeredFiniteElements->clear();
500 
501  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
502  {
503  EntFiniteElement_multiIndex::iterator miit = fe_ent_ptr->begin();
504  EntFiniteElement_multiIndex::iterator hi_miit = fe_ent_ptr->end();
505  // iterate all finite element entities in database
506  for (; miit != hi_miit; miit++) {
507  // if element is in problem
508  if (((*miit)->getId() & problem_ptr->getBitFEId()).any()) {
509  BitRefLevel prb_bit = problem_ptr->getBitRefLevel();
510  BitRefLevel prb_mask = problem_ptr->getMaskBitRefLevel();
511  BitRefLevel fe_bit = (*miit)->getBitRefLevel();
512  // if entity is not problem refinement level
513  if ((fe_bit & prb_mask) != fe_bit)
514  continue;
515  if ((fe_bit & prb_bit) != prb_bit)
516  continue;
517  // get dof uids for rows and columns
518  CHKERR(*miit)->getRowDofView(*dofs_field_ptr, dofs_rows);
519  if (!square_matrix) {
520  CHKERR(*miit)->getColDofView(*dofs_field_ptr, dofs_cols);
521  }
522  }
523  }
524  }
525 
526  // Add row dofs to problem
527  {
528  // zero rows
529  problem_ptr->nbDofsRow = 0;
530  problem_ptr->nbLocDofsRow = 0;
531  problem_ptr->nbGhostDofsRow = 0;
532 
533  // add dofs for rows
534  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
535  hi_miit;
536  hi_miit = dofs_rows.get<0>().end();
537 
538  int count_dofs = dofs_rows.get<1>().count(true);
539  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
540  boost::shared_ptr<std::vector<NumeredDofEntity>>(
541  new std::vector<NumeredDofEntity>());
542  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
543  dofs_array->reserve(count_dofs);
544  miit = dofs_rows.get<0>().begin();
545  for (; miit != hi_miit; miit++) {
546  if ((*miit)->getActive()) {
547  dofs_array->emplace_back(*miit);
548  dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
549  }
550  }
551  auto hint = problem_ptr->numeredDofsRows->end();
552  for (auto &v : *dofs_array) {
553  hint = problem_ptr->numeredDofsRows->emplace_hint(hint, dofs_array, &v);
554  }
555  }
556 
557  // Add col dofs to problem
558  if (!square_matrix) {
559  // zero cols
560  problem_ptr->nbDofsCol = 0;
561  problem_ptr->nbLocDofsCol = 0;
562  problem_ptr->nbGhostDofsCol = 0;
563 
564  // add dofs for cols
565  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
566  hi_miit;
567  hi_miit = dofs_cols.get<0>().end();
568 
569  int count_dofs = 0;
570  miit = dofs_cols.get<0>().begin();
571  for (; miit != hi_miit; miit++) {
572  if (!(*miit)->getActive()) {
573  continue;
574  }
575  count_dofs++;
576  }
577 
578  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
579  boost::shared_ptr<std::vector<NumeredDofEntity>>(
580  new std::vector<NumeredDofEntity>());
581  problem_ptr->getColDofsSequence()->push_back(dofs_array);
582  dofs_array->reserve(count_dofs);
583  miit = dofs_cols.get<0>().begin();
584  for (; miit != hi_miit; miit++) {
585  if (!(*miit)->getActive()) {
586  continue;
587  }
588  dofs_array->emplace_back(*miit);
589  dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
590  }
591  auto hint = problem_ptr->numeredDofsCols->end();
592  for (auto &v : *dofs_array) {
593  hint = problem_ptr->numeredDofsCols->emplace_hint(hint, dofs_array, &v);
594  }
595  } else {
596  problem_ptr->numeredDofsCols = problem_ptr->numeredDofsRows;
597  problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
598  problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
599  }
600 
601  // job done, some debugging and postprocessing
602  if (verb > QUIET) {
603  PetscSynchronizedPrintf(
604  m_field.get_comm(), "Problem %s Nb. rows %u Nb. cols %u\n",
605  problem_ptr->getName().c_str(), problem_ptr->numeredDofsRows->size(),
606  problem_ptr->numeredDofsCols->size());
607  }
608  if (verb > VERBOSE) {
609  EntFiniteElement_multiIndex::iterator miit = fe_ent_ptr->begin();
610  EntFiniteElement_multiIndex::iterator hi_miit = fe_ent_ptr->end();
611  std::ostringstream ss;
612  ss << "rank " << m_field.get_comm_rank() << " ";
613  ss << "FEs data for problem " << *problem_ptr << std::endl;
614  for (; miit != hi_miit; miit++) {
615  ss << "rank " << m_field.get_comm_rank() << " ";
616  ss << **miit << std::endl;
617  }
618  ss << "rank " << m_field.get_comm_rank() << " ";
619  ss << "FEs row dofs " << *problem_ptr << " Nb. row dof "
620  << problem_ptr->getNbDofsRow() << std::endl;
621  NumeredDofEntity_multiIndex::iterator miit_dd_row =
622  problem_ptr->numeredDofsRows->begin();
623  for (; miit_dd_row != problem_ptr->numeredDofsRows->end(); miit_dd_row++) {
624  ss << "rank " << m_field.get_comm_rank() << " ";
625  ss << **miit_dd_row << std::endl;
626  }
627  ss << "rank " << m_field.get_comm_rank() << " ";
628  ss << "FEs col dofs " << *problem_ptr << " Nb. col dof "
629  << problem_ptr->getNbDofsCol() << std::endl;
630  NumeredDofEntity_multiIndex::iterator miit_dd_col =
631  problem_ptr->numeredDofsCols->begin();
632  for (; miit_dd_col != problem_ptr->numeredDofsCols->end(); miit_dd_col++) {
633  ss << "rank " << m_field.get_comm_rank() << " ";
634  ss << **miit_dd_col << std::endl;
635  }
636  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
637  }
638 
639  if (verb > QUIET) {
640  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
641  }
642  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
643  // uses this function knows what
644  // he is doing
645 
646  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
647 
649 }
PetscLogEvent MOFEM_EVENT_ProblemsManager
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
virtual MoFEMErrorCode get_dofs(const DofEntity_multiIndex **dofs_ptr) const =0
Get dofs multi index.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, const UId &, &DofEntity::getGlobalUniqueId > >, ordered_non_unique< const_mem_fun< DofEntity, char, &DofEntity::getActive > > > > DofEntity_multiIndex_active_view
multi-index view on DofEntity activity
virtual int get_comm_rank() const =0
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< DofEntity, UId, &DofEntity::globalUId > >, ordered_non_unique< tag< Composite_Name_And_Ent_And_EntDofIdx_mi_tag >, composite_key< DofEntity, const_mem_fun< DofEntity::interface_type_Field, boost::string_ref, &DofEntity::getNameRef >, const_mem_fun< DofEntity, EntityHandle, &DofEntity::getEnt >, const_mem_fun< DofEntity, DofIdx, &DofEntity::getEntDofIdx > > >, ordered_non_unique< tag< Unique_Ent_mi_tag >, const_mem_fun< DofEntity, const UId &, &DofEntity::getEntGlobalUniqueId > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< DofEntity::interface_type_Field, boost::string_ref, &DofEntity::getNameRef > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< DofEntity, EntityHandle, &DofEntity::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< DofEntity, const_mem_fun< DofEntity::interface_type_Field, boost::string_ref, &DofEntity::getNameRef >, const_mem_fun< DofEntity, EntityHandle, &DofEntity::getEnt > > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< DofEntity, const_mem_fun< DofEntity::interface_type_Field, boost::string_ref, &DofEntity::getNameRef >, const_mem_fun< DofEntity::interface_type_RefEntity, EntityType, &DofEntity::getEntType > > > > > DofEntity_multiIndex
MultiIndex container keeps DofEntity.
virtual MoFEMErrorCode get_ents_finite_elements(const EntFiniteElement_multiIndex **fe_ent_ptr) const =0
Get entities finite elements multi-index.
#define CHKERR
Inline error check.
Definition: definitions.h:596
multi_index_container< boost::shared_ptr< EntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< EntFiniteElement, UId, &EntFiniteElement::globalUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef > >, ordered_non_unique< tag< BitFEId_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, BitFEId, &EntFiniteElement::getId >, LtBit< BitFEId > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_RefEntity, EntityType, &EntFiniteElement::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< EntFiniteElement, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > > > > > EntFiniteElement_multiIndex
MultiIndex container for EntFiniteElement.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem

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

652  {
653  MoFEM::Interface &m_field = cOre;
655 
656  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
657  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
658  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
659  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
660  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
661  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
662 
663  const Problem *problem_ptr;
664  CHKERR m_field.get_problem(name, &problem_ptr);
665  CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
666  square_matrix, verb);
667 
670 
672 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
#define CHKERR
Inline error check.
Definition: definitions.h:596
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,...
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

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

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

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

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

2884  {
2885  MoFEM::Interface &m_field = cOre;
2886  const Problem *problem_ptr;
2888  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2889  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2890  auto fit =
2891  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2892  .lower_bound(fe_name);
2893  auto hi_fe_it =
2894  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2895  .upper_bound(fe_name);
2896  std::vector<EntityHandle> fe_vec;
2897  fe_vec.reserve(std::distance(fit, hi_fe_it));
2898  for (; fit != hi_fe_it; fit++)
2899  fe_vec.push_back(fit->get()->getEnt());
2900  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2901  fe_vec.size());
2903 }
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:477
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

2908  {
2909  MoFEM::Interface &m_field = cOre;
2910  const Problem *problem_ptr;
2912  CHKERR m_field.get_problem(name, &problem_ptr);
2913  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2914  fe_name, layout);
2916 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

2141  {
2142 
2143  MoFEM::Interface &m_field = cOre;
2144  const Problem_multiIndex *problems_ptr;
2146 
2148  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
2149 
2150  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2151 
2152  // find p_miit
2153  CHKERR m_field.get_problems(&problems_ptr);
2154  ProblemByName &problems_by_name =
2155  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2156  ProblemByName::iterator p_miit = problems_by_name.find(name);
2157  if (p_miit == problems_by_name.end()) {
2158  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2159  "problem with name < %s > not defined (top tip check spelling)",
2160  name.c_str());
2161  }
2162  if (verb > QUIET)
2163  PetscPrintf(m_field.get_comm(),
2164  "Compose problem %s from rows of %s and columns of %s\n",
2165  p_miit->getName().c_str(), problem_for_rows.c_str(),
2166  problem_for_cols.c_str());
2167 
2168  // find p_miit_row
2169  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
2170  if (p_miit_row == problems_by_name.end()) {
2171  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2172  "problem with name < %s > not defined (top tip check spelling)",
2173  problem_for_rows.c_str());
2174  }
2175  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
2176  p_miit_row->numeredDofsRows;
2177 
2178  // find p_mit_col
2179  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
2180  if (p_miit_col == problems_by_name.end()) {
2181  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2182  "problem with name < %s > not defined (top tip check spelling)",
2183  problem_for_cols.c_str());
2184  }
2185  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
2186  p_miit_col->numeredDofsCols;
2187 
2188  bool copy[] = {copy_rows, copy_cols};
2189  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
2190  p_miit->numeredDofsRows, p_miit->numeredDofsCols};
2191 
2192  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
2193  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
2194  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
2195  dofs_col};
2196 
2197  for (int ss = 0; ss < 2; ss++) {
2198 
2199  // build indices
2200  *nb_local_dofs[ss] = 0;
2201  if (!copy[ss]) {
2202 
2203  // only copy indices which are belong to some elements if this problem
2204  std::vector<int> is_local, is_new;
2205 
2206  NumeredDofEntityByUId &dofs_by_uid =
2207  copied_dofs[ss]->get<Unique_mi_tag>();
2208  for (NumeredDofEntity_multiIndex::iterator dit =
2209  composed_dofs[ss]->begin();
2210  dit != composed_dofs[ss]->end(); dit++) {
2211  NumeredDofEntityByUId::iterator diit =
2212  dofs_by_uid.find((*dit)->getGlobalUniqueId());
2213  if (diit == dofs_by_uid.end()) {
2214  SETERRQ(
2216  "data inconsistency, could not find dof in composite problem");
2217  }
2218  int part_number = (*diit)->getPart(); // get part number
2219  int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
2220  bool success;
2221  success = composed_dofs[ss]->modify(
2222  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2223  petsc_global_dof));
2224  if (!success) {
2225  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2226  "modification unsuccessful");
2227  }
2228  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2229  success = composed_dofs[ss]->modify(
2230  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
2231  if (!success) {
2232  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2233  "modification unsuccessful");
2234  }
2235  is_local.push_back(petsc_global_dof);
2236  }
2237  }
2238 
2239  AO ao;
2240  CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
2241  NULL, &ao);
2242 
2243  // apply local to global mapping
2244  is_local.resize(0);
2245  for (NumeredDofEntity_multiIndex::iterator dit =
2246  composed_dofs[ss]->begin();
2247  dit != composed_dofs[ss]->end(); dit++) {
2248  is_local.push_back((*dit)->getPetscGlobalDofIdx());
2249  }
2250  CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2251  int idx2 = 0;
2252  for (NumeredDofEntity_multiIndex::iterator dit =
2253  composed_dofs[ss]->begin();
2254  dit != composed_dofs[ss]->end(); dit++) {
2255  int part_number = (*dit)->getPart(); // get part number
2256  int petsc_global_dof = is_local[idx2++];
2257  bool success;
2258  success = composed_dofs[ss]->modify(
2259  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2260  petsc_global_dof));
2261  if (!success) {
2262  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2263  "modification unsuccessful");
2264  }
2265  }
2266 
2267  CHKERR AODestroy(&ao);
2268 
2269  } else {
2270 
2271  for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2272  dit != copied_dofs[ss]->end(); dit++) {
2273  std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2274  p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2275  new NumeredDofEntity((*dit)->getDofEntityPtr())));
2276  if (p.second) {
2277  (*nb_dofs[ss])++;
2278  }
2279  int dof_idx = (*dit)->getDofIdx();
2280  int part_number = (*dit)->getPart(); // get part number
2281  int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2282  if (part_number == (unsigned int)m_field.get_comm_rank()) {
2283  const bool success = composed_dofs[ss]->modify(
2284  p.first, NumeredDofEntity_part_and_all_indices_change(
2285  part_number, dof_idx, petsc_global_dof,
2286  (*nb_local_dofs[ss])++));
2287  if (!success) {
2288  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2289  "modification unsuccessful");
2290  }
2291  } else {
2292  const bool success = composed_dofs[ss]->modify(
2293  p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
2294  part_number, dof_idx, petsc_global_dof));
2295  if (!success) {
2296  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2297  "modification unsuccessful");
2298  }
2299  }
2300  }
2301  }
2302  }
2303 
2304  CHKERR printPartitionedProblem(&*p_miit, verb);
2305  CHKERR debugPartitionedProblem(&*p_miit, verb);
2306 
2308 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
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 MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.

◆ partitionFiniteElements()

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

partition finite 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_check_approx_in_2d.cpp, hcurl_curl_operator.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 2417 of file ProblemsManager.cpp.

2420  {
2421  MoFEM::Interface &m_field = cOre;
2422  const Problem_multiIndex *problems_ptr;
2423  const EntFiniteElement_multiIndex *fe_ent_ptr;
2425 
2427  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2428 
2429  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2430  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2431 
2432  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2433  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2434  "adjacencies not build");
2435 
2437  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2438 
2440  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2441  "problem not partitioned");
2442 
2443  if (low_proc == -1)
2444  low_proc = m_field.get_comm_rank();
2445  if (hi_proc == -1)
2446  hi_proc = m_field.get_comm_rank();
2447 
2448  // Find pointer to problem of given name
2449  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2450  // Get pointer to problem data struture
2451  CHKERR m_field.get_problems(&problems_ptr);
2452  ProblemByName &problems =
2453  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2454  ProblemByName::iterator p_miit = problems.find(name);
2455  if (p_miit == problems.end()) {
2456  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2457  "problem < %s > not found (top tip: check spelling)",
2458  name.c_str());
2459  }
2460 
2461  // Get reference on finite elements multi-index on the problem
2462  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2463  *p_miit->numeredFiniteElements;
2464 
2465  // Clear all elements and data, build it again
2466  problem_finite_elements.clear();
2467 
2468  // Check if dofs and columns are the same, i.e. structurally symmetric problem
2469  bool do_cols_prob = true;
2470  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
2471  do_cols_prob = false;
2472  }
2473 
2474  // Loop over all elements in database and if right element is there add it
2475  // to problem finite element multi-index
2476  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
2477  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); efit++) {
2478 
2479  // if element is not part of problem
2480  if (((*efit)->getId() & p_miit->getBitFEId()).none())
2481  continue;
2482 
2483  BitRefLevel prb_bit = p_miit->getBitRefLevel();
2484  BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
2485  BitRefLevel fe_bit = (*efit)->getBitRefLevel();
2486  // if entity is not problem refinement level
2487  if ((fe_bit & prb_mask) != fe_bit)
2488  continue;
2489  if ((fe_bit & prb_bit) != prb_bit)
2490  continue;
2491 
2492  // create element
2493  boost::shared_ptr<NumeredEntFiniteElement> numered_fe(
2494  new NumeredEntFiniteElement(*efit));
2495 
2496  // check if rows and columns are the same on this element
2497  bool do_cols_fe = true;
2498  if ((numered_fe->sPtr->row_field_ents_view ==
2499  numered_fe->sPtr->col_field_ents_view) &&
2500  !do_cols_prob) {
2501  do_cols_fe = false;
2502  numered_fe->cols_dofs = numered_fe->rows_dofs;
2503  } else {
2504  // different dofs on rows and columns
2505  numered_fe->cols_dofs = boost::shared_ptr<FENumeredDofEntity_multiIndex>(
2507  }
2508  // get pointer to dofs multi-index on rows and columns
2509  auto rows_dofs = numered_fe->rows_dofs;
2510  auto cols_dofs = numered_fe->cols_dofs;
2511  // clear multi-indices
2512  rows_dofs->clear();
2513  if (do_cols_fe) {
2514  cols_dofs->clear();
2515  }
2518 
2519  // set partition to the element
2520  {
2521  if (part_from_moab) {
2522  // if partition is taken from moab partition
2523  int proc = (*efit)->getPartProc();
2524  if (proc == -1 && (*efit)->getEntType() == MBVERTEX) {
2525  proc = (*efit)->getOwnerProc();
2526  }
2527  NumeredEntFiniteElement_change_part(proc).operator()(numered_fe);
2528  } else {
2529  // count partition of the dofs in row, the larges dofs with given
2530  // partition is used to set partition of the element
2531  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows), rows_view,
2532  moab::Interface::UNION);
2533  std::vector<int> parts(m_field.get_comm_size(), 0);
2534  NumeredDofEntity_multiIndex_uid_view_ordered::iterator viit_rows;
2535  viit_rows = rows_view.begin();
2536  for (; viit_rows != rows_view.end(); viit_rows++) {
2537  parts[(*viit_rows)->pArt]++;
2538  }
2539  std::vector<int>::iterator pos =
2540  max_element(parts.begin(), parts.end());
2541  unsigned int max_part = std::distance(parts.begin(), pos);
2542  NumeredEntFiniteElement_change_part(max_part).operator()(numered_fe);
2543  }
2544  }
2545 
2546  // set dofs on rows and columns (if are different)
2547  if ((numered_fe->getPart() >= (unsigned int)low_proc) &&
2548  (numered_fe->getPart() <= (unsigned int)hi_proc)) {
2549 
2550  NumeredDofEntity_multiIndex_uid_view_ordered *dofs_view[] = {&rows_view,
2551  &cols_view};
2552  FENumeredDofEntity_multiIndex *fe_dofs[] = {rows_dofs.get(),
2553  cols_dofs.get()};
2554 
2555  for (int ss = 0; ss != (do_cols_fe ? 2 : 1); ss++) {
2556 
2557  if (ss == 0) {
2558  if (part_from_moab) {
2559  // get row_view
2560  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows),
2561  *dofs_view[ss],
2562  moab::Interface::UNION);
2563  }
2564  } else {
2565  // get cols_views
2566  CHKERR(*efit)->getColDofView(*(p_miit->numeredDofsCols),
2567  *dofs_view[ss], moab::Interface::UNION);
2568  }
2569 
2570  // Following reserve memory in sequences, only two allocations are here,
2571  // once for array of objects, next for array of shared pointers
2572 
2573  // reserve memory for field dofs
2574  boost::shared_ptr<std::vector<FENumeredDofEntity>> dofs_array(
2575  new std::vector<FENumeredDofEntity>());
2576 
2577  if (!ss) {
2578  numered_fe->getRowDofsSequence() = dofs_array;
2579  if (!do_cols_fe)
2580  numered_fe->getColDofsSequence() = dofs_array;
2581  } else
2582  numered_fe->getColDofsSequence() = dofs_array;
2583 
2584  auto vit = dofs_view[ss]->begin();
2585  auto hi_vit = dofs_view[ss]->end();
2586 
2587  dofs_array->reserve(std::distance(vit, hi_vit));
2588 
2589  // create elements objects
2590  for (; vit != hi_vit; vit++) {
2591  boost::shared_ptr<SideNumber> side_number_ptr;
2592  side_number_ptr = (*efit)->getSideNumberPtr((*vit)->getEnt());
2593  dofs_array->emplace_back(side_number_ptr, *vit);
2594  }
2595 
2596  // finally add DoFS to multi-indices
2597  auto hint = fe_dofs[ss]->end();
2598  for (auto &v : *dofs_array)
2599  hint = fe_dofs[ss]->emplace_hint(hint, dofs_array, &v);
2600  }
2601  }
2602  if (!numered_fe->sPtr->row_field_ents_view->empty() &&
2603  !numered_fe->sPtr->col_field_ents_view->empty()) {
2604 
2605  // Add element to the problem
2606  auto p = problem_finite_elements.insert(numered_fe);
2607  if (!p.second)
2608  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2609 
2610  if (verb >= VERY_VERBOSE) {
2611  std::ostringstream ss;
2612  ss << *p_miit << std::endl;
2613  ss << *p.first << std::endl;
2615  FENumeredDofEntityByUId::iterator miit =
2616  (*p.first)->rows_dofs->get<Unique_mi_tag>().begin();
2617  for (; miit != (*p.first)->rows_dofs->get<Unique_mi_tag>().end();
2618  miit++)
2619  ss << "rows: " << *(*miit) << std::endl;
2620  miit = (*p.first)->cols_dofs->get<Unique_mi_tag>().begin();
2621  for (; miit != (*p.first)->cols_dofs->get<Unique_mi_tag>().end();
2622  miit++)
2623  ss << "cols: " << *(*miit) << std::endl;
2624  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2625  }
2626  }
2627  }
2628 
2629  if (verb >= VERBOSE) {
2630  typedef NumeredEntFiniteElement_multiIndex::index<Part_mi_tag>::type
2631  NumeredEntFiniteElementPart;
2632  NumeredEntFiniteElementPart::iterator miit, hi_miit;
2633  miit = problem_finite_elements.get<Part_mi_tag>().lower_bound(
2634  m_field.get_comm_rank());
2635  hi_miit = problem_finite_elements.get<Part_mi_tag>().upper_bound(
2636  m_field.get_comm_rank());
2637  int count = std::distance(miit, hi_miit);
2638  std::ostringstream ss;
2639  ss << *p_miit;
2640  ss << " Nb. elems " << count << " on proc " << m_field.get_comm_rank()
2641  << std::endl;
2642  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2643  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2644  }
2645 
2648 }
virtual int get_comm_size() const =0
multi_index_container< boost::shared_ptr< FENumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, const UId &, &FENumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FENumeredDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FENumeredDofEntity::sideNumberPtr > > > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > > > > > FENumeredDofEntity_multiIndex
MultiIndex container keeps FENumeredDofEntity.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
FENumeredDofEntity_multiIndex::index< Unique_mi_tag >::type FENumeredDofEntityByUId
Dof entity multi-index by UId.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getGlobalUniqueId > >, 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_EntFiniteElement, 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_EntFiniteElement, 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 MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
virtual MoFEMErrorCode get_ents_finite_elements(const EntFiniteElement_multiIndex **fe_ent_ptr) const =0
Get entities finite elements multi-index.
#define CHKERR
Inline error check.
Definition: definitions.h:596
multi_index_container< boost::shared_ptr< EntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< EntFiniteElement, UId, &EntFiniteElement::globalUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef > >, ordered_non_unique< tag< BitFEId_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, BitFEId, &EntFiniteElement::getId >, LtBit< BitFEId > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_RefEntity, EntityType, &EntFiniteElement::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< EntFiniteElement, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > > > > > EntFiniteElement_multiIndex
MultiIndex container for EntFiniteElement.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0

◆ 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_check_approx_in_2d.cpp, hcurl_curl_operator.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 2650 of file ProblemsManager.cpp.

2651  {
2652  MoFEM::Interface &m_field = cOre;
2653  const Problem_multiIndex *problems_ptr;
2655 
2657  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2658  "partition of problem not build");
2660  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2661  "partitions finite elements not build");
2662 
2663  // find p_miit
2664  CHKERR m_field.get_problems(&problems_ptr);
2665 
2666  // get problem pointer
2667  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2668  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2669  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2670  name.c_str());
2671 
2672  // get reference to number of ghost dofs
2673  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2674  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2675  nb_row_ghost_dofs = 0;
2676  nb_col_ghost_dofs = 0;
2677 
2678  // do work if more than one processor
2679  if (m_field.get_comm_size() > 1) {
2680 
2682  ghost_idx_col_view;
2683 
2684  // get elements on this partition
2685  auto fe_range =
2686  p_miit->numeredFiniteElements->get<Part_mi_tag>().equal_range(
2687  m_field.get_comm_rank());
2688 
2689  // get dofs on elements which are not part of this partition
2690 
2691  // rows
2692  auto hint_r = ghost_idx_row_view.begin();
2693  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2694  for (auto &dof_ptr : *(*fe_ptr)->rows_dofs) {
2695  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2696  hint_r = ghost_idx_row_view.emplace_hint(
2697  hint_r, dof_ptr->getNumeredDofEntityPtr());
2698  }
2699  }
2700  }
2701 
2702  // columns
2703  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2704  auto hint_c = ghost_idx_col_view.begin();
2705  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2706  for (auto &dof_ptr : *(*fe_range.first)->cols_dofs) {
2707  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2708  hint_c = ghost_idx_col_view.emplace_hint(
2709  hint_c, dof_ptr->getNumeredDofEntityPtr());
2710  }
2711  }
2712  }
2713  }
2714 
2715  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2716  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2717 
2718  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2719  &ghost_idx_row_view, &ghost_idx_col_view};
2720  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2721  &p_miit->numeredDofsRows->get<Unique_mi_tag>(),
2722  &p_miit->numeredDofsCols->get<Unique_mi_tag>()};
2723 
2724  int loop_size = 2;
2725  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2726  loop_size = 1;
2727  }
2728 
2729  // set local ghost dofs indices
2730  for (int ss = 0; ss != loop_size; ++ss) {
2731  for (auto &gid : *ghost_idx_view[ss]) {
2732  NumeredDofEntityByUId::iterator dof =
2733  dof_by_uid_no_const[ss]->find(gid->getGlobalUniqueId());
2734  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2735  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2736  "inconsistent data, ghost dof already set");
2737  bool success = dof_by_uid_no_const[ss]->modify(
2738  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2739  if (PetscUnlikely(!success))
2740  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2741  "modification unsuccessful");
2742  (*nb_ghost_dofs[ss])++;
2743  }
2744  }
2745  if (loop_size == 1) {
2746  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2747  }
2748  }
2749 
2750  if (verb > QUIET) {
2751  std::ostringstream ss;
2752  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2753  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2754  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2755  << p_miit->getNbLocalDofsCol() << std::endl;
2756  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2757  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2758  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2759  << p_miit->getNbLocalDofsRow() << std::endl;
2760  if (verb > VERBOSE) {
2761  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2762  p_miit->numeredDofsCols->begin();
2763  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2764  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2765  continue;
2766  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2767  continue;
2768  ss << *(*miit_dd_col) << std::endl;
2769  }
2770  }
2771  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2772  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2773  }
2774 
2777 }
virtual int get_comm_size() const =0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
#define CHKERR
Inline error check.
Definition: definitions.h:596
int DofIdx
Index of DOF.
Definition: Types.hpp:29
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.

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

2781  {
2782  MoFEM::Interface &m_field = cOre;
2783  const Problem_multiIndex *problems_ptr;
2785 
2787  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2788  "partition of problem not build");
2790  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2791  "partitions finite elements not build");
2792 
2793  // get problem pointer
2794  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2795  // find p_miit
2796  CHKERR m_field.get_problems(&problems_ptr);
2797  ProblemsByName &problems_set =
2798  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2799  ProblemsByName::iterator p_miit = problems_set.find(name);
2800 
2801  // get reference to number of ghost dofs
2802  // get number of local dofs
2803  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2804  &(p_miit->nbGhostDofsCol)};
2805  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2806  for (int ss = 0; ss != 2; ++ss) {
2807  (*nb_ghost_dofs[ss]) = 0;
2808  }
2809 
2810  // do work if more than one processor
2811  if (m_field.get_comm_size() > 1) {
2812  // determine if rows on columns are different from dofs on rows
2813  int loop_size = 2;
2814  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2815  loop_size = 1;
2816  }
2817 
2818  typedef decltype(p_miit->numeredDofsRows) NumbDofTypeSharedPtr;
2819  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredDofsRows,
2820  p_miit->numeredDofsCols};
2821 
2822  // interate over dofs on rows and dofs on columns
2823  for (int ss = 0; ss != loop_size; ++ss) {
2824 
2825  // create dofs view by uid
2826  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2827 
2828  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2829  ghost_idx_view.reserve(std::distance(r.first, r.second));
2830  for (; r.first != r.second; ++r.first)
2831  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2832 
2833  auto cmp = [](auto a, auto b) {
2834  return (*a)->getGlobalUniqueId() < (*b)->getGlobalUniqueId();
2835  };
2836  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2837 
2838  // intare over dofs which have negative local index
2839  for (auto gid_it : ghost_idx_view) {
2840  bool success = numered_dofs[ss]->modify(
2841  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2842  if (!success)
2843  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2844  "modification unsuccessful");
2845  ++(*nb_ghost_dofs[ss]);
2846  }
2847  }
2848  if (loop_size == 1) {
2849  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2850  }
2851  }
2852 
2853  if (verb > QUIET) {
2854  std::ostringstream ss;
2855  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2856  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2857  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2858  << p_miit->getNbLocalDofsRow() << std::endl;
2859  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2860  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2861  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2862  << p_miit->getNbLocalDofsCol() << std::endl;
2863  if (verb > VERBOSE) {
2864  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2865  p_miit->numeredDofsCols->begin();
2866  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2867  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2868  continue;
2869  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2870  continue;
2871  ss << *(*miit_dd_col) << std::endl;
2872  }
2873  }
2874  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2875  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2876  }
2877 
2880 }
virtual int get_comm_size() const =0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
#define CHKERR
Inline error check.
Definition: definitions.h:596
int DofIdx
Index of DOF.
Definition: Types.hpp:29
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0

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

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

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

1954  {
1955  MoFEM::Interface &m_field = cOre;
1956  const Problem_multiIndex *problems_ptr;
1958 
1959  if (!(cOre.getBuildMoFEM() & (1 << 0)))
1960  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1961  if (!(cOre.getBuildMoFEM() & (1 << 1)))
1962  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1963  if (!(cOre.getBuildMoFEM() & (1 << 2)))
1964  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1965  "entFEAdjacencies not build");
1966  if (!(cOre.getBuildMoFEM() & (1 << 3)))
1967  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1968 
1969  if (verb > QUIET)
1970  PetscPrintf(m_field.get_comm(), "Partition problem %s\n", name.c_str());
1971 
1972  typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
1973  NumeredDofEntitysByIdx;
1974  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
1975 
1976  // Find problem pointer by name
1977  CHKERR m_field.get_problems(&problems_ptr);
1978  auto &problems_set =
1979  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1980  auto p_miit = problems_set.find(name);
1981  if (p_miit == problems_set.end()) {
1982  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1983  "problem with name %s not defined (top tip check spelling)",
1984  name.c_str());
1985  }
1986  int nb_dofs_row = p_miit->getNbDofsRow();
1987 
1988  Mat Adj;
1989  CHKERR m_field.getInterface<MatrixManager>()
1990  ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1991 
1992  int m, n;
1993  CHKERR MatGetSize(Adj, &m, &n);
1994  if (verb > VERY_VERBOSE)
1995  MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1996 
1997  // partitioning
1998  MatPartitioning part;
1999  IS is;
2000  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
2001  CHKERR MatPartitioningSetAdjacency(part, Adj);
2002  CHKERR MatPartitioningSetFromOptions(part);
2003  CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
2004  CHKERR MatPartitioningApply(part, &is);
2005  if (verb > VERY_VERBOSE)
2006  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
2007 
2008  // gather
2009  IS is_gather, is_num, is_gather_num;
2010  CHKERR ISAllGather(is, &is_gather);
2011  CHKERR ISPartitioningToNumbering(is, &is_num);
2012  CHKERR ISAllGather(is_num, &is_gather_num);
2013  const int *part_number, *petsc_idx;
2014  CHKERR ISGetIndices(is_gather, &part_number);
2015  CHKERR ISGetIndices(is_gather_num, &petsc_idx);
2016  int size_is_num, size_is_gather;
2017  CHKERR ISGetSize(is_gather, &size_is_gather);
2018  if (size_is_gather != (int)nb_dofs_row)
2019  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
2020  size_is_gather, nb_dofs_row);
2021 
2022  CHKERR ISGetSize(is_num, &size_is_num);
2023  if (size_is_num != (int)nb_dofs_row)
2024  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
2025  size_is_num, nb_dofs_row);
2026 
2027  bool square_matrix = false;
2028  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols)
2029  square_matrix = true;
2030 
2031  if (!square_matrix) {
2032  // FIXME: This is for back compatibility, if deprecate interface function
2033  // build interfaces is removed, this part of the code will be obsolete
2034  auto mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().begin();
2035  auto hi_mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().end();
2036  auto mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().begin();
2037  auto hi_mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().end();
2038  for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
2039  if (mit_col == hi_mit_col) {
2040  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2041  "check finite element definition, nb. of rows is not equal to "
2042  "number for columns");
2043  }
2044  if (mit_row->get()->getGlobalUniqueId() !=
2045  mit_col->get()->getGlobalUniqueId()) {
2046  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2047  "check finite element definition, nb. of rows is not equal to "
2048  "number for columns");
2049  }
2050  }
2051  }
2052 
2053  if (verb > VERBOSE)
2054  PetscPrintf(m_field.get_comm(), "\tloop problem dofs");
2055 
2056  // Set petsc global indices
2057  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
2058  p_miit->numeredDofsRows->get<Idx_mi_tag>());
2059  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2060  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2061  nb_row_local_dofs = 0;
2062  nb_row_ghost_dofs = 0;
2063 
2064  for (auto miit_dofs_row = dofs_row_by_idx_no_const.begin();
2065  miit_dofs_row != dofs_row_by_idx_no_const.end(); miit_dofs_row++) {
2066  const int part = part_number[(*miit_dofs_row)->dofIdx];
2067  if (part == (unsigned int)m_field.get_comm_rank()) {
2068  const bool success = dofs_row_by_idx_no_const.modify(
2069  miit_dofs_row,
2070  NumeredDofEntity_part_and_indices_change(
2071  part, petsc_idx[(*miit_dofs_row)->dofIdx], nb_row_local_dofs++));
2072  if (!success) {
2073  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2074  "modification unsuccessful");
2075  }
2076  } else {
2077  const bool success = dofs_row_by_idx_no_const.modify(
2078  miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
2079  part, petsc_idx[(*miit_dofs_row)->dofIdx]));
2080  if (!success) {
2081  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2082  "modification unsuccessful");
2083  }
2084  }
2085  }
2086 
2087  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2088  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2089  if (square_matrix) {
2090  nb_col_local_dofs = nb_row_local_dofs;
2091  nb_col_ghost_dofs = nb_row_ghost_dofs;
2092  } else {
2093  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2094  const_cast<NumeredDofEntitysByIdx &>(
2095  p_miit->numeredDofsCols->get<Idx_mi_tag>());
2096  nb_col_local_dofs = 0;
2097  nb_col_ghost_dofs = 0;
2098  for (auto miit_dofs_col = dofs_col_by_idx_no_const.begin();
2099  miit_dofs_col != dofs_col_by_idx_no_const.end(); miit_dofs_col++) {
2100  const int part = part_number[(*miit_dofs_col)->dofIdx];
2101  if (part == (unsigned int)m_field.get_comm_rank()) {
2102  const bool success = dofs_col_by_idx_no_const.modify(
2103  miit_dofs_col, NumeredDofEntity_part_and_indices_change(
2104  part, petsc_idx[(*miit_dofs_col)->dofIdx],
2105  nb_col_local_dofs++));
2106  if (!success) {
2107  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2108  "modification unsuccessful");
2109  }
2110  } else {
2111  const bool success = dofs_col_by_idx_no_const.modify(
2112  miit_dofs_col, NumeredDofEntity_part_and_glob_idx_change(
2113  part, petsc_idx[(*miit_dofs_col)->dofIdx]));
2114  if (!success) {
2115  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2116  "modification unsuccessful");
2117  }
2118  }
2119  }
2120  }
2121 
2122  if (verb > VERBOSE)
2123  PetscPrintf(m_field.get_comm(), " <- done\n");
2124 
2125  CHKERR ISRestoreIndices(is_gather, &part_number);
2126  CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
2127  CHKERR ISDestroy(&is_num);
2128  CHKERR ISDestroy(&is_gather_num);
2129  CHKERR ISDestroy(&is_gather);
2130  CHKERR ISDestroy(&is);
2131  CHKERR MatPartitioningDestroy(&part);
2132  CHKERR MatDestroy(&Adj);
2133  CHKERR printPartitionedProblem(&*p_miit, verb);
2134 
2135  cOre.getBuildMoFEM() |= 1 << 4;
2137 }
virtual int get_comm_size() const =0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
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
virtual MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0

◆ 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_check_approx_in_2d.cpp, hcurl_curl_operator.cpp, hcurl_divergence_operator_2d.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, and quad_polynomial_approximation.cpp.

Definition at line 1805 of file ProblemsManager.cpp.

1806  {
1807 
1808  MoFEM::Interface &m_field = cOre;
1809  const Problem_multiIndex *problems_ptr;
1812  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1813  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1814  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1815  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1816  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1818  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1819  if (verb > 0) {
1820  PetscPrintf(m_field.get_comm(), "Simple partition problem %s\n",
1821  name.c_str());
1822  }
1823  CHKERR m_field.get_problems(&problems_ptr);
1824  // find p_miit
1825  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1826  ProblemByName &problems_set =
1827  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1828  ProblemByName::iterator p_miit = problems_set.find(name);
1829  if (p_miit == problems_set.end()) {
1830  SETERRQ1(PETSC_COMM_SELF, 1,
1831  "problem < %s > is not found (top tip: check spelling)",
1832  name.c_str());
1833  }
1834  typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1835  Idx_mi_tag>::type NumeredDofEntitysByIdx;
1836  NumeredDofEntitysByIdx &dofs_row_by_idx =
1837  p_miit->numeredDofsRows->get<Idx_mi_tag>();
1838  NumeredDofEntitysByIdx &dofs_col_by_idx =
1839  p_miit->numeredDofsCols->get<Idx_mi_tag>();
1840  boost::multi_index::index<NumeredDofEntity_multiIndex,
1841  Idx_mi_tag>::type::iterator miit_row,
1842  hi_miit_row;
1843  boost::multi_index::index<NumeredDofEntity_multiIndex,
1844  Idx_mi_tag>::type::iterator miit_col,
1845  hi_miit_col;
1846  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1847  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1848  nb_row_local_dofs = 0;
1849  nb_row_ghost_dofs = 0;
1850  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1851  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1852  nb_col_local_dofs = 0;
1853  nb_col_ghost_dofs = 0;
1854 
1855  bool square_matrix = false;
1856  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
1857  square_matrix = true;
1858  }
1859 
1860  // get row range of local indices
1861  PetscLayout layout_row;
1862  const int *ranges_row;
1863 
1864  DofIdx nb_dofs_row = dofs_row_by_idx.size();
1865  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1866  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1867  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1868  CHKERR PetscLayoutSetUp(layout_row);
1869  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1870  // get col range of local indices
1871  PetscLayout layout_col;
1872  const int *ranges_col;
1873  if (!square_matrix) {
1874  DofIdx nb_dofs_col = dofs_col_by_idx.size();
1875  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1876  CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1877  CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1878  CHKERR PetscLayoutSetUp(layout_col);
1879  CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1880  }
1881  for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1882  part++) {
1883  miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1884  hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1885  if (std::distance(miit_row, hi_miit_row) !=
1886  ranges_row[part + 1] - ranges_row[part]) {
1887  SETERRQ4(
1888  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1889  "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1890  "rstart (%d != %d - %d = %d) ",
1891  std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1892  ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1893  }
1894  // loop rows
1895  for (; miit_row != hi_miit_row; miit_row++) {
1896  bool success = dofs_row_by_idx.modify(
1897  miit_row,
1898  NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1899  if (!success)
1900  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1901  "modification unsuccessful");
1902  if (part == (unsigned int)m_field.get_comm_rank()) {
1903  success = dofs_row_by_idx.modify(
1904  miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1905  if (!success)
1906  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1907  "modification unsuccessful");
1908  }
1909  }
1910  if (!square_matrix) {
1911  miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1912  hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1913  if (std::distance(miit_col, hi_miit_col) !=
1914  ranges_col[part + 1] - ranges_col[part]) {
1915  SETERRQ4(
1916  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1917  "data inconsistency, std::distance(miit_col,hi_miit_col) != rend - "
1918  "rstart (%d != %d - %d = %d) ",
1919  std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1920  ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1921  }
1922  // loop cols
1923  for (; miit_col != hi_miit_col; miit_col++) {
1924  bool success = dofs_col_by_idx.modify(
1925  miit_col, NumeredDofEntity_part_and_glob_idx_change(
1926  part, (*miit_col)->dofIdx));
1927  if (!success)
1928  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1929  "modification unsuccessful");
1930  if (part == (unsigned int)m_field.get_comm_rank()) {
1931  success = dofs_col_by_idx.modify(
1932  miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1933  if (!success)
1934  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1935  "modification unsuccessful");
1936  }
1937  }
1938  }
1939  }
1940  CHKERR PetscLayoutDestroy(&layout_row);
1941  if (!square_matrix) {
1942  CHKERR PetscLayoutDestroy(&layout_col);
1943  }
1944  if (square_matrix) {
1945  nb_col_local_dofs = nb_row_local_dofs;
1946  nb_col_ghost_dofs = nb_row_ghost_dofs;
1947  }
1948  CHKERR printPartitionedProblem(&*p_miit, verb);
1951 }
virtual int get_comm_size() const =0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:179
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode get_problems(const Problem_multiIndex **problems_ptr) const =0
Get pointer to problems multi-index.
multi_index_container< Problem, indexed_by< ordered_unique< tag< Meshset_mi_tag >, member< Problem, EntityHandle, &Problem::meshset > >, hashed_unique< tag< BitProblemId_mi_tag >, const_mem_fun< Problem, BitProblemId, &Problem::getId >, HashBit< BitProblemId >, EqBit< BitProblemId > >, hashed_unique< tag< Problem_mi_tag >, const_mem_fun< Problem, std::string, &Problem::getName > > > > Problem_multiIndex
MultiIndex for entities for Problem.
#define CHKERR
Inline error check.
Definition: definitions.h:596
int DofIdx
Index of DOF.
Definition: Types.hpp:29
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Unique_Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getEntGlobalUniqueId > >, 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< FieldName_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef > >, 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 > >, ordered_non_unique< tag< Composite_Name_And_Ent_And_EntDofIdx_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, DofIdx, &NumeredDofEntity::getEntDofIdx > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > >, ordered_non_unique< tag< Composite_Name_Ent_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.

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

2921  {
2922 
2923  MoFEM::Interface &m_field = cOre;
2925 
2926  const Problem *prb_ptr;
2927  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2928 
2929  decltype(prb_ptr->numeredDofsRows) numered_dofs[2] = {
2930  prb_ptr->numeredDofsRows, nullptr};
2931  if (prb_ptr->numeredDofsRows != prb_ptr->numeredDofsCols)
2932  numered_dofs[1] = prb_ptr->numeredDofsCols;
2933 
2934  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2935  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2936  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2937 
2938  for (int s = 0; s != 2; ++s)
2939  if (numered_dofs[s]) {
2940 
2941  typedef multi_index_container<
2942 
2943  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2944 
2945  >
2946  NumeredDofEntity_it_view_multiIndex;
2947 
2948  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2949 
2950  // Set -1 to global and local dofs indices
2951  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2952  ++pit) {
2953  auto lo =
2954  numered_dofs[s]
2955  ->get<Composite_Name_And_Ent_And_EntDofIdx_mi_tag>()
2956  .lower_bound(boost::make_tuple(field_name, pit->first, 0));
2957  auto hi = numered_dofs[s]
2958  ->get<Composite_Name_And_Ent_And_EntDofIdx_mi_tag>()
2959  .lower_bound(boost::make_tuple(field_name, pit->second,
2961  for (; lo != hi; ++lo)
2962  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2963  (*lo)->getDofCoeffIdx() <= hi_coeff)
2964  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2965  }
2966 
2967  if (verb >= VERY_NOISY) {
2968  for (auto &dof : dofs_it_view) {
2969  std::ostringstream ss;
2970  ss << **dof;
2971  PetscSynchronizedPrintf(m_field.get_comm(), "%s\n", ss.str().c_str());
2972  }
2973  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2974  }
2975 
2976  // set negative index
2977  auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2978  for (auto dit : dofs_it_view) {
2979  bool success = numered_dofs[s]->modify(dit, mod);
2980  if (!success)
2981  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2982  "can not set negative indices");
2983  }
2984 
2985  // create weak view
2986  std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2987  dosf_weak_view.reserve(dofs_it_view.size());
2988  for (auto dit : dofs_it_view)
2989  dosf_weak_view.push_back(*dit);
2990 
2991  if (verb >= NOISY)
2992  PetscSynchronizedPrintf(
2993  m_field.get_comm(),
2994  "Number of DOFs in multi-index %d and to delete %d\n",
2995  numered_dofs[s]->size(), dofs_it_view.size());
2996 
2997  // erase dofs from problem
2998  for (auto weak_dit : dosf_weak_view)
2999  if (auto dit = weak_dit.lock()) {
3000  numered_dofs[s]->erase(dit->getGlobalUniqueId());
3001  }
3002 
3003  if (verb >= NOISY)
3004  PetscSynchronizedPrintf(
3005  m_field.get_comm(),
3006  "Number of DOFs in multi-index after delete %d\n",
3007  numered_dofs[s]->size());
3008 
3009  // get current number of ghost dofs
3010  int nb_local_dofs = 0;
3011  int nb_ghost_dofs = 0;
3012  for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
3013  dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
3014  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
3015  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
3016  ++nb_local_dofs;
3017  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
3018  ++nb_ghost_dofs;
3019  else
3020  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Imposible case");
3021  }
3022 
3023  // get indices
3024  auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
3025  const int nb_dofs = numered_dofs[s]->size();
3026  indices.clear();
3027  indices.reserve(nb_dofs);
3028  for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
3029  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3030  bool add = true;
3031  if (only_local) {
3032  if ((*dit)->getPetscLocalDofIdx() < 0 ||
3033  (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
3034  add = false;
3035  }
3036  }
3037  if (add)
3038  indices.push_back(decltype(tag)::get_index(dit));
3039  }
3040  };
3041 
3042  auto get_indices_by_uid = [&](auto tag, auto &indices) {
3043  const int nb_dofs = numered_dofs[s]->size();
3044  indices.clear();
3045  indices.reserve(nb_dofs);
3046  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3047  ++dit)
3048  indices.push_back(decltype(tag)::get_index(dit));
3049  };
3050 
3051  auto concatenate_dofs = [&](auto tag, auto &indices,
3052  const auto local_only) {
3054  get_indices_by_tag(tag, indices, local_only);
3055 
3056  AO ao;
3057  if (local_only)
3058  CHKERR AOCreateMapping(m_field.get_comm(), indices.size(),
3059  &*indices.begin(), PETSC_NULL, &ao);
3060  else
3061  CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
3062  &*indices.begin(), PETSC_NULL, &ao);
3063 
3064  get_indices_by_uid(tag, indices);
3065  CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3066  CHKERR AODestroy(&ao);
3068  };
3069 
3070  // set indices index
3071  auto set_concatinated_indices = [&]() {
3072  std::vector<int> global_indices;
3073  std::vector<int> local_indices;
3075  CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3076  CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3077  auto gi = global_indices.begin();
3078  auto li = local_indices.begin();
3079  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3080  ++dit) {
3081  auto mod = NumeredDofEntity_part_and_all_indices_change(
3082  (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3083  bool success = numered_dofs[s]->modify(dit, mod);
3084  if (!success)
3085  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3086  "can not set negative indices");
3087  ++gi;
3088  ++li;
3089  }
3091  };
3092  CHKERR set_concatinated_indices();
3093 
3094  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3095  m_field.get_comm());
3096  *(local_nbdof_ptr[s]) = nb_local_dofs;
3097  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3098 
3099  if (debug)
3100  for (auto dof : (*numered_dofs[s])) {
3101  if (dof->getPetscGlobalDofIdx() < 0) {
3102  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3103  "Negative global idx");
3104  }
3105  if (dof->getPetscLocalDofIdx() < 0) {
3106  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3107  "Negative local idx");
3108  }
3109  }
3110 
3111  } else {
3112 
3113  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3114  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3115  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3116  }
3117 
3118  if (verb > QUIET) {
3119  PetscSynchronizedPrintf(
3120  m_field.get_comm(),
3121  "removed ents on rank %d from problem %s dofs [ %d / %d local, %d / %d "
3122  "ghost, %d / %d global]\n",
3123  m_field.get_comm_rank(), prb_ptr->getName().c_str(),
3124  prb_ptr->getNbLocalDofsRow(), prb_ptr->getNbLocalDofsCol(),
3125  prb_ptr->getNbGhostDofsRow(), prb_ptr->getNbGhostDofsCol(),
3126  prb_ptr->getNbDofsRow(), prb_ptr->getNbDofsCol());
3127  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
3128  }
3129 
3131 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:298
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0