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