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