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