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