v0.8.23
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 1448 of file ProblemsManager.cpp.

1451  {
1453  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1454  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1455  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1456  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1457  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1458  MoFEM::Interface &m_field = cOre;
1459  const Problem_multiIndex *problems_ptr;
1461 
1462  CHKERR m_field.clear_problem(out_name);
1463  CHKERR m_field.get_problems(&problems_ptr);
1464  // get reference to all problems
1465  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1466  ProblemByName &problems_by_name =
1467  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1468 
1469  // Get iterators to out problem, i.e. build problem
1470  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1471  if (out_problem_it == problems_by_name.end()) {
1472  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1473  "problem with name < %s > not defined (top tip check spelling)",
1474  out_name.c_str());
1475  }
1476  // Make data structure for composed-problem data
1477  out_problem_it->composedProblemsData =
1478  boost::make_shared<ComposedProblemsData>();
1479  boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1480  out_problem_it->getComposedProblemsData();
1481 
1482  const std::vector<std::string> *add_prb[] = {&add_row_problems,
1483  &add_col_problems};
1484  std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1485  &cmp_prb_data->colProblemsAdd};
1486  std::vector<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1487 
1488  // Get local indices counter
1489  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1490  &out_problem_it->nbLocDofsCol};
1491  // Get global indices counter
1492  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1493 
1494  // Set number of ghost nodes to zero
1495  {
1496  out_problem_it->nbDofsRow = 0;
1497  out_problem_it->nbDofsCol = 0;
1498  out_problem_it->nbLocDofsRow = 0;
1499  out_problem_it->nbLocDofsCol = 0;
1500  out_problem_it->nbGhostDofsRow = 0;
1501  out_problem_it->nbGhostDofsCol = 0;
1502  }
1503  int nb_dofs_reserve[] = {0, 0};
1504 
1505  // Loop over rows and columns in the main problem and sub-problems
1506  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1507  // cerr << "SS " << ss << endl;
1508  // cerr << add_prb[ss]->size() << endl;
1509  // cerr << add_prb_ptr[ss]->size() << endl;
1510  add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1511  add_prb_is[ss]->reserve(add_prb[ss]->size());
1512  for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1513  vit != add_prb[ss]->end(); vit++) {
1514  ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1515  if (prb_it == problems_by_name.end()) {
1516  SETERRQ1(
1517  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1518  "problem with name < %s > not defined (top tip check spelling)",
1519  vit->c_str());
1520  }
1521  add_prb_ptr[ss]->push_back(&*prb_it);
1522  // set number of dofs on rows and columns
1523  if (ss == 0) {
1524  // row
1525  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1526  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1527  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredDofsRows->size();
1528  } else {
1529  // column
1530  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1531  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1532  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredDofsCols->size();
1533  }
1534  }
1535  }
1536  // if squre problem, rows and columns are the same
1537  if (square_matrix) {
1538  add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1539  add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1540  out_problem_it->numeredDofsCols = out_problem_it->numeredDofsRows;
1541  *nb_dofs[1] = *nb_dofs[0];
1542  *nb_local_dofs[1] = *nb_local_dofs[0];
1543  }
1544 
1545  // reserve memory for dofs
1546  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1547  // Reserve memory
1548  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1549  dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1550  dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1551  if (!ss)
1552  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1553  else
1554  out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1555  }
1556 
1557  // Push back DOFs
1558  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1559  NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1560  dit,
1561  hi_dit;
1562  int shift_glob = 0;
1563  int shift_loc = 0;
1564  for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1565  PetscInt *dofs_out_idx_ptr;
1566  int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1567  CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1568  if (ss == 0) {
1569  dit = (*add_prb_ptr[ss])[pp]
1570  ->numeredDofsRows->get<PetscGlobalIdx_mi_tag>()
1571  .begin();
1572  hi_dit = (*add_prb_ptr[ss])[pp]
1573  ->numeredDofsRows->get<PetscGlobalIdx_mi_tag>()
1574  .end();
1575  } else {
1576  dit = (*add_prb_ptr[ss])[pp]
1577  ->numeredDofsCols->get<PetscGlobalIdx_mi_tag>()
1578  .begin();
1579  hi_dit = (*add_prb_ptr[ss])[pp]
1580  ->numeredDofsCols->get<PetscGlobalIdx_mi_tag>()
1581  .end();
1582  }
1583  int is_nb = 0;
1584  for (; dit != hi_dit; dit++) {
1585  BitRefLevel prb_bit = out_problem_it->getBitRefLevel();
1586  BitRefLevel prb_mask = out_problem_it->getMaskBitRefLevel();
1587  BitRefLevel dof_bit = dit->get()->getBitRefLevel();
1588  if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1589  continue;
1590  const int rank = m_field.get_comm_rank();
1591  const int part = dit->get()->getPart();
1592  const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1593  const int loc_idx =
1594  (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1595  : -1;
1596  dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1597  glob_idx, loc_idx, part);
1598  if (part == rank) {
1599  dofs_out_idx_ptr[is_nb++] = glob_idx;
1600  }
1601  }
1602  if (is_nb > nb_local_dofs) {
1603  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1604  "Data inconsistency");
1605  }
1606  IS is;
1607  CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1608  PETSC_OWN_POINTER, &is);
1609  (*add_prb_is[ss]).push_back(is);
1610  if (ss == 0) {
1611  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1612  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1613  } else {
1614  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1615  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1616  }
1617  if (square_matrix) {
1618  (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1619  (*add_prb_is[1]).push_back(is);
1620  CHKERR PetscObjectReference((PetscObject)is);
1621  }
1622  }
1623  }
1624 
1625  if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1626  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1627  }
1628 
1629  // Insert DOFs to problem multi-index
1630  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1631  auto hint = (ss == 0) ? out_problem_it->numeredDofsRows->end()
1632  : out_problem_it->numeredDofsCols->end();
1633  for (auto &v : *dofs_array[ss])
1634  hint = (ss == 0) ? out_problem_it->numeredDofsRows->emplace_hint(
1635  hint, dofs_array[ss], &v)
1636  : out_problem_it->numeredDofsCols->emplace_hint(
1637  hint, dofs_array[ss], &v);
1638  }
1639 
1640  // Compress DOFs
1641  *nb_dofs[0] = 0;
1642  *nb_dofs[1] = 0;
1643  *nb_local_dofs[0] = 0;
1644  *nb_local_dofs[1] = 0;
1645  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1646 
1647  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1648  if (ss == 0) {
1649  dofs_ptr = out_problem_it->numeredDofsRows;
1650  } else {
1651  dofs_ptr = out_problem_it->numeredDofsCols;
1652  }
1653  NumeredDofEntityByUId::iterator dit, hi_dit;
1654  dit = dofs_ptr->get<Unique_mi_tag>().begin();
1655  hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1656  std::vector<int> idx;
1657  idx.reserve(std::distance(dit, hi_dit));
1658  // set dofs in order entity and dof number on entity
1659  for (; dit != hi_dit; dit++) {
1660  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1661  bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1662  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1663  if (!success) {
1664  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1665  "modification unsuccessful");
1666  }
1667  idx.push_back(dit->get()->getPetscGlobalDofIdx());
1668  } else {
1669  if (dit->get()->getPetscLocalDofIdx() != -1) {
1670  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1671  "local index should be negative");
1672  }
1673  }
1674  }
1675  if (square_matrix) {
1676  *nb_local_dofs[1] = *nb_local_dofs[0];
1677  }
1678 
1679  // set new dofs mapping
1680  IS is;
1681  CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1682  PETSC_USE_POINTER, &is);
1683  CHKERR ISGetSize(is, nb_dofs[ss]);
1684  if (square_matrix) {
1685  *nb_dofs[1] = *nb_dofs[0];
1686  }
1687  // {
1688  // PetscSynchronizedPrintf(m_field.get_comm(),"nb dofs %d %d
1689  // %d\n",*nb_local_dofs[0],idx.size(),*nb_dofs[0]);
1690  // PetscSynchronizedFlush(m_field.get_comm(),PETSC_STDOUT);
1691  // }
1692  AO ao;
1693  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1694  for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1695  CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1696 
1697  // Set DOFs numeration
1698  {
1699  std::vector<int> idx_new;
1700  idx_new.reserve(dofs_ptr->size());
1701  for (NumeredDofEntityByUId::iterator dit =
1702  dofs_ptr->get<Unique_mi_tag>().begin();
1703  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1704  idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1705  }
1706  // set new global dofs numeration
1707  IS is_new;
1708  CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1709  &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1710  CHKERR AOApplicationToPetscIS(ao, is_new);
1711  // set global indices to multi-index
1712  std::vector<int>::iterator vit = idx_new.begin();
1713  for (NumeredDofEntityByUId::iterator dit =
1714  dofs_ptr->get<Unique_mi_tag>().begin();
1715  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1716  bool success =
1717  dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1718  dit->get()->getPart(), *(vit++)));
1719  if (!success) {
1720  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1721  "modification unsuccessful");
1722  }
1723  }
1724  CHKERR ISDestroy(&is_new);
1725  }
1726  CHKERR ISDestroy(&is);
1727  CHKERR AODestroy(&ao);
1728  }
1729 
1730  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1731  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1732 
1733  // Inidcate that porble has been build
1736 
1738 }
virtual MoFEMErrorCode clear_problem(const std::string &name, int verb=DEFAULT_VERBOSITY)=0
clear problem
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
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:406
virtual MPI_Comm & get_comm() const =0

