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