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

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

459  {
460 
461  MoFEM::Interface &m_field = cOre;
463  if (!(cOre.getBuildMoFEM() & (1 << 0)))
464  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
465  if (!(cOre.getBuildMoFEM() & (1 << 1)))
466  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
467  if (!(cOre.getBuildMoFEM() & (1 << 2)))
468  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
469  const Problem *problem_ptr;
470  CHKERR m_field.get_problem(name, &problem_ptr);
471  CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
472  cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
473  // function knows what he is doing
475 }
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:481
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:600
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411

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

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

651  {
652  MoFEM::Interface &m_field = cOre;
654 
655  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
656  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
657  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
658  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
659  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
660  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
661 
662  const Problem *problem_ptr;
663  CHKERR m_field.get_problem(name, &problem_ptr);
664  CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
665  square_matrix, verb);
666 
669 
671 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
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:600
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:411

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

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

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

2883  {
2884  MoFEM::Interface &m_field = cOre;
2885  const Problem *problem_ptr;
2887  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2888  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2889  auto fit =
2890  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2891  .lower_bound(fe_name);
2892  auto hi_fe_it =
2893  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2894  .upper_bound(fe_name);
2895  std::vector<EntityHandle> fe_vec;
2896  fe_vec.reserve(std::distance(fit, hi_fe_it));
2897  for (; fit != hi_fe_it; fit++)
2898  fe_vec.push_back(fit->get()->getEnt());
2899  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2900  fe_vec.size());
2902 }
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:481
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:600
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411

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

2907  {
2908  MoFEM::Interface &m_field = cOre;
2909  const Problem *problem_ptr;
2911  CHKERR m_field.get_problem(name, &problem_ptr);
2912  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2913  fe_name, layout);
2915 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
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:600
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411

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

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