◆ 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, forces_and_sources_testing_edge_element.cpp, and forces_and_sources_testing_flat_prism_element.cpp.

Definition at line 394 of file ProblemsManager.cpp.

396  {
397 
398  MoFEM::Interface &m_field = cOre;
400  if (!(cOre.getBuildMoFEM() & (1 << 0)))
401  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
402  if (!(cOre.getBuildMoFEM() & (1 << 1)))
403  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
404  if (!(cOre.getBuildMoFEM() & (1 << 2)))
405  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
406  const Problem *problem_ptr;
407  CHKERR m_field.get_problem(name, &problem_ptr);
408  CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
409  cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
410  // function knows what he is doing
412 }
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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

416  {
417  MoFEM::Interface &m_field = cOre;
418  const EntFiniteElement_multiIndex *fe_ent_ptr;
419  const DofEntity_multiIndex *dofs_field_ptr;
421  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
422 
423  // Note: Only allowed changes on problem_ptr structure which not influence
424  // multi-index.
425 
426  if (problem_ptr->getBitRefLevel().none()) {
427  SETERRQ1(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
428  problem_ptr->getName().c_str());
429  }
430  CHKERR m_field.clear_problem(problem_ptr->getName());
431  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
432  CHKERR m_field.get_dofs(&dofs_field_ptr);
433 
434  // zero finite elements
435  problem_ptr->numeredFiniteElements->clear();
436 
437  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
438  {
439  EntFiniteElement_multiIndex::iterator miit = fe_ent_ptr->begin();
440  EntFiniteElement_multiIndex::iterator hi_miit = fe_ent_ptr->end();
441  // iterate all finite element entities in database
442  for (; miit != hi_miit; miit++) {
443  // if element is in problem
444  if (((*miit)->getId() & problem_ptr->getBitFEId()).any()) {
445  BitRefLevel prb_bit = problem_ptr->getBitRefLevel();
446  BitRefLevel prb_mask = problem_ptr->getMaskBitRefLevel();
447  BitRefLevel fe_bit = (*miit)->getBitRefLevel();
448  // if entity is not problem refinement level
449  if ((fe_bit & prb_mask) != fe_bit)
450  continue;
451  if ((fe_bit & prb_bit) != prb_bit)
452  continue;
453  // get dof uids for rows and columns
454  CHKERR(*miit)->getRowDofView(*dofs_field_ptr, dofs_rows);
455  if (!square_matrix) {
456  CHKERR(*miit)->getColDofView(*dofs_field_ptr, dofs_cols);
457  }
458  }
459  }
460  }
461 
462  // Add row dofs to problem
463  {
464  // zero rows
465  problem_ptr->nbDofsRow = 0;
466  problem_ptr->nbLocDofsRow = 0;
467  problem_ptr->nbGhostDofsRow = 0;
468 
469  // add dofs for rows
470  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
471  hi_miit;
472  hi_miit = dofs_rows.get<0>().end();
473 
474  int count_dofs = dofs_rows.get<1>().count(true);
475  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
476  boost::shared_ptr<std::vector<NumeredDofEntity>>(
477  new std::vector<NumeredDofEntity>());
478  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
479  dofs_array->reserve(count_dofs);
480  miit = dofs_rows.get<0>().begin();
481  for (; miit != hi_miit; miit++) {
482  if ((*miit)->getActive()) {
483  dofs_array->emplace_back(*miit);
484  dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
485  }
486  }
487  auto hint = problem_ptr->numeredDofsRows->end();
488  for (auto &v : *dofs_array) {
489  hint = problem_ptr->numeredDofsRows->emplace_hint(hint, dofs_array, &v);
490  }
491  }
492 
493  // Add col dofs to problem
494  if (!square_matrix) {
495  // zero cols
496  problem_ptr->nbDofsCol = 0;
497  problem_ptr->nbLocDofsCol = 0;
498  problem_ptr->nbGhostDofsCol = 0;
499 
500  // add dofs for cols
501  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
502  hi_miit;
503  hi_miit = dofs_cols.get<0>().end();
504 
505  int count_dofs = 0;
506  miit = dofs_cols.get<0>().begin();
507  for (; miit != hi_miit; miit++) {
508  if (!(*miit)->getActive()) {
509  continue;
510  }
511  count_dofs++;
512  }
513 
514  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
515  boost::shared_ptr<std::vector<NumeredDofEntity>>(
516  new std::vector<NumeredDofEntity>());
517  problem_ptr->getColDofsSequence()->push_back(dofs_array);
518  dofs_array->reserve(count_dofs);
519  miit = dofs_cols.get<0>().begin();
520  for (; miit != hi_miit; miit++) {
521  if (!(*miit)->getActive()) {
522  continue;
523  }
524  dofs_array->emplace_back(*miit);
525  dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
526  }
527  auto hint = problem_ptr->numeredDofsCols->end();
528  for (auto &v : *dofs_array) {
529  hint = problem_ptr->numeredDofsCols->emplace_hint(hint, dofs_array, &v);
530  }
531  } else {
532  problem_ptr->numeredDofsCols = problem_ptr->numeredDofsRows;
533  problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
534  problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
535  }
536 
537  // job done, some debugging and postprocessing
538  if (verb > QUIET) {
539  PetscSynchronizedPrintf(
540  m_field.get_comm(), "Problem %s Nb. rows %u Nb. cols %u\n",
541  problem_ptr->getName().c_str(), problem_ptr->numeredDofsRows->size(),
542  problem_ptr->numeredDofsCols->size());
543  }
544  if (verb > VERBOSE) {
545  EntFiniteElement_multiIndex::iterator miit = fe_ent_ptr->begin();
546  EntFiniteElement_multiIndex::iterator hi_miit = fe_ent_ptr->end();
547  std::ostringstream ss;
548  ss << "rank " << m_field.get_comm_rank() << " ";
549  ss << "FEs data for problem " << *problem_ptr << std::endl;
550  for (; miit != hi_miit; miit++) {
551  ss << "rank " << m_field.get_comm_rank() << " ";
552  ss << **miit << std::endl;
553  }
554  ss << "rank " << m_field.get_comm_rank() << " ";
555  ss << "FEs row dofs " << *problem_ptr << " Nb. row dof "
556  << problem_ptr->getNbDofsRow() << std::endl;
557  NumeredDofEntity_multiIndex::iterator miit_dd_row =
558  problem_ptr->numeredDofsRows->begin();
559  for (; miit_dd_row != problem_ptr->numeredDofsRows->end(); miit_dd_row++) {
560  ss << "rank " << m_field.get_comm_rank() << " ";
561  ss << **miit_dd_row << std::endl;
562  }
563  ss << "rank " << m_field.get_comm_rank() << " ";
564  ss << "FEs col dofs " << *problem_ptr << " Nb. col dof "
565  << problem_ptr->getNbDofsCol() << std::endl;
566  NumeredDofEntity_multiIndex::iterator miit_dd_col =
567  problem_ptr->numeredDofsCols->begin();
568  for (; miit_dd_col != problem_ptr->numeredDofsCols->end(); miit_dd_col++) {
569  ss << "rank " << m_field.get_comm_rank() << " ";
570  ss << **miit_dd_col << std::endl;
571  }
572  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
573  }
574 
575  if (verb > QUIET) {
576  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
577  }
578  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
579  // uses this function knows what
580  // he is doing
581 
582  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
583 
585 }
virtual MoFEMErrorCode clear_problem(const std::string &name, int verb=DEFAULT_VERBOSITY)=0
clear problem
PetscLogEvent MOFEM_EVENT_ProblemsManager
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
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:406
virtual MPI_Comm & get_comm() const =0

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

