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