2419  {
2420  MoFEM::Interface &m_field = cOre;
2421  const Problem_multiIndex *problems_ptr;
2422  const EntFiniteElement_multiIndex *fe_ent_ptr;
2424 
2426  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2427 
2428  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2429  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2430 
2431  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2432  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2433  "adjacencies not build");
2434 
2436  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2437 
2439  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2440  "problem not partitioned");
2441 
2442  if (low_proc == -1)
2443  low_proc = m_field.get_comm_rank();
2444  if (hi_proc == -1)
2445  hi_proc = m_field.get_comm_rank();
2446 
2447  // Find pointer to problem of given name
2449  // Get pointer to problem data struture
2450  CHKERR m_field.get_problems(&problems_ptr);
2451  ProblemByName &problems =
2452  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2453  ProblemByName::iterator p_miit = problems.find(name);
2454  if (p_miit == problems.end()) {
2455  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2456  "problem < %s > not found (top tip: check spelling)",
2457  name.c_str());
2458  }
2459 
2460  // Get reference on finite elements multi-index on the problem
2461  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2462  *p_miit->numeredFiniteElements;
2463 
2464  // Clear all elements and data, build it again
2465  problem_finite_elements.clear();
2466 
2467  // Check if dofs and columns are the same, i.e. structurally symmetric problem
2468  bool do_cols_prob = true;
2469  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
2470  do_cols_prob = false;
2471  }
2472 
2473  // Loop over all elements in database and if right element is there add it
2474  // to problem finite element multi-index
2475  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
2476  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); efit++) {
2477 
2478  // if element is not part of problem
2479  if (((*efit)->getId() & p_miit->getBitFEId()).none())
2480  continue;
2481 
2482  BitRefLevel prb_bit = p_miit->getBitRefLevel();
2483  BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
2484  BitRefLevel fe_bit = (*efit)->getBitRefLevel();
2485  // if entity is not problem refinement level
2486  if ((fe_bit & prb_mask) != fe_bit)
2487  continue;
2488  if ((fe_bit & prb_bit) != prb_bit)
2489  continue;
2490 
2491  // create element
2492  boost::shared_ptr<NumeredEntFiniteElement> numered_fe(
2493  new NumeredEntFiniteElement(*efit));
2494 
2495  // check if rows and columns are the same on this element
2496  bool do_cols_fe = true;
2497  if ((numered_fe->sPtr->row_field_ents_view ==
2498  numered_fe->sPtr->col_field_ents_view) &&
2499  !do_cols_prob) {
2500  do_cols_fe = false;
2501  numered_fe->cols_dofs = numered_fe->rows_dofs;
2502  } else {
2503  // different dofs on rows and columns
2504  numered_fe->cols_dofs = boost::shared_ptr<FENumeredDofEntity_multiIndex>(
2506  }
2507  // get pointer to dofs multi-index on rows and columns
2508  auto rows_dofs = numered_fe->rows_dofs;
2509  auto cols_dofs = numered_fe->cols_dofs;
2510  // clear multi-indices
2511  rows_dofs->clear();
2512  if (do_cols_fe) {
2513  cols_dofs->clear();
2514  }
2517 
2518  // set partition to the element
2519  {
2520  if (part_from_moab) {
2521  // if partition is taken from moab partition
2522  int proc = (*efit)->getPartProc();
2523  if (proc == -1 && (*efit)->getEntType() == MBVERTEX) {
2524  proc = (*efit)->getOwnerProc();
2525  }
2526  NumeredEntFiniteElement_change_part(proc).operator()(numered_fe);
2527  } else {
2528  // count partition of the dofs in row, the larges dofs with given
2529  // partition is used to set partition of the element
2530  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows), rows_view,
2531  moab::Interface::UNION);
2532  std::vector<int> parts(m_field.get_comm_size(), 0);
2533  NumeredDofEntity_multiIndex_uid_view_ordered::iterator viit_rows;
2534  viit_rows = rows_view.begin();
2535  for (; viit_rows != rows_view.end(); viit_rows++) {
2536  parts[(*viit_rows)->pArt]++;
2537  }
2538  std::vector<int>::iterator pos =
2539  max_element(parts.begin(), parts.end());
2540  unsigned int max_part = std::distance(parts.begin(), pos);
2541  NumeredEntFiniteElement_change_part(max_part).operator()(numered_fe);
2542  }
2543  }
2544 
2545  // set dofs on rows and columns (if are different)
2546  if ((numered_fe->getPart() >= (unsigned int)low_proc) &&
2547  (numered_fe->getPart() <= (unsigned int)hi_proc)) {
2548 
2549  NumeredDofEntity_multiIndex_uid_view_ordered *dofs_view[] = {&rows_view,
2550  &cols_view};
2551  FENumeredDofEntity_multiIndex *fe_dofs[] = {rows_dofs.get(),
2552  cols_dofs.get()};
2553 
2554  for (int ss = 0; ss != (do_cols_fe ? 2 : 1); ss++) {
2555 
2556  if (ss == 0) {
2557  if (part_from_moab) {
2558  // get row_view
2559  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows),
2560  *dofs_view[ss],
2561  moab::Interface::UNION);
2562  }
2563  } else {
2564  // get cols_views
2565  CHKERR(*efit)->getColDofView(*(p_miit->numeredDofsCols),
2566  *dofs_view[ss], moab::Interface::UNION);
2567  }
2568 
2569  // Following reserve memory in sequences, only two allocations are here,
2570  // once for array of objects, next for array of shared pointers
2571 
2572  // reserve memory for field dofs
2573  boost::shared_ptr<std::vector<FENumeredDofEntity>> dofs_array(
2574  new std::vector<FENumeredDofEntity>());
2575 
2576  if (!ss) {
2577  numered_fe->getRowDofsSequence() = dofs_array;
2578  if (!do_cols_fe)
2579  numered_fe->getColDofsSequence() = dofs_array;
2580  } else
2581  numered_fe->getColDofsSequence() = dofs_array;
2582 
2583  auto vit = dofs_view[ss]->begin();
2584  auto hi_vit = dofs_view[ss]->end();
2585 
2586  dofs_array->reserve(std::distance(vit, hi_vit));
2587 
2588  // create elements objects
2589  for (; vit != hi_vit; vit++) {
2590  boost::shared_ptr<SideNumber> side_number_ptr;
2591  side_number_ptr = (*efit)->getSideNumberPtr((*vit)->getEnt());
2592  dofs_array->emplace_back(side_number_ptr, *vit);
2593  }
2594 
2595  // finally add DoFS to multi-indices
2596  auto hint = fe_dofs[ss]->end();
2597  for (auto &v : *dofs_array)
2598  hint = fe_dofs[ss]->emplace_hint(hint, dofs_array, &v);
2599  }
2600  }
2601  if (!numered_fe->sPtr->row_field_ents_view->empty() &&
2602  !numered_fe->sPtr->col_field_ents_view->empty()) {
2603 
2604  // Add element to the problem
2605  auto p = problem_finite_elements.insert(numered_fe);
2606  if (!p.second)
2607  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2608 
2609  if (verb >= VERY_VERBOSE) {
2610  std::ostringstream ss;
2611  ss << *p_miit << std::endl;
2612  ss << *p.first << std::endl;
2614  FENumeredDofEntityByUId::iterator miit =
2615  (*p.first)->rows_dofs->get<Unique_mi_tag>().begin();
2616  for (; miit != (*p.first)->rows_dofs->get<Unique_mi_tag>().end();
2617  miit++)
2618  ss << "rows: " << *(*miit) << std::endl;
2619  miit = (*p.first)->cols_dofs->get<Unique_mi_tag>().begin();
2620  for (; miit != (*p.first)->cols_dofs->get<Unique_mi_tag>().end();
2621  miit++)
2622  ss << "cols: " << *(*miit) << std::endl;
2623  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2624  }
2625  }
2626  }
2627 
2628  if (verb >= VERBOSE) {
2630  NumeredEntFiniteElementPart;
2631  NumeredEntFiniteElementPart::iterator miit, hi_miit;
2632  miit = problem_finite_elements.get<Part_mi_tag>().lower_bound(
2633  m_field.get_comm_rank());
2634  hi_miit = problem_finite_elements.get<Part_mi_tag>().upper_bound(
2635  m_field.get_comm_rank());
2636  int count = std::distance(miit, hi_miit);
2637  std::ostringstream ss;
2638  ss << *p_miit;
2639  ss << " Nb. elems " << count << " on proc " << m_field.get_comm_rank()
2640  << std::endl;
2641  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2642  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2643  }
2644 
2647 }
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:481
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:600
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:411
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_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 2649 of file ProblemsManager.cpp.