588  {
589  MoFEM::Interface &m_field = cOre;
591 
592  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
593  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
594  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
595  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
596  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
597  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
598 
599  const Problem *problem_ptr;
600  CHKERR m_field.get_problem(name, &problem_ptr);
601  CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
602  square_matrix, verb);
603 
606 
608 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
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:406

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

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

1196  {
1197  MoFEM::Interface &m_field = cOre;
1198  const Problem_multiIndex *problems_ptr;
1200 
1201  CHKERR m_field.clear_problem(out_name);
1202  CHKERR m_field.get_problems(&problems_ptr);
1203 
1204  // get reference to all problems
1205  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1206  ProblemByName &problems_by_name =
1207  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1208 
1209  // get iterators to out problem, i.e. build problem
1210  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1211  if (out_problem_it == problems_by_name.end()) {
1212  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1213  "subproblem with name < %s > not defined (top tip check spelling)",
1214  out_name.c_str());
1215  }
1216  // get iterator to main problem, i.e. out problem is subproblem of main
1217  // problem
1218  ProblemByName::iterator main_problem_it = problems_by_name.find(main_problem);
1219  if (main_problem_it == problems_by_name.end()) {
1220  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1221  "problem of subproblem with name < %s > not defined (top tip "
1222  "check spelling)",
1223  main_problem.c_str());
1224  }
1225 
1226  // get dofs for row & columns for out problem,
1227  boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1228  out_problem_it->numeredDofsRows, out_problem_it->numeredDofsCols};
1229  // get dofs for row & columns for main problem
1230  boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1231  main_problem_it->numeredDofsRows, main_problem_it->numeredDofsCols};
1232  // get local indices counter
1233  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1234  &out_problem_it->nbLocDofsCol};
1235  // get global indices counter
1236  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1237 
1238  // set number of ghost nodes to zero
1239  {
1240  out_problem_it->nbGhostDofsRow = 0;
1241  out_problem_it->nbGhostDofsCol = 0;
1242  }
1243 
1244  // put rows & columns field names in array
1245  std::vector<std::string> fields[] = {fields_row, fields_col};
1246  const map<std::string, std::pair<EntityType, EntityType>> *entityMap[] = {
1247  entityMapRow, entityMapCol};
1248 
1249  // make data structure fos sub-problem data
1250  out_problem_it->subProblemData =
1251  boost::make_shared<Problem::SubProblemData>();
1252 
1253  // Loop over rows and columns
1254  for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1255 
1256  // reset dofs and columns counters
1257  (*nb_local_dofs[ss]) = 0;
1258  (*nb_dofs[ss]) = 0;
1259  // clear arrays
1260  out_problem_dofs[ss]->clear();
1261 
1262  // If DOFs are cleared clear finite elements too.
1263  out_problem_it->numeredFiniteElements->clear();
1264 
1265  // get dofs by field name and insert them in out problem multi-indices
1266  for (auto field : fields[ss]) {
1267 
1268  // Following reserve memory in sequences, only two allocations are here,
1269  // once for array of objects, next for array of shared pointers
1270 
1271  // aliased sequence of pointer is killed with element
1272  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1273  boost::make_shared<std::vector<NumeredDofEntity>>();
1274  // reserve memory for field dofs
1275  if (!ss)
1276  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1277  else
1278  out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1279 
1280  // create elements objects
1281  auto dit =
1282  main_problem_dofs[ss]->get<FieldName_mi_tag>().lower_bound(field);
1283  auto hi_dit =
1284  main_problem_dofs[ss]->get<FieldName_mi_tag>().upper_bound(field);
1285 
1286  auto add_dit_to_dofs_array = [&](auto &dit) {
1287  dofs_array->emplace_back(
1288  dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1289  dit->get()->getPetscGlobalDofIdx(),
1290  dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1291  };
1292 
1293  if (entityMap[ss]) {
1294  auto mit = entityMap[ss]->find(field);
1295  if (mit != entityMap[ss]->end()) {
1296  EntityType lo_type = mit->second.first;
1297  EntityType hi_type = mit->second.second;
1298  int count = 0;
1299  for (auto diit = dit; diit != hi_dit; ++diit) {
1300  EntityType ent_type = (*diit)->getEntType();
1301  if (ent_type >= lo_type && ent_type <= hi_type)
1302  ++count;
1303  }
1304  dofs_array->reserve(count);
1305  for (; dit != hi_dit; ++dit) {
1306  EntityType ent_type = (*dit)->getEntType();
1307  if (ent_type >= lo_type && ent_type <= hi_type)
1308  add_dit_to_dofs_array(dit);
1309  }
1310  } else {
1311  dofs_array->reserve(std::distance(dit, hi_dit));
1312  for (; dit != hi_dit; dit++)
1313  add_dit_to_dofs_array(dit);
1314  }
1315  } else {
1316  dofs_array->reserve(std::distance(dit, hi_dit));
1317  for (; dit != hi_dit; dit++)
1318  add_dit_to_dofs_array(dit);
1319  }
1320 
1321  // fill multi-index
1322  auto hint = out_problem_dofs[ss]->end();
1323  for (auto &v : *dofs_array)
1324  hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1325  }
1326  // Set local indexes
1327  {
1328  auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1329  auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1330  for (; dit != hi_dit; dit++) {
1331  int idx = -1; // if dof is not part of partition, set local index to -1
1332  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1333  idx = (*nb_local_dofs[ss])++;
1334  }
1335  bool success = out_problem_dofs[ss]->modify(
1336  out_problem_dofs[ss]->project<0>(dit),
1337  NumeredDofEntity_local_idx_change(idx));
1338  if (!success) {
1339  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1340  "operation unsuccessful");
1341  }
1342  };
1343  }
1344  // Set global indexes, compress global indices
1345  {
1346  auto dit =
1347  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1348  auto hi_dit =
1349  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1350  out_problem_dofs[ss]->size());
1351  const int nb = std::distance(dit, hi_dit);
1352  // get main problem global indices
1353  std::vector<int> main_indices(nb);
1354  for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1355  *it = dit->get()->getPetscGlobalDofIdx();
1356  }
1357  // create is with global dofs
1358  IS is;
1359  CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1360  PETSC_USE_POINTER, &is);
1361  // create map form main problem global indices to out problem global
1362  // indices
1363  AO ao;
1364  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1365  if (ss == 0) {
1366  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->rowIs));
1367  // CHKERR ISSort(out_problem_it->getSubData()->rowIs);
1368  out_problem_it->getSubData()->rowMap = ao;
1369  CHKERR PetscObjectReference((PetscObject)ao);
1370  } else {
1371  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->colIs));
1372  // CHKERR ISSort(out_problem_it->getSubData()->colIs);
1373  out_problem_it->getSubData()->colMap = ao;
1374  CHKERR PetscObjectReference((PetscObject)ao);
1375  }
1376  CHKERR AOApplicationToPetscIS(ao, is);
1377  // set global number of DOFs
1378  CHKERR ISGetSize(is, nb_dofs[ss]);
1379  CHKERR ISDestroy(&is);
1380  // set out problem global indices after applying map
1381  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1382  for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1383  dit++, it++) {
1384  bool success = out_problem_dofs[ss]->modify(
1385  out_problem_dofs[ss]->project<0>(dit),
1386  NumeredDofEntity_part_and_all_indices_change(
1387  dit->get()->getPart(), *it, *it,
1388  dit->get()->getPetscLocalDofIdx()));
1389  if (!success) {
1390  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1391  "operation unsuccessful");
1392  }
1393  }
1394  // set global indices to nodes not on this part
1395  {
1396  NumeredDofEntityByLocalIdx::iterator dit =
1397  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1398  NumeredDofEntityByLocalIdx::iterator hi_dit =
1399  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1400  const int nb = std::distance(dit, hi_dit);
1401  std::vector<int> main_indices_non_local(nb);
1402  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1403  dit++, it++) {
1404  *it = dit->get()->getPetscGlobalDofIdx();
1405  }
1406  IS is;
1407  CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1408  &*main_indices_non_local.begin(),
1409  PETSC_USE_POINTER, &is);
1410  CHKERR AOApplicationToPetscIS(ao, is);
1411  CHKERR ISDestroy(&is);
1412  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1413  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1414  dit++, it++) {
1415  bool success = out_problem_dofs[ss]->modify(
1416  out_problem_dofs[ss]->project<0>(dit),
1417  NumeredDofEntity_part_and_all_indices_change(
1418  dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1419  dit->get()->getPetscLocalDofIdx()));
1420  if (!success) {
1421  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1422  "operation unsuccessful");
1423  }
1424  }
1425  }
1426  CHKERR AODestroy(&ao);
1427  }
1428  }
1429 
1430  if (square_matrix) {
1431  out_problem_it->numeredDofsCols = out_problem_it->numeredDofsRows;
1432  out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1433  out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1434  out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1435  out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1436  CHKERR PetscObjectReference(
1437  (PetscObject)out_problem_it->getSubData()->rowIs);
1438  CHKERR PetscObjectReference(
1439  (PetscObject)out_problem_it->getSubData()->rowMap);
1440  }
1441 
1442  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1443  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1444 
1446 }
virtual MoFEMErrorCode clear_problem(const std::string &name, int verb=DEFAULT_VERBOSITY)=0
clear problem
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
virtual MPI_Comm & get_comm() const =0

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

