v0.14.0
ProblemsManager.cpp
Go to the documentation of this file.
1 /** \file ProblemsManager.cpp
2  * \brief Managing complexities for problem
3  * \ingroup mofem_problems_manager
4  */
5 
6 namespace MoFEM {
7 
8 #define ProblemManagerFunctionBegin \
9  MoFEMFunctionBegin; \
10  MOFEM_LOG_CHANNEL("WORLD"); \
11  MOFEM_LOG_CHANNEL("SYNC"); \
12  MOFEM_LOG_FUNCTION(); \
13  MOFEM_LOG_TAG("SYNC", "ProblemsManager"); \
14  MOFEM_LOG_TAG("WORLD", "ProblemsManager")
15 
16 struct IdxDataType {
17  IdxDataType(const UId uid, const int dof) {
18  bcopy(&uid, dAta, 4 * sizeof(int));
19  dAta[4] = dof;
20  }
21 
22 private:
23  int dAta[5];
24 };
25 
27  IdxDataTypePtr(const int *ptr) : pTr(ptr) {}
28  inline int getDofIdx() const {
29  int global_dof = pTr[4];
30  return global_dof;
31  }
32  inline UId getUId() const {
33  unsigned int b0, b1, b2, b3;
34  bcopy(&pTr[0], &b0, sizeof(int));
35  bcopy(&pTr[1], &b1, sizeof(int));
36  bcopy(&pTr[2], &b2, sizeof(int));
37  bcopy(&pTr[3], &b3, sizeof(int));
38  UId uid = static_cast<UId>(b0) | static_cast<UId>(b1) << 8 * sizeof(int) |
39  static_cast<UId>(b2) << 16 * sizeof(int) |
40  static_cast<UId>(b3) << 24 * sizeof(int);
41  return uid;
42  }
43 
44 private:
45  const int *pTr;
46 };
47 
49 ProblemsManager::query_interface(boost::typeindex::type_index type_index,
50  UnknownInterface **iface) const {
51  *iface = const_cast<ProblemsManager *>(this);
52  return 0;
53 }
54 
56  : cOre(const_cast<MoFEM::Core &>(core)),
57  buildProblemFromFields(PETSC_FALSE),
58  synchroniseProblemEntities(PETSC_FALSE) {
59  PetscLogEventRegister("ProblemsManager", 0, &MOFEM_EVENT_ProblemsManager);
60 }
61 
63  MoFEM::Interface &m_field = cOre;
65  CHKERR PetscOptionsBegin(m_field.get_comm(), "", "Problem manager", "none");
66  {
67  CHKERR PetscOptionsBool(
68  "-problem_build_from_fields",
69  "Add DOFs to problem directly from fields not through DOFs on elements",
71  }
72  ierr = PetscOptionsEnd();
73  CHKERRG(ierr);
75 }
76 
78  const Range &ents, const int dim, const int adj_dim, const int n_parts,
79  Tag *th_vertex_weights, Tag *th_edge_weights, Tag *th_part_weights,
80  int verb, const bool debug) {
81  return static_cast<MoFEM::Interface &>(cOre)
82  .getInterface<CommInterface>()
83  ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
84  th_edge_weights, th_part_weights, verb, debug);
85 }
86 
88  const bool square_matrix,
89  int verb) {
90 
91  MoFEM::Interface &m_field = cOre;
93  if (!(cOre.getBuildMoFEM() & (1 << 0)))
94  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
95  if (!(cOre.getBuildMoFEM() & (1 << 1)))
96  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
97  if (!(cOre.getBuildMoFEM() & (1 << 2)))
98  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
99  const Problem *problem_ptr;
100  CHKERR m_field.get_problem(name, &problem_ptr);
101  CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
102  cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
103  // function knows what he is doing
105 }
106 
108  const bool square_matrix,
109  int verb) {
110  MoFEM::Interface &m_field = cOre;
111  auto dofs_field_ptr = m_field.get_dofs();
112  auto ents_field_ptr = m_field.get_field_ents();
113  auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
115  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
116 
117  // Note: Only allowed changes on problem_ptr structure which not influence
118  // multi-index.
119 
120  if (problem_ptr->getBitRefLevel().none()) {
121  SETERRQ1(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
122  problem_ptr->getName().c_str());
123  }
124  CHKERR m_field.clear_problem(problem_ptr->getName());
125 
126  // zero finite elements
127  problem_ptr->numeredFiniteElementsPtr->clear();
128 
129  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
130  auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
132  for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
133 
134  const auto uid = (*it)->getLocalUniqueId();
135 
136  auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
137  for (auto lo = r.first; lo != r.second; ++lo) {
138 
139  if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
140  std::array<bool, 2> row_col = {false, false};
141 
142  const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
143  const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
144  const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
145 
146  // if entity is not problem refinement level
147  if ((fe_bit & prb_mask) != fe_bit)
148  continue;
149  if ((fe_bit & prb_bit).none())
150  continue;
151 
152  auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
153  auto dit = nb_dofs->lower_bound(uid);
154  decltype(dit) hi_dit;
155  if (dit != nb_dofs->end()) {
156  hi_dit = nb_dofs->upper_bound(
157  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
158  view.insert(dit, hi_dit);
159  row_col[rc] = true;
160  }
161  };
162 
163  if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
164  add_to_view(dofs_field_ptr, dofs_rows, ROW);
165 
166  if (!square_matrix)
167  if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
168  add_to_view(dofs_field_ptr, dofs_cols, COL);
169 
170  if (square_matrix && row_col[ROW])
171  break;
172  else if (row_col[ROW] && row_col[COL])
173  break;
174  }
175  }
176  }
178  };
179 
180  CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
181 
182  // Add row dofs to problem
183  {
184  // zero rows
185  problem_ptr->nbDofsRow = 0;
186  problem_ptr->nbLocDofsRow = 0;
187  problem_ptr->nbGhostDofsRow = 0;
188 
189  // add dofs for rows
190  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
191  hi_miit;
192  hi_miit = dofs_rows.get<0>().end();
193 
194  int count_dofs = dofs_rows.get<1>().count(true);
195  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
196  boost::shared_ptr<std::vector<NumeredDofEntity>>(
197  new std::vector<NumeredDofEntity>());
198  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
199  dofs_array->reserve(count_dofs);
200  miit = dofs_rows.get<0>().begin();
201  for (; miit != hi_miit; miit++) {
202  if ((*miit)->getActive()) {
203  dofs_array->emplace_back(*miit);
204  dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
205  }
206  }
207  auto hint = problem_ptr->numeredRowDofsPtr->end();
208  for (auto &v : *dofs_array) {
209  hint = problem_ptr->numeredRowDofsPtr->emplace_hint(hint, dofs_array, &v);
210  }
211  }
212 
213  // Add col dofs to problem
214  if (!square_matrix) {
215  // zero cols
216  problem_ptr->nbDofsCol = 0;
217  problem_ptr->nbLocDofsCol = 0;
218  problem_ptr->nbGhostDofsCol = 0;
219 
220  // add dofs for cols
221  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
222  hi_miit;
223  hi_miit = dofs_cols.get<0>().end();
224 
225  int count_dofs = 0;
226  miit = dofs_cols.get<0>().begin();
227  for (; miit != hi_miit; miit++) {
228  if (!(*miit)->getActive()) {
229  continue;
230  }
231  count_dofs++;
232  }
233 
234  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
235  boost::shared_ptr<std::vector<NumeredDofEntity>>(
236  new std::vector<NumeredDofEntity>());
237  problem_ptr->getColDofsSequence()->push_back(dofs_array);
238  dofs_array->reserve(count_dofs);
239  miit = dofs_cols.get<0>().begin();
240  for (; miit != hi_miit; miit++) {
241  if (!(*miit)->getActive()) {
242  continue;
243  }
244  dofs_array->emplace_back(*miit);
245  dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
246  }
247  auto hint = problem_ptr->numeredColDofsPtr->end();
248  for (auto &v : *dofs_array) {
249  hint = problem_ptr->numeredColDofsPtr->emplace_hint(hint, dofs_array, &v);
250  }
251  } else {
252  problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
253  problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
254  problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
255  }
256 
257  // job done, some debugging and postprocessing
258  if (verb >= VERBOSE) {
259  MOFEM_LOG("SYNC", Sev::verbose)
260  << problem_ptr->getName() << " Nb. local dofs "
261  << problem_ptr->numeredRowDofsPtr->size() << " by "
262  << problem_ptr->numeredColDofsPtr->size();
263  MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
264  }
265 
266  if (verb >= NOISY) {
267  MOFEM_LOG("SYNC", Sev::noisy)
268  << "FEs row dofs " << *problem_ptr << " Nb. row dof "
269  << problem_ptr->getNbDofsRow();
270  for (auto &miit : *problem_ptr->numeredRowDofsPtr)
271  MOFEM_LOG("SYNC", Sev::noisy) << *miit;
272 
273  MOFEM_LOG("SYNC", Sev::noisy)
274  << "FEs col dofs " << *problem_ptr << " Nb. col dof "
275  << problem_ptr->getNbDofsCol();
276  for (auto &miit : *problem_ptr->numeredColDofsPtr)
277  MOFEM_LOG("SYNC", Sev::noisy) << *miit;
278  MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
279  }
280 
281  cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
282  // uses this function knows
283  // what he is doing
284 
285  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
286 
288 }
289 
291  const std::string name, const bool square_matrix, int verb) {
292  MoFEM::Interface &m_field = cOre;
294 
295  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
296  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
297  if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
298  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
299  if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
300  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
301 
302  const Problem *problem_ptr;
303  CHKERR m_field.get_problem(name, &problem_ptr);
304  CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
305  square_matrix, verb);
306 
309 
311 }
312 
314  Problem *problem_ptr, const bool square_matrix, int verb) {
315  MoFEM::Interface &m_field = cOre;
316  auto fields_ptr = m_field.get_fields();
317  auto fe_ptr = m_field.get_finite_elements();
318  auto dofs_field_ptr = m_field.get_dofs();
320  PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
321 
322  // clear data structures
323  CHKERR m_field.clear_problem(problem_ptr->getName());
324 
325  CHKERR getOptions();
326 
327  if (problem_ptr->getBitRefLevel().none())
328  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
329  "problem <%s> refinement level not set",
330  problem_ptr->getName().c_str());
331 
332  int loop_size = 2;
333  if (square_matrix) {
334  loop_size = 1;
335  problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
336  } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
337  problem_ptr->numeredColDofsPtr =
338  boost::shared_ptr<NumeredDofEntity_multiIndex>(
340  }
341 
342  const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
343  const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
344 
345  // // get rows and cols dofs view based on data on elements
346  DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
347 
348  // Add DOFs to problem by visiting all elements and adding DOFs from
349  // elements to the problem
350  if (buildProblemFromFields == PETSC_FALSE) {
351 
352  auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
353  auto ents_field_ptr = m_field.get_field_ents();
354  auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
356  for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
357  ++it) {
358 
359  const auto uid = (*it)->getLocalUniqueId();
360 
361  auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
362  for (auto lo = r.first; lo != r.second; ++lo) {
363 
364  if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
365  std::array<bool, 2> row_col = {false, false};
366 
367  const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
368 
369  // if entity is not problem refinement level
370  if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
371 
372  auto add_to_view = [&](auto dofs, auto &view, auto rc) {
373  auto dit = dofs->lower_bound(uid);
374  decltype(dit) hi_dit;
375  if (dit != dofs->end()) {
376  hi_dit = dofs->upper_bound(
377  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
378  view.insert(dit, hi_dit);
379  row_col[rc] = true;
380  }
381  };
382 
383  if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
384  add_to_view(dofs_field_ptr, dofs_rows, ROW);
385 
386  if (!square_matrix)
387  if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
388  add_to_view(dofs_field_ptr, dofs_cols, COL);
389 
390  if (square_matrix && row_col[ROW])
391  break;
392  else if (row_col[ROW] && row_col[COL])
393  break;
394  }
395  }
396  }
397  }
399  };
400 
401  CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
402  }
403 
404  // Add DOFS to the problem by searching all the filedes, and adding to
405  // problem owned or shared DOFs
406  if (buildProblemFromFields == PETSC_TRUE) {
407  // Get fields IDs on elements
408  BitFieldId fields_ids_row, fields_ids_col;
409  for (auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
410  if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
411  fields_ids_row |= fit->get()->getBitFieldIdRow();
412  fields_ids_col |= fit->get()->getBitFieldIdCol();
413  }
414  }
415  // Get fields DOFs
416  for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
417  if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
418 
419  auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
420  FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
421  auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
422  FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
423 
424  for (; dit != hi_dit; dit++) {
425 
426  const int owner_proc = dit->get()->getOwnerProc();
427  if (owner_proc != m_field.get_comm_rank()) {
428  const unsigned char pstatus = dit->get()->getPStatus();
429  if (pstatus == 0) {
430  continue;
431  }
432  }
433 
434  const auto &dof_bit = (*dit)->getBitRefLevel();
435  // if entity is not problem refinement level
436  if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
437 
438  if ((fit->get()->getId() & fields_ids_row).any()) {
439  dofs_rows.insert(*dit);
440  }
441  if (!square_matrix) {
442  if ((fit->get()->getId() & fields_ids_col).any()) {
443  dofs_cols.insert(*dit);
444  }
445  }
446  }
447  }
448  }
449  }
450  }
451 
453  // Get fields IDs on elements
454  BitFieldId fields_ids_row, fields_ids_col;
455  BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
456  for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
457  fit != fe_ptr->end(); fit++) {
458  if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
459  fields_ids_row |= fit->get()->getBitFieldIdRow();
460  fields_ids_col |= fit->get()->getBitFieldIdCol();
461  }
462  }
463 
464  DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
465  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
466  DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
467  hi_miit;
468  miit = dofs_ptr[ss]->get<1>().lower_bound(1);
469  hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
470  Range ents_to_synchronise;
471  for (; miit != hi_miit; ++miit) {
472  if (miit->get()->getEntDofIdx() != 0)
473  continue;
474  ents_to_synchronise.insert(miit->get()->getEnt());
475  }
476  Range tmp_ents = ents_to_synchronise;
477  CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
478  ents_to_synchronise, nullptr, verb);
479  ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
480  for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
481  if ((fit->get()->getId() & *fields_ids[ss]).any()) {
482  const auto bit_number = (*fit)->getBitNumber();
483  for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
484  pit != ents_to_synchronise.pair_end(); ++pit) {
485  const auto f = pit->first;
486  const auto s = pit->second;
487  const auto lo_uid =
489  const auto hi_uid =
491 
492  auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
493  auto hi_dit =
494  dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
495 
496  dofs_ptr[ss]->insert(dit, hi_dit);
497  }
498  }
499  }
500  }
501  }
502 
503  // add dofs for rows and cols and set ownership
504  DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
505  boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
506  problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
507  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
508  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
509  &problem_ptr->nbLocDofsCol};
510  int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
511  &problem_ptr->nbGhostDofsCol};
512  for (int ss = 0; ss < 2; ss++) {
513  *(nbdof_ptr[ss]) = 0;
514  *(local_nbdof_ptr[ss]) = 0;
515  *(ghost_nbdof_ptr[ss]) = 0;
516  }
517 
518  // Loop over dofs on rows and columns and add to multi-indices in dofs
519  // problem structure, set partition for each dof
520  int nb_local_dofs[] = {0, 0};
521  for (int ss = 0; ss < loop_size; ss++) {
522  DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
523  hi_miit;
524  miit = dofs_ptr[ss]->get<1>().lower_bound(1);
525  hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
526  for (; miit != hi_miit; miit++) {
527  int owner_proc = (*miit)->getOwnerProc();
528  if (owner_proc == m_field.get_comm_rank()) {
529  nb_local_dofs[ss]++;
530  }
531  }
532  }
533 
534  // get layout
535  int start_ranges[2], end_ranges[2];
536  for (int ss = 0; ss != loop_size; ss++) {
537  PetscLayout layout;
538  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
539  CHKERR PetscLayoutSetBlockSize(layout, 1);
540  CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
541  CHKERR PetscLayoutSetUp(layout);
542  CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
543  CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
544  CHKERR PetscLayoutDestroy(&layout);
545  }
546  if (square_matrix) {
547  nbdof_ptr[1] = nbdof_ptr[0];
548  nb_local_dofs[1] = nb_local_dofs[0];
549  start_ranges[1] = start_ranges[0];
550  end_ranges[1] = end_ranges[0];
551  }
552 
553  // if(sizeof(UId) != SIZEOFUID) {
554  // SETERRQ2(
555  // PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,
556  // "check size of UId, size of UId is %u != %u",
557  // sizeof(UId),SIZEOFUID
558  // );
559  // }
560 
561  // set local and global indices on own dofs
562  const size_t idx_data_type_size = sizeof(IdxDataType);
563  const size_t data_block_size = idx_data_type_size / sizeof(int);
564 
565  if (sizeof(IdxDataType) % sizeof(int)) {
566  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
567  }
568  std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
569  m_field.get_comm_size()),
570  ids_data_packed_cols(m_field.get_comm_size());
571 
572  // Loop over dofs on this processor and prepare those dofs to send on
573  // another proc
574  for (int ss = 0; ss != loop_size; ++ss) {
575 
576  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
577  hi_miit;
578  hi_miit = dofs_ptr[ss]->get<0>().end();
579 
580  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
581  boost::shared_ptr<std::vector<NumeredDofEntity>>(
582  new std::vector<NumeredDofEntity>());
583  int nb_dofs_to_add = 0;
584  miit = dofs_ptr[ss]->get<0>().begin();
585  for (; miit != hi_miit; ++miit) {
586  // Only set global idx for dofs on this processor part
587  if (!(miit->get()->getActive()))
588  continue;
589  ++nb_dofs_to_add;
590  }
591  dofs_array->reserve(nb_dofs_to_add);
592  if (ss == 0) {
593  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
594  } else {
595  problem_ptr->getColDofsSequence()->push_back(dofs_array);
596  }
597 
598  int &local_idx = *local_nbdof_ptr[ss];
599  miit = dofs_ptr[ss]->get<0>().begin();
600  for (; miit != hi_miit; ++miit) {
601 
602  // Only set global idx for dofs on this processor part
603  if (!(miit->get()->getActive()))
604  continue;
605 
606  dofs_array->emplace_back(*miit);
607 
608  int owner_proc = dofs_array->back().getOwnerProc();
609  if (owner_proc < 0) {
610  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
611  "data inconsistency");
612  }
613 
614  if (owner_proc != m_field.get_comm_rank()) {
615  dofs_array->back().pArt = owner_proc;
616  dofs_array->back().dofIdx = -1;
617  dofs_array->back().petscGloablDofIdx = -1;
618  dofs_array->back().petscLocalDofIdx = -1;
619  } else {
620 
621  // Set part and indexes
622  int glob_idx = start_ranges[ss] + local_idx;
623  dofs_array->back().pArt = owner_proc;
624  dofs_array->back().dofIdx = glob_idx;
625  dofs_array->back().petscGloablDofIdx = glob_idx;
626  dofs_array->back().petscLocalDofIdx = local_idx;
627  local_idx++;
628 
629  unsigned char pstatus = dofs_array->back().getPStatus();
630  // check id dof is shared, if that is a case global idx added to data
631  // structure and is sended to other processors
632  if (pstatus > 0) {
633 
634  for (int proc = 0;
635  proc < MAX_SHARING_PROCS &&
636  -1 != dofs_array->back().getSharingProcsPtr()[proc];
637  proc++) {
638  // make it different for rows and columns when needed
639  if (ss == 0) {
640  ids_data_packed_rows[dofs_array->back()
641  .getSharingProcsPtr()[proc]]
642  .emplace_back(dofs_array->back().getGlobalUniqueId(),
643  glob_idx);
644  } else {
645  ids_data_packed_cols[dofs_array->back()
646  .getSharingProcsPtr()[proc]]
647  .emplace_back(dofs_array->back().getGlobalUniqueId(),
648  glob_idx);
649  }
650  if (!(pstatus & PSTATUS_MULTISHARED)) {
651  break;
652  }
653  }
654  }
655  }
656  }
657 
658  auto hint = numered_dofs_ptr[ss]->end();
659  for (auto &v : *dofs_array)
660  hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
661  }
662  if (square_matrix) {
663  local_nbdof_ptr[1] = local_nbdof_ptr[0];
664  }
665 
666  int nsends_rows = 0, nsends_cols = 0;
667  // Non zero lengths[i] represent a message to i of length lengths[i].
668  std::vector<int> lengths_rows(m_field.get_comm_size()),
669  lengths_cols(m_field.get_comm_size());
670  lengths_rows.clear();
671  lengths_cols.clear();
672  for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
673  lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
674  lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
675  if (!ids_data_packed_rows[proc].empty())
676  nsends_rows++;
677  if (!ids_data_packed_cols[proc].empty())
678  nsends_cols++;
679  }
680 
681  MPI_Status *status;
682  CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
683 
684  // Do rows
685  int nrecvs_rows; // number of messages received
686  int *onodes_rows; // list of node-ids from which messages are expected
687  int *olengths_rows; // corresponding message lengths
688  int **rbuf_row; // must bee freed by user
689 
690  // make sure it is a PETSc comm
691  MPI_Comm comm;
692  CHKERR PetscCommDuplicate(m_field.get_comm(), &comm, NULL);
693 
694  {
695 
696  // rows
697 
698  // Computes the number of messages a node expects to receive
699  CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
700  &nrecvs_rows);
701  // std::cerr << nrecvs_rows << std::endl;
702 
703  // Computes info about messages that a MPI-node will receive, including
704  // (from-id,length) pairs for each message.
705  CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
706  &lengths_rows[0], &onodes_rows,
707  &olengths_rows);
708 
709  // Gets a unique new tag from a PETSc communicator. All processors that
710  // share the communicator MUST call this routine EXACTLY the same number
711  // of times. This tag should only be used with the current objects
712  // communicator; do NOT use it with any other MPI communicator.
713  int tag_row;
714  CHKERR PetscCommGetNewTag(comm, &tag_row);
715 
716  // Allocate a buffer sufficient to hold messages of size specified in
717  // olengths. And post Irecvs on these buffers using node info from onodes
718  MPI_Request *r_waits_row; // must bee freed by user
719  // rbuf has a pointers to messeges. It has size of of nrecvs (number of
720  // messages) +1. In the first index a block is allocated,
721  // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
722 
723  CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
724  olengths_rows, &rbuf_row, &r_waits_row);
725  CHKERR PetscFree(onodes_rows);
726 
727  MPI_Request *s_waits_row; // status of sens messages
728  CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
729 
730  // Send messeges
731  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
732  if (!lengths_rows[proc])
733  continue; // no message to send to this proc
734  CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
735  lengths_rows[proc], // message length
736  MPIU_INT, proc, // to proc
737  tag_row, comm, s_waits_row + kk);
738  kk++;
739  }
740 
741  if (nrecvs_rows) {
742  CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
743  }
744  if (nsends_rows) {
745  CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
746  }
747 
748  CHKERR PetscFree(r_waits_row);
749  CHKERR PetscFree(s_waits_row);
750  }
751 
752  // cols
753  int nrecvs_cols = nrecvs_rows;
754  int *olengths_cols = olengths_rows;
755  PetscInt **rbuf_col = rbuf_row;
756  if (!square_matrix) {
757 
758  // Computes the number of messages a node expects to receive
759  CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
760  &nrecvs_cols);
761 
762  // Computes info about messages that a MPI-node will receive, including
763  // (from-id,length) pairs for each message.
764  int *onodes_cols;
765  CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
766  &lengths_cols[0], &onodes_cols,
767  &olengths_cols);
768 
769  // Gets a unique new tag from a PETSc communicator.
770  int tag_col;
771  CHKERR PetscCommGetNewTag(comm, &tag_col);
772 
773  // Allocate a buffer sufficient to hold messages of size specified in
774  // olengths. And post Irecvs on these buffers using node info from onodes
775  MPI_Request *r_waits_col; // must bee freed by user
776  CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
777  olengths_cols, &rbuf_col, &r_waits_col);
778  CHKERR PetscFree(onodes_cols);
779 
780  MPI_Request *s_waits_col; // status of sens messages
781  CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
782 
783  // Send messages
784  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
785  if (!lengths_cols[proc])
786  continue; // no message to send to this proc
787  CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
788  lengths_cols[proc], // message length
789  MPIU_INT, proc, // to proc
790  tag_col, comm, s_waits_col + kk);
791  kk++;
792  }
793 
794  if (nrecvs_cols) {
795  CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
796  }
797  if (nsends_cols) {
798  CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
799  }
800 
801  CHKERR PetscFree(r_waits_col);
802  CHKERR PetscFree(s_waits_col);
803  }
804 
805  CHKERR PetscCommDestroy(&comm);
806  CHKERR PetscFree(status);
807 
808  DofEntity_multiIndex_global_uid_view dofs_glob_uid_view;
809  auto hint = dofs_glob_uid_view.begin();
810  for (auto dof : *m_field.get_dofs())
811  dofs_glob_uid_view.emplace_hint(hint, dof);
812 
813  // set values received from other processors
814  for (int ss = 0; ss != loop_size; ++ss) {
815 
816  int nrecvs;
817  int *olengths;
818  int **data_procs;
819  if (ss == 0) {
820  nrecvs = nrecvs_rows;
821  olengths = olengths_rows;
822  data_procs = rbuf_row;
823  } else {
824  nrecvs = nrecvs_cols;
825  olengths = olengths_cols;
826  data_procs = rbuf_col;
827  }
828 
829  UId uid;
830  for (int kk = 0; kk != nrecvs; ++kk) {
831  int len = olengths[kk];
832  int *data_from_proc = data_procs[kk];
833  for (int dd = 0; dd < len; dd += data_block_size) {
834  uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
835  auto ddit = dofs_glob_uid_view.find(uid);
836 
837  if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
838 
839 #ifndef NDEBUG
840  auto ents_field_ptr = m_field.get_field_ents();
841  MOFEM_LOG("SELF", Sev::error)
842  << "DOF is shared or multishared between processors. For example "
843  "if order of field on given entity is set in inconsistently, "
844  "has different value on two processor, error such as this is "
845  "triggered";
846 
847  MOFEM_LOG("SELF", Sev::error) << "UId " << uid << " is not found";
848  const auto owner_proc = FieldEntity::getOwnerFromUniqueId(uid);
849  MOFEM_LOG("SELF", Sev::error)
850  << "Problematic UId owner proc is " << owner_proc;
851  const auto uid_handle = FieldEntity::getHandleFromUniqueId(uid);
852  MOFEM_LOG("SELF", Sev::error)
853  << "Problematic UId entity owning processor handle is "
854  << uid_handle << " entity type "
855  << moab::CN::EntityTypeName(type_from_handle(uid_handle));
856  const auto uid_bit_number =
858  MOFEM_LOG("SELF", Sev::error)
859  << "Problematic UId field is "
860  << m_field.get_field_name(
861  field_bit_from_bit_number(uid_bit_number));
862 
864  field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
865  auto fe_it = field_view.find(FieldEntity::getGlobalUniqueIdCalculate(
866  owner_proc, uid_bit_number, uid_handle));
867  if (fe_it == field_view.end()) {
868  MOFEM_LOG("SELF", Sev::error)
869  << "Also, no field entity in database for given global UId";
870  } else {
871  MOFEM_LOG("SELF", Sev::error) << "Field entity in databse exist "
872  "(but have no DOF wih give UId";
873  MOFEM_LOG("SELF", Sev::error) << **fe_it;
874 
875  // Save file with missing entity
876  auto error_file_name =
877  "error_with_missing_entity_" +
878  boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
879  ".vtk";
880  MOFEM_LOG("SELF", Sev::error)
881  << "Look to file < " << error_file_name
882  << " > it contains entity with missing DOF.";
883 
884  auto tmp_msh = get_temp_meshset_ptr(m_field.get_moab());
885  const auto local_fe_ent = (*fe_it)->getEnt();
886  CHKERR m_field.get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
887  CHKERR m_field.get_moab().write_file(error_file_name.c_str(), "VTK",
888  "", tmp_msh->get_ptr(), 1);
889  }
890 #endif
891 
892  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
893  "DOF with global UId not found (Compile code in Debug to "
894  "learn more about problem");
895  }
896 
897  auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
898 
899  if (dit != numered_dofs_ptr[ss]->end()) {
900 
901  int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
902 #ifndef NDEBUG
903  if (PetscUnlikely(global_idx < 0))
904  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
905  "received negative dof");
906 #endif
907  bool success;
908  success = numered_dofs_ptr[ss]->modify(
909  dit, NumeredDofEntity_mofem_index_change(global_idx));
910 
911 #ifndef NDEBUG
912  if (PetscUnlikely(!success))
913  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
914  "modification unsuccessful");
915 #endif
916  success = numered_dofs_ptr[ss]->modify(
917  dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
918  global_idx));
919 #ifndef NDEBUG
920  if (PetscUnlikely(!success))
921  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
922  "modification unsuccessful");
923 #endif
924  };
925 
926 #ifndef NDEBUG
927  if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
928 
929  // Dof is shared on this processor, however there is no element
930  // which have this dof. If DOF is not shared and received from other
931  // processor, but not marked as a shared on other that means that is
932  // data inconstancy and error should be thrown.
933 
934  std::ostringstream zz;
935  zz << **ddit << std::endl;
936  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
937  "data inconsistency, dofs is not shared, but received "
938  "from other proc\n"
939  "%s",
940  zz.str().c_str());
941  }
942 #endif
943  }
944  }
945  }
946 
947  if (square_matrix) {
948  (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
949  (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
950  }
951 
952  CHKERR PetscFree(olengths_rows);
953  CHKERR PetscFree(rbuf_row[0]);
954  CHKERR PetscFree(rbuf_row);
955  if (!square_matrix) {
956  CHKERR PetscFree(olengths_cols);
957  CHKERR PetscFree(rbuf_col[0]);
958  CHKERR PetscFree(rbuf_col);
959  }
960 
961  if (square_matrix) {
962  if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
963  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
964  "data inconsistency for square_matrix %d!=%d",
965  numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
966  }
967  if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
968  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
969  "data inconsistency for square_matrix");
970  }
971  }
972 
973  CHKERR printPartitionedProblem(problem_ptr, verb);
974  CHKERR debugPartitionedProblem(problem_ptr, verb);
975 
976  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
977 
979 }
980 
982  const std::string out_name,
983 
984  const std::vector<std::string> &fields_row,
985  const std::vector<std::string> &fields_col,
986 
987  const std::string main_problem, const bool square_matrix,
988 
989  const map<std::string, boost::shared_ptr<Range>> *entityMapRow,
990  const map<std::string, boost::shared_ptr<Range>> *entityMapCol,
991 
992  int verb) {
993  MoFEM::Interface &m_field = cOre;
994  auto problems_ptr = m_field.get_problems();
996 
997  CHKERR m_field.clear_problem(out_name);
998 
999  // get reference to all problems
1000  using ProblemByName = decltype(problems_ptr->get<Problem_mi_tag>());
1001  auto &problems_by_name =
1002  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1003 
1004  // get iterators to out problem, i.e. build problem
1005  auto out_problem_it = problems_by_name.find(out_name);
1006  if (out_problem_it == problems_by_name.end()) {
1007  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1008  "subproblem with name < %s > not defined (top tip check spelling)",
1009  out_name.c_str());
1010  }
1011  // get iterator to main problem, i.e. out problem is subproblem of main
1012  // problem
1013  auto main_problem_it = problems_by_name.find(main_problem);
1014  if (main_problem_it == problems_by_name.end()) {
1015  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1016  "problem of subproblem with name < %s > not defined (top tip "
1017  "check spelling)",
1018  main_problem.c_str());
1019  }
1020 
1021  // get dofs for row & columns for out problem,
1022  boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1023  out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1024  // get dofs for row & columns for main problem
1025  boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1026  main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1027  // get local indices counter
1028  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1029  &out_problem_it->nbLocDofsCol};
1030  // get global indices counter
1031  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1032 
1033  // set number of ghost nodes to zero
1034  {
1035  out_problem_it->nbGhostDofsRow = 0;
1036  out_problem_it->nbGhostDofsCol = 0;
1037  }
1038 
1039  // put rows & columns field names in array
1040  std::vector<std::string> fields[] = {fields_row, fields_col};
1041  const map<std::string, boost::shared_ptr<Range>> *entityMap[] = {
1042  entityMapRow, entityMapCol};
1043 
1044  // make data structure fos sub-problem data
1045  out_problem_it->subProblemData =
1046  boost::make_shared<Problem::SubProblemData>();
1047 
1048  // Loop over rows and columns
1049  for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1050 
1051  // reset dofs and columns counters
1052  (*nb_local_dofs[ss]) = 0;
1053  (*nb_dofs[ss]) = 0;
1054  // clear arrays
1055  out_problem_dofs[ss]->clear();
1056 
1057  // If DOFs are cleared clear finite elements too.
1058  out_problem_it->numeredFiniteElementsPtr->clear();
1059 
1060  // get dofs by field name and insert them in out problem multi-indices
1061  for (auto field : fields[ss]) {
1062 
1063  // Following reserve memory in sequences, only two allocations are here,
1064  // once for array of objects, next for array of shared pointers
1065 
1066  // aliased sequence of pointer is killed with element
1067  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1068  boost::make_shared<std::vector<NumeredDofEntity>>();
1069  // reserve memory for field dofs
1070  if (!ss)
1071  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1072  else
1073  out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1074 
1075  // create elements objects
1076  auto bit_number = m_field.get_field_bit_number(field);
1077 
1078  auto add_dit_to_dofs_array = [&](auto &dit) {
1079  if (dit->get()->getPetscGlobalDofIdx() >= 0)
1080  dofs_array->emplace_back(
1081  dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1082  dit->get()->getPetscGlobalDofIdx(),
1083  dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1084  };
1085 
1086  auto get_dafult_dof_range = [&]() {
1087  auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1088  FieldEntity::getLoBitNumberUId(bit_number));
1089  auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1090  FieldEntity::getHiBitNumberUId(bit_number));
1091  return std::make_pair(dit, hi_dit);
1092  };
1093 
1094  auto check = [&](auto &field) -> boost::shared_ptr<Range> {
1095  if (entityMap[ss]) {
1096  auto mit = entityMap[ss]->find(field);
1097  if (mit != entityMap[ss]->end()) {
1098  return mit->second;
1099  } else {
1100  return nullptr;
1101  }
1102  } else {
1103  return nullptr;
1104  }
1105  };
1106 
1107  if (auto r_ptr = check(field)) {
1108  for (auto p = r_ptr->pair_begin(); p != r_ptr->pair_end(); ++p) {
1109  const auto lo_ent = p->first;
1110  const auto hi_ent = p->second;
1111  auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1112  DofEntity::getLoFieldEntityUId(bit_number, lo_ent));
1113  auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1114  DofEntity::getHiFieldEntityUId(bit_number, hi_ent));
1115  dofs_array->reserve(std::distance(dit, hi_dit));
1116  for (; dit != hi_dit; dit++) {
1117  add_dit_to_dofs_array(dit);
1118  }
1119  }
1120  } else {
1121  auto [dit, hi_dit] = get_dafult_dof_range();
1122  dofs_array->reserve(std::distance(dit, hi_dit));
1123  for (; dit != hi_dit; dit++)
1124  add_dit_to_dofs_array(dit);
1125  }
1126 
1127  // fill multi-index
1128  auto hint = out_problem_dofs[ss]->end();
1129  for (auto &v : *dofs_array)
1130  hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1131  }
1132  // Set local indexes
1133  {
1134  auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1135  auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1136  for (; dit != hi_dit; dit++) {
1137  int idx = -1; // if dof is not part of partition, set local index to -1
1138  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1139  idx = (*nb_local_dofs[ss])++;
1140  }
1141  bool success = out_problem_dofs[ss]->modify(
1142  out_problem_dofs[ss]->project<0>(dit),
1144  if (!success) {
1145  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1146  "operation unsuccessful");
1147  }
1148  };
1149  }
1150  // Set global indexes, compress global indices
1151  {
1152  auto dit =
1153  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1154  auto hi_dit =
1155  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1156  out_problem_dofs[ss]->size());
1157  const int nb = std::distance(dit, hi_dit);
1158  // get main problem global indices
1159  std::vector<int> main_indices(nb);
1160  for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1161  *it = dit->get()->getPetscGlobalDofIdx();
1162  }
1163  // create is with global dofs
1164  IS is;
1165  CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1166  PETSC_USE_POINTER, &is);
1167  // create map form main problem global indices to out problem global
1168  // indices
1169  AO ao;
1170  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1171  if (ss == 0) {
1172  IS is_dup;
1173  CHKERR ISDuplicate(is, &is_dup);
1174  out_problem_it->getSubData()->rowIs = SmartPetscObj<IS>(is_dup, false);
1175  out_problem_it->getSubData()->rowMap = SmartPetscObj<AO>(ao, true);
1176  } else {
1177  IS is_dup;
1178  CHKERR ISDuplicate(is, &is_dup);
1179  out_problem_it->getSubData()->colIs = SmartPetscObj<IS>(is_dup, false);
1180  out_problem_it->getSubData()->colMap = SmartPetscObj<AO>(ao, true);
1181  }
1182  CHKERR AOApplicationToPetscIS(ao, is);
1183  // set global number of DOFs
1184  CHKERR ISGetSize(is, nb_dofs[ss]);
1185  CHKERR ISDestroy(&is);
1186  // set out problem global indices after applying map
1187  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1188  for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1189  dit++, it++) {
1190  bool success = out_problem_dofs[ss]->modify(
1191  out_problem_dofs[ss]->project<0>(dit),
1193  dit->get()->getPart(), *it, *it,
1194  dit->get()->getPetscLocalDofIdx()));
1195  if (!success) {
1196  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1197  "operation unsuccessful");
1198  }
1199  }
1200  // set global indices to nodes not on this part
1201  {
1202  NumeredDofEntityByLocalIdx::iterator dit =
1203  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1204  NumeredDofEntityByLocalIdx::iterator hi_dit =
1205  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1206  const int nb = std::distance(dit, hi_dit);
1207  std::vector<int> main_indices_non_local(nb);
1208  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1209  dit++, it++) {
1210  *it = dit->get()->getPetscGlobalDofIdx();
1211  }
1212  IS is;
1213  CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1214  &*main_indices_non_local.begin(),
1215  PETSC_USE_POINTER, &is);
1216  CHKERR AOApplicationToPetscIS(ao, is);
1217  CHKERR ISDestroy(&is);
1218  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1219  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1220  dit++, it++) {
1221  bool success = out_problem_dofs[ss]->modify(
1222  out_problem_dofs[ss]->project<0>(dit),
1224  dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1225  dit->get()->getPetscLocalDofIdx()));
1226  if (!success) {
1227  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1228  "operation unsuccessful");
1229  }
1230  }
1231  }
1232  CHKERR AODestroy(&ao);
1233  }
1234  }
1235 
1236  if (square_matrix) {
1237  out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1238  out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1239  out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1240  out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1241  out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1242  }
1243 
1244  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1245  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1246 
1248 }
1249 
1251  const std::string out_name, const std::vector<std::string> add_row_problems,
1252  const std::vector<std::string> add_col_problems, const bool square_matrix,
1253  int verb) {
1255  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1256  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1257  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1258  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1259  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1260  MoFEM::Interface &m_field = cOre;
1261  auto problems_ptr = m_field.get_problems();
1263 
1264  CHKERR m_field.clear_problem(out_name);
1265  // get reference to all problems
1267  ProblemByName &problems_by_name =
1268  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1269 
1270  // Get iterators to out problem, i.e. build problem
1271  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1272  if (out_problem_it == problems_by_name.end()) {
1273  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1274  "problem with name < %s > not defined (top tip check spelling)",
1275  out_name.c_str());
1276  }
1277  // Make data structure for composed-problem data
1278  out_problem_it->composedProblemsData =
1279  boost::make_shared<ComposedProblemsData>();
1280  boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1281  out_problem_it->getComposedProblemsData();
1282 
1283  const std::vector<std::string> *add_prb[] = {&add_row_problems,
1284  &add_col_problems};
1285  std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1286  &cmp_prb_data->colProblemsAdd};
1287  std::vector<SmartPetscObj<IS>> *add_prb_is[] = {&cmp_prb_data->rowIs,
1288  &cmp_prb_data->colIs};
1289 
1290  // Get local indices counter
1291  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1292  &out_problem_it->nbLocDofsCol};
1293  // Get global indices counter
1294  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1295 
1296  // Set number of ghost nodes to zero
1297  {
1298  out_problem_it->nbDofsRow = 0;
1299  out_problem_it->nbDofsCol = 0;
1300  out_problem_it->nbLocDofsRow = 0;
1301  out_problem_it->nbLocDofsCol = 0;
1302  out_problem_it->nbGhostDofsRow = 0;
1303  out_problem_it->nbGhostDofsCol = 0;
1304  }
1305  int nb_dofs_reserve[] = {0, 0};
1306 
1307  // Loop over rows and columns in the main problem and sub-problems
1308  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1309  add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1310  add_prb_is[ss]->reserve(add_prb[ss]->size());
1311  for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1312  vit != add_prb[ss]->end(); vit++) {
1313  ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1314  if (prb_it == problems_by_name.end()) {
1315  SETERRQ1(
1316  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1317  "problem with name < %s > not defined (top tip check spelling)",
1318  vit->c_str());
1319  }
1320  add_prb_ptr[ss]->push_back(&*prb_it);
1321  // set number of dofs on rows and columns
1322  if (ss == 0) {
1323  // row
1324  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1325  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1326  nb_dofs_reserve[ss] +=
1327  add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1328  } else {
1329  // column
1330  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1331  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1332  nb_dofs_reserve[ss] +=
1333  add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1334  }
1335  }
1336  }
1337  // if squre problem, rows and columns are the same
1338  if (square_matrix) {
1339  add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1340  add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1341  out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1342  *nb_dofs[1] = *nb_dofs[0];
1343  *nb_local_dofs[1] = *nb_local_dofs[0];
1344  }
1345 
1346  // reserve memory for dofs
1347  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1348  // Reserve memory
1349  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1350  dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1351  dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1352  if (!ss)
1353  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1354  else
1355  out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1356  }
1357 
1358  // Push back DOFs
1359  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1360  NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1361  dit,
1362  hi_dit;
1363  int shift_glob = 0;
1364  int shift_loc = 0;
1365  for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1366  PetscInt *dofs_out_idx_ptr;
1367  int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1368  CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1369  if (ss == 0) {
1370  dit = (*add_prb_ptr[ss])[pp]
1371  ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1372  .begin();
1373  hi_dit = (*add_prb_ptr[ss])[pp]
1374  ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1375  .end();
1376  } else {
1377  dit = (*add_prb_ptr[ss])[pp]
1378  ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1379  .begin();
1380  hi_dit = (*add_prb_ptr[ss])[pp]
1381  ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1382  .end();
1383  }
1384  int is_nb = 0;
1385  for (; dit != hi_dit; dit++) {
1386  const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
1387  const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
1388  const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
1389  if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1390  continue;
1391  const int rank = m_field.get_comm_rank();
1392  const int part = dit->get()->getPart();
1393  const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1394  const int loc_idx =
1395  (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1396  : -1;
1397  dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1398  glob_idx, loc_idx, part);
1399  if (part == rank) {
1400  dofs_out_idx_ptr[is_nb++] = glob_idx;
1401  }
1402  }
1403  if (is_nb > nb_local_dofs) {
1404  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1405  "Data inconsistency");
1406  }
1407  IS is;
1408  CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1409  PETSC_OWN_POINTER, &is);
1410  auto smart_is = SmartPetscObj<IS>(is);
1411  (*add_prb_is[ss]).push_back(smart_is);
1412  if (ss == 0) {
1413  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1414  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1415  } else {
1416  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1417  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1418  }
1419  if (square_matrix) {
1420  (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1421  (*add_prb_is[1]).push_back(smart_is);
1422  }
1423  }
1424  }
1425 
1426  if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1427  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1428  }
1429 
1430  // Insert DOFs to problem multi-index
1431  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1432  auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1433  : out_problem_it->numeredColDofsPtr->end();
1434  for (auto &v : *dofs_array[ss])
1435  hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1436  hint, dofs_array[ss], &v)
1437  : out_problem_it->numeredColDofsPtr->emplace_hint(
1438  hint, dofs_array[ss], &v);
1439  }
1440 
1441  // Compress DOFs
1442  *nb_dofs[0] = 0;
1443  *nb_dofs[1] = 0;
1444  *nb_local_dofs[0] = 0;
1445  *nb_local_dofs[1] = 0;
1446  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1447 
1448  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1449  if (ss == 0) {
1450  dofs_ptr = out_problem_it->numeredRowDofsPtr;
1451  } else {
1452  dofs_ptr = out_problem_it->numeredColDofsPtr;
1453  }
1454  NumeredDofEntityByUId::iterator dit, hi_dit;
1455  dit = dofs_ptr->get<Unique_mi_tag>().begin();
1456  hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1457  std::vector<int> idx;
1458  idx.reserve(std::distance(dit, hi_dit));
1459  // set dofs in order entity and dof number on entity
1460  for (; dit != hi_dit; dit++) {
1461  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1462  bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1463  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1464  if (!success) {
1465  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1466  "modification unsuccessful");
1467  }
1468  idx.push_back(dit->get()->getPetscGlobalDofIdx());
1469  } else {
1470  if (dit->get()->getPetscLocalDofIdx() != -1) {
1471  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1472  "local index should be negative");
1473  }
1474  }
1475  }
1476  if (square_matrix) {
1477  *nb_local_dofs[1] = *nb_local_dofs[0];
1478  }
1479 
1480  // set new dofs mapping
1481  IS is;
1482  CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1483  PETSC_USE_POINTER, &is);
1484  CHKERR ISGetSize(is, nb_dofs[ss]);
1485  if (square_matrix) {
1486  *nb_dofs[1] = *nb_dofs[0];
1487  }
1488 
1489  AO ao;
1490  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1491  for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1492  CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1493 
1494  // Set DOFs numeration
1495  {
1496  std::vector<int> idx_new;
1497  idx_new.reserve(dofs_ptr->size());
1498  for (NumeredDofEntityByUId::iterator dit =
1499  dofs_ptr->get<Unique_mi_tag>().begin();
1500  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1501  idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1502  }
1503  // set new global dofs numeration
1504  IS is_new;
1505  CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1506  &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1507  CHKERR AOApplicationToPetscIS(ao, is_new);
1508  // set global indices to multi-index
1509  std::vector<int>::iterator vit = idx_new.begin();
1510  for (NumeredDofEntityByUId::iterator dit =
1511  dofs_ptr->get<Unique_mi_tag>().begin();
1512  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1513  bool success =
1514  dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1515  dit->get()->getPart(), *(vit++)));
1516  if (!success) {
1517  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1518  "modification unsuccessful");
1519  }
1520  }
1521  CHKERR ISDestroy(&is_new);
1522  }
1523  CHKERR ISDestroy(&is);
1524  CHKERR AODestroy(&ao);
1525  }
1526 
1527  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1528  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1529 
1530  // Inidcate that porble has been build
1533 
1535 }
1536 
1538  int verb) {
1539 
1540  MoFEM::Interface &m_field = cOre;
1541  auto problems_ptr = m_field.get_problems();
1543 
1545  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1546  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1547  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1548  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1549  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1551  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1552  MOFEM_LOG("WORLD", Sev::verbose) << "Simple partition problem " << name;
1553 
1554  // find p_miit
1556  ProblemByName &problems_set =
1557  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1558  ProblemByName::iterator p_miit = problems_set.find(name);
1559  if (p_miit == problems_set.end()) {
1560  SETERRQ1(PETSC_COMM_SELF, 1,
1561  "problem < %s > is not found (top tip: check spelling)",
1562  name.c_str());
1563  }
1564  typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1565  Idx_mi_tag>::type NumeredDofEntitysByIdx;
1566  NumeredDofEntitysByIdx &dofs_row_by_idx =
1567  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>();
1568  NumeredDofEntitysByIdx &dofs_col_by_idx =
1569  p_miit->numeredColDofsPtr->get<Idx_mi_tag>();
1570  boost::multi_index::index<NumeredDofEntity_multiIndex,
1571  Idx_mi_tag>::type::iterator miit_row,
1572  hi_miit_row;
1573  boost::multi_index::index<NumeredDofEntity_multiIndex,
1574  Idx_mi_tag>::type::iterator miit_col,
1575  hi_miit_col;
1576  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1577  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1578  nb_row_local_dofs = 0;
1579  nb_row_ghost_dofs = 0;
1580  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1581  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1582  nb_col_local_dofs = 0;
1583  nb_col_ghost_dofs = 0;
1584 
1585  bool square_matrix = false;
1586  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1587  square_matrix = true;
1588  }
1589 
1590  // get row range of local indices
1591  PetscLayout layout_row;
1592  const int *ranges_row;
1593 
1594  DofIdx nb_dofs_row = dofs_row_by_idx.size();
1595  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1596  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1597  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1598  CHKERR PetscLayoutSetUp(layout_row);
1599  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1600  // get col range of local indices
1601  PetscLayout layout_col;
1602  const int *ranges_col;
1603  if (!square_matrix) {
1604  DofIdx nb_dofs_col = dofs_col_by_idx.size();
1605  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1606  CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1607  CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1608  CHKERR PetscLayoutSetUp(layout_col);
1609  CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1610  }
1611  for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1612  part++) {
1613  miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1614  hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1615  if (std::distance(miit_row, hi_miit_row) !=
1616  ranges_row[part + 1] - ranges_row[part]) {
1617  SETERRQ4(
1618  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1619  "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1620  "rstart (%d != %d - %d = %d) ",
1621  std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1622  ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1623  }
1624  // loop rows
1625  for (; miit_row != hi_miit_row; miit_row++) {
1626  bool success = dofs_row_by_idx.modify(
1627  miit_row,
1628  NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1629  if (!success)
1630  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1631  "modification unsuccessful");
1632  if (part == (unsigned int)m_field.get_comm_rank()) {
1633  success = dofs_row_by_idx.modify(
1634  miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1635  if (!success)
1636  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1637  "modification unsuccessful");
1638  }
1639  }
1640  if (!square_matrix) {
1641  miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1642  hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1643  if (std::distance(miit_col, hi_miit_col) !=
1644  ranges_col[part + 1] - ranges_col[part]) {
1645  SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1646  "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1647  "rend - "
1648  "rstart (%d != %d - %d = %d) ",
1649  std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1650  ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1651  }
1652  // loop cols
1653  for (; miit_col != hi_miit_col; miit_col++) {
1654  bool success = dofs_col_by_idx.modify(
1656  part, (*miit_col)->dofIdx));
1657  if (!success)
1658  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1659  "modification unsuccessful");
1660  if (part == (unsigned int)m_field.get_comm_rank()) {
1661  success = dofs_col_by_idx.modify(
1662  miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1663  if (!success)
1664  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1665  "modification unsuccessful");
1666  }
1667  }
1668  }
1669  }
1670  CHKERR PetscLayoutDestroy(&layout_row);
1671  if (!square_matrix) {
1672  CHKERR PetscLayoutDestroy(&layout_col);
1673  }
1674  if (square_matrix) {
1675  nb_col_local_dofs = nb_row_local_dofs;
1676  nb_col_ghost_dofs = nb_row_ghost_dofs;
1677  }
1678  CHKERR printPartitionedProblem(&*p_miit, verb);
1681 }
1682 
1684  int verb) {
1685  MoFEM::Interface &m_field = cOre;
1686  auto problems_ptr = m_field.get_problems();
1688 
1689  MOFEM_LOG("WORLD", Sev::noisy) << "Partition problem " << name;
1690 
1691  using NumeredDofEntitysByIdx =
1693  using ProblemsByName = Problem_multiIndex::index<Problem_mi_tag>::type;
1694 
1695  // Find problem pointer by name
1696  auto &problems_set =
1697  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1698  auto p_miit = problems_set.find(name);
1699  if (p_miit == problems_set.end()) {
1700  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1701  "problem with name %s not defined (top tip check spelling)",
1702  name.c_str());
1703  }
1704  int nb_dofs_row = p_miit->getNbDofsRow();
1705 
1706  if (m_field.get_comm_size() != 1) {
1707 
1708  if (!(cOre.getBuildMoFEM() & (1 << 0)))
1709  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1710  if (!(cOre.getBuildMoFEM() & (1 << 1)))
1711  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1712  if (!(cOre.getBuildMoFEM() & (1 << 2)))
1713  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1714  "entFEAdjacencies not build");
1715  if (!(cOre.getBuildMoFEM() & (1 << 3)))
1716  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1717 
1718  Mat Adj;
1719  CHKERR m_field.getInterface<MatrixManager>()
1720  ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1721 
1722  int m, n;
1723  CHKERR MatGetSize(Adj, &m, &n);
1724  if (verb > VERY_VERBOSE)
1725  MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1726 
1727  // partitioning
1728  MatPartitioning part;
1729  IS is;
1730  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
1731  CHKERR MatPartitioningSetAdjacency(part, Adj);
1732  CHKERR MatPartitioningSetFromOptions(part);
1733  CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
1734  CHKERR MatPartitioningApply(part, &is);
1735  if (verb > VERY_VERBOSE)
1736  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
1737 
1738  // gather
1739  IS is_gather, is_num, is_gather_num;
1740  CHKERR ISAllGather(is, &is_gather);
1741  CHKERR ISPartitioningToNumbering(is, &is_num);
1742  CHKERR ISAllGather(is_num, &is_gather_num);
1743  const int *part_number, *petsc_idx;
1744  CHKERR ISGetIndices(is_gather, &part_number);
1745  CHKERR ISGetIndices(is_gather_num, &petsc_idx);
1746  int size_is_num, size_is_gather;
1747  CHKERR ISGetSize(is_gather, &size_is_gather);
1748  if (size_is_gather != (int)nb_dofs_row)
1749  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1750  "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
1751 
1752  CHKERR ISGetSize(is_num, &size_is_num);
1753  if (size_is_num != (int)nb_dofs_row)
1754  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1755  "data inconsistency %d != %d", size_is_num, nb_dofs_row);
1756 
1757  bool square_matrix = false;
1758  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1759  square_matrix = true;
1760 
1761  // if (!square_matrix) {
1762  // // FIXME: This is for back compatibility, if deprecate interface
1763  // function
1764  // // build interfaces is removed, this part of the code will be obsolete
1765  // auto mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().begin();
1766  // auto hi_mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().end();
1767  // auto mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().begin();
1768  // auto hi_mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().end();
1769  // for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
1770  // if (mit_col == hi_mit_col) {
1771  // SETERRQ(
1772  // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1773  // "check finite element definition, nb. of rows is not equal to "
1774  // "number for columns");
1775  // }
1776  // if (mit_row->get()->getLocalUniqueId() !=
1777  // mit_col->get()->getLocalUniqueId()) {
1778  // SETERRQ(
1779  // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1780  // "check finite element definition, nb. of rows is not equal to "
1781  // "number for columns");
1782  // }
1783  // }
1784  // }
1785 
1786  auto number_dofs = [&](auto &dofs_idx, auto &counter) {
1788  for (auto miit_dofs_row = dofs_idx.begin();
1789  miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
1790  const int part = part_number[(*miit_dofs_row)->dofIdx];
1791  if (part == (unsigned int)m_field.get_comm_rank()) {
1792  const bool success = dofs_idx.modify(
1793  miit_dofs_row,
1795  part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
1796  if (!success) {
1797  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1798  "modification unsuccessful");
1799  }
1800  } else {
1801  const bool success = dofs_idx.modify(
1803  part, petsc_idx[(*miit_dofs_row)->dofIdx]));
1804  if (!success) {
1805  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1806  "modification unsuccessful");
1807  }
1808  }
1809  }
1811  };
1812 
1813  // Set petsc global indices
1814  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1815  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1816  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1817  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1818  nb_row_local_dofs = 0;
1819  nb_row_ghost_dofs = 0;
1820 
1821  CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1822 
1823  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1824  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1825  if (square_matrix) {
1826  nb_col_local_dofs = nb_row_local_dofs;
1827  nb_col_ghost_dofs = nb_row_ghost_dofs;
1828  } else {
1829  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1830  const_cast<NumeredDofEntitysByIdx &>(
1831  p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1832  nb_col_local_dofs = 0;
1833  nb_col_ghost_dofs = 0;
1834  CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1835  }
1836 
1837  CHKERR ISRestoreIndices(is_gather, &part_number);
1838  CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
1839  CHKERR ISDestroy(&is_num);
1840  CHKERR ISDestroy(&is_gather_num);
1841  CHKERR ISDestroy(&is_gather);
1842  CHKERR ISDestroy(&is);
1843  CHKERR MatPartitioningDestroy(&part);
1844  CHKERR MatDestroy(&Adj);
1845  CHKERR printPartitionedProblem(&*p_miit, verb);
1846  } else {
1847 
1848  auto number_dofs = [&](auto &dof_idx, auto &counter) {
1850  for (auto miit_dofs_row = dof_idx.begin(); miit_dofs_row != dof_idx.end();
1851  miit_dofs_row++) {
1852  const bool success = dof_idx.modify(
1853  miit_dofs_row,
1854  NumeredDofEntity_part_and_indices_change(0, counter, counter));
1855  ++counter;
1856  if (!success) {
1857  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1858  "modification unsuccessful");
1859  }
1860  }
1862  };
1863 
1864  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1865  p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1866  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1867  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1868  nb_row_local_dofs = 0;
1869  nb_row_ghost_dofs = 0;
1870 
1871  CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1872 
1873  bool square_matrix = false;
1874  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1875  square_matrix = true;
1876 
1877  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1878  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1879  if (square_matrix) {
1880  nb_col_local_dofs = nb_row_local_dofs;
1881  nb_col_ghost_dofs = nb_row_ghost_dofs;
1882  } else {
1883  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1884  const_cast<NumeredDofEntitysByIdx &>(
1885  p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1886  nb_col_local_dofs = 0;
1887  nb_col_ghost_dofs = 0;
1888  CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1889  }
1890  }
1891 
1892  cOre.getBuildMoFEM() |= 1 << 4;
1894 }
1895 
1897  const std::string name, const std::string problem_for_rows, bool copy_rows,
1898  const std::string problem_for_cols, bool copy_cols, int verb) {
1899  MoFEM::Interface &m_field = cOre;
1900  auto problems_ptr = m_field.get_problems();
1902 
1904  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
1905 
1907 
1908  // find p_miit
1909  ProblemByName &problems_by_name =
1910  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1911  ProblemByName::iterator p_miit = problems_by_name.find(name);
1912  if (p_miit == problems_by_name.end()) {
1913  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
1914  "problem with name < %s > not defined (top tip check spelling)",
1915  name.c_str());
1916  }
1917  if (verb > QUIET)
1918  MOFEM_LOG("WORLD", Sev::inform)
1919  << p_miit->getName() << " from rows of " << problem_for_rows
1920  << " and columns of " << problem_for_cols;
1921 
1922  // find p_miit_row
1923  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
1924  if (p_miit_row == problems_by_name.end()) {
1925  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1926  "problem with name < %s > not defined (top tip check spelling)",
1927  problem_for_rows.c_str());
1928  }
1929  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
1930  p_miit_row->numeredRowDofsPtr;
1931 
1932  // find p_mit_col
1933  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
1934  if (p_miit_col == problems_by_name.end()) {
1935  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1936  "problem with name < %s > not defined (top tip check spelling)",
1937  problem_for_cols.c_str());
1938  }
1939  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
1940  p_miit_col->numeredColDofsPtr;
1941 
1942  bool copy[] = {copy_rows, copy_cols};
1943  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
1944  p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
1945 
1946  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
1947  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
1948  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
1949  dofs_col};
1950 
1951  for (int ss = 0; ss < 2; ss++) {
1952 
1953  // build indices
1954  *nb_local_dofs[ss] = 0;
1955  if (!copy[ss]) {
1956 
1957  // only copy indices which are belong to some elements if this problem
1958  std::vector<int> is_local, is_new;
1959 
1960  NumeredDofEntityByUId &dofs_by_uid =
1961  copied_dofs[ss]->get<Unique_mi_tag>();
1962  for (NumeredDofEntity_multiIndex::iterator dit =
1963  composed_dofs[ss]->begin();
1964  dit != composed_dofs[ss]->end(); dit++) {
1965  NumeredDofEntityByUId::iterator diit =
1966  dofs_by_uid.find((*dit)->getLocalUniqueId());
1967  if (diit == dofs_by_uid.end()) {
1968  SETERRQ(
1970  "data inconsistency, could not find dof in composite problem");
1971  }
1972  int part_number = (*diit)->getPart(); // get part number
1973  int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
1974  bool success;
1975  success = composed_dofs[ss]->modify(
1977  petsc_global_dof));
1978  if (!success) {
1979  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1980  "modification unsuccessful");
1981  }
1982  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
1983  success = composed_dofs[ss]->modify(
1984  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1985  if (!success) {
1986  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1987  "modification unsuccessful");
1988  }
1989  is_local.push_back(petsc_global_dof);
1990  }
1991  }
1992 
1993  AO ao;
1994  CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
1995  NULL, &ao);
1996 
1997  // apply local to global mapping
1998  is_local.resize(0);
1999  for (NumeredDofEntity_multiIndex::iterator dit =
2000  composed_dofs[ss]->begin();
2001  dit != composed_dofs[ss]->end(); dit++) {
2002  is_local.push_back((*dit)->getPetscGlobalDofIdx());
2003  }
2004  CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2005  int idx2 = 0;
2006  for (NumeredDofEntity_multiIndex::iterator dit =
2007  composed_dofs[ss]->begin();
2008  dit != composed_dofs[ss]->end(); dit++) {
2009  int part_number = (*dit)->getPart(); // get part number
2010  int petsc_global_dof = is_local[idx2++];
2011  bool success;
2012  success = composed_dofs[ss]->modify(
2014  petsc_global_dof));
2015  if (!success) {
2016  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2017  "modification unsuccessful");
2018  }
2019  }
2020 
2021  CHKERR AODestroy(&ao);
2022 
2023  } else {
2024 
2025  for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2026  dit != copied_dofs[ss]->end(); dit++) {
2027  std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2028  p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2029  new NumeredDofEntity((*dit)->getDofEntityPtr())));
2030  if (p.second) {
2031  (*nb_dofs[ss])++;
2032  }
2033  int dof_idx = (*dit)->getDofIdx();
2034  int part_number = (*dit)->getPart(); // get part number
2035  int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2036  if (part_number == (unsigned int)m_field.get_comm_rank()) {
2037  const bool success = composed_dofs[ss]->modify(
2039  part_number, dof_idx, petsc_global_dof,
2040  (*nb_local_dofs[ss])++));
2041  if (!success) {
2042  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2043  "modification unsuccessful");
2044  }
2045  } else {
2046  const bool success = composed_dofs[ss]->modify(
2048  part_number, dof_idx, petsc_global_dof));
2049  if (!success) {
2050  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2051  "modification unsuccessful");
2052  }
2053  }
2054  }
2055  }
2056  }
2057 
2058  CHKERR printPartitionedProblem(&*p_miit, verb);
2059  CHKERR debugPartitionedProblem(&*p_miit, verb);
2060 
2062 }
2063 
2065 ProblemsManager::printPartitionedProblem(const Problem *problem_ptr, int verb) {
2066  MoFEM::Interface &m_field = cOre;
2068 
2069  if (verb > QUIET) {
2070 
2071  MOFEM_LOG("SYNC", Sev::inform)
2072  << problem_ptr->getName() << " Nb. local dof "
2073  << problem_ptr->getNbLocalDofsRow() << " by "
2074  << problem_ptr->getNbLocalDofsCol() << " nb global dofs "
2075  << problem_ptr->getNbDofsRow() << " by " << problem_ptr->getNbDofsCol();
2076 
2077  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2078  }
2079 
2081 }
2082 
2084 ProblemsManager::debugPartitionedProblem(const Problem *problem_ptr, int verb) {
2085  MoFEM::Interface &m_field = cOre;
2087 
2088  auto save_ent = [](moab::Interface &moab, const std::string name,
2089  const EntityHandle ent) {
2091  EntityHandle out_meshset;
2092  CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
2093  CHKERR moab.add_entities(out_meshset, &ent, 1);
2094  CHKERR moab.write_file(name.c_str(), "VTK", "", &out_meshset, 1);
2095  CHKERR moab.delete_entities(&out_meshset, 1);
2097  };
2098 
2099  if (debug > 0) {
2100 
2102  NumeredDofEntitysByIdx;
2103  NumeredDofEntitysByIdx::iterator dit, hi_dit;
2104  const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
2105  &(problem_ptr->numeredRowDofsPtr->get<Idx_mi_tag>()),
2106  &(problem_ptr->numeredColDofsPtr->get<Idx_mi_tag>())};
2107 
2108  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
2109  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
2110  &problem_ptr->nbLocDofsCol};
2111 
2112  for (int ss = 0; ss < 2; ss++) {
2113 
2114  dit = numered_dofs_ptr[ss]->begin();
2115  hi_dit = numered_dofs_ptr[ss]->end();
2116  for (; dit != hi_dit; dit++) {
2117  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2118  if ((*dit)->getPetscLocalDofIdx() < 0) {
2119  std::ostringstream zz;
2120  zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2121  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2122  "local dof index for %d (0-row, 1-col) not set, i.e. has "
2123  "negative value\n %s",
2124  ss, zz.str().c_str());
2125  }
2126  if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
2127  std::ostringstream zz;
2128  zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2129  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2130  "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
2131  zz.str().c_str());
2132  }
2133  } else {
2134  if ((*dit)->getPetscGlobalDofIdx() < 0) {
2135 
2136  const EntityHandle ent = (*dit)->getEnt();
2137  CHKERR save_ent(
2138  m_field.get_moab(),
2139  "debug_part" +
2140  boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
2141  "_negative_global_index.vtk",
2142  ent);
2143 
2144  std::ostringstream zz;
2145  zz << "rank " << m_field.get_comm_rank() << " "
2146  << dit->get()->getBitRefLevel() << " " << **dit;
2147  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2148  "global dof index for %d (0-row, 1-col) row not set, i.e. "
2149  "has negative value\n %s",
2150  ss, zz.str().c_str());
2151  }
2152  if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
2153  std::ostringstream zz;
2154  zz << "rank " << m_field.get_comm_rank() << " nb_dofs "
2155  << *nbdof_ptr[ss] << " " << **dit;
2156  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2157  "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
2158  zz.str().c_str());
2159  }
2160  }
2161  }
2162  }
2163  }
2165 }
2166 
2168  bool part_from_moab,
2169  int low_proc,
2170  int hi_proc, int verb) {
2171  MoFEM::Interface &m_field = cOre;
2172  auto problems_ptr = m_field.get_problems();
2173  auto fe_ent_ptr = m_field.get_ents_finite_elements();
2175 
2177  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2178 
2179  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2180  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2181 
2182  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2183  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2184  "adjacencies not build");
2185 
2187  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2188 
2190  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2191  "problem not partitioned");
2192 
2193  if (low_proc == -1)
2194  low_proc = m_field.get_comm_rank();
2195  if (hi_proc == -1)
2196  hi_proc = m_field.get_comm_rank();
2197 
2198  // Find pointer to problem of given name
2200  auto &problems =
2201  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2202  ProblemByName::iterator p_miit = problems.find(name);
2203  if (p_miit == problems.end()) {
2204  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2205  "problem < %s > not found (top tip: check spelling)",
2206  name.c_str());
2207  }
2208 
2209  // Get reference on finite elements multi-index on the problem
2210  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2211  *p_miit->numeredFiniteElementsPtr;
2212 
2213  // Clear all elements and data, build it again
2214  problem_finite_elements.clear();
2215 
2216  // Check if dofs and columns are the same, i.e. structurally symmetric
2217  // problem
2218  bool do_cols_prob = true;
2219  if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2220  do_cols_prob = false;
2221  }
2222 
2223  auto get_good_elems = [&]() {
2224  auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
2225  good_elems.reserve(fe_ent_ptr->size());
2226 
2227  const auto prb_bit = p_miit->getBitRefLevel();
2228  const auto prb_mask = p_miit->getBitRefLevelMask();
2229 
2230  // Loop over all elements in database and if right element is there add it
2231  // to problem finite element multi-index
2232  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2233 
2234  // if element is not part of problem
2235  if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2236 
2237  const auto &fe_bit = (*efit)->getBitRefLevel();
2238 
2239  // if entity is not problem refinement level
2240  if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit).any())
2241  good_elems.emplace_back(efit);
2242  }
2243  }
2244 
2245  return good_elems;
2246  };
2247 
2248  auto good_elems = get_good_elems();
2249 
2250  auto numbered_good_elems_ptr =
2251  boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2252  numbered_good_elems_ptr->reserve(good_elems.size());
2253  for (auto &efit : good_elems)
2254  numbered_good_elems_ptr->emplace_back(NumeredEntFiniteElement(*efit));
2255 
2256  if (!do_cols_prob) {
2257  for (auto &fe : *numbered_good_elems_ptr) {
2258  if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2259  fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2260  }
2261  }
2262  }
2263 
2264  if (part_from_moab) {
2265  for (auto &fe : *numbered_good_elems_ptr) {
2266  // if partition is taken from moab partition
2267  int proc = fe.getPartProc();
2268  if (proc == -1 && fe.getEntType() == MBVERTEX)
2269  proc = fe.getOwnerProc();
2270  fe.part = proc;
2271  }
2272  }
2273 
2274  for (auto &fe : *numbered_good_elems_ptr) {
2275 
2277  CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2278 
2279  if (!part_from_moab) {
2280  std::vector<int> parts(m_field.get_comm_size(), 0);
2281  for (auto &dof_ptr : rows_view)
2282  parts[dof_ptr->pArt]++;
2283  std::vector<int>::iterator pos = max_element(parts.begin(), parts.end());
2284  const auto max_part = std::distance(parts.begin(), pos);
2285  fe.part = max_part;
2286  }
2287  }
2288 
2289  for (auto &fe : *numbered_good_elems_ptr) {
2290 
2291  auto check_fields_and_dofs = [&]() {
2292  if (!part_from_moab) {
2293  if (fe.getBitFieldIdRow().none() && m_field.get_comm_size() == 0) {
2294  MOFEM_LOG("WORLD", Sev::warning)
2295  << "At least one field has to be added to element row to "
2296  "determine partition of finite element. Check element " +
2297  boost::lexical_cast<std::string>(fe.getName());
2298  }
2299  }
2300 
2301  return true;
2302  };
2303 
2304  if (check_fields_and_dofs()) {
2305  // Add element to the problem
2306  auto p = problem_finite_elements.insert(
2307  boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
2308  &fe));
2309  if (!p.second)
2310  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2311  }
2312  }
2313 
2314  if (verb >= VERBOSE) {
2315  auto elements_on_rank =
2316  problem_finite_elements.get<Part_mi_tag>().equal_range(
2317  m_field.get_comm_rank());
2318  MOFEM_LOG("SYNC", Sev::verbose)
2319  << p_miit->getName() << " nb. elems "
2320  << std::distance(elements_on_rank.first, elements_on_rank.second);
2321  auto fe_ptr = m_field.get_finite_elements();
2322  for (auto &fe : *fe_ptr) {
2323  auto e_range =
2324  problem_finite_elements.get<Composite_Name_And_Part_mi_tag>()
2325  .equal_range(
2326  boost::make_tuple(fe->getName(), m_field.get_comm_rank()));
2327  MOFEM_LOG("SYNC", Sev::noisy)
2328  << "Element " << fe->getName() << " nb. elems "
2329  << std::distance(e_range.first, e_range.second);
2330  }
2331 
2332  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
2333  }
2334 
2337 }
2338 
2340  int verb) {
2341  MoFEM::Interface &m_field = cOre;
2342  auto problems_ptr = m_field.get_problems();
2344 
2346  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2347  "partition of problem not build");
2349  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2350  "partitions finite elements not build");
2351 
2352  // get problem pointer
2353  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2354  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2355  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2356  name.c_str());
2357 
2358  // get reference to number of ghost dofs
2359  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2360  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2361  nb_row_ghost_dofs = 0;
2362  nb_col_ghost_dofs = 0;
2363 
2364  // do work if more than one processor
2365  if (m_field.get_comm_size() > 1) {
2366 
2368  ghost_idx_col_view;
2369 
2370  // get elements on this partition
2371  auto fe_range =
2372  p_miit->numeredFiniteElementsPtr->get<Part_mi_tag>().equal_range(
2373  m_field.get_comm_rank());
2374 
2375  // get dofs on elements which are not part of this partition
2376 
2377  struct Inserter {
2378  using Vec = std::vector<boost::shared_ptr<NumeredDofEntity>>;
2379  using It = Vec::iterator;
2380  It operator()(Vec &dofs_view, It &hint,
2381  boost::shared_ptr<NumeredDofEntity> &&dof) {
2382  dofs_view.emplace_back(dof);
2383  return dofs_view.end();
2384  }
2385  };
2386 
2387  // rows
2388  std::vector<boost::shared_ptr<NumeredDofEntity>> fe_vec_view;
2389  auto hint_r = ghost_idx_row_view.begin();
2390  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2391 
2392  fe_vec_view.clear();
2393  CHKERR EntFiniteElement::getDofView((*fe_ptr)->getRowFieldEnts(),
2394  *(p_miit->getNumeredRowDofsPtr()),
2395  fe_vec_view, Inserter());
2396 
2397  for (auto &dof_ptr : fe_vec_view) {
2398  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2399  hint_r = ghost_idx_row_view.emplace_hint(hint_r, dof_ptr);
2400  }
2401  }
2402  }
2403 
2404  // columns
2405  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2406 
2407  auto hint_c = ghost_idx_col_view.begin();
2408  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2409 
2410  fe_vec_view.clear();
2411  CHKERR EntFiniteElement::getDofView((*fe_ptr)->getColFieldEnts(),
2412  *(p_miit->getNumeredColDofsPtr()),
2413  fe_vec_view, Inserter());
2414 
2415  for (auto &dof_ptr : fe_vec_view) {
2416  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2417  hint_c = ghost_idx_col_view.emplace_hint(hint_c, dof_ptr);
2418  }
2419  }
2420  }
2421  }
2422 
2423  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2424  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2425 
2426  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2427  &ghost_idx_row_view, &ghost_idx_col_view};
2428  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2429  &p_miit->numeredRowDofsPtr->get<Unique_mi_tag>(),
2430  &p_miit->numeredColDofsPtr->get<Unique_mi_tag>()};
2431 
2432  int loop_size = 2;
2433  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2434  loop_size = 1;
2435  }
2436 
2437  // set local ghost dofs indices
2438  for (int ss = 0; ss != loop_size; ++ss) {
2439  for (auto &gid : *ghost_idx_view[ss]) {
2440  NumeredDofEntityByUId::iterator dof =
2441  dof_by_uid_no_const[ss]->find(gid->getLocalUniqueId());
2442  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2443  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2444  "inconsistent data, ghost dof already set");
2445  bool success = dof_by_uid_no_const[ss]->modify(
2446  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2447  if (PetscUnlikely(!success))
2448  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2449  "modification unsuccessful");
2450  (*nb_ghost_dofs[ss])++;
2451  }
2452  }
2453  if (loop_size == 1) {
2454  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2455  }
2456  }
2457 
2458  if (verb > QUIET) {
2459  MOFEM_LOG("SYNC", Sev::inform)
2460  << " FEs ghost dofs on problem " << p_miit->getName()
2461  << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2462  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2463  << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2464 
2465  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2466  }
2467 
2470 }
2471 
2474  int verb) {
2475  MoFEM::Interface &m_field = cOre;
2476  auto problems_ptr = m_field.get_problems();
2478 
2480  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2481  "partition of problem not build");
2483  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2484  "partitions finite elements not build");
2485 
2486  // get problem pointer
2487  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2488  // find p_miit
2489  ProblemsByName &problems_set =
2490  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2491  ProblemsByName::iterator p_miit = problems_set.find(name);
2492 
2493  // get reference to number of ghost dofs
2494  // get number of local dofs
2495  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2496  &(p_miit->nbGhostDofsCol)};
2497  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2498  for (int ss = 0; ss != 2; ++ss) {
2499  (*nb_ghost_dofs[ss]) = 0;
2500  }
2501 
2502  // do work if more than one processor
2503  if (m_field.get_comm_size() > 1) {
2504  // determine if rows on columns are different from dofs on rows
2505  int loop_size = 2;
2506  if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2507  loop_size = 1;
2508  }
2509 
2510  typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2511  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2512  p_miit->numeredColDofsPtr};
2513 
2514  // iterate over dofs on rows and dofs on columns
2515  for (int ss = 0; ss != loop_size; ++ss) {
2516 
2517  // create dofs view by uid
2518  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2519 
2520  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2521  ghost_idx_view.reserve(std::distance(r.first, r.second));
2522  for (; r.first != r.second; ++r.first)
2523  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2524 
2525  auto cmp = [](auto a, auto b) {
2526  return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2527  };
2528  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2529 
2530  // iterate over dofs which have negative local index
2531  for (auto gid_it : ghost_idx_view) {
2532  bool success = numered_dofs[ss]->modify(
2533  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2534  if (!success)
2535  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2536  "modification unsuccessful");
2537  ++(*nb_ghost_dofs[ss]);
2538  }
2539  }
2540  if (loop_size == 1) {
2541  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2542  }
2543  }
2544 
2545  if (verb > QUIET) {
2546  MOFEM_LOG("SYNC", Sev::inform)
2547  << " FEs ghost dofs on problem " << p_miit->getName()
2548  << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2549  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2550  << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2551 
2552  MOFEM_LOG_SYNCHRONISE(m_field.get_comm())
2553  }
2554 
2557 }
2558 
2560  const std::string &fe_name,
2561  EntityHandle *meshset) const {
2562  MoFEM::Interface &m_field = cOre;
2563  const Problem *problem_ptr;
2564  const FiniteElement_multiIndex *fes_ptr;
2566 
2567  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2568  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2569  CHKERR m_field.get_finite_elements(&fes_ptr);
2570 
2571  auto fe_miit = fes_ptr->get<FiniteElement_name_mi_tag>().find(fe_name);
2572  if (fe_miit != fes_ptr->get<FiniteElement_name_mi_tag>().end()) {
2573  auto fit =
2574  problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
2576  0, (*fe_miit)->getFEUId()));
2577  auto hi_fe_it =
2578  problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
2580  get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
2581  std::vector<EntityHandle> fe_vec;
2582  fe_vec.reserve(std::distance(fit, hi_fe_it));
2583  for (; fit != hi_fe_it; fit++)
2584  fe_vec.push_back(fit->get()->getEnt());
2585  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2586  fe_vec.size());
2587  }
2588 
2590 }
2591 
2594  const std::string &fe_name,
2595  PetscLayout *layout) const {
2596  MoFEM::Interface &m_field = cOre;
2597  const Problem *problem_ptr;
2599  CHKERR m_field.get_problem(name, &problem_ptr);
2600  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2601  fe_name, layout);
2603 }
2604 
2606  const std::string problem_name, const std::string field_name,
2607  const Range ents, const int lo_coeff, const int hi_coeff,
2608  const int lo_order, const int hi_order, int verb, const bool debug) {
2609 
2610  MoFEM::Interface &m_field = cOre;
2612 
2613  const Problem *prb_ptr;
2614  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2615 
2616  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2617  prb_ptr->numeredRowDofsPtr, nullptr};
2618  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2619  numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2620 
2621  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2622  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2623  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2624 
2625  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2626  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2627  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2628  // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2629  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2630  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2631 
2632  for (int s = 0; s != 2; ++s)
2633  if (numered_dofs[s]) {
2634 
2635  typedef multi_index_container<
2636 
2637  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2638 
2639  >
2640  NumeredDofEntity_it_view_multiIndex;
2641 
2642  const auto bit_number = m_field.get_field_bit_number(field_name);
2643  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2644 
2645  // Set -1 to global and local dofs indices
2646  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2647  ++pit) {
2648  auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2649  DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2650  auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2651  DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2652 
2653  for (; lo != hi; ++lo)
2654  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2655  (*lo)->getDofCoeffIdx() <= hi_coeff &&
2656  (*lo)->getDofOrder() >= lo_order &&
2657  (*lo)->getDofOrder() <= hi_order)
2658  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2659  }
2660 
2661  if (verb > QUIET) {
2662  for (auto &dof : dofs_it_view)
2663  MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2664  MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
2665  }
2666 
2667  // create weak view
2668  std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
2669  dofs_weak_view.reserve(dofs_it_view.size());
2670  for (auto dit : dofs_it_view)
2671  dofs_weak_view.push_back(*dit);
2672 
2673  if (verb >= NOISY)
2674  MOFEM_LOG_C("SYNC", Sev::noisy,
2675  "Number of DOFs in multi-index %d and to delete %d\n",
2676  numered_dofs[s]->size(), dofs_it_view.size());
2677 
2678  // erase dofs from problem
2679  for (auto weak_dit : dofs_weak_view)
2680  if (auto dit = weak_dit.lock()) {
2681  numered_dofs[s]->erase(dit->getLocalUniqueId());
2682  }
2683 
2684  if (verb >= NOISY)
2685  MOFEM_LOG_C("SYNC", Sev::noisy,
2686  "Number of DOFs in multi-index after delete %d\n",
2687  numered_dofs[s]->size());
2688 
2689  // get current number of ghost dofs
2690  int nb_local_dofs = 0;
2691  int nb_ghost_dofs = 0;
2692  for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2693  dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2694  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2695  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2696  ++nb_local_dofs;
2697  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2698  ++nb_ghost_dofs;
2699  else
2700  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2701  "Impossible case. You could run problem on no distributed "
2702  "mesh. That is not implemented. Dof local index is %d",
2703  (*dit)->getPetscLocalDofIdx());
2704  }
2705 
2706  // get indices
2707  auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2708  const int nb_dofs = numered_dofs[s]->size();
2709  indices.clear();
2710  indices.reserve(nb_dofs);
2711  for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2712  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2713  bool add = true;
2714  if (only_local) {
2715  if ((*dit)->getPetscLocalDofIdx() < 0 ||
2716  (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2717  add = false;
2718  }
2719  }
2720  if (add)
2721  indices.push_back(decltype(tag)::get_index(dit));
2722  }
2723  };
2724 
2725  auto get_indices_by_uid = [&](auto tag, auto &indices) {
2726  const int nb_dofs = numered_dofs[s]->size();
2727  indices.clear();
2728  indices.reserve(nb_dofs);
2729  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2730  ++dit)
2731  indices.push_back(decltype(tag)::get_index(dit));
2732  };
2733 
2734  auto get_sub_ao = [&](auto sub_data) {
2735  if (s == 0) {
2736  return sub_data->getSmartRowMap();
2737  } else {
2738  return sub_data->getSmartColMap();
2739  }
2740  };
2741 
2742  auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
2743  if (s == 0) {
2744  sub_data->rowIs = is;
2745  sub_data->rowMap = ao;
2746  sub_data->colIs = is; // if square problem col is equal to row
2747  sub_data->colMap = ao;
2748  } else {
2749  sub_data->colIs = is;
2750  sub_data->colMap = ao;
2751  }
2752  };
2753 
2754  auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
2755  if (s == 0) {
2756  if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
2757  sub_data->colIs = sub_data->getSmartRowIs();
2758  sub_data->colMap = sub_data->getSmartRowMap();
2759  }
2760  }
2761  };
2762 
2763  auto concatenate_dofs = [&](auto tag, auto &indices,
2764  const auto local_only) {
2766  get_indices_by_tag(tag, indices, local_only);
2767 
2768  SmartPetscObj<AO> ao;
2769  // Create AO from app indices (i.e. old), to pestc indices (new after
2770  // remove)
2771  if (local_only)
2772  ao = createAOMapping(m_field.get_comm(), indices.size(),
2773  &*indices.begin(), PETSC_NULL);
2774  else
2775  ao = createAOMapping(PETSC_COMM_SELF, indices.size(),
2776  &*indices.begin(), PETSC_NULL);
2777 
2778  // Set mapping to sub dm data
2779  if (local_only) {
2780  if (auto sub_data = prb_ptr->getSubData()) {
2781  // create is and then map it to main problem of sub-problem
2782  auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2783  &*indices.begin(), PETSC_COPY_VALUES);
2784  // get old app, i.e. oroginal befor sub indices, and ao, from app,
2785  // to petsc sub indices.
2786  auto sub_ao = get_sub_ao(sub_data);
2787  CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
2788  sub_ao = createAOMappingIS(sub_is, PETSC_NULL);
2789  // set new sub ao
2790  set_sub_is_and_ao(sub_data, sub_is, sub_ao);
2791  apply_symmetry(sub_data);
2792  } else {
2793  // create sub data
2794  prb_ptr->getSubData() =
2795  boost::make_shared<Problem::SubProblemData>();
2796  auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2797  &*indices.begin(), PETSC_COPY_VALUES);
2798  // set sub is ao
2799  set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
2800  apply_symmetry(prb_ptr->getSubData());
2801  }
2802  }
2803 
2804  get_indices_by_uid(tag, indices);
2805  CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
2806 
2808  };
2809 
2810  // set indices index
2811  auto set_concatenated_indices = [&]() {
2812  std::vector<int> global_indices;
2813  std::vector<int> local_indices;
2815  CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
2816  CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
2817  auto gi = global_indices.begin();
2818  auto li = local_indices.begin();
2819  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2820  ++dit) {
2822  (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
2823  bool success = numered_dofs[s]->modify(dit, mod);
2824  if (!success)
2825  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2826  "can not set negative indices");
2827  ++gi;
2828  ++li;
2829  }
2831  };
2832  CHKERR set_concatenated_indices();
2833 
2834  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2835  m_field.get_comm());
2836  *(local_nbdof_ptr[s]) = nb_local_dofs;
2837  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
2838 
2839  if (debug)
2840  for (auto dof : (*numered_dofs[s])) {
2841  if (dof->getPetscGlobalDofIdx() < 0) {
2842  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2843  "Negative global idx");
2844  }
2845  if (dof->getPetscLocalDofIdx() < 0) {
2846  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2847  "Negative local idx");
2848  }
2849  }
2850 
2851  } else {
2852 
2853  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2854  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2855  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
2856  }
2857 
2858  if (verb > QUIET) {
2859  MOFEM_LOG_C(
2860  "WORLD", Sev::inform,
2861  "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
2862  prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
2863  prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
2864  MOFEM_LOG_C("SYNC", Sev::verbose,
2865  "Removed DOFs from problem %s dofs [ %d / %d "
2866  "(before %d / %d) local, %d / %d (before %d / %d)]",
2867  prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
2868  prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
2869  nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
2870  prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
2871  nb_init_ghost_col_dofs);
2872  MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2873  }
2874 
2876 }
2877 
2879  const std::string problem_name, const std::string field_name,
2880  const Range ents, const int lo_coeff, const int hi_coeff,
2881  const int lo_order, const int hi_order, int verb, const bool debug) {
2882 
2883  MoFEM::Interface &m_field = cOre;
2885 
2886  const Problem *prb_ptr;
2887  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2888 
2889  decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2890  prb_ptr->numeredRowDofsPtr, nullptr};
2891  if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2892  numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2893 
2894  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2895  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2896  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2897 
2898  const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2899  const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2900  const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2901  // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2902  const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2903  const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2904 
2905  const std::array<int, 2> nb_init_dofs = {nb_init_row_dofs, nb_init_col_dofs};
2906 
2907  for (int s = 0; s != 2; ++s)
2908  if (numered_dofs[s]) {
2909 
2910  typedef multi_index_container<
2911 
2912  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2913 
2914  >
2915  NumeredDofEntity_it_view_multiIndex;
2916 
2917  const auto bit_number = m_field.get_field_bit_number(field_name);
2918  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2919 
2920  // Set -1 to global and local dofs indices
2921  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2922  ++pit) {
2923  auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2924  DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2925  auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2926  DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2927 
2928  for (; lo != hi; ++lo)
2929  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2930  (*lo)->getDofCoeffIdx() <= hi_coeff &&
2931  (*lo)->getDofOrder() >= lo_order &&
2932  (*lo)->getDofOrder() <= hi_order)
2933  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2934  }
2935 
2936  if (verb > QUIET) {
2937  for (auto &dof : dofs_it_view)
2938  MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2939  }
2940 
2941  // set negative index
2942  auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2943  for (auto dit : dofs_it_view) {
2944  bool success = numered_dofs[s]->modify(dit, mod);
2945  if (!success)
2946  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2947  "can not set negative indices");
2948  }
2949 
2950  // create weak view
2951  std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2952  dosf_weak_view.reserve(dofs_it_view.size());
2953  for (auto dit : dofs_it_view)
2954  dosf_weak_view.push_back(*dit);
2955 
2956  if (verb >= NOISY)
2957  MOFEM_LOG_C("SYNC", Sev::noisy,
2958  "Number of DOFs in multi-index %d and to delete %d\n",
2959  numered_dofs[s]->size(), dofs_it_view.size());
2960 
2961  // erase dofs from problem
2962  for (auto weak_dit : dosf_weak_view)
2963  if (auto dit = weak_dit.lock()) {
2964  numered_dofs[s]->erase(dit->getLocalUniqueId());
2965  }
2966 
2967  if (verb >= NOISY)
2968  MOFEM_LOG_C("SYNC", Sev::noisy,
2969  "Number of DOFs in multi-index after delete %d\n",
2970  numered_dofs[s]->size());
2971 
2972  // get current number of ghost dofs
2973  int nb_global_dof = 0;
2974  int nb_local_dofs = 0;
2975  int nb_ghost_dofs = 0;
2976 
2977  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2978  ++dit) {
2979 
2980  if ((*dit)->getDofIdx() >= 0) {
2981 
2982  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2983  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2984  ++nb_local_dofs;
2985  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2986  ++nb_ghost_dofs;
2987 
2988  ++nb_global_dof;
2989  }
2990  }
2991 
2992  if (debug) {
2993  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2994  m_field.get_comm());
2995  if (*(nbdof_ptr[s]) != nb_global_dof)
2996  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2997  "Number of local DOFs do not add up %d != %d",
2998  *(nbdof_ptr[s]), nb_global_dof);
2999  }
3000 
3001  *(nbdof_ptr[s]) = nb_global_dof;
3002  *(local_nbdof_ptr[s]) = nb_local_dofs;
3003  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3004 
3005  // get indices
3006  auto get_indices_by_tag = [&](auto tag) {
3007  std::vector<int> indices;
3008  indices.resize(nb_init_dofs[s], -1);
3009  for (auto dit = numered_dofs[s]->get<Idx_mi_tag>().lower_bound(0);
3010  dit != numered_dofs[s]->get<Idx_mi_tag>().end(); ++dit) {
3011  indices[(*dit)->getDofIdx()] = decltype(tag)::get_index(dit);
3012  }
3013  return indices;
3014  };
3015 
3016  auto renumber = [&](auto tag, auto &indices) {
3018  int idx = 0;
3019  for (auto dit = numered_dofs[s]->get<decltype(tag)>().lower_bound(0);
3020  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3021  indices[(*dit)->getDofIdx()] = idx++;
3022  }
3024  };
3025 
3026  auto get_sub_ao = [&](auto sub_data) {
3027  if (s == 0) {
3028  return sub_data->getSmartRowMap();
3029  } else {
3030  return sub_data->getSmartColMap();
3031  }
3032  };
3033 
3034  auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
3035  if (s == 0) {
3036  sub_data->rowIs = is;
3037  sub_data->rowMap = ao;
3038  sub_data->colIs = is;
3039  sub_data->colMap = ao;
3040  } else {
3041  sub_data->colIs = is;
3042  sub_data->colMap = ao;
3043  }
3044  };
3045 
3046  auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
3047  if (s == 0) {
3048  if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
3049  sub_data->colIs = sub_data->getSmartRowIs();
3050  sub_data->colMap = sub_data->getSmartRowMap();
3051  }
3052  }
3053  };
3054 
3055  auto set_sub_data = [&](auto &indices) {
3057  if (auto sub_data = prb_ptr->getSubData()) {
3058  // create is and then map it to main problem of sub-problem
3059  auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3060  &*indices.begin(), PETSC_COPY_VALUES);
3061  // get old app, i.e. oroginal befor sub indices, and ao, from
3062  // app, to petsc sub indices.
3063  auto sub_ao = get_sub_ao(sub_data);
3064  CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
3065  sub_ao = createAOMappingIS(sub_is, PETSC_NULL);
3066  // set new sub ao
3067  set_sub_is_and_ao(sub_data, sub_is, sub_ao);
3068  apply_symmetry(sub_data);
3069  } else {
3070  prb_ptr->getSubData() = boost::make_shared<Problem::SubProblemData>();
3071  auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3072  &*indices.begin(), PETSC_COPY_VALUES);
3073  auto sub_ao = createAOMappingIS(sub_is, PETSC_NULL);
3074  // set sub is ao
3075  set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, sub_ao);
3076  apply_symmetry(prb_ptr->getSubData());
3077  }
3079  };
3080 
3081  auto global_indices = get_indices_by_tag(PetscGlobalIdx_mi_tag());
3082  auto local_indices = get_indices_by_tag(PetscLocalIdx_mi_tag());
3083  CHKERR set_sub_data(global_indices);
3084  CHKERR renumber(PetscGlobalIdx_mi_tag(), global_indices);
3085  CHKERR renumber(PetscLocalIdx_mi_tag(), local_indices);
3086 
3087  int i = 0;
3088  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3089  ++dit) {
3090  auto idx = (*dit)->getDofIdx();
3091  if (idx >= 0) {
3093  (*dit)->getPart(), i++, global_indices[idx], local_indices[idx]);
3094  bool success = numered_dofs[s]->modify(dit, mod);
3095  if (!success)
3096  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3097  "can not set negative indices");
3098  } else {
3100  (*dit)->getPart(), -1, -1, -1);
3101  bool success = numered_dofs[s]->modify(dit, mod);
3102  if (!success)
3103  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3104  "can not set negative indices");
3105  }
3106  };
3107 
3108  if (debug) {
3109  for (auto dof : (*numered_dofs[s])) {
3110  if (dof->getDofIdx() >= 0 && dof->getPetscGlobalDofIdx() < 0) {
3111  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3112  "Negative global idx");
3113  }
3114  }
3115  }
3116 
3117  } else {
3118 
3119  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3120  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3121  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3122  }
3123 
3124  if (verb >= NOISY)
3125  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
3126 
3127  if (verb > QUIET) {
3128  MOFEM_LOG_C(
3129  "WORLD", Sev::inform,
3130  "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
3131  prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
3132  prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
3133  MOFEM_LOG_C("SYNC", Sev::verbose,
3134  "Removed DOFs from problem %s dofs [ %d / %d "
3135  "(before %d / %d) local, %d / %d (before %d / %d)]",
3136  prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
3137  prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
3138  nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
3139  prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
3140  nb_init_ghost_col_dofs);
3141  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
3142  }
3144 }
3145 
3147  const std::string problem_name, const std::string field_name,
3148  const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
3149  Range *ents_ptr, const int lo_coeff, const int hi_coeff, const int lo_order,
3150  const int hi_order, int verb, const bool debug) {
3151  MoFEM::Interface &m_field = cOre;
3153 
3154  auto bit_manager = m_field.getInterface<BitRefManager>();
3155 
3156  Range ents;
3157  if (ents_ptr) {
3158  ents = *ents_ptr;
3159  CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3160  ents, verb);
3161  } else {
3162  CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3163  verb);
3164  }
3165 
3166  CHKERR removeDofsOnEntities(problem_name, field_name, ents, lo_coeff,
3167  hi_coeff, lo_order, hi_order, verb, debug);
3168 
3170 }
3171 
3173  const std::string problem_name, const std::string field_name,
3174  const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
3175  Range *ents_ptr, const int lo_coeff, const int hi_coeff, const int lo_order,
3176  const int hi_order, int verb, const bool debug) {
3177  MoFEM::Interface &m_field = cOre;
3179 
3180  auto bit_manager = m_field.getInterface<BitRefManager>();
3181 
3182  Range ents;
3183  if (ents_ptr) {
3184  ents = *ents_ptr;
3185  CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3186  ents, verb);
3187  } else {
3188  CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3189  verb);
3190  }
3191 
3193  lo_coeff, hi_coeff, lo_order,
3194  hi_order, verb, debug);
3195 
3197 }
3198 
3200 ProblemsManager::markDofs(const std::string problem_name, RowColData rc,
3201  const enum MarkOP op, const Range ents,
3202  std::vector<unsigned char> &marker) const {
3203 
3204  Interface &m_field = cOre;
3205  const Problem *problem_ptr;
3207  CHKERR m_field.get_problem(problem_name, &problem_ptr);
3208  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3209  switch (rc) {
3210  case ROW:
3211  dofs = problem_ptr->getNumeredRowDofsPtr();
3212  break;
3213  case COL:
3214  dofs = problem_ptr->getNumeredColDofsPtr();
3215  default:
3216  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
3217  }
3218  marker.resize(dofs->size(), 0);
3219  std::vector<unsigned char> marker_tmp;
3220 
3221  switch (op) {
3222  case MarkOP::OR:
3223  for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
3224  auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
3225  auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
3226  for (; lo != hi; ++lo)
3227  marker[(*lo)->getPetscLocalDofIdx()] |= 1;
3228  }
3229  break;
3230  case MarkOP::AND:
3231  marker_tmp.resize(dofs->size(), 0);
3232  for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
3233  auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
3234  auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
3235  for (; lo != hi; ++lo)
3236  marker_tmp[(*lo)->getPetscLocalDofIdx()] = 1;
3237  }
3238  for (int i = 0; i != marker.size(); ++i) {
3239  marker[i] &= marker_tmp[i];
3240  }
3241  break;
3242  }
3243 
3245 }
3246 
3248  const std::string problem_name, RowColData rc, const std::string field_name,
3249  const int lo, const int hi, const enum ProblemsManager::MarkOP op,
3250  const unsigned char c, std::vector<unsigned char> &marker) const {
3251 
3252  Interface &m_field = cOre;
3253  const Problem *problem_ptr;
3255  CHKERR m_field.get_problem(problem_name, &problem_ptr);
3256  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3257  switch (rc) {
3258  case ROW:
3259  dofs = problem_ptr->getNumeredRowDofsPtr();
3260  break;
3261  case COL:
3262  dofs = problem_ptr->getNumeredColDofsPtr();
3263  default:
3264  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
3265  }
3266  marker.resize(dofs->size(), 0);
3267 
3268  auto dof_lo = dofs->get<Unique_mi_tag>().lower_bound(
3270  auto dof_hi = dofs->get<Unique_mi_tag>().upper_bound(
3272 
3273  auto marker_ref = [marker](auto &it) -> unsigned int & {
3274  return marker[(*it)->getPetscLocalDofIdx()];
3275  };
3276 
3277  switch (op) {
3278  case MarkOP::OR:
3279  for (; dof_lo != dof_hi; ++dof_lo)
3280  if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3281  (*dof_lo)->getDofCoeffIdx() <= hi)
3282  marker[(*dof_lo)->getPetscLocalDofIdx()] |= c;
3283  break;
3284  case MarkOP::AND:
3285  for (; dof_lo != dof_hi; ++dof_lo)
3286  if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3287  (*dof_lo)->getDofCoeffIdx() <= hi)
3288  marker[(*dof_lo)->getPetscLocalDofIdx()] &= c;
3289  break;
3290  }
3291 
3293 }
3294 
3296 ProblemsManager::addFieldToEmptyFieldBlocks(const std::string problem_name,
3297  const std::string row_field,
3298  const std::string col_field) const {
3299 
3300  Interface &m_field = cOre;
3302 
3303  const auto problem_ptr = m_field.get_problem(problem_name);
3304  auto get_field_id = [&](const std::string field_name) {
3305  return m_field.get_field_structure(field_name)->getId();
3306  };
3307  const auto row_id = get_field_id(row_field);
3308  const auto col_id = get_field_id(col_field);
3309 
3310  problem_ptr->addFieldToEmptyFieldBlocks(BlockFieldPair(row_id, col_id));
3311 
3313 }
3314 
3315 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
MoFEM::ProblemsManager::buildProblemFromFields
PetscBool buildProblemFromFields
Definition: ProblemsManager.hpp:30
MoFEM::IdxDataTypePtr::getUId
UId getUId() const
Definition: ProblemsManager.cpp:32
MoFEM::CoreInterface::get_field_name
virtual std::string get_field_name(const BitFieldId id) const =0
get field name from id
MoFEM::FieldEntity::getGlobalUniqueIdCalculate
UId getGlobalUniqueIdCalculate() const
Calculate global UId.
Definition: FieldEntsMultiIndices.hpp:216
MoFEM::CoreInterface::get_ents_finite_elements
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:21
MoFEM::Problem::nbGhostDofsRow
DofIdx nbGhostDofsRow
Number of ghost DOFs in row.
Definition: ProblemsMultiIndices.hpp:69
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
MoFEM::ProblemsManager::MarkOP
MarkOP
Definition: ProblemsManager.hpp:365
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
ProblemManagerFunctionBegin
#define ProblemManagerFunctionBegin
Definition: ProblemsManager.cpp:8
MoFEM::NumeredDofEntity_part_and_mofem_glob_idx_change
Definition: DofsMultiIndices.hpp:549
MoFEM::field_bit_from_bit_number
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1908
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
MoFEM::createISGeneral
auto createISGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
Creates a data structure for an index set containing a list of integers.
Definition: PetscSmartObj.hpp:287
MoFEM::Problem::getNumberOfElementsByNameAndPart
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
Definition: ProblemsMultiIndices.cpp:98
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::Problem::getNbDofsCol
DofIdx getNbDofsCol() const
Definition: ProblemsMultiIndices.hpp:377
NOISY
@ NOISY
Definition: definitions.h:211
MoFEM::IdxDataTypePtr::getDofIdx
int getDofIdx() const
Definition: ProblemsManager.cpp:28
MoFEM::Problem_mi_tag
Definition: TagMultiIndices.hpp:70
MoFEM::ProblemsManager::partitionMesh
DEPRECATED MoFEMErrorCode partitionMesh(const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
Set partition tag to each finite element in the problem.
Definition: ProblemsManager.cpp:77
MoFEM::ProblemsManager::inheritPartition
MoFEMErrorCode 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
Definition: ProblemsManager.cpp:1896
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::NumeredDofEntity_part_and_indices_change
Definition: DofsMultiIndices.hpp:564
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::NumeredDofEntity_part_and_glob_idx_change
Definition: DofsMultiIndices.hpp:537
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::PetscLocalIdx_mi_tag
Definition: TagMultiIndices.hpp:45
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::NumeredDofEntity_part_and_all_indices_change
Definition: DofsMultiIndices.hpp:610
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEM::Composite_Name_And_Part_mi_tag
Definition: TagMultiIndices.hpp:87
MoFEM::Problem::getBitFEId
BitFEId getBitFEId() const
Definition: ProblemsMultiIndices.cpp:45
MoFEM::CoreTmp< 0 >::BUILD_ADJ
@ BUILD_ADJ
Definition: Core.hpp:214
MoFEM::Problem::getRowDofsSequence
auto & getRowDofsSequence() const
Get reference to sequence data numbered dof container.
Definition: ProblemsMultiIndices.hpp:506
MoFEM::CoreInterface::get_finite_elements
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
MoFEM::Problem::getNbGhostDofsCol
DofIdx getNbGhostDofsCol() const
Definition: ProblemsMultiIndices.hpp:381
MoFEM::Problem::getNbGhostDofsRow
DofIdx getNbGhostDofsRow() const
Definition: ProblemsMultiIndices.hpp:380
MoFEM::ProblemsManager::printPartitionedProblem
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
Definition: ProblemsManager.cpp:2065
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MoFEM::ProblemsManager::debugPartitionedProblem
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
Definition: ProblemsManager.cpp:2084
MoFEM::ProblemsManager::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: ProblemsManager.cpp:49
MoFEM::DofEntity_multiIndex_active_view
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getLocalUniqueId > >, ordered_non_unique< const_mem_fun< DofEntity, char, &DofEntity::getActive > > > > DofEntity_multiIndex_active_view
multi-index view on DofEntity activity
Definition: DofsMultiIndices.hpp:370
sdf.r
int r
Definition: sdf.py:8
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::CoreInterface::get_dofs
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
NumeredDofEntity_multiIndex
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
Definition: DofsMultiIndices.hpp:469
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::NumeredDofEntityByUId
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.
Definition: DofsMultiIndices.hpp:476
MoFEM::IdxDataTypePtr::pTr
const int * pTr
Definition: ProblemsManager.cpp:45
MoFEM::Field::getId
const BitFieldId & getId() const
Get unique field id.
Definition: FieldMultiIndices.hpp:128
VERBOSE
@ VERBOSE
Definition: definitions.h:209
MoFEM::createAOMapping
auto createAOMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[])
Creates an application mapping using two integer arrays.
Definition: PetscSmartObj.hpp:338
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::ProblemsManager::cOre
MoFEM::Core & cOre
Definition: ProblemsManager.hpp:26
MoFEM::ProblemsManager::removeDofsOnEntitiesNotDistributed
MoFEMErrorCode removeDofsOnEntitiesNotDistributed(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, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2878
MoFEM::NumeredDofEntity
keeps information about indexed dofs for the problem
Definition: DofsMultiIndices.hpp:226
MoFEM::NumeredDofEntity_local_idx_change
Definition: DofsMultiIndices.hpp:584
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::Part_mi_tag
Definition: TagMultiIndices.hpp:53
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh
MoFEMErrorCode partitionGhostDofsOnDistributedMesh(const std::string name, int verb=VERBOSE)
determine ghost nodes on distributed meshes
Definition: ProblemsManager.cpp:2473
MoFEM::Problem::numeredColDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
Definition: ProblemsMultiIndices.hpp:75
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::FieldEntity::getHandleFromUniqueId
static auto getHandleFromUniqueId(const UId uid)
Get the Handle From Unique Id.
Definition: FieldEntsMultiIndices.hpp:192
MoFEM::Types::UId
uint128_t UId
Unique Id.
Definition: Types.hpp:31
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1864
RowColData
RowColData
RowColData.
Definition: definitions.h:123
MoFEM::ProblemsManager::removeDofsOnEntities
MoFEMErrorCode 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, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2605
convert.type
type
Definition: convert.py:64
MoFEM::Problem::getNbDofsRow
DofIdx getNbDofsRow() const
Definition: ProblemsMultiIndices.hpp:376
MoFEM::FieldEntity::getLocalUniqueIdCalculate
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Definition: FieldEntsMultiIndices.hpp:136
MoFEM::EntFiniteElement::getLocalUniqueIdCalculate
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
Definition: FEMultiIndices.hpp:528
MoFEM::CoreTmp< 0 >::BUILD_FIELD
@ BUILD_FIELD
Definition: Core.hpp:212
MoFEM::IdxDataType
Definition: ProblemsManager.cpp:16
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition: ProblemsMultiIndices.hpp:119
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
COL
@ COL
Definition: definitions.h:123
MoFEM::Problem::nbDofsCol
DofIdx nbDofsCol
Global number of DOFs in col.
Definition: ProblemsMultiIndices.hpp:66
MoFEM::ProblemsManager::addFieldToEmptyFieldBlocks
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string problem_name, const std::string row_field, const std::string col_field) const
add empty block to problem
Definition: ProblemsManager.cpp:3296
MoFEM::NumeredDofEntity_multiIndex_uid_view_ordered
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
Definition: DofsMultiIndices.hpp:503
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::ProblemsManager::getOptions
MoFEMErrorCode getOptions()
Definition: ProblemsManager.cpp:62
MoFEM::ProblemsManager::markDofs
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Definition: ProblemsManager.cpp:3200
MoFEM::CoreInterface::get_field_ents
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
MoFEM::NumeredDofEntity_mofem_index_change
Definition: DofsMultiIndices.hpp:597
MoFEM::CoreTmp< 0 >::BUILD_PROBLEM
@ BUILD_PROBLEM
Definition: Core.hpp:215
MoFEM::Problem::getColDofsSequence
auto & getColDofsSequence() const
Get reference to sequence data numbered dof container.
Definition: ProblemsMultiIndices.hpp:522
MoFEM::BlockFieldPair
Problem::BlockFieldPair BlockFieldPair
Definition: ProblemsMultiIndices.hpp:562
MoFEM::Problem::nbLocDofsRow
DofIdx nbLocDofsRow
Local number of DOFs in row.
Definition: ProblemsMultiIndices.hpp:67
MoFEM::IdxDataTypePtr
Definition: ProblemsManager.cpp:26
MoFEM::CoreTmp< 0 >::PARTITION_PROBLEM
@ PARTITION_PROBLEM
Definition: Core.hpp:216
MoFEM::ProblemsManager::getProblemElementsLayout
MoFEMErrorCode getProblemElementsLayout(const std::string name, const std::string &fe_name, PetscLayout *layout) const
Get layout of elements in the problem.
Definition: ProblemsManager.cpp:2593
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1876
MoFEM::PetscGlobalIdx_mi_tag
Definition: TagMultiIndices.hpp:38
MoFEM::FieldEntity_multiIndex_global_uid_view
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FieldEntity, UId, &FieldEntity::getGlobalUniqueId > > > > FieldEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
Definition: FieldEntsMultiIndices.hpp:444
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:378
MoFEM::CoreInterface::get_problems
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::ProblemsManager::partitionFiniteElements
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
Definition: ProblemsManager.cpp:2167
MoFEM::Problem::numeredFiniteElementsPtr
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
Definition: ProblemsMultiIndices.hpp:77
MoFEM::DofEntity_multiIndex_global_uid_view
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getGlobalUniqueId > > > > DofEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
Definition: DofsMultiIndices.hpp:349
MoFEM::DofEntity::getLoFieldEntityUId
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:69
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
NumeredEntFiniteElement_multiIndex
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >> >> NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
Definition: FEMultiIndices.hpp:830
MoFEM::IdxDataTypePtr::IdxDataTypePtr
IdxDataTypePtr(const int *ptr)
Definition: ProblemsManager.cpp:27
MoFEM::ProblemsManager::buildComposedProblem
MoFEMErrorCode buildComposedProblem(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
Definition: ProblemsManager.cpp:1250
MoFEM::Problem::nbGhostDofsCol
DofIdx nbGhostDofsCol
Number of ghost DOFs in col.
Definition: ProblemsMultiIndices.hpp:70
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
convert.n
n
Definition: convert.py:82
MoFEM::CoreInterface::get_fields
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::CoreTmp< 0 >::PARTITION_FE
@ PARTITION_FE
Definition: Core.hpp:217
Range
MoFEM::CoreTmp< 0 >::getBuildMoFEM
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition: Core.hpp:225
MoFEM::IdxDataType::IdxDataType
IdxDataType(const UId uid, const int dof)
Definition: ProblemsManager.cpp:17
MoFEM::Problem::nbLocDofsCol
DofIdx nbLocDofsCol
Local number of DOFs in colIs.
Definition: ProblemsMultiIndices.hpp:68
MoFEM::ProblemsManager::MOFEM_EVENT_ProblemsManager
PetscLogEvent MOFEM_EVENT_ProblemsManager
Definition: ProblemsManager.hpp:427
FTensor::dd
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::IdxDataType::dAta
int dAta[5]
Definition: ProblemsManager.cpp:23
MoFEM::EntFiniteElement::getDofView
static MoFEMErrorCode getDofView(const FE_ENTS &fe_ents_view, const MOFEM_DOFS &mofem_dofs, MOFEM_DOFS_VIEW &dofs_view, INSERTER &&inserter)
Definition: FEMultiIndices.hpp:583
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::DofEntity::getHiFieldEntityUId
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:75
MoFEM::ProblemsManager::partitionSimpleProblem
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
Definition: ProblemsManager.cpp:1537
MoFEM::ProblemsManager::synchroniseProblemEntities
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
Definition: ProblemsManager.hpp:33
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::FiniteElement_name_mi_tag
Definition: TagMultiIndices.hpp:26
MoFEM::CoreInterface::get_ents_elements_adjacency
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::Problem::nbDofsRow
DofIdx nbDofsRow
Global number of DOFs in row.
Definition: ProblemsMultiIndices.hpp:65
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ProblemsManager::buildProblemOnDistributedMesh
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
Definition: ProblemsManager.cpp:290
MoFEM::Idx_mi_tag
Definition: TagMultiIndices.hpp:31
MoFEM::CoreInterface::clear_problem
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
MoFEM::Problem::getBitRefLevelMask
BitRefLevel getBitRefLevelMask() const
Definition: ProblemsMultiIndices.hpp:384
MoFEM::ProblemsManager::ProblemsManager
ProblemsManager(const MoFEM::Core &core)
Definition: ProblemsManager.cpp:55
MoFEM::Problem::getBitRefLevel
BitRefLevel getBitRefLevel() const
Definition: ProblemsMultiIndices.hpp:383
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
MoFEM::Problem::getNbLocalDofsCol
DofIdx getNbLocalDofsCol() const
Definition: ProblemsMultiIndices.hpp:379
MoFEM::CoreTmp< 0 >::BUILD_FE
@ BUILD_FE
Definition: Core.hpp:213
MoFEM::Problem::getNumeredColDofsPtr
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
Definition: ProblemsMultiIndices.hpp:87
FiniteElement_multiIndex
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.
Definition: FEMultiIndices.hpp:849
MoFEM::NumeredEntFiniteElement
Partitioned (Indexed) Finite Element in Problem.
Definition: FEMultiIndices.hpp:728
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::ProblemsManager::getFEMeshset
MoFEMErrorCode getFEMeshset(const std::string prb_name, const std::string &fe_name, EntityHandle *meshset) const
create add entities of finite element in the problem
Definition: ProblemsManager.cpp:2559
marker
auto marker
set bit to marker
Definition: hanging_node_approx.cpp:82
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::CoreTmp< 0 >::PARTITION_GHOST_DOFS
@ PARTITION_GHOST_DOFS
Definition: Core.hpp:218
MoFEM::SmartPetscObj< IS >
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
MoFEM::ProblemsManager::modifyMarkDofs
MoFEMErrorCode modifyMarkDofs(const std::string problem_name, RowColData rc, const std::string field_name, const int lo, const int hi, const enum MarkOP op, const unsigned char c, std::vector< unsigned char > &marker) const
Mark DOFs.
Definition: ProblemsManager.cpp:3247
MoFEM::FieldEntity::getFieldBitNumberFromUniqueId
static auto getFieldBitNumberFromUniqueId(const UId uid)
Get the Field Bit Number From Unique Id.
Definition: FieldEntsMultiIndices.hpp:205
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
VERY_VERBOSE
@ VERY_VERBOSE
Definition: definitions.h:210
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
MoFEM::ProblemsManager::buildSubProblem
MoFEMErrorCode 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, boost::shared_ptr< Range >> *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range >> *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
Definition: ProblemsManager.cpp:981
MoFEM::Types::DofIdx
int DofIdx
Index of DOF.
Definition: Types.hpp:18
MoFEM::FieldEntity::getOwnerFromUniqueId
static auto getOwnerFromUniqueId(const UId uid)
Definition: FieldEntsMultiIndices.hpp:177