2650  {
2651  MoFEM::Interface &m_field = cOre;
2652  const Problem_multiIndex *problems_ptr;
2654 
2656  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2657  "partition of problem not build");
2659  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2660  "partitions finite elements not build");
2661 
2662  // find p_miit
2663  CHKERR m_field.get_problems(&problems_ptr);
2664 
2665  // get problem pointer
2666  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2667  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2668  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2669  name.c_str());
2670 
2671  // get reference to number of ghost dofs
2672  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2673  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2674  nb_row_ghost_dofs = 0;
2675  nb_col_ghost_dofs = 0;
2676 
2677  // do work if more than one processor
2678  if (m_field.get_comm_size() > 1) {
2679 
2681  ghost_idx_col_view;
2682 
2683  // get elements on this partition
2684  auto fe_range =
2685  p_miit->numeredFiniteElements->get<Part_mi_tag>().equal_range(
2686  m_field.get_comm_rank());
2687 
2688  // get dofs on elements which are not part of this partition
2689 
2690  // rows
2691  auto hint_r = ghost_idx_row_view.begin();
2692  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2693  for (auto &dof_ptr : *(*fe_ptr)->rows_dofs) {
2694  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2695  hint_r = ghost_idx_row_view.emplace_hint(
2696  hint_r, dof_ptr->getNumeredDofEntityPtr());
2697  }
2698  }
2699  }
2700 
2701  // columns
2702  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2703  auto hint_c = ghost_idx_col_view.begin();
2704  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2705  for (auto &dof_ptr : *(*fe_range.first)->cols_dofs) {
2706  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2707  hint_c = ghost_idx_col_view.emplace_hint(
2708  hint_c, dof_ptr->getNumeredDofEntityPtr());
2709  }
2710  }
2711  }
2712  }
2713 
2714  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2715  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2716 
2717  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2718  &ghost_idx_row_view, &ghost_idx_col_view};
2719  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2720  &p_miit->numeredDofsRows->get<Unique_mi_tag>(),
2721  &p_miit->numeredDofsCols->get<Unique_mi_tag>()};
2722 
2723  int loop_size = 2;
2724  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2725  loop_size = 1;
2726  }
2727 
2728  // set local ghost dofs indices
2729  for (int ss = 0; ss != loop_size; ++ss) {
2730  for (auto &gid : *ghost_idx_view[ss]) {
2731  NumeredDofEntityByUId::iterator dof =
2732  dof_by_uid_no_const[ss]->find(gid->getGlobalUniqueId());
2733  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2734  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2735  "inconsistent data, ghost dof already set");
2736  bool success = dof_by_uid_no_const[ss]->modify(
2737  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2738  if (PetscUnlikely(!success))
2739  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2740  "modification unsuccessful");
2741  (*nb_ghost_dofs[ss])++;
2742  }
2743  }
2744  if (loop_size == 1) {
2745  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2746  }
2747  }
2748 
2749  if (verb > QUIET) {
2750  std::ostringstream ss;
2751  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2752  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2753  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2754  << p_miit->getNbLocalDofsCol() << std::endl;
2755  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2756  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2757  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2758  << p_miit->getNbLocalDofsRow() << std::endl;
2759  if (verb > VERBOSE) {
2760  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2761  p_miit->numeredDofsCols->begin();
2762  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2763  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2764  continue;
2765  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2766  continue;
2767  ss << *(*miit_dd_col) << std::endl;
2768  }
2769  }
2770  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2771  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2772  }
2773 
2776 }
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:481
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:600
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:411
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 2779 of file ProblemsManager.cpp.