2819  {
2820  MoFEM::Interface &m_field = cOre;
2821  const Problem *problem_ptr;
2823  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2824  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2825  auto fit =
2826  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2827  .lower_bound(fe_name);
2828  auto hi_fe_it =
2829  problem_ptr->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
2830  .upper_bound(fe_name);
2831  std::vector<EntityHandle> fe_vec;
2832  fe_vec.reserve(std::distance(fit, hi_fe_it));
2833  for (; fit != hi_fe_it; fit++)
2834  fe_vec.push_back(fit->get()->getEnt());
2835  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2836  fe_vec.size());
2838 }
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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

2843  {
2844  MoFEM::Interface &m_field = cOre;
2845  const Problem *problem_ptr;
2847  CHKERR m_field.get_problem(name, &problem_ptr);
2848  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2849  fe_name, layout);
2851 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

2076  {
2077 
2078  MoFEM::Interface &m_field = cOre;
2079  const Problem_multiIndex *problems_ptr;
2081 
2083  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
2084 
2085  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2086 
2087  // find p_miit
2088  CHKERR m_field.get_problems(&problems_ptr);
2089  ProblemByName &problems_by_name =
2090  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2091  ProblemByName::iterator p_miit = problems_by_name.find(name);
2092  if (p_miit == problems_by_name.end()) {
2093  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2094  "problem with name < %s > not defined (top tip check spelling)",
2095  name.c_str());
2096  }
2097  if (verb > QUIET)
2098  PetscPrintf(m_field.get_comm(),
2099  "Compose problem %s from rows of %s and columns of %s\n",
2100  p_miit->getName().c_str(), problem_for_rows.c_str(),
2101  problem_for_cols.c_str());
2102 
2103  // find p_miit_row
2104  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
2105  if (p_miit_row == problems_by_name.end()) {
2106  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2107  "problem with name < %s > not defined (top tip check spelling)",
2108  problem_for_rows.c_str());
2109  }
2110  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
2111  p_miit_row->numeredDofsRows;
2112 
2113  // find p_mit_col
2114  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
2115  if (p_miit_col == problems_by_name.end()) {
2116  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2117  "problem with name < %s > not defined (top tip check spelling)",
2118  problem_for_cols.c_str());
2119  }
2120  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
2121  p_miit_col->numeredDofsCols;
2122 
2123  bool copy[] = {copy_rows, copy_cols};
2124  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
2125  p_miit->numeredDofsRows, p_miit->numeredDofsCols};
2126 
2127  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
2128  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
2129  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
2130  dofs_col};
2131 
2132  for (int ss = 0; ss < 2; ss++) {
2133 
2134  // build indices
2135  *nb_local_dofs[ss] = 0;
2136  if (!copy[ss]) {
2137 
2138  // only copy indices which are belong to some elements if this problem
2139  std::vector<int> is_local, is_new;
2140 
2141  NumeredDofEntityByUId &dofs_by_uid =
2142  copied_dofs[ss]->get<Unique_mi_tag>();
2143  for (NumeredDofEntity_multiIndex::iterator dit =
2144  composed_dofs[ss]->begin();
2145  dit != composed_dofs[ss]->end(); dit++) {
2146  NumeredDofEntityByUId::iterator diit =
2147  dofs_by_uid.find((*dit)->getGlobalUniqueId());
2148  if (diit == dofs_by_uid.end()) {
2149  SETERRQ(
2151  "data inconsistency, could not find dof in composite problem");
2152  }
2153  int part_number = (*diit)->getPart(); // get part number
2154  int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
2155  bool success;
2156  success = composed_dofs[ss]->modify(
2157  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2158  petsc_global_dof));
2159  if (!success) {
2160  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2161  "modification unsuccessful");
2162  }
2163  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2164  success = composed_dofs[ss]->modify(
2165  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
2166  if (!success) {
2167  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2168  "modification unsuccessful");
2169  }
2170  is_local.push_back(petsc_global_dof);
2171  }
2172  }
2173 
2174  AO ao;
2175  CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
2176  NULL, &ao);
2177 
2178  // apply local to global mapping
2179  is_local.resize(0);
2180  for (NumeredDofEntity_multiIndex::iterator dit =
2181  composed_dofs[ss]->begin();
2182  dit != composed_dofs[ss]->end(); dit++) {
2183  is_local.push_back((*dit)->getPetscGlobalDofIdx());
2184  }
2185  CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2186  int idx2 = 0;
2187  for (NumeredDofEntity_multiIndex::iterator dit =
2188  composed_dofs[ss]->begin();
2189  dit != composed_dofs[ss]->end(); dit++) {
2190  int part_number = (*dit)->getPart(); // get part number
2191  int petsc_global_dof = is_local[idx2++];
2192  bool success;
2193  success = composed_dofs[ss]->modify(
2194  dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2195  petsc_global_dof));
2196  if (!success) {
2197  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2198  "modification unsuccessful");
2199  }
2200  }
2201 
2202  CHKERR AODestroy(&ao);
2203 
2204  } else {
2205 
2206  for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2207  dit != copied_dofs[ss]->end(); dit++) {
2208  std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2209  p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2210  new NumeredDofEntity((*dit)->getDofEntityPtr())));
2211  if (p.second) {
2212  (*nb_dofs[ss])++;
2213  }
2214  int dof_idx = (*dit)->getDofIdx();
2215  int part_number = (*dit)->getPart(); // get part number
2216  int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2217  if (part_number == (unsigned int)m_field.get_comm_rank()) {
2218  const bool success = composed_dofs[ss]->modify(
2219  p.first, NumeredDofEntity_part_and_all_indices_change(
2220  part_number, dof_idx, petsc_global_dof,
2221  (*nb_local_dofs[ss])++));
2222  if (!success) {
2223  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2224  "modification unsuccessful");
2225  }
2226  } else {
2227  const bool success = composed_dofs[ss]->modify(
2228  p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
2229  part_number, dof_idx, petsc_global_dof));
2230  if (!success) {
2231  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2232  "modification unsuccessful");
2233  }
2234  }
2235  }
2236  }
2237  }
2238 
2239  CHKERR printPartitionedProblem(&*p_miit, verb);
2240  CHKERR debugPartitionedProblem(&*p_miit, verb);
2241 
2243 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
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, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, and remove_entities_from_problem.cpp.

Definition at line 2352 of file ProblemsManager.cpp.

