v0.9.0
ProblemsManager.cpp
Go to the documentation of this file.
1 /** \file ProblemsManager.cpp
2  * \brief Managing complexities for problem
3  * \ingroup mofem_problems_manager
4  */
5 
6 /* 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.getInterface<CommInterface>()->synchroniseEntities(
815  ents_to_synchronise, verb);
816  ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
817  for (Field_multiIndex::iterator fit = fields_ptr->begin();
818  fit != fields_ptr->end(); fit++) {
819  if ((fit->get()->getId() & *fields_ids[ss]).any()) {
820  for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
821  pit != ents_to_synchronise.pair_end(); ++pit) {
822  const EntityHandle f = pit->first;
823  const EntityHandle s = pit->second;
824  DofEntity_multiIndex::index<
825  Composite_Name_And_Ent_mi_tag>::type::iterator dit,
826  hi_dit;
827  dit = dofs_field_ptr->get<Composite_Name_And_Ent_mi_tag>()
828  .lower_bound(boost::make_tuple(fit->get()->getName(), f));
829  hi_dit =
830  dofs_field_ptr->get<Composite_Name_And_Ent_mi_tag>()
831  .upper_bound(boost::make_tuple(fit->get()->getName(), s));
832  dofs_ptr[ss]->insert(dit, hi_dit);
833  }
834  }
835  }
836  }
837  }
838 
839  // add dofs for rows and cols and set ownership
840  DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
841  boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
842  problem_ptr->numeredDofsRows, problem_ptr->numeredDofsCols};
843  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
844  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
845  &problem_ptr->nbLocDofsCol};
846  int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
847  &problem_ptr->nbGhostDofsCol};
848  for (int ss = 0; ss < 2; ss++) {
849  *(nbdof_ptr[ss]) = 0;
850  *(local_nbdof_ptr[ss]) = 0;
851  *(ghost_nbdof_ptr[ss]) = 0;
852  }
853 
854  // Loop over dofs on rows and columns and add to multi-indices in dofs problem
855  // structure, set partition for each dof
856  int nb_local_dofs[] = {0, 0};
857  for (int ss = 0; ss < loop_size; ss++) {
858  DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
859  hi_miit;
860  miit = dofs_ptr[ss]->get<1>().lower_bound(1);
861  hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
862  for (; miit != hi_miit; miit++) {
863  int owner_proc = (*miit)->getOwnerProc();
864  if (owner_proc == m_field.get_comm_rank()) {
865  nb_local_dofs[ss]++;
866  }
867  }
868  }
869 
870  // get layout
871  int start_ranges[2], end_ranges[2];
872  for (int ss = 0; ss != loop_size; ss++) {
873  PetscLayout layout;
874  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
875  CHKERR PetscLayoutSetBlockSize(layout, 1);
876  CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
877  CHKERR PetscLayoutSetUp(layout);
878  CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
879  CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
880  CHKERR PetscLayoutDestroy(&layout);
881  }
882  if (square_matrix) {
883  nbdof_ptr[1] = nbdof_ptr[0];
884  nb_local_dofs[1] = nb_local_dofs[0];
885  start_ranges[1] = start_ranges[0];
886  end_ranges[1] = end_ranges[0];
887  }
888 
889  // if(sizeof(UId) != SIZEOFUID) {
890  // SETERRQ2(
891  // PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,
892  // "check size of UId, size of UId is %u != %u",
893  // sizeof(UId),SIZEOFUID
894  // );
895  // }
896 
897  // set local and global indices on own dofs
898  const size_t idx_data_type_size = sizeof(IdxDataType);
899  const size_t data_block_size = idx_data_type_size / sizeof(int);
900 
901  if (sizeof(IdxDataType) % sizeof(int)) {
902  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
903  }
904  std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
905  m_field.get_comm_size()),
906  ids_data_packed_cols(m_field.get_comm_size());
907 
908  // Loop over dofs on this processor and prepare those dofs to send on another
909  // proc
910  for (int ss = 0; ss != loop_size; ++ss) {
911 
912  DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
913  hi_miit;
914  hi_miit = dofs_ptr[ss]->get<0>().end();
915 
916  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
917  boost::shared_ptr<std::vector<NumeredDofEntity>>(
918  new std::vector<NumeredDofEntity>());
919  int nb_dofs_to_add = 0;
920  miit = dofs_ptr[ss]->get<0>().begin();
921  for (; miit != hi_miit; ++miit) {
922  // Only set global idx for dofs on this processor part
923  if (!(miit->get()->getActive()))
924  continue;
925  ++nb_dofs_to_add;
926  }
927  dofs_array->reserve(nb_dofs_to_add);
928  if (ss == 0) {
929  problem_ptr->getRowDofsSequence()->push_back(dofs_array);
930  } else {
931  problem_ptr->getColDofsSequence()->push_back(dofs_array);
932  }
933 
934  int &local_idx = *local_nbdof_ptr[ss];
935  miit = dofs_ptr[ss]->get<0>().begin();
936  for (; miit != hi_miit; ++miit) {
937 
938  // Only set global idx for dofs on this processor part
939  if (!(miit->get()->getActive()))
940  continue;
941 
942  dofs_array->emplace_back(*miit);
943 
944  int owner_proc = dofs_array->back().getOwnerProc();
945  if (owner_proc < 0) {
946  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
947  "data inconsistency");
948  }
949 
950  if (owner_proc != m_field.get_comm_rank()) {
951  dofs_array->back().pArt = owner_proc;
952  dofs_array->back().dofIdx = -1;
953  dofs_array->back().petscGloablDofIdx = -1;
954  dofs_array->back().petscLocalDofIdx = -1;
955  } else {
956 
957  // Set part and indexes
958  int glob_idx = start_ranges[ss] + local_idx;
959  dofs_array->back().pArt = owner_proc;
960  dofs_array->back().dofIdx = glob_idx;
961  dofs_array->back().petscGloablDofIdx = glob_idx;
962  dofs_array->back().petscLocalDofIdx = local_idx;
963  local_idx++;
964 
965  unsigned char pstatus = dofs_array->back().getPStatus();
966  // check id dof is shared, if that is a case global idx added to data
967  // structure and is sended to other processors
968  if (pstatus > 0) {
969 
970  for (int proc = 0;
971  proc < MAX_SHARING_PROCS &&
972  -1 != dofs_array->back().getSharingProcsPtr()[proc];
973  proc++) {
974  // make it different for rows and columns when needed
975  if (ss == 0) {
976  ids_data_packed_rows[dofs_array->back()
977  .getSharingProcsPtr()[proc]]
978  .emplace_back(dofs_array->back().getGlobalUniqueId(),
979  glob_idx);
980  } else {
981  ids_data_packed_cols[dofs_array->back()
982  .getSharingProcsPtr()[proc]]
983  .emplace_back(dofs_array->back().getGlobalUniqueId(),
984  glob_idx);
985  }
986  if (!(pstatus & PSTATUS_MULTISHARED)) {
987  break;
988  }
989  }
990  }
991  }
992  }
993 
994  auto hint = numered_dofs_ptr[ss]->end();
995  for (auto &v : *dofs_array)
996  hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
997  }
998  if (square_matrix) {
999  local_nbdof_ptr[1] = local_nbdof_ptr[0];
1000  }
1001 
1002  int nsends_rows = 0, nsends_cols = 0;
1003  // Non zero lengths[i] represent a message to i of length lengths[i].
1004  std::vector<int> lengths_rows(m_field.get_comm_size()),
1005  lengths_cols(m_field.get_comm_size());
1006  lengths_rows.clear();
1007  lengths_cols.clear();
1008  for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
1009  lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
1010  lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
1011  if (!ids_data_packed_rows[proc].empty())
1012  nsends_rows++;
1013  if (!ids_data_packed_cols[proc].empty())
1014  nsends_cols++;
1015  }
1016 
1017  MPI_Status *status;
1018  CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
1019 
1020  // Do rows
1021  int nrecvs_rows; // number of messages received
1022  int *onodes_rows; // list of node-ids from which messages are expected
1023  int *olengths_rows; // corresponding message lengths
1024  int **rbuf_row; // must bee freed by user
1025 
1026  {
1027 
1028  // make sure it is a PETSc comm
1029  CHKERR PetscCommDuplicate(m_field.get_comm(), &m_field.get_comm(), NULL);
1030 
1031  // rows
1032 
1033  // Computes the number of messages a node expects to receive
1034  CHKERR PetscGatherNumberOfMessages(m_field.get_comm(), NULL,
1035  &lengths_rows[0], &nrecvs_rows);
1036  // std::cerr << nrecvs_rows << std::endl;
1037 
1038  // Computes info about messages that a MPI-node will receive, including
1039  // (from-id,length) pairs for each message.
1040  CHKERR PetscGatherMessageLengths(m_field.get_comm(), nsends_rows,
1041  nrecvs_rows, &lengths_rows[0],
1042  &onodes_rows, &olengths_rows);
1043 
1044  // Gets a unique new tag from a PETSc communicator. All processors that
1045  // share the communicator MUST call this routine EXACTLY the same number of
1046  // times. This tag should only be used with the current objects
1047  // communicator; do NOT use it with any other MPI communicator.
1048  int tag_row;
1049  CHKERR PetscCommGetNewTag(m_field.get_comm(), &tag_row);
1050 
1051  // Allocate a buffer sufficient to hold messages of size specified in
1052  // olengths. And post Irecvs on these buffers using node info from onodes
1053  MPI_Request *r_waits_row; // must bee freed by user
1054  // rbuf has a pointers to messeges. It has size of of nrecvs (number of
1055  // messages) +1. In the first index a block is allocated,
1056  // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
1057 
1058  CHKERR PetscPostIrecvInt(m_field.get_comm(), tag_row, nrecvs_rows,
1059  onodes_rows, olengths_rows, &rbuf_row,
1060  &r_waits_row);
1061  CHKERR PetscFree(onodes_rows);
1062 
1063  MPI_Request *s_waits_row; // status of sens messages
1064  CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
1065 
1066  // Send messeges
1067  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
1068  if (!lengths_rows[proc])
1069  continue; // no message to send to this proc
1070  CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
1071  lengths_rows[proc], // message length
1072  MPIU_INT, proc, // to proc
1073  tag_row, m_field.get_comm(), s_waits_row + kk);
1074  kk++;
1075  }
1076 
1077  if (nrecvs_rows) {
1078  CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
1079  }
1080  if (nsends_rows) {
1081  CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
1082  }
1083 
1084  CHKERR PetscFree(r_waits_row);
1085  CHKERR PetscFree(s_waits_row);
1086  }
1087 
1088  // cols
1089  int nrecvs_cols = nrecvs_rows;
1090  int *olengths_cols = olengths_rows;
1091  PetscInt **rbuf_col = rbuf_row;
1092  if (!square_matrix) {
1093 
1094  // Computes the number of messages a node expects to receive
1095  CHKERR PetscGatherNumberOfMessages(m_field.get_comm(), NULL,
1096  &lengths_cols[0], &nrecvs_cols);
1097 
1098  // Computes info about messages that a MPI-node will receive, including
1099  // (from-id,length) pairs for each message.
1100  int *onodes_cols;
1101  CHKERR PetscGatherMessageLengths(m_field.get_comm(), nsends_cols,
1102  nrecvs_cols, &lengths_cols[0],
1103  &onodes_cols, &olengths_cols);
1104 
1105  // Gets a unique new tag from a PETSc communicator.
1106  int tag_col;
1107  CHKERR PetscCommGetNewTag(m_field.get_comm(), &tag_col);
1108 
1109  // Allocate a buffer sufficient to hold messages of size specified in
1110  // olengths. And post Irecvs on these buffers using node info from onodes
1111  MPI_Request *r_waits_col; // must bee freed by user
1112  CHKERR PetscPostIrecvInt(m_field.get_comm(), tag_col, nrecvs_cols,
1113  onodes_cols, olengths_cols, &rbuf_col,
1114  &r_waits_col);
1115  CHKERR PetscFree(onodes_cols);
1116 
1117  MPI_Request *s_waits_col; // status of sens messages
1118  CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
1119 
1120  // Send messeges
1121  for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
1122  if (!lengths_cols[proc])
1123  continue; // no message to send to this proc
1124  CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
1125  lengths_cols[proc], // message length
1126  MPIU_INT, proc, // to proc
1127  tag_col, m_field.get_comm(), s_waits_col + kk);
1128  kk++;
1129  }
1130 
1131  if (nrecvs_cols) {
1132  CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
1133  }
1134  if (nsends_cols) {
1135  CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
1136  }
1137 
1138  CHKERR PetscFree(r_waits_col);
1139  CHKERR PetscFree(s_waits_col);
1140  }
1141 
1142  CHKERR PetscFree(status);
1143 
1144  // set values received from other processors
1145  for (int ss = 0; ss != loop_size; ++ss) {
1146 
1147  int nrecvs;
1148  int *olengths;
1149  int **data_procs;
1150  if (ss == 0) {
1151  nrecvs = nrecvs_rows;
1152  olengths = olengths_rows;
1153  data_procs = rbuf_row;
1154  } else {
1155  nrecvs = nrecvs_cols;
1156  olengths = olengths_cols;
1157  data_procs = rbuf_col;
1158  }
1159 
1160  const DofEntity_multiIndex *dofs_ptr;
1161  CHKERR m_field.get_dofs(&dofs_ptr);
1162 
1163  UId uid;
1164  NumeredDofEntity_multiIndex::iterator dit;
1165  for (int kk = 0; kk != nrecvs; ++kk) {
1166  int len = olengths[kk];
1167  int *data_from_proc = data_procs[kk];
1168  for (int dd = 0; dd < len; dd += data_block_size) {
1169  uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
1170  dit = numered_dofs_ptr[ss]->find(uid);
1171  if (dit == numered_dofs_ptr[ss]->end()) {
1172  // Dof is shared to this processor, however there is no element which
1173  // have this dof
1174  // continue;
1175  DofEntity_multiIndex::iterator ddit = dofs_ptr->find(uid);
1176  if (ddit != dofs_ptr->end()) {
1177  unsigned char pstatus = ddit->get()->getPStatus();
1178  if (pstatus > 0) {
1179  continue;
1180  } else {
1181  std::ostringstream zz;
1182  zz << **ddit << std::endl;
1183  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1184  "data inconsistency, dofs is not shared, but received "
1185  "from other proc\n"
1186  "%s",
1187  zz.str().c_str());
1188  }
1189  } else {
1190  std::ostringstream zz;
1191  zz << uid << std::endl;
1192  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1193  "no such dof %s in mofem database", zz.str().c_str());
1194  }
1195  std::ostringstream zz;
1196  zz << uid << std::endl;
1197  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1198  "dof %s not found", zz.str().c_str());
1199  }
1200  int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
1201  if (global_idx < 0) {
1202  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1203  "received negative dof");
1204  }
1205  bool success;
1206  success = numered_dofs_ptr[ss]->modify(
1207  dit, NumeredDofEntity_mofem_index_change(global_idx));
1208  if (!success)
1209  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1210  "modification unsuccessful");
1211  success = numered_dofs_ptr[ss]->modify(
1212  dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
1213  global_idx));
1214  if (!success)
1215  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1216  "modification unsuccessful");
1217  }
1218  }
1219  }
1220 
1221  if (square_matrix) {
1222  (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
1223  (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
1224  }
1225 
1226  CHKERR PetscFree(olengths_rows);
1227  CHKERR PetscFree(rbuf_row[0]);
1228  CHKERR PetscFree(rbuf_row);
1229  if (!square_matrix) {
1230  CHKERR PetscFree(olengths_cols);
1231  CHKERR PetscFree(rbuf_col[0]);
1232  CHKERR PetscFree(rbuf_col);
1233  }
1234 
1235  if (square_matrix) {
1236  if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
1237  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1238  "data inconsistency for square_matrix %d!=%d",
1239  numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
1240  }
1241  if (problem_ptr->numeredDofsRows != problem_ptr->numeredDofsCols) {
1242  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1243  "data inconsistency for square_matrix");
1244  }
1245  }
1246 
1247  CHKERR printPartitionedProblem(problem_ptr, verb);
1248  CHKERR debugPartitionedProblem(problem_ptr, verb);
1249 
1250  PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
1251 
1253 }
1254 
1256  const std::string out_name, const std::vector<std::string> &fields_row,
1257  const std::vector<std::string> &fields_col, const std::string main_problem,
1258  const bool square_matrix,
1259  const map<std::string, std::pair<EntityType, EntityType>> *entityMapRow,
1260  const map<std::string, std::pair<EntityType, EntityType>> *entityMapCol,
1261  int verb) {
1262  MoFEM::Interface &m_field = cOre;
1263  const Problem_multiIndex *problems_ptr;
1265 
1266  CHKERR m_field.clear_problem(out_name);
1267  CHKERR m_field.get_problems(&problems_ptr);
1268 
1269  // get reference to all problems
1270  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1271  ProblemByName &problems_by_name =
1272  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1273 
1274  // get iterators to out problem, i.e. build problem
1275  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1276  if (out_problem_it == problems_by_name.end()) {
1277  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1278  "subproblem with name < %s > not defined (top tip check spelling)",
1279  out_name.c_str());
1280  }
1281  // get iterator to main problem, i.e. out problem is subproblem of main
1282  // problem
1283  ProblemByName::iterator main_problem_it = problems_by_name.find(main_problem);
1284  if (main_problem_it == problems_by_name.end()) {
1285  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1286  "problem of subproblem with name < %s > not defined (top tip "
1287  "check spelling)",
1288  main_problem.c_str());
1289  }
1290 
1291  // get dofs for row & columns for out problem,
1292  boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1293  out_problem_it->numeredDofsRows, out_problem_it->numeredDofsCols};
1294  // get dofs for row & columns for main problem
1295  boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1296  main_problem_it->numeredDofsRows, main_problem_it->numeredDofsCols};
1297  // get local indices counter
1298  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1299  &out_problem_it->nbLocDofsCol};
1300  // get global indices counter
1301  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1302 
1303  // set number of ghost nodes to zero
1304  {
1305  out_problem_it->nbGhostDofsRow = 0;
1306  out_problem_it->nbGhostDofsCol = 0;
1307  }
1308 
1309  // put rows & columns field names in array
1310  std::vector<std::string> fields[] = {fields_row, fields_col};
1311  const map<std::string, std::pair<EntityType, EntityType>> *entityMap[] = {
1312  entityMapRow, entityMapCol};
1313 
1314  // make data structure fos sub-problem data
1315  out_problem_it->subProblemData =
1316  boost::make_shared<Problem::SubProblemData>();
1317 
1318  // Loop over rows and columns
1319  for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1320 
1321  // reset dofs and columns counters
1322  (*nb_local_dofs[ss]) = 0;
1323  (*nb_dofs[ss]) = 0;
1324  // clear arrays
1325  out_problem_dofs[ss]->clear();
1326 
1327  // If DOFs are cleared clear finite elements too.
1328  out_problem_it->numeredFiniteElements->clear();
1329 
1330  // get dofs by field name and insert them in out problem multi-indices
1331  for (auto field : fields[ss]) {
1332 
1333  // Following reserve memory in sequences, only two allocations are here,
1334  // once for array of objects, next for array of shared pointers
1335 
1336  // aliased sequence of pointer is killed with element
1337  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1338  boost::make_shared<std::vector<NumeredDofEntity>>();
1339  // reserve memory for field dofs
1340  if (!ss)
1341  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1342  else
1343  out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1344 
1345  // create elements objects
1346  auto dit =
1347  main_problem_dofs[ss]->get<FieldName_mi_tag>().lower_bound(field);
1348  auto hi_dit =
1349  main_problem_dofs[ss]->get<FieldName_mi_tag>().upper_bound(field);
1350 
1351  auto add_dit_to_dofs_array = [&](auto &dit) {
1352  dofs_array->emplace_back(
1353  dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1354  dit->get()->getPetscGlobalDofIdx(),
1355  dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1356  };
1357 
1358  if (entityMap[ss]) {
1359  auto mit = entityMap[ss]->find(field);
1360  if (mit != entityMap[ss]->end()) {
1361  EntityType lo_type = mit->second.first;
1362  EntityType hi_type = mit->second.second;
1363  int count = 0;
1364  for (auto diit = dit; diit != hi_dit; ++diit) {
1365  EntityType ent_type = (*diit)->getEntType();
1366  if (ent_type >= lo_type && ent_type <= hi_type)
1367  ++count;
1368  }
1369  dofs_array->reserve(count);
1370  for (; dit != hi_dit; ++dit) {
1371  EntityType ent_type = (*dit)->getEntType();
1372  if (ent_type >= lo_type && ent_type <= hi_type)
1373  add_dit_to_dofs_array(dit);
1374  }
1375  } else {
1376  dofs_array->reserve(std::distance(dit, hi_dit));
1377  for (; dit != hi_dit; dit++)
1378  add_dit_to_dofs_array(dit);
1379  }
1380  } else {
1381  dofs_array->reserve(std::distance(dit, hi_dit));
1382  for (; dit != hi_dit; dit++)
1383  add_dit_to_dofs_array(dit);
1384  }
1385 
1386  // fill multi-index
1387  auto hint = out_problem_dofs[ss]->end();
1388  for (auto &v : *dofs_array)
1389  hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1390  }
1391  // Set local indexes
1392  {
1393  auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1394  auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1395  for (; dit != hi_dit; dit++) {
1396  int idx = -1; // if dof is not part of partition, set local index to -1
1397  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1398  idx = (*nb_local_dofs[ss])++;
1399  }
1400  bool success = out_problem_dofs[ss]->modify(
1401  out_problem_dofs[ss]->project<0>(dit),
1403  if (!success) {
1404  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1405  "operation unsuccessful");
1406  }
1407  };
1408  }
1409  // Set global indexes, compress global indices
1410  {
1411  auto dit =
1412  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1413  auto hi_dit =
1414  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1415  out_problem_dofs[ss]->size());
1416  const int nb = std::distance(dit, hi_dit);
1417  // get main problem global indices
1418  std::vector<int> main_indices(nb);
1419  for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1420  *it = dit->get()->getPetscGlobalDofIdx();
1421  }
1422  // create is with global dofs
1423  IS is;
1424  CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1425  PETSC_USE_POINTER, &is);
1426  // create map form main problem global indices to out problem global
1427  // indices
1428  AO ao;
1429  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1430  if (ss == 0) {
1431  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->rowIs));
1432  // CHKERR ISSort(out_problem_it->getSubData()->rowIs);
1433  out_problem_it->getSubData()->rowMap = ao;
1434  CHKERR PetscObjectReference((PetscObject)ao);
1435  } else {
1436  CHKERR ISDuplicate(is, &(out_problem_it->getSubData()->colIs));
1437  // CHKERR ISSort(out_problem_it->getSubData()->colIs);
1438  out_problem_it->getSubData()->colMap = ao;
1439  CHKERR PetscObjectReference((PetscObject)ao);
1440  }
1441  CHKERR AOApplicationToPetscIS(ao, is);
1442  // set global number of DOFs
1443  CHKERR ISGetSize(is, nb_dofs[ss]);
1444  CHKERR ISDestroy(&is);
1445  // set out problem global indices after applying map
1446  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1447  for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1448  dit++, it++) {
1449  bool success = out_problem_dofs[ss]->modify(
1450  out_problem_dofs[ss]->project<0>(dit),
1452  dit->get()->getPart(), *it, *it,
1453  dit->get()->getPetscLocalDofIdx()));
1454  if (!success) {
1455  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1456  "operation unsuccessful");
1457  }
1458  }
1459  // set global indices to nodes not on this part
1460  {
1461  NumeredDofEntityByLocalIdx::iterator dit =
1462  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1463  NumeredDofEntityByLocalIdx::iterator hi_dit =
1464  out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1465  const int nb = std::distance(dit, hi_dit);
1466  std::vector<int> main_indices_non_local(nb);
1467  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1468  dit++, it++) {
1469  *it = dit->get()->getPetscGlobalDofIdx();
1470  }
1471  IS is;
1472  CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1473  &*main_indices_non_local.begin(),
1474  PETSC_USE_POINTER, &is);
1475  CHKERR AOApplicationToPetscIS(ao, is);
1476  CHKERR ISDestroy(&is);
1477  dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1478  for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1479  dit++, it++) {
1480  bool success = out_problem_dofs[ss]->modify(
1481  out_problem_dofs[ss]->project<0>(dit),
1483  dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1484  dit->get()->getPetscLocalDofIdx()));
1485  if (!success) {
1486  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1487  "operation unsuccessful");
1488  }
1489  }
1490  }
1491  CHKERR AODestroy(&ao);
1492  }
1493  }
1494 
1495  if (square_matrix) {
1496  out_problem_it->numeredDofsCols = out_problem_it->numeredDofsRows;
1497  out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1498  out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1499  out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1500  out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1501  CHKERR PetscObjectReference(
1502  (PetscObject)out_problem_it->getSubData()->rowIs);
1503  CHKERR PetscObjectReference(
1504  (PetscObject)out_problem_it->getSubData()->rowMap);
1505  }
1506 
1507  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1508  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1509 
1511 }
1512 
1514  const std::string out_name, const std::vector<std::string> add_row_problems,
1515  const std::vector<std::string> add_col_problems, const bool square_matrix,
1516  int verb) {
1518  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1519  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1520  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1521  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1522  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1523  MoFEM::Interface &m_field = cOre;
1524  const Problem_multiIndex *problems_ptr;
1526 
1527  CHKERR m_field.clear_problem(out_name);
1528  CHKERR m_field.get_problems(&problems_ptr);
1529  // get reference to all problems
1530  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1531  ProblemByName &problems_by_name =
1532  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1533 
1534  // Get iterators to out problem, i.e. build problem
1535  ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1536  if (out_problem_it == problems_by_name.end()) {
1537  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1538  "problem with name < %s > not defined (top tip check spelling)",
1539  out_name.c_str());
1540  }
1541  // Make data structure for composed-problem data
1542  out_problem_it->composedProblemsData =
1543  boost::make_shared<ComposedProblemsData>();
1544  boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1545  out_problem_it->getComposedProblemsData();
1546 
1547  const std::vector<std::string> *add_prb[] = {&add_row_problems,
1548  &add_col_problems};
1549  std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1550  &cmp_prb_data->colProblemsAdd};
1551  std::vector<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1552 
1553  // Get local indices counter
1554  int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1555  &out_problem_it->nbLocDofsCol};
1556  // Get global indices counter
1557  int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1558 
1559  // Set number of ghost nodes to zero
1560  {
1561  out_problem_it->nbDofsRow = 0;
1562  out_problem_it->nbDofsCol = 0;
1563  out_problem_it->nbLocDofsRow = 0;
1564  out_problem_it->nbLocDofsCol = 0;
1565  out_problem_it->nbGhostDofsRow = 0;
1566  out_problem_it->nbGhostDofsCol = 0;
1567  }
1568  int nb_dofs_reserve[] = {0, 0};
1569 
1570  // Loop over rows and columns in the main problem and sub-problems
1571  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1572  // cerr << "SS " << ss << endl;
1573  // cerr << add_prb[ss]->size() << endl;
1574  // cerr << add_prb_ptr[ss]->size() << endl;
1575  add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1576  add_prb_is[ss]->reserve(add_prb[ss]->size());
1577  for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1578  vit != add_prb[ss]->end(); vit++) {
1579  ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1580  if (prb_it == problems_by_name.end()) {
1581  SETERRQ1(
1582  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1583  "problem with name < %s > not defined (top tip check spelling)",
1584  vit->c_str());
1585  }
1586  add_prb_ptr[ss]->push_back(&*prb_it);
1587  // set number of dofs on rows and columns
1588  if (ss == 0) {
1589  // row
1590  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1591  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1592  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredDofsRows->size();
1593  } else {
1594  // column
1595  *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1596  *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1597  nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredDofsCols->size();
1598  }
1599  }
1600  }
1601  // if squre problem, rows and columns are the same
1602  if (square_matrix) {
1603  add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1604  add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1605  out_problem_it->numeredDofsCols = out_problem_it->numeredDofsRows;
1606  *nb_dofs[1] = *nb_dofs[0];
1607  *nb_local_dofs[1] = *nb_local_dofs[0];
1608  }
1609 
1610  // reserve memory for dofs
1611  boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1612  // Reserve memory
1613  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1614  dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1615  dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1616  if (!ss)
1617  out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1618  else
1619  out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1620  }
1621 
1622  // Push back DOFs
1623  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1624  NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1625  dit,
1626  hi_dit;
1627  int shift_glob = 0;
1628  int shift_loc = 0;
1629  for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1630  PetscInt *dofs_out_idx_ptr;
1631  int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1632  CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1633  if (ss == 0) {
1634  dit = (*add_prb_ptr[ss])[pp]
1635  ->numeredDofsRows->get<PetscGlobalIdx_mi_tag>()
1636  .begin();
1637  hi_dit = (*add_prb_ptr[ss])[pp]
1638  ->numeredDofsRows->get<PetscGlobalIdx_mi_tag>()
1639  .end();
1640  } else {
1641  dit = (*add_prb_ptr[ss])[pp]
1642  ->numeredDofsCols->get<PetscGlobalIdx_mi_tag>()
1643  .begin();
1644  hi_dit = (*add_prb_ptr[ss])[pp]
1645  ->numeredDofsCols->get<PetscGlobalIdx_mi_tag>()
1646  .end();
1647  }
1648  int is_nb = 0;
1649  for (; dit != hi_dit; dit++) {
1650  BitRefLevel prb_bit = out_problem_it->getBitRefLevel();
1651  BitRefLevel prb_mask = out_problem_it->getMaskBitRefLevel();
1652  BitRefLevel dof_bit = dit->get()->getBitRefLevel();
1653  if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1654  continue;
1655  const int rank = m_field.get_comm_rank();
1656  const int part = dit->get()->getPart();
1657  const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1658  const int loc_idx =
1659  (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1660  : -1;
1661  dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1662  glob_idx, loc_idx, part);
1663  if (part == rank) {
1664  dofs_out_idx_ptr[is_nb++] = glob_idx;
1665  }
1666  }
1667  if (is_nb > nb_local_dofs) {
1668  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1669  "Data inconsistency");
1670  }
1671  IS is;
1672  CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1673  PETSC_OWN_POINTER, &is);
1674  (*add_prb_is[ss]).push_back(is);
1675  if (ss == 0) {
1676  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1677  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1678  } else {
1679  shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1680  shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1681  }
1682  if (square_matrix) {
1683  (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1684  (*add_prb_is[1]).push_back(is);
1685  CHKERR PetscObjectReference((PetscObject)is);
1686  }
1687  }
1688  }
1689 
1690  if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1691  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1692  }
1693 
1694  // Insert DOFs to problem multi-index
1695  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1696  auto hint = (ss == 0) ? out_problem_it->numeredDofsRows->end()
1697  : out_problem_it->numeredDofsCols->end();
1698  for (auto &v : *dofs_array[ss])
1699  hint = (ss == 0) ? out_problem_it->numeredDofsRows->emplace_hint(
1700  hint, dofs_array[ss], &v)
1701  : out_problem_it->numeredDofsCols->emplace_hint(
1702  hint, dofs_array[ss], &v);
1703  }
1704 
1705  // Compress DOFs
1706  *nb_dofs[0] = 0;
1707  *nb_dofs[1] = 0;
1708  *nb_local_dofs[0] = 0;
1709  *nb_local_dofs[1] = 0;
1710  for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1711 
1712  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1713  if (ss == 0) {
1714  dofs_ptr = out_problem_it->numeredDofsRows;
1715  } else {
1716  dofs_ptr = out_problem_it->numeredDofsCols;
1717  }
1718  NumeredDofEntityByUId::iterator dit, hi_dit;
1719  dit = dofs_ptr->get<Unique_mi_tag>().begin();
1720  hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1721  std::vector<int> idx;
1722  idx.reserve(std::distance(dit, hi_dit));
1723  // set dofs in order entity and dof number on entity
1724  for (; dit != hi_dit; dit++) {
1725  if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1726  bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1727  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1728  if (!success) {
1729  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1730  "modification unsuccessful");
1731  }
1732  idx.push_back(dit->get()->getPetscGlobalDofIdx());
1733  } else {
1734  if (dit->get()->getPetscLocalDofIdx() != -1) {
1735  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1736  "local index should be negative");
1737  }
1738  }
1739  }
1740  if (square_matrix) {
1741  *nb_local_dofs[1] = *nb_local_dofs[0];
1742  }
1743 
1744  // set new dofs mapping
1745  IS is;
1746  CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1747  PETSC_USE_POINTER, &is);
1748  CHKERR ISGetSize(is, nb_dofs[ss]);
1749  if (square_matrix) {
1750  *nb_dofs[1] = *nb_dofs[0];
1751  }
1752  // {
1753  // PetscSynchronizedPrintf(m_field.get_comm(),"nb dofs %d %d
1754  // %d\n",*nb_local_dofs[0],idx.size(),*nb_dofs[0]);
1755  // PetscSynchronizedFlush(m_field.get_comm(),PETSC_STDOUT);
1756  // }
1757  AO ao;
1758  CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1759  for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1760  CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1761 
1762  // Set DOFs numeration
1763  {
1764  std::vector<int> idx_new;
1765  idx_new.reserve(dofs_ptr->size());
1766  for (NumeredDofEntityByUId::iterator dit =
1767  dofs_ptr->get<Unique_mi_tag>().begin();
1768  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1769  idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1770  }
1771  // set new global dofs numeration
1772  IS is_new;
1773  CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1774  &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1775  CHKERR AOApplicationToPetscIS(ao, is_new);
1776  // set global indices to multi-index
1777  std::vector<int>::iterator vit = idx_new.begin();
1778  for (NumeredDofEntityByUId::iterator dit =
1779  dofs_ptr->get<Unique_mi_tag>().begin();
1780  dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1781  bool success =
1782  dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1783  dit->get()->getPart(), *(vit++)));
1784  if (!success) {
1785  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1786  "modification unsuccessful");
1787  }
1788  }
1789  CHKERR ISDestroy(&is_new);
1790  }
1791  CHKERR ISDestroy(&is);
1792  CHKERR AODestroy(&ao);
1793  }
1794 
1795  CHKERR printPartitionedProblem(&*out_problem_it, verb);
1796  CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1797 
1798  // Inidcate that porble has been build
1801 
1803 }
1804 
1806  int verb) {
1807 
1808  MoFEM::Interface &m_field = cOre;
1809  const Problem_multiIndex *problems_ptr;
1812  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1813  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1814  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1815  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1816  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1818  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1819  if (verb > 0) {
1820  PetscPrintf(m_field.get_comm(), "Simple partition problem %s\n",
1821  name.c_str());
1822  }
1823  CHKERR m_field.get_problems(&problems_ptr);
1824  // find p_miit
1825  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1826  ProblemByName &problems_set =
1827  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1828  ProblemByName::iterator p_miit = problems_set.find(name);
1829  if (p_miit == problems_set.end()) {
1830  SETERRQ1(PETSC_COMM_SELF, 1,
1831  "problem < %s > is not found (top tip: check spelling)",
1832  name.c_str());
1833  }
1834  typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1835  Idx_mi_tag>::type NumeredDofEntitysByIdx;
1836  NumeredDofEntitysByIdx &dofs_row_by_idx =
1837  p_miit->numeredDofsRows->get<Idx_mi_tag>();
1838  NumeredDofEntitysByIdx &dofs_col_by_idx =
1839  p_miit->numeredDofsCols->get<Idx_mi_tag>();
1840  boost::multi_index::index<NumeredDofEntity_multiIndex,
1841  Idx_mi_tag>::type::iterator miit_row,
1842  hi_miit_row;
1843  boost::multi_index::index<NumeredDofEntity_multiIndex,
1844  Idx_mi_tag>::type::iterator miit_col,
1845  hi_miit_col;
1846  DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1847  DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1848  nb_row_local_dofs = 0;
1849  nb_row_ghost_dofs = 0;
1850  DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1851  DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1852  nb_col_local_dofs = 0;
1853  nb_col_ghost_dofs = 0;
1854 
1855  bool square_matrix = false;
1856  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
1857  square_matrix = true;
1858  }
1859 
1860  // get row range of local indices
1861  PetscLayout layout_row;
1862  const int *ranges_row;
1863 
1864  DofIdx nb_dofs_row = dofs_row_by_idx.size();
1865  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1866  CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1867  CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1868  CHKERR PetscLayoutSetUp(layout_row);
1869  CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1870  // get col range of local indices
1871  PetscLayout layout_col;
1872  const int *ranges_col;
1873  if (!square_matrix) {
1874  DofIdx nb_dofs_col = dofs_col_by_idx.size();
1875  CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1876  CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1877  CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1878  CHKERR PetscLayoutSetUp(layout_col);
1879  CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1880  }
1881  for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1882  part++) {
1883  miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1884  hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1885  if (std::distance(miit_row, hi_miit_row) !=
1886  ranges_row[part + 1] - ranges_row[part]) {
1887  SETERRQ4(
1888  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1889  "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1890  "rstart (%d != %d - %d = %d) ",
1891  std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1892  ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1893  }
1894  // loop rows
1895  for (; miit_row != hi_miit_row; miit_row++) {
1896  bool success = dofs_row_by_idx.modify(
1897  miit_row,
1898  NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1899  if (!success)
1900  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1901  "modification unsuccessful");
1902  if (part == (unsigned int)m_field.get_comm_rank()) {
1903  success = dofs_row_by_idx.modify(
1904  miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1905  if (!success)
1906  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1907  "modification unsuccessful");
1908  }
1909  }
1910  if (!square_matrix) {
1911  miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1912  hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1913  if (std::distance(miit_col, hi_miit_col) !=
1914  ranges_col[part + 1] - ranges_col[part]) {
1915  SETERRQ4(
1916  PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1917  "data inconsistency, std::distance(miit_col,hi_miit_col) != rend - "
1918  "rstart (%d != %d - %d = %d) ",
1919  std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1920  ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1921  }
1922  // loop cols
1923  for (; miit_col != hi_miit_col; miit_col++) {
1924  bool success = dofs_col_by_idx.modify(
1926  part, (*miit_col)->dofIdx));
1927  if (!success)
1928  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1929  "modification unsuccessful");
1930  if (part == (unsigned int)m_field.get_comm_rank()) {
1931  success = dofs_col_by_idx.modify(
1932  miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1933  if (!success)
1934  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1935  "modification unsuccessful");
1936  }
1937  }
1938  }
1939  }
1940  CHKERR PetscLayoutDestroy(&layout_row);
1941  if (!square_matrix) {
1942  CHKERR PetscLayoutDestroy(&layout_col);
1943  }
1944  if (square_matrix) {
1945  nb_col_local_dofs = nb_row_local_dofs;
1946  nb_col_ghost_dofs = nb_row_ghost_dofs;
1947  }
1948  CHKERR printPartitionedProblem(&*p_miit, verb);
1951 }
1952 
1954  int verb) {
1955  MoFEM::Interface &m_field = cOre;
1956  const Problem_multiIndex *problems_ptr;
1958 
1959  if (!(cOre.getBuildMoFEM() & (1 << 0)))
1960  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1961  if (!(cOre.getBuildMoFEM() & (1 << 1)))
1962  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1963  if (!(cOre.getBuildMoFEM() & (1 << 2)))
1964  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1965  "entFEAdjacencies not build");
1966  if (!(cOre.getBuildMoFEM() & (1 << 3)))
1967  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1968 
1969  if (verb > QUIET)
1970  PetscPrintf(m_field.get_comm(), "Partition problem %s\n", name.c_str());
1971 
1972  typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
1973  NumeredDofEntitysByIdx;
1974  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
1975 
1976  // Find problem pointer by name
1977  CHKERR m_field.get_problems(&problems_ptr);
1978  auto &problems_set =
1979  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1980  auto p_miit = problems_set.find(name);
1981  if (p_miit == problems_set.end()) {
1982  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1983  "problem with name %s not defined (top tip check spelling)",
1984  name.c_str());
1985  }
1986  int nb_dofs_row = p_miit->getNbDofsRow();
1987 
1988  Mat Adj;
1989  CHKERR m_field.getInterface<MatrixManager>()
1990  ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1991 
1992  int m, n;
1993  CHKERR MatGetSize(Adj, &m, &n);
1994  if (verb > VERY_VERBOSE)
1995  MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1996 
1997  // partitioning
1998  MatPartitioning part;
1999  IS is;
2000  CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
2001  CHKERR MatPartitioningSetAdjacency(part, Adj);
2002  CHKERR MatPartitioningSetFromOptions(part);
2003  CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
2004  CHKERR MatPartitioningApply(part, &is);
2005  if (verb > VERY_VERBOSE)
2006  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
2007 
2008  // gather
2009  IS is_gather, is_num, is_gather_num;
2010  CHKERR ISAllGather(is, &is_gather);
2011  CHKERR ISPartitioningToNumbering(is, &is_num);
2012  CHKERR ISAllGather(is_num, &is_gather_num);
2013  const int *part_number, *petsc_idx;
2014  CHKERR ISGetIndices(is_gather, &part_number);
2015  CHKERR ISGetIndices(is_gather_num, &petsc_idx);
2016  int size_is_num, size_is_gather;
2017  CHKERR ISGetSize(is_gather, &size_is_gather);
2018  if (size_is_gather != (int)nb_dofs_row)
2019  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
2020  size_is_gather, nb_dofs_row);
2021 
2022  CHKERR ISGetSize(is_num, &size_is_num);
2023  if (size_is_num != (int)nb_dofs_row)
2024  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "data inconsistency %d != %d",
2025  size_is_num, nb_dofs_row);
2026 
2027  bool square_matrix = false;
2028  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols)
2029  square_matrix = true;
2030 
2031  if (!square_matrix) {
2032  // FIXME: This is for back compatibility, if deprecate interface function
2033  // build interfaces is removed, this part of the code will be obsolete
2034  auto mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().begin();
2035  auto hi_mit_row = p_miit->numeredDofsRows->get<Idx_mi_tag>().end();
2036  auto mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().begin();
2037  auto hi_mit_col = p_miit->numeredDofsCols->get<Idx_mi_tag>().end();
2038  for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
2039  if (mit_col == hi_mit_col) {
2040  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2041  "check finite element definition, nb. of rows is not equal to "
2042  "number for columns");
2043  }
2044  if (mit_row->get()->getGlobalUniqueId() !=
2045  mit_col->get()->getGlobalUniqueId()) {
2046  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2047  "check finite element definition, nb. of rows is not equal to "
2048  "number for columns");
2049  }
2050  }
2051  }
2052 
2053  if (verb > VERBOSE)
2054  PetscPrintf(m_field.get_comm(), "\tloop problem dofs");
2055 
2056  // Set petsc global indices
2057  auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
2058  p_miit->numeredDofsRows->get<Idx_mi_tag>());
2059  int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2060  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2061  nb_row_local_dofs = 0;
2062  nb_row_ghost_dofs = 0;
2063 
2064  for (auto miit_dofs_row = dofs_row_by_idx_no_const.begin();
2065  miit_dofs_row != dofs_row_by_idx_no_const.end(); miit_dofs_row++) {
2066  const int part = part_number[(*miit_dofs_row)->dofIdx];
2067  if (part == (unsigned int)m_field.get_comm_rank()) {
2068  const bool success = dofs_row_by_idx_no_const.modify(
2069  miit_dofs_row,
2071  part, petsc_idx[(*miit_dofs_row)->dofIdx], nb_row_local_dofs++));
2072  if (!success) {
2073  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2074  "modification unsuccessful");
2075  }
2076  } else {
2077  const bool success = dofs_row_by_idx_no_const.modify(
2079  part, petsc_idx[(*miit_dofs_row)->dofIdx]));
2080  if (!success) {
2081  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2082  "modification unsuccessful");
2083  }
2084  }
2085  }
2086 
2087  int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2088  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2089  if (square_matrix) {
2090  nb_col_local_dofs = nb_row_local_dofs;
2091  nb_col_ghost_dofs = nb_row_ghost_dofs;
2092  } else {
2093  NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2094  const_cast<NumeredDofEntitysByIdx &>(
2095  p_miit->numeredDofsCols->get<Idx_mi_tag>());
2096  nb_col_local_dofs = 0;
2097  nb_col_ghost_dofs = 0;
2098  for (auto miit_dofs_col = dofs_col_by_idx_no_const.begin();
2099  miit_dofs_col != dofs_col_by_idx_no_const.end(); miit_dofs_col++) {
2100  const int part = part_number[(*miit_dofs_col)->dofIdx];
2101  if (part == (unsigned int)m_field.get_comm_rank()) {
2102  const bool success = dofs_col_by_idx_no_const.modify(
2104  part, petsc_idx[(*miit_dofs_col)->dofIdx],
2105  nb_col_local_dofs++));
2106  if (!success) {
2107  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2108  "modification unsuccessful");
2109  }
2110  } else {
2111  const bool success = dofs_col_by_idx_no_const.modify(
2113  part, petsc_idx[(*miit_dofs_col)->dofIdx]));
2114  if (!success) {
2115  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2116  "modification unsuccessful");
2117  }
2118  }
2119  }
2120  }
2121 
2122  if (verb > VERBOSE)
2123  PetscPrintf(m_field.get_comm(), " <- done\n");
2124 
2125  CHKERR ISRestoreIndices(is_gather, &part_number);
2126  CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
2127  CHKERR ISDestroy(&is_num);
2128  CHKERR ISDestroy(&is_gather_num);
2129  CHKERR ISDestroy(&is_gather);
2130  CHKERR ISDestroy(&is);
2131  CHKERR MatPartitioningDestroy(&part);
2132  CHKERR MatDestroy(&Adj);
2133  CHKERR printPartitionedProblem(&*p_miit, verb);
2134 
2135  cOre.getBuildMoFEM() |= 1 << 4;
2137 }
2138 
2140  const std::string name, const std::string problem_for_rows, bool copy_rows,
2141  const std::string problem_for_cols, bool copy_cols, int verb) {
2142 
2143  MoFEM::Interface &m_field = cOre;
2144  const Problem_multiIndex *problems_ptr;
2146 
2148  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
2149 
2150  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2151 
2152  // find p_miit
2153  CHKERR m_field.get_problems(&problems_ptr);
2154  ProblemByName &problems_by_name =
2155  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2156  ProblemByName::iterator p_miit = problems_by_name.find(name);
2157  if (p_miit == problems_by_name.end()) {
2158  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2159  "problem with name < %s > not defined (top tip check spelling)",
2160  name.c_str());
2161  }
2162  if (verb > QUIET)
2163  PetscPrintf(m_field.get_comm(),
2164  "Compose problem %s from rows of %s and columns of %s\n",
2165  p_miit->getName().c_str(), problem_for_rows.c_str(),
2166  problem_for_cols.c_str());
2167 
2168  // find p_miit_row
2169  ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
2170  if (p_miit_row == problems_by_name.end()) {
2171  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2172  "problem with name < %s > not defined (top tip check spelling)",
2173  problem_for_rows.c_str());
2174  }
2175  const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
2176  p_miit_row->numeredDofsRows;
2177 
2178  // find p_mit_col
2179  ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
2180  if (p_miit_col == problems_by_name.end()) {
2181  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2182  "problem with name < %s > not defined (top tip check spelling)",
2183  problem_for_cols.c_str());
2184  }
2185  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
2186  p_miit_col->numeredDofsCols;
2187 
2188  bool copy[] = {copy_rows, copy_cols};
2189  boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
2190  p_miit->numeredDofsRows, p_miit->numeredDofsCols};
2191 
2192  int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
2193  int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
2194  boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
2195  dofs_col};
2196 
2197  for (int ss = 0; ss < 2; ss++) {
2198 
2199  // build indices
2200  *nb_local_dofs[ss] = 0;
2201  if (!copy[ss]) {
2202 
2203  // only copy indices which are belong to some elements if this problem
2204  std::vector<int> is_local, is_new;
2205 
2206  NumeredDofEntityByUId &dofs_by_uid =
2207  copied_dofs[ss]->get<Unique_mi_tag>();
2208  for (NumeredDofEntity_multiIndex::iterator dit =
2209  composed_dofs[ss]->begin();
2210  dit != composed_dofs[ss]->end(); dit++) {
2211  NumeredDofEntityByUId::iterator diit =
2212  dofs_by_uid.find((*dit)->getGlobalUniqueId());
2213  if (diit == dofs_by_uid.end()) {
2214  SETERRQ(
2216  "data inconsistency, could not find dof in composite problem");
2217  }
2218  int part_number = (*diit)->getPart(); // get part number
2219  int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
2220  bool success;
2221  success = composed_dofs[ss]->modify(
2223  petsc_global_dof));
2224  if (!success) {
2225  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2226  "modification unsuccessful");
2227  }
2228  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2229  success = composed_dofs[ss]->modify(
2230  dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
2231  if (!success) {
2232  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2233  "modification unsuccessful");
2234  }
2235  is_local.push_back(petsc_global_dof);
2236  }
2237  }
2238 
2239  AO ao;
2240  CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
2241  NULL, &ao);
2242 
2243  // apply local to global mapping
2244  is_local.resize(0);
2245  for (NumeredDofEntity_multiIndex::iterator dit =
2246  composed_dofs[ss]->begin();
2247  dit != composed_dofs[ss]->end(); dit++) {
2248  is_local.push_back((*dit)->getPetscGlobalDofIdx());
2249  }
2250  CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2251  int idx2 = 0;
2252  for (NumeredDofEntity_multiIndex::iterator dit =
2253  composed_dofs[ss]->begin();
2254  dit != composed_dofs[ss]->end(); dit++) {
2255  int part_number = (*dit)->getPart(); // get part number
2256  int petsc_global_dof = is_local[idx2++];
2257  bool success;
2258  success = composed_dofs[ss]->modify(
2260  petsc_global_dof));
2261  if (!success) {
2262  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2263  "modification unsuccessful");
2264  }
2265  }
2266 
2267  CHKERR AODestroy(&ao);
2268 
2269  } else {
2270 
2271  for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2272  dit != copied_dofs[ss]->end(); dit++) {
2273  std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2274  p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2275  new NumeredDofEntity((*dit)->getDofEntityPtr())));
2276  if (p.second) {
2277  (*nb_dofs[ss])++;
2278  }
2279  int dof_idx = (*dit)->getDofIdx();
2280  int part_number = (*dit)->getPart(); // get part number
2281  int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2282  if (part_number == (unsigned int)m_field.get_comm_rank()) {
2283  const bool success = composed_dofs[ss]->modify(
2285  part_number, dof_idx, petsc_global_dof,
2286  (*nb_local_dofs[ss])++));
2287  if (!success) {
2288  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2289  "modification unsuccessful");
2290  }
2291  } else {
2292  const bool success = composed_dofs[ss]->modify(
2294  part_number, dof_idx, petsc_global_dof));
2295  if (!success) {
2296  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2297  "modification unsuccessful");
2298  }
2299  }
2300  }
2301  }
2302  }
2303 
2304  CHKERR printPartitionedProblem(&*p_miit, verb);
2305  CHKERR debugPartitionedProblem(&*p_miit, verb);
2306 
2308 }
2309 
2311 ProblemsManager::printPartitionedProblem(const Problem *problem_ptr, int verb) {
2312  MoFEM::Interface &m_field = cOre;
2314  if (verb > 0) {
2315  std::ostringstream ss;
2316  ss << "partition_problem: rank = " << m_field.get_comm_rank()
2317  << " FEs row ghost dofs " << *problem_ptr << " Nb. local dof "
2318  << problem_ptr->getNbLocalDofsRow() << " nb global row dofs "
2319  << problem_ptr->getNbDofsRow() << std::endl;
2320  ss << "partition_problem: rank = " << m_field.get_comm_rank()
2321  << " FEs col ghost dofs " << *problem_ptr << " Nb. local dof "
2322  << problem_ptr->getNbLocalDofsCol() << " nb global col dofs "
2323  << problem_ptr->getNbDofsCol() << std::endl;
2324  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2325  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2326  }
2327  if (verb > 1) {
2328  // FIXME mess if printed on more than one processors
2329  // std::ostringstream ss;
2330  std::cout << "rank = " << m_field.get_comm_rank() << " FEs row dofs "
2331  << *problem_ptr << " Nb. row dof " << problem_ptr->getNbDofsRow()
2332  << " Nb. local dof " << problem_ptr->getNbLocalDofsRow()
2333  << std::endl;
2334  NumeredDofEntity_multiIndex::iterator miit_dd_row =
2335  problem_ptr->numeredDofsRows->begin();
2336  for (; miit_dd_row != problem_ptr->numeredDofsRows->end(); miit_dd_row++) {
2337  std::cout << **miit_dd_row << std::endl;
2338  }
2339  std::cout << "rank = " << m_field.get_comm_rank() << " FEs col dofs "
2340  << *problem_ptr << " Nb. col dof " << problem_ptr->getNbDofsCol()
2341  << " Nb. local dof " << problem_ptr->getNbLocalDofsCol()
2342  << std::endl;
2343  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2344  problem_ptr->numeredDofsCols->begin();
2345  for (; miit_dd_col != problem_ptr->numeredDofsCols->end(); miit_dd_col++) {
2346  std::cout << **miit_dd_col << std::endl;
2347  }
2348  // PetscSynchronizedPrintf(comm,ss.str().c_str());
2349  // PetscSynchronizedFlush(comm,PETSC_STDOUT);
2350  }
2352 }
2353 
2355 ProblemsManager::debugPartitionedProblem(const Problem *problem_ptr, int verb) {
2356  MoFEM::Interface &m_field = cOre;
2358  if (debug > 0) {
2359 
2360  typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
2361  NumeredDofEntitysByIdx;
2362  NumeredDofEntitysByIdx::iterator dit, hi_dit;
2363  const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
2364  &(problem_ptr->numeredDofsRows->get<Idx_mi_tag>()),
2365  &(problem_ptr->numeredDofsCols->get<Idx_mi_tag>())};
2366 
2367  int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
2368  int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
2369  &problem_ptr->nbLocDofsCol};
2370 
2371  for (int ss = 0; ss < 2; ss++) {
2372 
2373  dit = numered_dofs_ptr[ss]->begin();
2374  hi_dit = numered_dofs_ptr[ss]->end();
2375  for (; dit != hi_dit; dit++) {
2376  if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2377  if ((*dit)->getPetscLocalDofIdx() < 0) {
2378  std::ostringstream zz;
2379  zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2380  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
2381  "local dof index for %d (0-row, 1-col) not set, i.e. has "
2382  "negative value\n %s",
2383  ss, zz.str().c_str());
2384  }
2385  if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
2386  std::ostringstream zz;
2387  zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2388  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
2389  "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
2390  zz.str().c_str());
2391  }
2392  } else {
2393  if ((*dit)->getPetscGlobalDofIdx() < 0) {
2394  std::ostringstream zz;
2395  zz << "rank " << m_field.get_comm_rank() << " "
2396  << dit->get()->getBitRefLevel() << " " << **dit;
2397  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
2398  "global dof index for %d (0-row, 1-col) row not set, i.e. "
2399  "has negative value\n %s",
2400  ss, zz.str().c_str());
2401  }
2402  if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
2403  std::ostringstream zz;
2404  zz << "rank " << m_field.get_comm_rank() << " nb_dofs "
2405  << *nbdof_ptr[ss] << " " << **dit;
2406  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
2407  "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
2408  zz.str().c_str());
2409  }
2410  }
2411  }
2412  }
2413  }
2415 }
2416 
2418  bool part_from_moab,
2419  int low_proc,
2420  int hi_proc, int verb) {
2421  MoFEM::Interface &m_field = cOre;
2422  const Problem_multiIndex *problems_ptr;
2423  const EntFiniteElement_multiIndex *fe_ent_ptr;
2425 
2427  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2428 
2429  if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2430  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2431 
2432  if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2433  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2434  "adjacencies not build");
2435 
2437  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2438 
2440  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2441  "problem not partitioned");
2442 
2443  if (low_proc == -1)
2444  low_proc = m_field.get_comm_rank();
2445  if (hi_proc == -1)
2446  hi_proc = m_field.get_comm_rank();
2447 
2448  // Find pointer to problem of given name
2449  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2450  // Get pointer to problem data struture
2451  CHKERR m_field.get_problems(&problems_ptr);
2452  ProblemByName &problems =
2453  const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2454  ProblemByName::iterator p_miit = problems.find(name);
2455  if (p_miit == problems.end()) {
2456  SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2457  "problem < %s > not found (top tip: check spelling)",
2458  name.c_str());
2459  }
2460 
2461  // Get reference on finite elements multi-index on the problem
2462  NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2463  *p_miit->numeredFiniteElements;
2464 
2465  // Clear all elements and data, build it again
2466  problem_finite_elements.clear();
2467 
2468  // Check if dofs and columns are the same, i.e. structurally symmetric problem
2469  bool do_cols_prob = true;
2470  if (p_miit->numeredDofsRows == p_miit->numeredDofsCols) {
2471  do_cols_prob = false;
2472  }
2473 
2474  // Loop over all elements in database and if right element is there add it
2475  // to problem finite element multi-index
2476  CHKERR m_field.get_ents_finite_elements(&fe_ent_ptr);
2477  for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); efit++) {
2478 
2479  // if element is not part of problem
2480  if (((*efit)->getId() & p_miit->getBitFEId()).none())
2481  continue;
2482 
2483  BitRefLevel prb_bit = p_miit->getBitRefLevel();
2484  BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
2485  BitRefLevel fe_bit = (*efit)->getBitRefLevel();
2486  // if entity is not problem refinement level
2487  if ((fe_bit & prb_mask) != fe_bit)
2488  continue;
2489  if ((fe_bit & prb_bit) != prb_bit)
2490  continue;
2491 
2492  // create element
2493  boost::shared_ptr<NumeredEntFiniteElement> numered_fe(
2494  new NumeredEntFiniteElement(*efit));
2495 
2496  // check if rows and columns are the same on this element
2497  bool do_cols_fe = true;
2498  if ((numered_fe->sPtr->row_field_ents_view ==
2499  numered_fe->sPtr->col_field_ents_view) &&
2500  !do_cols_prob) {
2501  do_cols_fe = false;
2502  numered_fe->cols_dofs = numered_fe->rows_dofs;
2503  } else {
2504  // different dofs on rows and columns
2505  numered_fe->cols_dofs = boost::shared_ptr<FENumeredDofEntity_multiIndex>(
2507  }
2508  // get pointer to dofs multi-index on rows and columns
2509  auto rows_dofs = numered_fe->rows_dofs;
2510  auto cols_dofs = numered_fe->cols_dofs;
2511  // clear multi-indices
2512  rows_dofs->clear();
2513  if (do_cols_fe) {
2514  cols_dofs->clear();
2515  }
2518 
2519  // set partition to the element
2520  {
2521  if (part_from_moab) {
2522  // if partition is taken from moab partition
2523  int proc = (*efit)->getPartProc();
2524  if (proc == -1 && (*efit)->getEntType() == MBVERTEX) {
2525  proc = (*efit)->getOwnerProc();
2526  }
2527  NumeredEntFiniteElement_change_part(proc).operator()(numered_fe);
2528  } else {
2529  // count partition of the dofs in row, the larges dofs with given
2530  // partition is used to set partition of the element
2531  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows), rows_view,
2532  moab::Interface::UNION);
2533  std::vector<int> parts(m_field.get_comm_size(), 0);
2534  NumeredDofEntity_multiIndex_uid_view_ordered::iterator viit_rows;
2535  viit_rows = rows_view.begin();
2536  for (; viit_rows != rows_view.end(); viit_rows++) {
2537  parts[(*viit_rows)->pArt]++;
2538  }
2539  std::vector<int>::iterator pos =
2540  max_element(parts.begin(), parts.end());
2541  unsigned int max_part = std::distance(parts.begin(), pos);
2542  NumeredEntFiniteElement_change_part(max_part).operator()(numered_fe);
2543  }
2544  }
2545 
2546  // set dofs on rows and columns (if are different)
2547  if ((numered_fe->getPart() >= (unsigned int)low_proc) &&
2548  (numered_fe->getPart() <= (unsigned int)hi_proc)) {
2549 
2550  NumeredDofEntity_multiIndex_uid_view_ordered *dofs_view[] = {&rows_view,
2551  &cols_view};
2552  FENumeredDofEntity_multiIndex *fe_dofs[] = {rows_dofs.get(),
2553  cols_dofs.get()};
2554 
2555  for (int ss = 0; ss != (do_cols_fe ? 2 : 1); ss++) {
2556 
2557  if (ss == 0) {
2558  if (part_from_moab) {
2559  // get row_view
2560  CHKERR(*efit)->getRowDofView(*(p_miit->numeredDofsRows),
2561  *dofs_view[ss],
2562  moab::Interface::UNION);
2563  }
2564  } else {
2565  // get cols_views
2566  CHKERR(*efit)->getColDofView(*(p_miit->numeredDofsCols),
2567  *dofs_view[ss], moab::Interface::UNION);
2568  }
2569 
2570  // Following reserve memory in sequences, only two allocations are here,
2571  // once for array of objects, next for array of shared pointers
2572 
2573  // reserve memory for field dofs
2574  boost::shared_ptr<std::vector<FENumeredDofEntity>> dofs_array(
2575  new std::vector<FENumeredDofEntity>());
2576 
2577  if (!ss) {
2578  numered_fe->getRowDofsSequence() = dofs_array;
2579  if (!do_cols_fe)
2580  numered_fe->getColDofsSequence() = dofs_array;
2581  } else
2582  numered_fe->getColDofsSequence() = dofs_array;
2583 
2584  auto vit = dofs_view[ss]->begin();
2585  auto hi_vit = dofs_view[ss]->end();
2586 
2587  dofs_array->reserve(std::distance(vit, hi_vit));
2588 
2589  // create elements objects
2590  for (; vit != hi_vit; vit++) {
2591  boost::shared_ptr<SideNumber> side_number_ptr;
2592  side_number_ptr = (*efit)->getSideNumberPtr((*vit)->getEnt());
2593  dofs_array->emplace_back(side_number_ptr, *vit);
2594  }
2595 
2596  // finally add DoFS to multi-indices
2597  auto hint = fe_dofs[ss]->end();
2598  for (auto &v : *dofs_array)
2599  hint = fe_dofs[ss]->emplace_hint(hint, dofs_array, &v);
2600  }
2601  }
2602  if (!numered_fe->sPtr->row_field_ents_view->empty() &&
2603  !numered_fe->sPtr->col_field_ents_view->empty()) {
2604 
2605  // Add element to the problem
2606  auto p = problem_finite_elements.insert(numered_fe);
2607  if (!p.second)
2608  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2609 
2610  if (verb >= VERY_VERBOSE) {
2611  std::ostringstream ss;
2612  ss << *p_miit << std::endl;
2613  ss << *p.first << std::endl;
2615  FENumeredDofEntityByUId::iterator miit =
2616  (*p.first)->rows_dofs->get<Unique_mi_tag>().begin();
2617  for (; miit != (*p.first)->rows_dofs->get<Unique_mi_tag>().end();
2618  miit++)
2619  ss << "rows: " << *(*miit) << std::endl;
2620  miit = (*p.first)->cols_dofs->get<Unique_mi_tag>().begin();
2621  for (; miit != (*p.first)->cols_dofs->get<Unique_mi_tag>().end();
2622  miit++)
2623  ss << "cols: " << *(*miit) << std::endl;
2624  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2625  }
2626  }
2627  }
2628 
2629  if (verb >= VERBOSE) {
2630  typedef NumeredEntFiniteElement_multiIndex::index<Part_mi_tag>::type
2631  NumeredEntFiniteElementPart;
2632  NumeredEntFiniteElementPart::iterator miit, hi_miit;
2633  miit = problem_finite_elements.get<Part_mi_tag>().lower_bound(
2634  m_field.get_comm_rank());
2635  hi_miit = problem_finite_elements.get<Part_mi_tag>().upper_bound(
2636  m_field.get_comm_rank());
2637  int count = std::distance(miit, hi_miit);
2638  std::ostringstream ss;
2639  ss << *p_miit;
2640  ss << " Nb. elems " << count << " on proc " << m_field.get_comm_rank()
2641  << std::endl;
2642  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2643  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2644  }
2645 
2648 }
2649 
2651  int verb) {
2652  MoFEM::Interface &m_field = cOre;
2653  const Problem_multiIndex *problems_ptr;
2655 
2657  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2658  "partition of problem not build");
2660  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2661  "partitions finite elements not build");
2662 
2663  // find p_miit
2664  CHKERR m_field.get_problems(&problems_ptr);
2665 
2666  // get problem pointer
2667  auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2668  if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2669  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2670  name.c_str());
2671 
2672  // get reference to number of ghost dofs
2673  int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2674  int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2675  nb_row_ghost_dofs = 0;
2676  nb_col_ghost_dofs = 0;
2677 
2678  // do work if more than one processor
2679  if (m_field.get_comm_size() > 1) {
2680 
2682  ghost_idx_col_view;
2683 
2684  // get elements on this partition
2685  auto fe_range =
2686  p_miit->numeredFiniteElements->get<Part_mi_tag>().equal_range(
2687  m_field.get_comm_rank());
2688 
2689  // get dofs on elements which are not part of this partition
2690 
2691  // rows
2692  auto hint_r = ghost_idx_row_view.begin();
2693  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2694  for (auto &dof_ptr : *(*fe_ptr)->rows_dofs) {
2695  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2696  hint_r = ghost_idx_row_view.emplace_hint(
2697  hint_r, dof_ptr->getNumeredDofEntityPtr());
2698  }
2699  }
2700  }
2701 
2702  // columns
2703  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2704  auto hint_c = ghost_idx_col_view.begin();
2705  for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2706  for (auto &dof_ptr : *(*fe_range.first)->cols_dofs) {
2707  if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2708  hint_c = ghost_idx_col_view.emplace_hint(
2709  hint_c, dof_ptr->getNumeredDofEntityPtr());
2710  }
2711  }
2712  }
2713  }
2714 
2715  int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2716  int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2717 
2718  NumeredDofEntity_multiIndex_uid_view_ordered *ghost_idx_view[2] = {
2719  &ghost_idx_row_view, &ghost_idx_col_view};
2720  NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2721  &p_miit->numeredDofsRows->get<Unique_mi_tag>(),
2722  &p_miit->numeredDofsCols->get<Unique_mi_tag>()};
2723 
2724  int loop_size = 2;
2725  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2726  loop_size = 1;
2727  }
2728 
2729  // set local ghost dofs indices
2730  for (int ss = 0; ss != loop_size; ++ss) {
2731  for (auto &gid : *ghost_idx_view[ss]) {
2732  NumeredDofEntityByUId::iterator dof =
2733  dof_by_uid_no_const[ss]->find(gid->getGlobalUniqueId());
2734  if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2735  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2736  "inconsistent data, ghost dof already set");
2737  bool success = dof_by_uid_no_const[ss]->modify(
2738  dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2739  if (PetscUnlikely(!success))
2740  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2741  "modification unsuccessful");
2742  (*nb_ghost_dofs[ss])++;
2743  }
2744  }
2745  if (loop_size == 1) {
2746  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2747  }
2748  }
2749 
2750  if (verb > QUIET) {
2751  std::ostringstream ss;
2752  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2753  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2754  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2755  << p_miit->getNbLocalDofsCol() << std::endl;
2756  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2757  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2758  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2759  << p_miit->getNbLocalDofsRow() << std::endl;
2760  if (verb > VERBOSE) {
2761  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2762  p_miit->numeredDofsCols->begin();
2763  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2764  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2765  continue;
2766  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2767  continue;
2768  ss << *(*miit_dd_col) << std::endl;
2769  }
2770  }
2771  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2772  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2773  }
2774 
2777 }
2778 
2781  int verb) {
2782  MoFEM::Interface &m_field = cOre;
2783  const Problem_multiIndex *problems_ptr;
2785 
2787  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2788  "partition of problem not build");
2790  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2791  "partitions finite elements not build");
2792 
2793  // get problem pointer
2794  typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2795  // find p_miit
2796  CHKERR m_field.get_problems(&problems_ptr);
2797  ProblemsByName &problems_set =
2798  const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2799  ProblemsByName::iterator p_miit = problems_set.find(name);
2800 
2801  // get reference to number of ghost dofs
2802  // get number of local dofs
2803  DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2804  &(p_miit->nbGhostDofsCol)};
2805  DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2806  for (int ss = 0; ss != 2; ++ss) {
2807  (*nb_ghost_dofs[ss]) = 0;
2808  }
2809 
2810  // do work if more than one processor
2811  if (m_field.get_comm_size() > 1) {
2812  // determine if rows on columns are different from dofs on rows
2813  int loop_size = 2;
2814  if (p_miit->numeredDofsCols == p_miit->numeredDofsRows) {
2815  loop_size = 1;
2816  }
2817 
2818  typedef decltype(p_miit->numeredDofsRows) NumbDofTypeSharedPtr;
2819  NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredDofsRows,
2820  p_miit->numeredDofsCols};
2821 
2822  // interate over dofs on rows and dofs on columns
2823  for (int ss = 0; ss != loop_size; ++ss) {
2824 
2825  // create dofs view by uid
2826  auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2827 
2828  std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2829  ghost_idx_view.reserve(std::distance(r.first, r.second));
2830  for (; r.first != r.second; ++r.first)
2831  ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2832 
2833  auto cmp = [](auto a, auto b) {
2834  return (*a)->getGlobalUniqueId() < (*b)->getGlobalUniqueId();
2835  };
2836  sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2837 
2838  // intare over dofs which have negative local index
2839  for (auto gid_it : ghost_idx_view) {
2840  bool success = numered_dofs[ss]->modify(
2841  gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2842  if (!success)
2843  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2844  "modification unsuccessful");
2845  ++(*nb_ghost_dofs[ss]);
2846  }
2847  }
2848  if (loop_size == 1) {
2849  (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2850  }
2851  }
2852 
2853  if (verb > QUIET) {
2854  std::ostringstream ss;
2855  ss << "partition_ghost_row_dofs: rank = " << m_field.get_comm_rank()
2856  << " FEs row ghost dofs " << *p_miit << " Nb. row ghost dof "
2857  << p_miit->getNbGhostDofsRow() << " Nb. local dof "
2858  << p_miit->getNbLocalDofsRow() << std::endl;
2859  ss << "partition_ghost_col_dofs: rank = " << m_field.get_comm_rank()
2860  << " FEs col ghost dofs " << *p_miit << " Nb. col ghost dof "
2861  << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2862  << p_miit->getNbLocalDofsCol() << std::endl;
2863  if (verb > VERBOSE) {
2864  NumeredDofEntity_multiIndex::iterator miit_dd_col =
2865  p_miit->numeredDofsCols->begin();
2866  for (; miit_dd_col != p_miit->numeredDofsCols->end(); miit_dd_col++) {
2867  if ((*miit_dd_col)->pArt == (unsigned int)m_field.get_comm_rank())
2868  continue;
2869  if ((*miit_dd_col)->petscLocalDofIdx == (DofIdx)-1)
2870  continue;
2871  ss << *(*miit_dd_col) << std::endl;
2872  }
2873  }
2874  PetscSynchronizedPrintf(m_field.get_comm(), ss.str().c_str());
2875  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2876  }
2877 
2880 }
2881 
2883  const std::string fe_name,
2884  EntityHandle *meshset) const {
2885  MoFEM::Interface &m_field = cOre;
2886  const Problem *problem_ptr;
2888  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2889  CHKERR m_field.get_problem(prb_name, &problem_ptr);
2890  auto fit =
2892  .lower_bound(fe_name);
2893  auto hi_fe_it =
2895  .upper_bound(fe_name);
2896  std::vector<EntityHandle> fe_vec;
2897  fe_vec.reserve(std::distance(fit, hi_fe_it));
2898  for (; fit != hi_fe_it; fit++)
2899  fe_vec.push_back(fit->get()->getEnt());
2900  CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2901  fe_vec.size());
2903 }
2904 
2907  const std::string fe_name,
2908  PetscLayout *layout) const {
2909  MoFEM::Interface &m_field = cOre;
2910  const Problem *problem_ptr;
2912  CHKERR m_field.get_problem(name, &problem_ptr);
2913  CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2914  fe_name, layout);
2916 }
2917 
2919  const std::string problem_name, const std::string field_name,
2920  const Range ents, const int lo_coeff, const int hi_coeff, int verb,
2921  const bool debug) {
2922 
2923  MoFEM::Interface &m_field = cOre;
2925 
2926  const Problem *prb_ptr;
2927  CHKERR m_field.get_problem(problem_name, &prb_ptr);
2928 
2929  decltype(prb_ptr->numeredDofsRows) numered_dofs[2] = {
2930  prb_ptr->numeredDofsRows, nullptr};
2931  if (prb_ptr->numeredDofsRows != prb_ptr->numeredDofsCols)
2932  numered_dofs[1] = prb_ptr->numeredDofsCols;
2933 
2934  int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2935  int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2936  int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2937 
2938  for (int s = 0; s != 2; ++s)
2939  if (numered_dofs[s]) {
2940 
2941  typedef multi_index_container<
2942 
2943  NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2944 
2945  >
2946  NumeredDofEntity_it_view_multiIndex;
2947 
2948  NumeredDofEntity_it_view_multiIndex dofs_it_view;
2949 
2950  // Set -1 to global and local dofs indices
2951  for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2952  ++pit) {
2953  auto lo =
2954  numered_dofs[s]
2956  .lower_bound(boost::make_tuple(field_name, pit->first, 0));
2957  auto hi = numered_dofs[s]
2959  .lower_bound(boost::make_tuple(field_name, pit->second,
2961  for (; lo != hi; ++lo)
2962  if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2963  (*lo)->getDofCoeffIdx() <= hi_coeff)
2964  dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2965  }
2966 
2967  if (verb >= VERY_NOISY) {
2968  for (auto &dof : dofs_it_view) {
2969  std::ostringstream ss;
2970  ss << **dof;
2971  PetscSynchronizedPrintf(m_field.get_comm(), "%s\n", ss.str().c_str());
2972  }
2973  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
2974  }
2975 
2976  // set negative index
2977  auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2978  for (auto dit : dofs_it_view) {
2979  bool success = numered_dofs[s]->modify(dit, mod);
2980  if (!success)
2981  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2982  "can not set negative indices");
2983  }
2984 
2985  // create weak view
2986  std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2987  dosf_weak_view.reserve(dofs_it_view.size());
2988  for (auto dit : dofs_it_view)
2989  dosf_weak_view.push_back(*dit);
2990 
2991  if (verb >= NOISY)
2992  PetscSynchronizedPrintf(
2993  m_field.get_comm(),
2994  "Number of DOFs in multi-index %d and to delete %d\n",
2995  numered_dofs[s]->size(), dofs_it_view.size());
2996 
2997  // erase dofs from problem
2998  for (auto weak_dit : dosf_weak_view)
2999  if (auto dit = weak_dit.lock()) {
3000  numered_dofs[s]->erase(dit->getGlobalUniqueId());
3001  }
3002 
3003  if (verb >= NOISY)
3004  PetscSynchronizedPrintf(
3005  m_field.get_comm(),
3006  "Number of DOFs in multi-index after delete %d\n",
3007  numered_dofs[s]->size());
3008 
3009  // get current number of ghost dofs
3010  int nb_local_dofs = 0;
3011  int nb_ghost_dofs = 0;
3012  for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
3013  dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
3014  if ((*dit)->getPetscLocalDofIdx() >= 0 &&
3015  (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
3016  ++nb_local_dofs;
3017  else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
3018  ++nb_ghost_dofs;
3019  else
3020  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Imposible case");
3021  }
3022 
3023  // get indices
3024  auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
3025  const int nb_dofs = numered_dofs[s]->size();
3026  indices.clear();
3027  indices.reserve(nb_dofs);
3028  for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
3029  dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3030  bool add = true;
3031  if (only_local) {
3032  if ((*dit)->getPetscLocalDofIdx() < 0 ||
3033  (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
3034  add = false;
3035  }
3036  }
3037  if (add)
3038  indices.push_back(decltype(tag)::get_index(dit));
3039  }
3040  };
3041 
3042  auto get_indices_by_uid = [&](auto tag, auto &indices) {
3043  const int nb_dofs = numered_dofs[s]->size();
3044  indices.clear();
3045  indices.reserve(nb_dofs);
3046  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3047  ++dit)
3048  indices.push_back(decltype(tag)::get_index(dit));
3049  };
3050 
3051  auto concatenate_dofs = [&](auto tag, auto &indices,
3052  const auto local_only) {
3054  get_indices_by_tag(tag, indices, local_only);
3055 
3056  AO ao;
3057  if (local_only)
3058  CHKERR AOCreateMapping(m_field.get_comm(), indices.size(),
3059  &*indices.begin(), PETSC_NULL, &ao);
3060  else
3061  CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
3062  &*indices.begin(), PETSC_NULL, &ao);
3063 
3064  get_indices_by_uid(tag, indices);
3065  CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3066  CHKERR AODestroy(&ao);
3068  };
3069 
3070  // set indices index
3071  auto set_concatinated_indices = [&]() {
3072  std::vector<int> global_indices;
3073  std::vector<int> local_indices;
3075  CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3076  CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3077  auto gi = global_indices.begin();
3078  auto li = local_indices.begin();
3079  for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3080  ++dit) {
3082  (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3083  bool success = numered_dofs[s]->modify(dit, mod);
3084  if (!success)
3085  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3086  "can not set negative indices");
3087  ++gi;
3088  ++li;
3089  }
3091  };
3092  CHKERR set_concatinated_indices();
3093 
3094  MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3095  m_field.get_comm());
3096  *(local_nbdof_ptr[s]) = nb_local_dofs;
3097  *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3098 
3099  if (debug)
3100  for (auto dof : (*numered_dofs[s])) {
3101  if (dof->getPetscGlobalDofIdx() < 0) {
3102  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3103  "Negative global idx");
3104  }
3105  if (dof->getPetscLocalDofIdx() < 0) {
3106  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3107  "Negative local idx");
3108  }
3109  }
3110 
3111  } else {
3112 
3113  *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3114  *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3115  *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3116  }
3117 
3118  if (verb > QUIET) {
3119  PetscSynchronizedPrintf(
3120  m_field.get_comm(),
3121  "removed ents on rank %d from problem %s dofs [ %d / %d local, %d / %d "
3122  "ghost, %d / %d global]\n",
3123  m_field.get_comm_rank(), prb_ptr->getName().c_str(),
3124  prb_ptr->getNbLocalDofsRow(), prb_ptr->getNbLocalDofsCol(),
3125  prb_ptr->getNbGhostDofsRow(), prb_ptr->getNbGhostDofsCol(),
3126  prb_ptr->getNbDofsRow(), prb_ptr->getNbDofsCol());
3127  PetscSynchronizedFlush(m_field.get_comm(), PETSC_STDOUT);
3128  }
3129 
3131 }
3132 
3133 MoFEMErrorCode ProblemsManager::markDofs(const std::string problem_name,
3134  RowColData rc, const Range ents,
3135  std::vector<bool> &marker) {
3136 
3137  Interface &m_field = cOre;
3138  const Problem *problem_ptr;
3140  CHKERR m_field.get_problem(problem_name, &problem_ptr);
3141  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3142  switch (rc) {
3143  case ROW:
3144  dofs = problem_ptr->getNumeredDofsRows();
3145  break;
3146  case COL:
3147  dofs = problem_ptr->getNumeredDofsCols();
3148  default:
3149  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Should be row or column");
3150  }
3151  marker.resize(dofs->size());
3152  marker.clear();
3153  for(auto p = ents.pair_begin();p!=ents.pair_end(); ++p) {
3154  auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
3155  auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
3156  for (; lo != hi; ++lo)
3157  marker[(*lo)->getPetscLocalDofIdx()] = true;
3158  }
3160 }
3161 
3162 } // MOFEM namespace
DofIdx getNbDofsRow() const
MoFEM interface unique ID.
Managing BitRefLevels.
BitFEId getBitFEId() const
static const MOFEMuuid IDD_MOFEMProblemsManager
virtual moab::Interface & get_moab()=0
PetscLogEvent MOFEM_EVENT_ProblemsManager
DofIdx nbLocDofsCol
Local number of DOFs in colIs.
DofIdx getNbGhostDofsCol() const
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
DofIdx nbDofsCol
Global number of DOFs in col.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
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:477
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:544
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:508
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
RowColData
RowColData.
Definition: definitions.h:186
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
const boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredDofsCols() const
get access to numeredDofsCols storing DOFs on cols
virtual MoFEMErrorCode get_dofs(const DofEntity_multiIndex **dofs_ptr) const =0
Get dofs multi index.
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
virtual MoFEMErrorCode get_finite_elements(const FiniteElement_multiIndex **fe_ptr) const =0
Get finite elements multi-index.
DofIdx nbGhostDofsRow
Number of ghost DOFs in row.
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:596
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.
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:298
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()
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const Range ents, std::vector< bool > &marker)
Create vector with marked indices.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
std::string getName() const
virtual MoFEMErrorCode get_fields(const Field_multiIndex **fields_ptr) const =0
Get fields multi-index from database.
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElements
store finite elements
DofIdx getNbLocalDofsCol() const
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredDofsCols
store DOFs on columns for this problem
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
uint128_t UId
Unique Id.
Definition: Types.hpp:41
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
Matrix manager is used to build and partition problems.
const boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredDofsRows() const
get access to numeredDofsRows storing DOFs on rows
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Unique_Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getEntGlobalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_And_EntDofIdx_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, DofIdx, &NumeredDofEntity::getEntDofIdx > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > >, ordered_non_unique< tag< Composite_Name_Ent_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
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