2780  {
2781  MoFEM::Interface &m_field = cOre;
2782  const Problem_multiIndex *problems_ptr;
2784 
2786  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2787  "partition of problem not build");
2789  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2790  "partitions finite elements not build");
2791 
2792  // get problem pointer
2793  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2794  // find p_miit
2795  CHKERR m_field.get_problems(&problems_ptr);
2796  ProblemsByName &problems_set =
2797  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2798  ProblemsByName::iterator p_miit = problems_set.find(name);
2799 
2800  // get reference to number of ghost dofs
2801  // get number of local dofs
2802  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2803  &(p_miit->nbGhostDofsCol)};
2804  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2805  for (int ss = 0; ss != 2; ++ss) {
2806  (*nb_ghost_dofs[ss]) = 0;
2807  }
2808 
2809  // do work if more than one processor
2810  if (m_field.get_comm_size() > 1) {
2811  // determine if rows on columns are different from dofs on rows
2812  int loop_size = 2;
2813  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2814  loop_size = 1;
2815  }
2816 
2817  typedef decltype(p_miit->numeredDofsRows) NumbDofTypeSharedPtr;
2818  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredDofsRows,
2819  p_miit->numeredDofsCols};
2820 
2821  // interate over dofs on rows and dofs on columns
2822  for (int ss = 0; ss != loop_size; ++ss) {
2823 
2824  // create dofs view by uid
2825  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2826 
2827  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2828  ghost_idx_view.reserve(std::distance(r.first, r.second));
2829  for (; r.first != r.second; ++r.first)
2830  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2831 
2832  auto cmp = [](auto a, auto b) {
2833  return (*a)->getGlobalUniqueId() < (*b)->getGlobalUniqueId();
2834  };
2835  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2836 
2837  // intare over dofs which have negative local index
2838  for (auto gid_it : ghost_idx_view) {
2839  bool success = numered_dofs[ss]->modify(
2840  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2841  if (!success)
2842  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2843  "modification unsuccessful");
2844  ++(*nb_ghost_dofs[ss]);
2845  }
2846  }
2847  if (loop_size == 1) {
2848  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2849  }
2850  }
2851 
2852  if (verb > QUIET) {
2853  std::ostringstream ss;
2854  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2855  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2856  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2857  << p_miit->getNbLocalDofsRow() << std::endl;
2858  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2859  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2860  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2861  << p_miit->getNbLocalDofsCol() << std::endl;
2862  if (verb > VERBOSE) {
2863  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2864  p_miit->numeredDofsCols->begin();
2865  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2866  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2867  continue;
2868  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2869  continue;
2870  ss << *(*miit_dd_col) << std::endl;
2871  }
2872  }
2873  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2874  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2875  }
2876 
2879 }
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:481
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:600
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:411
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 96 of file ProblemsManager.cpp.

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

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

Definition at line 1804 of file ProblemsManager.cpp.

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

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