2355  {
2356  MoFEM::Interface &m_field = cOre;
2357  const Problem_multiIndex *problems_ptr;
2358  const EntFiniteElement_multiIndex *fe_ent_ptr;
2360 
2362  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2363 
2364  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2365  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2366 
2367  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2368  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2369  "adjacencies not build");
2370 
2372  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2373 
2375  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2376  "problem not partitioned");
2377 
2378  if (low_proc == -1)
2379  low_proc = m_field.get_comm_rank();
2380  if (hi_proc == -1)
2381  hi_proc = m_field.get_comm_rank();
2382 
2383  // Find pointer to problem of given name
2384  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2385  // Get pointer to problem data struture
2386  CHKERR m_field.get_problems(&problems_ptr);
2387  ProblemByName &problems =
2388  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2389  ProblemByName::iterator p_miit = problems.find(name);
2390  if (p_miit == problems.end()) {
2391  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2392  "problem < %s > not found (top tip: check spelling)",
2393  name.c_str());
2394  }
2395 
2396  // Get reference on finite elements multi-index on the problem
2397  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2398  *p_miit->numeredFiniteElements;
2399 
2400  // Clear all elements and data, build it again
2401  problem_finite_elements.clear();
2402 
2403  // Check if dofs and columns are the same, i.e. structurally symmetric problem
2404  bool do_cols_prob = true;
2405  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
2406  do_cols_prob = false;
2407  }
2408 
2409  // Loop over all elements in database and if right element is there add it
2410  // to problem finite element multi-index
2411  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
2412  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); efit++) {
2413 
2414  // if element is not part of problem
2415  if (((*efit)->getId() & p_miit->getBitFEId()).none())
2416  continue;
2417 
2418  BitRefLevel prb_bit = p_miit->getBitRefLevel();
2419  BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
2420  BitRefLevel fe_bit = (*efit)->getBitRefLevel();
2421  // if entity is not problem refinement level
2422  if ((fe_bit & prb_mask) != fe_bit)
2423  continue;
2424  if ((fe_bit & prb_bit) != prb_bit)
2425  continue;
2426 
2427  // create element
2428  boost::shared_ptr<NumeredEntFiniteElement> numered_fe(
2429  new NumeredEntFiniteElement(*efit));
2430 
2431  // check if rows and columns are the same on this element
2432  bool do_cols_fe = true;
2433  if ((numered_fe->sPtr->row_field_ents_view ==
2434  numered_fe->sPtr->col_field_ents_view) &&
2435  !do_cols_prob) {
2436  do_cols_fe = false;
2437  numered_fe->cols_dofs = numered_fe->rows_dofs;
2438  } else {
2439  // different dofs on rows and columns
2440  numered_fe->cols_dofs = boost::shared_ptr<FENumeredDofEntity_multiIndex>(
2442  }
2443  // get pointer to dofs multi-index on rows and columns
2444  auto rows_dofs = numered_fe->rows_dofs;
2445  auto cols_dofs = numered_fe->cols_dofs;
2446  // clear multi-indices
2447  rows_dofs->clear();
2448  if (do_cols_fe) {
2449  cols_dofs->clear();
2450  }
2453 
2454  // set partition to the element
2455  {
2456  if (part_from_moab) {
2457  // if partition is taken from moab partition
2458  int proc = (*efit)->getPartProc();
2459  if (proc == -1 && (*efit)->getEntType() == MBVERTEX) {
2460  proc = (*efit)->getOwnerProc();
2461  }
2462  NumeredEntFiniteElement_change_part(proc).operator()(numered_fe);
2463  } else {
2464  // count partition of the dofs in row, the larges dofs with given
2465  // partition is used to set partition of the element
2466  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows), rows_view,
2467  moab::Interface::UNION);
2468  std::vector<int> parts(m_field.get_comm_size(), 0);
2469  NumeredDofEntity_multiIndex_uid_view_ordered::iterator viit_rows;
2470  viit_rows = rows_view.begin();
2471  for (; viit_rows != rows_view.end(); viit_rows++) {
2472  parts[(*viit_rows)->pArt]++;
2473  }
2474  std::vector<int>::iterator pos =
2475  max_element(parts.begin(), parts.end());
2476  unsigned int max_part = std::distance(parts.begin(), pos);
2477  NumeredEntFiniteElement_change_part(max_part).operator()(numered_fe);
2478  }
2479  }
2480 
2481  // set dofs on rows and columns (if are different)
2482  if ((numered_fe->getPart() >= (unsigned int)low_proc) &&
2483  (numered_fe->getPart() <= (unsigned int)hi_proc)) {
2484 
2485  NumeredDofEntity_multiIndex_uid_view_ordered *dofs_view[] = {&rows_view,
2486  &cols_view};
2487  FENumeredDofEntity_multiIndex *fe_dofs[] = {rows_dofs.get(),
2488  cols_dofs.get()};
2489 
2490  for (int ss = 0; ss != (do_cols_fe ? 2 : 1); ss++) {
2491 
2492  if (ss == 0) {
2493  if (part_from_moab) {
2494  // get row_view
2495  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows),
2496  *dofs_view[ss],
2497  moab::Interface::UNION);
2498  }
2499  } else {
2500  // get cols_views
2501  CHKERR(*efit)->getColDofView(*(p_miit->numeredDofsCols),
2502  *dofs_view[ss], moab::Interface::UNION);
2503  }
2504 
2505  // Following reserve memory in sequences, only two allocations are here,
2506  // once for array of objects, next for array of shared pointers
2507 
2508  // reserve memory for field dofs
2509  boost::shared_ptr<std::vector<FENumeredDofEntity>> dofs_array(
2510  new std::vector<FENumeredDofEntity>());
2511 
2512  if (!ss) {
2513  numered_fe->getRowDofsSequence() = dofs_array;
2514  if (!do_cols_fe)
2515  numered_fe->getColDofsSequence() = dofs_array;
2516  } else
2517  numered_fe->getColDofsSequence() = dofs_array;
2518 
2519  auto vit = dofs_view[ss]->begin();
2520  auto hi_vit = dofs_view[ss]->end();
2521 
2522  dofs_array->reserve(std::distance(vit, hi_vit));
2523 
2524  // create elements objects
2525  for (; vit != hi_vit; vit++) {
2526  boost::shared_ptr<SideNumber> side_number_ptr;
2527  side_number_ptr = (*efit)->getSideNumberPtr((*vit)->getEnt());
2528  dofs_array->emplace_back(side_number_ptr, *vit);
2529  }
2530 
2531  // finally add DoFS to multi-indices
2532  auto hint = fe_dofs[ss]->end();
2533  for (auto &v : *dofs_array)
2534  hint = fe_dofs[ss]->emplace_hint(hint, dofs_array, &v);
2535  }
2536  }
2537  if (!numered_fe->sPtr->row_field_ents_view->empty() &&
2538  !numered_fe->sPtr->col_field_ents_view->empty()) {
2539 
2540  // Add element to the problem
2541  auto p = problem_finite_elements.insert(numered_fe);
2542  if (!p.second)
2543  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2544 
2545  if (verb >= VERY_VERBOSE) {
2546  std::ostringstream ss;
2547  ss << *p_miit << std::endl;
2548  ss << *p.first << std::endl;
2550  FENumeredDofEntityByUId::iterator miit =
2551  (*p.first)->rows_dofs->get<Unique_mi_tag>().begin();
2552  for (; miit != (*p.first)->rows_dofs->get<Unique_mi_tag>().end();
2553  miit++)
2554  ss << "rows: " << *(*miit) << std::endl;
2555  miit = (*p.first)->cols_dofs->get<Unique_mi_tag>().begin();
2556  for (; miit != (*p.first)->cols_dofs->get<Unique_mi_tag>().end();
2557  miit++)
2558  ss << "cols: " << *(*miit) << std::endl;
2559  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2560  }
2561  }
2562  }
2563 
2564  if (verb >= VERBOSE) {
2565  typedef NumeredEntFiniteElement_multiIndex::index<Part_mi_tag>::type
2566  NumeredEntFiniteElementPart;
2567  NumeredEntFiniteElementPart::iterator miit, hi_miit;
2568  miit = problem_finite_elements.get<Part_mi_tag>().lower_bound(
2569  m_field.get_comm_rank());
2570  hi_miit = problem_finite_elements.get<Part_mi_tag>().upper_bound(
2571  m_field.get_comm_rank());
2572  int count = std::distance(miit, hi_miit);
2573  std::ostringstream ss;
2574  ss << *p_miit;
2575  ss << " Nb. elems " << count << " on proc " << m_field.get_comm_rank()
2576  << std::endl;
2577  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2578  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2579  }
2580 
2583 }
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:476
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:595
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:406
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, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, and remove_entities_from_problem.cpp.

Definition at line 2585 of file ProblemsManager.cpp.

2586  {
2587  MoFEM::Interface &m_field = cOre;
2588  const Problem_multiIndex *problems_ptr;
2590 
2592  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2593  "partition of problem not build");
2595  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2596  "partitions finite elements not build");
2597 
2598  // find p_miit
2599  CHKERR m_field.get_problems(&problems_ptr);
2600 
2601  // get problem pointer
2602  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2603  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2604  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2605  name.c_str());
2606 
2607  // get reference to number of ghost dofs
2608  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2609  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2610  nb_row_ghost_dofs = 0;
2611  nb_col_ghost_dofs = 0;
2612 
2613  // do work if more than one processor
2614  if (m_field.get_comm_size() > 1) {
2615 
2617  ghost_idx_col_view;
2618 
2619  // get elements on this partition
2620  auto fe_range =
2621  p_miit->numeredFiniteElements->get<Part_mi_tag>().equal_range(
2622  m_field.get_comm_rank());
2623 
2624  // get dofs on elements which are not part of this partition
2625 
2626  // rows
2627  auto hint_r = ghost_idx_row_view.begin();
2628  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2629  for (auto &dof_ptr : *(*fe_ptr)->rows_dofs) {
2630  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2631  hint_r = ghost_idx_row_view.emplace_hint(
2632  hint_r, dof_ptr->getNumeredDofEntityPtr());
2633  }
2634  }
2635  }
2636 
2637  // columns
2638  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2639  auto hint_c = ghost_idx_col_view.begin();
2640  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2641  for (auto &dof_ptr : *(*fe_range.first)->cols_dofs) {
2642  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2643  hint_c = ghost_idx_col_view.emplace_hint(
2644  hint_c, dof_ptr->getNumeredDofEntityPtr());
2645  }
2646  }
2647  }
2648  }
2649 
2650  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2651  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2652 
2653  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2654  &ghost_idx_row_view, &ghost_idx_col_view};
2655  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2656  &p_miit->numeredDofsRows->get<Unique_mi_tag>(),
2657  &p_miit->numeredDofsCols->get<Unique_mi_tag>()};
2658 
2659  int loop_size = 2;
2660  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2661  loop_size = 1;
2662  }
2663 
2664  // set local ghost dofs indices
2665  for (int ss = 0; ss != loop_size; ++ss) {
2666  for (auto &gid : *ghost_idx_view[ss]) {
2667  NumeredDofEntityByUId::iterator dof =
2668  dof_by_uid_no_const[ss]->find(gid->getGlobalUniqueId());
2669  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2670  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2671  "inconsistent data, ghost dof already set");
2672  bool success = dof_by_uid_no_const[ss]->modify(
2673  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2674  if (PetscUnlikely(!success))
2675  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2676  "modification unsuccessful");
2677  (*nb_ghost_dofs[ss])++;
2678  }
2679  }
2680  if (loop_size == 1) {
2681  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2682  }
2683  }
2684 
2685  if (verb > QUIET) {
2686  std::ostringstream ss;
2687  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2688  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2689  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2690  << p_miit->getNbLocalDofsCol() << std::endl;
2691  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2692  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2693  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2694  << p_miit->getNbLocalDofsRow() << std::endl;
2695  if (verb > VERBOSE) {
2696  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2697  p_miit->numeredDofsCols->begin();
2698  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2699  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2700  continue;
2701  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2702  continue;
2703  ss << *(*miit_dd_col) << std::endl;
2704  }
2705  }
2706  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2707  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2708  }
2709 
2712 }
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:476
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:595
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:406
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 2715 of file ProblemsManager.cpp.

2716  {
2717  MoFEM::Interface &m_field = cOre;
2718  const Problem_multiIndex *problems_ptr;
2720 
2722  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2723  "partition of problem not build");
2725  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2726  "partitions finite elements not build");
2727 
2728  // get problem pointer
2729  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2730  // find p_miit
2731  CHKERR m_field.get_problems(&problems_ptr);
2732  ProblemsByName &problems_set =
2733  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2734  ProblemsByName::iterator p_miit = problems_set.find(name);
2735 
2736  // get reference to number of ghost dofs
2737  // get number of local dofs
2738  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2739  &(p_miit->nbGhostDofsCol)};
2740  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2741  for (int ss = 0; ss != 2; ++ss) {
2742  (*nb_ghost_dofs[ss]) = 0;
2743  }
2744 
2745  // do work if more than one processor
2746  if (m_field.get_comm_size() > 1) {
2747  // determine if rows on columns are different from dofs on rows
2748  int loop_size = 2;
2749  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2750  loop_size = 1;
2751  }
2752 
2753  typedef decltype(p_miit->numeredDofsRows) NumbDofTypeSharedPtr;
2754  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredDofsRows,
2755  p_miit->numeredDofsCols};
2756 
2757  // interate over dofs on rows and dofs on columns
2758  for (int ss = 0; ss != loop_size; ++ss) {
2759 
2760  // create dofs view by uid
2761  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2762 
2763  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2764  ghost_idx_view.reserve(std::distance(r.first, r.second));
2765  for (; r.first != r.second; ++r.first)
2766  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2767 
2768  auto cmp = [](auto a, auto b) {
2769  return (*a)->getGlobalUniqueId() < (*b)->getGlobalUniqueId();
2770  };
2771  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2772 
2773  // intare over dofs which have negative local index
2774  for (auto gid_it : ghost_idx_view) {
2775  bool success = numered_dofs[ss]->modify(
2776  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2777  if (!success)
2778  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2779  "modification unsuccessful");
2780  ++(*nb_ghost_dofs[ss]);
2781  }
2782  }
2783  if (loop_size == 1) {
2784  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2785  }
2786  }
2787 
2788  if (verb > QUIET) {
2789  std::ostringstream ss;
2790  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2791  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2792  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2793  << p_miit->getNbLocalDofsRow() << std::endl;
2794  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2795  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2796  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2797  << p_miit->getNbLocalDofsCol() << std::endl;
2798  if (verb > VERBOSE) {
2799  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2800  p_miit->numeredDofsCols->begin();
2801  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2802  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2803  continue;
2804  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2805  continue;
2806  ss << *(*miit_dd_col) << std::endl;
2807  }
2808  }
2809  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2810  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2811  }
2812 
2815 }
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:476
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:595
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:406
virtual MPI_Comm & get_comm() const =0

◆ partitionMesh()

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

Set partition tag to each finite element in the problemThis will use one of the mesh partitioning programs available from PETSc See http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPartitioningType.html

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

Definition at line 97 of file ProblemsManager.cpp.

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

Definition at line 1888 of file ProblemsManager.cpp.

1889  {
1890  MoFEM::Interface &m_field = cOre;
1891  const Problem_multiIndex *problems_ptr;
1893 
1894  if (!(cOre.getBuildMoFEM() & (1 << 0)))
1895  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1896  if (!(cOre.getBuildMoFEM() & (1 << 1)))
1897  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1898  if (!(cOre.getBuildMoFEM() & (1 << 2)))
1899  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1900  "entFEAdjacencies not build");
1901  if (!(cOre.getBuildMoFEM() & (1 << 3)))
1902  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1903 
1904  if (verb > QUIET)
1905  PetscPrintf(m_field.get_comm(), "Partition problem %s\n", name.c_str());
1906 
1907  typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
1908  NumeredDofEntitysByIdx;
1909  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
1910 
1911  // Find problem pointer by name
1912  CHKERR m_field.get_problems(&problems_ptr);
1913  auto &problems_set =
1914  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1915  auto p_miit = problems_set.find(name);
1916  if (p_miit == problems_set.end()) {
1917  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1918  "problem with name %s not defined (top tip check spelling)",
1919  name.c_str());
1920  }
1921  int nb_dofs_row = p_miit->getNbDofsRow();
1922 
1923  Mat Adj;
1924  CHKERR m_field.getInterface<MatrixManager>()
1925  ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1926 
1927  int m, n;
1928  CHKERR MatGetSize(Adj, &m, &n);
1929  if (verb > VERY_VERBOSE)
1930  MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1931 
1932  // partitioning
1933  MatPartitioning part;
1934  IS is;
1935  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
1936  CHKERR MatPartitioningSetAdjacency(part, Adj);
1937  CHKERR MatPartitioningSetFromOptions(part);
1938  CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
1939  CHKERR MatPartitioningApply(part, &is);
1940  if (verb > VERY_VERBOSE)
1941  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
1942 
1943  // gather
1944  IS is_gather, is_num, is_gather_num;
1945  CHKERR ISAllGather(is, &is_gather);
1946  CHKERR ISPartitioningToNumbering(is, &is_num);
1947  CHKERR ISAllGather(is_num, &is_gather_num);
1948  const int *part_number, *petsc_idx;
1949  CHKERR ISGetIndices(is_gather, &part_number);
1950  CHKERR ISGetIndices(is_gather_num, &petsc_idx);
1951  int size_is_num, size_is_gather;
1952  CHKERR ISGetSize(is_gather, &size_is_gather);
1953  if (size_is_gather != (int)nb_dofs_row)
1954  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
1955  size_is_gather, nb_dofs_row);
1956 
1957  CHKERR ISGetSize(is_num, &size_is_num);
1958  if (size_is_num != (int)nb_dofs_row)
1959  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
1960  size_is_num, nb_dofs_row);
1961 
1962  bool square_matrix = false;
1963  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols)
1964  square_matrix = true;
1965 
1966  if (!square_matrix) {
1967  // FIXME: This is for back compatibility, if deprecate interface function
1968  // build interfaces is removed, this part of the code will be obsolete
1969  auto mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().begin();
1970  auto hi_mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().end();
1971  auto mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().begin();
1972  auto hi_mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().end();
1973  for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
1974  if (mit_col == hi_mit_col) {
1975  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1976  "check finite element definition, nb. of rows is not equal to "
1977  "number for columns");
1978  }
1979  if (mit_row->get()->getGlobalUniqueId() !=
1980  mit_col->get()->getGlobalUniqueId()) {
1981  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1982  "check finite element definition, nb. of rows is not equal to "
1983  "number for columns");
1984  }
1985  }
1986  }
1987 
1988  if (verb > VERBOSE)
1989  PetscPrintf(m_field.get_comm(), "\tloop problem dofs");
1990 
1991  // Set petsc global indices
1992  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1993  p_miit->numeredDofsRows->get<Idx_mi_tag>());
1994  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1995  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1996  nb_row_local_dofs = 0;
1997  nb_row_ghost_dofs = 0;
1998 
1999  for (auto miit_dofs_row = dofs_row_by_idx_no_const.begin();
2000  miit_dofs_row != dofs_row_by_idx_no_const.end(); miit_dofs_row++) {
2001  const int part = part_number[(*miit_dofs_row)->dofIdx];
2002  if (part == (unsigned int)m_field.get_comm_rank()) {
2003  const bool success = dofs_row_by_idx_no_const.modify(
2004  miit_dofs_row,
2005  NumeredDofEntity_part_and_indices_change(
2006  part, petsc_idx[(*miit_dofs_row)->dofIdx], nb_row_local_dofs++));
2007  if (!success) {
2008  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2009  "modification unsuccessful");
2010  }
2011  } else {
2012  const bool success = dofs_row_by_idx_no_const.modify(
2013  miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
2014  part, petsc_idx[(*miit_dofs_row)->dofIdx]));
2015  if (!success) {
2016  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2017  "modification unsuccessful");
2018  }
2019  }
2020  }
2021 
2022  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2023  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2024  if (square_matrix) {
2025  nb_col_local_dofs = nb_row_local_dofs;
2026  nb_col_ghost_dofs = nb_row_ghost_dofs;
2027  } else {
2028  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2029  const_cast<NumeredDofEntitysByIdx &>(
2030  p_miit->numeredDofsCols->get<Idx_mi_tag>());
2031  nb_col_local_dofs = 0;
2032  nb_col_ghost_dofs = 0;
2033  for (auto miit_dofs_col = dofs_col_by_idx_no_const.begin();
2034  miit_dofs_col != dofs_col_by_idx_no_const.end(); miit_dofs_col++) {
2035  const int part = part_number[(*miit_dofs_col)->dofIdx];
2036  if (part == (unsigned int)m_field.get_comm_rank()) {
2037  const bool success = dofs_col_by_idx_no_const.modify(
2038  miit_dofs_col, NumeredDofEntity_part_and_indices_change(
2039  part, petsc_idx[(*miit_dofs_col)->dofIdx],
2040  nb_col_local_dofs++));
2041  if (!success) {
2042  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2043  "modification unsuccessful");
2044  }
2045  } else {
2046  const bool success = dofs_col_by_idx_no_const.modify(
2047  miit_dofs_col, NumeredDofEntity_part_and_glob_idx_change(
2048  part, petsc_idx[(*miit_dofs_col)->dofIdx]));
2049  if (!success) {
2050  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2051  "modification unsuccessful");
2052  }
2053  }
2054  }
2055  }
2056 
2057  if (verb > VERBOSE)
2058  PetscPrintf(m_field.get_comm(), " <- done\n");
2059 
2060  CHKERR ISRestoreIndices(is_gather, &part_number);
2061  CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
2062  CHKERR ISDestroy(&is_num);
2063  CHKERR ISDestroy(&is_gather_num);
2064  CHKERR ISDestroy(&is_gather);
2065  CHKERR ISDestroy(&is);
2066  CHKERR MatPartitioningDestroy(&part);
2067  CHKERR MatDestroy(&Adj);
2068  CHKERR printPartitionedProblem(&*p_miit, verb);
2069 
2070  cOre.getBuildMoFEM() |= 1 << 4;
2072 }
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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
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, and forces_and_sources_testing_flat_prism_element.cpp.

Definition at line 1740 of file ProblemsManager.cpp.

1741  {
1742 
1743  MoFEM::Interface &m_field = cOre;
1744  const Problem_multiIndex *problems_ptr;
1747  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1748  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1749  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1750  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1751  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1753  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1754  if (verb > 0) {
1755  PetscPrintf(m_field.get_comm(), "Simple partition problem %s\n",
1756  name.c_str());
1757  }
1758  CHKERR m_field.get_problems(&problems_ptr);
1759  // find p_miit
1760  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1761  ProblemByName &problems_set =
1762  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1763  ProblemByName::iterator p_miit = problems_set.find(name);
1764  if (p_miit == problems_set.end()) {
1765  SETERRQ1(PETSC_COMM_SELF, 1,
1766  "problem < %s > is not found (top tip: check spelling)",
1767  name.c_str());
1768  }
1769  typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1770  Idx_mi_tag>::type NumeredDofEntitysByIdx;
1771  NumeredDofEntitysByIdx &dofs_row_by_idx =
1772  p_miit->numeredDofsRows->get<Idx_mi_tag>();
1773  NumeredDofEntitysByIdx &dofs_col_by_idx =
1774  p_miit->numeredDofsCols->get<Idx_mi_tag>();
1775  boost::multi_index::index<NumeredDofEntity_multiIndex,
1776  Idx_mi_tag>::type::iterator miit_row,
1777  hi_miit_row;
1778  boost::multi_index::index<NumeredDofEntity_multiIndex,
1779  Idx_mi_tag>::type::iterator miit_col,
1780  hi_miit_col;
1781  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1782  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1783  nb_row_local_dofs = 0;
1784  nb_row_ghost_dofs = 0;
1785  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1786  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1787  nb_col_local_dofs = 0;
1788  nb_col_ghost_dofs = 0;
1789 
1790  bool square_matrix = false;
1791  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
1792  square_matrix = true;
1793  }
1794 
1795  // get row range of local indices
1796  PetscLayout layout_row;
1797  const int *ranges_row;
1798 
1799  DofIdx nb_dofs_row = dofs_row_by_idx.size();
1800  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1801  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1802  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1803  CHKERR PetscLayoutSetUp(layout_row);
1804  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1805  // get col range of local indices
1806  PetscLayout layout_col;
1807  const int *ranges_col;
1808  if (!square_matrix) {
1809  DofIdx nb_dofs_col = dofs_col_by_idx.size();
1810  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1811  CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1812  CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1813  CHKERR PetscLayoutSetUp(layout_col);
1814  CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1815  }
1816  for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1817  part++) {
1818  miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1819  hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1820  if (std::distance(miit_row, hi_miit_row) !=
1821  ranges_row[part + 1] - ranges_row[part]) {
1822  SETERRQ4(
1823  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1824  "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1825  "rstart (%d != %d - %d = %d) ",
1826  std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1827  ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1828  }
1829  // loop rows
1830  for (; miit_row != hi_miit_row; miit_row++) {
1831  bool success = dofs_row_by_idx.modify(
1832  miit_row,
1833  NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1834  if (!success)
1835  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1836  "modification unsuccessful");
1837  if (part == (unsigned int)m_field.get_comm_rank()) {
1838  success = dofs_row_by_idx.modify(
1839  miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1840  if (!success)
1841  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1842  "modification unsuccessful");
1843  }
1844  }
1845  if (!square_matrix) {
1846  miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1847  hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1848  if (std::distance(miit_col, hi_miit_col) !=
1849  ranges_col[part + 1] - ranges_col[part]) {
1850  SETERRQ4(
1851  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1852  "data inconsistency, std::distance(miit_col,hi_miit_col) != rend - "
1853  "rstart (%d != %d - %d = %d) ",
1854  std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1855  ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1856  }
1857  // loop cols
1858  for (; miit_col != hi_miit_col; miit_col++) {
1859  bool success = dofs_col_by_idx.modify(
1860  miit_col, NumeredDofEntity_part_and_glob_idx_change(
1861  part, (*miit_col)->dofIdx));
1862  if (!success)
1863  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1864  "modification unsuccessful");
1865  if (part == (unsigned int)m_field.get_comm_rank()) {
1866  success = dofs_col_by_idx.modify(
1867  miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1868  if (!success)
1869  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1870  "modification unsuccessful");
1871  }
1872  }
1873  }
1874  }
1875  CHKERR PetscLayoutDestroy(&layout_row);
1876  if (!square_matrix) {
1877  CHKERR PetscLayoutDestroy(&layout_col);
1878  }
1879  if (square_matrix) {
1880  nb_col_local_dofs = nb_row_local_dofs;
1881  nb_col_ghost_dofs = nb_row_ghost_dofs;
1882  }
1883  CHKERR printPartitionedProblem(&*p_miit, verb);
1886 }
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:476
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:595
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:406
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 2853 of file ProblemsManager.cpp.

2856  {
2857 
2858  MoFEM::Interface &m_field = cOre;
2860 
2861  const Problem *prb_ptr;
2862  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2863 
2864  decltype(prb_ptr->numeredDofsRows) numered_dofs[2] = {
2865  prb_ptr->numeredDofsRows, nullptr};
2866  if (prb_ptr->numeredDofsRows != prb_ptr->numeredDofsCols)
2867  numered_dofs[1] = prb_ptr->numeredDofsCols;
2868 
2869  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2870  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2871  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2872 
2873  for (int s = 0; s != 2; ++s)
2874  if (numered_dofs[s]) {
2875 
2876  typedef multi_index_container<
2877 
2878  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2879 
2880  >
2881  NumeredDofEntity_it_view_multiIndex;
2882 
2883  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2884 
2885  // Set -1 to global and local dofs indices
2886  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2887  ++pit) {
2888  auto lo =
2889  numered_dofs[s]
2890  ->get<Composite_Name_And_Ent_And_EntDofIdx_mi_tag>()
2891  .lower_bound(boost::make_tuple(field_name, pit->first, 0));
2892  auto hi = numered_dofs[s]
2893  ->get<Composite_Name_And_Ent_And_EntDofIdx_mi_tag>()
2894  .lower_bound(boost::make_tuple(field_name, pit->second,
2896  for (; lo != hi; ++lo)
2897  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2898  (*lo)->getDofCoeffIdx() <= hi_coeff)
2899  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2900  }
2901 
2902  if (verb >= VERY_NOISY) {
2903  for (auto &dof : dofs_it_view) {
2904  std::ostringstream ss;
2905  ss << **dof;
2906  PetscSynchronizedPrintf(m_field.get_comm(), "%s\n", ss.str().c_str());
2907  }
2908  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2909  }
2910 
2911  // set negative index
2912  auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2913  for (auto dit : dofs_it_view) {
2914  bool success = numered_dofs[s]->modify(dit, mod);
2915  if (!success)
2916  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2917  "can not set negative indices");
2918  }
2919 
2920  // create weak view
2921  std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2922  dosf_weak_view.reserve(dofs_it_view.size());
2923  for (auto dit : dofs_it_view)
2924  dosf_weak_view.push_back(*dit);
2925 
2926  if (verb >= NOISY)
2927  PetscSynchronizedPrintf(
2928  m_field.get_comm(),
2929  "Number of DOFs in multi-index %d and to delete %d\n",
2930  numered_dofs[s]->size(), dofs_it_view.size());
2931 
2932  // erase dofs from problem
2933  for (auto weak_dit : dosf_weak_view)
2934  if (auto dit = weak_dit.lock()) {
2935  numered_dofs[s]->erase(dit->getGlobalUniqueId());
2936  }
2937 
2938  if (verb >= NOISY)
2939  PetscSynchronizedPrintf(
2940  m_field.get_comm(),
2941  "Number of DOFs in multi-index after delete %d\n",
2942  numered_dofs[s]->size());
2943 
2944  // get current number of ghost dofs
2945  int nb_local_dofs = 0;
2946  int nb_ghost_dofs = 0;
2947  for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2948  dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2949  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2950  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2951  ++nb_local_dofs;
2952  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2953  ++nb_ghost_dofs;
2954  else
2955  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Imposible case");
2956  }
2957 
2958  // get indices
2959  auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2960  const int nb_dofs = numered_dofs[s]->size();
2961  indices.clear();
2962  indices.reserve(nb_dofs);
2963  for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2964  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2965  bool add = true;
2966  if (only_local) {
2967  if ((*dit)->getPetscLocalDofIdx() < 0 ||
2968  (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2969  add = false;
2970  }
2971  }
2972  if (add)
2973  indices.push_back(decltype(tag)::get_index(dit));
2974  }
2975  };
2976 
2977  auto get_indices_by_uid = [&](auto tag, auto &indices) {
2978  const int nb_dofs = numered_dofs[s]->size();
2979  indices.clear();
2980  indices.reserve(nb_dofs);
2981  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2982  ++dit) {
2983  const int idx = decltype(tag)::get_index(dit);
2984  indices.push_back(decltype(tag)::get_index(dit));
2985  }
2986  };
2987 
2988  auto concatenate_dofs = [&](auto tag, auto &indices,
2989  const auto local_only) {
2991  get_indices_by_tag(tag, indices, local_only);
2992 
2993  AO ao;
2994  if (local_only)
2995  CHKERR AOCreateMapping(m_field.get_comm(), indices.size(),
2996  &*indices.begin(), PETSC_NULL, &ao);
2997  else
2998  CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
2999  &*indices.begin(), PETSC_NULL, &ao);
3000 
3001  get_indices_by_uid(tag, indices);
3002  CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3003  CHKERR AODestroy(&ao);
3005  };
3006 
3007  // set indices index
3008  auto set_concatinated_indices = [&]() {
3009  std::vector<int> global_indices;
3010  std::vector<int> local_indices;
3012  CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3013  CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3014  auto gi = global_indices.begin();
3015  auto li = local_indices.begin();
3016  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3017  ++dit) {
3018  auto mod = NumeredDofEntity_part_and_all_indices_change(
3019  (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3020  bool success = numered_dofs[s]->modify(dit, mod);
3021  if (!success)
3022  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3023  "can not set negative indices");
3024  ++gi;
3025  ++li;
3026  }
3028  };
3029  CHKERR set_concatinated_indices();
3030 
3031  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3032  m_field.get_comm());
3033  *(local_nbdof_ptr[s]) = nb_local_dofs;
3034  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3035 
3036  if (debug)
3037  for (auto dof : (*numered_dofs[s])) {
3038  if (dof->getPetscGlobalDofIdx() < 0) {
3039  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3040  "Negative global idx");
3041  }
3042  if (dof->getPetscLocalDofIdx() < 0) {
3043  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3044  "Negative local idx");
3045  }
3046  }
3047 
3048  } else {
3049 
3050  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3051  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3052  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3053  }
3054 
3055  if (verb > QUIET) {
3056  PetscSynchronizedPrintf(
3057  m_field.get_comm(),
3058  "removed ents on rank %d from problem %s dofs [ %d / %d local, %d / %d "
3059  "ghost, %d / %d global]\n",
3060  m_field.get_comm_rank(), prb_ptr->getName().c_str(),
3061  prb_ptr->getNbLocalDofsRow(), prb_ptr->getNbLocalDofsCol(),
3062  prb_ptr->getNbGhostDofsRow(), prb_ptr->getNbGhostDofsCol(),
3063  prb_ptr->getNbDofsRow(), prb_ptr->getNbDofsCol());
3064  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
3065  }
3066 
3068 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:297
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
virtual MPI_Comm & get_comm() const =0