v0.15.0
Loading...
Searching...
No Matches
Classes | Namespaces | Macros
ProblemsManager.cpp File Reference

Managing complexities for problem. More...

Go to the source code of this file.

Classes

struct  MoFEM::IdxDataType
 
struct  MoFEM::IdxDataTypePtr
 

Namespaces

namespace  MoFEM
 implementation of Data Operators for Forces and Sources
 

Macros

#define ProblemManagerFunctionBegin
 

Detailed Description

Managing complexities for problem.

Definition in file ProblemsManager.cpp.

Macro Definition Documentation

◆ ProblemManagerFunctionBegin

#define ProblemManagerFunctionBegin
Value:
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_CHANNEL("SYNC"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("SYNC", "ProblemsManager"); \
MOFEM_LOG_TAG("WORLD", "ProblemsManager")
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

Definition at line 8 of file ProblemsManager.cpp.

15 {
16 IdxDataType(const UId uid, const int dof) {
17 bcopy(&uid, dAta, 4 * sizeof(int));
18 dAta[4] = dof;
19 }
20
21private:
22 int dAta[5];
23};
24
25struct IdxDataTypePtr {
26 IdxDataTypePtr(const int *ptr) : pTr(ptr) {}
27 inline int getDofIdx() const {
28 int global_dof = pTr[4];
29 return global_dof;
30 }
31 inline UId getUId() const {
32 unsigned int b0, b1, b2, b3;
33 bcopy(&pTr[0], &b0, sizeof(int));
34 bcopy(&pTr[1], &b1, sizeof(int));
35 bcopy(&pTr[2], &b2, sizeof(int));
36 bcopy(&pTr[3], &b3, sizeof(int));
37 UId uid = static_cast<UId>(b0) | static_cast<UId>(b1) << 8 * sizeof(int) |
38 static_cast<UId>(b2) << 16 * sizeof(int) |
39 static_cast<UId>(b3) << 24 * sizeof(int);
40 return uid;
41 }
42
43private:
44 const int *pTr;
45};
46
48ProblemsManager::query_interface(boost::typeindex::type_index type_index,
49 UnknownInterface **iface) const {
50 *iface = const_cast<ProblemsManager *>(this);
51 return 0;
52}
53
54ProblemsManager::ProblemsManager(const MoFEM::Core &core)
55 : cOre(const_cast<MoFEM::Core &>(core)),
56 buildProblemFromFields(PETSC_FALSE),
57 synchroniseProblemEntities(PETSC_FALSE) {
58 PetscLogEventRegister("ProblemsManager", 0, &MOFEM_EVENT_ProblemsManager);
59}
60
61MoFEMErrorCode ProblemsManager::getOptions() {
62 MoFEM::Interface &m_field = cOre;
64 PetscOptionsBegin(m_field.get_comm(), "", "Problem manager", "none");
65 {
66 CHKERR PetscOptionsBool(
67 "-problem_build_from_fields",
68 "Add DOFs to problem directly from fields not through DOFs on elements",
69 "", buildProblemFromFields, &buildProblemFromFields, NULL);
70 }
71 PetscOptionsEnd();
73}
74
75MoFEMErrorCode ProblemsManager::partitionMesh(
76 const Range &ents, const int dim, const int adj_dim, const int n_parts,
77 Tag *th_vertex_weights, Tag *th_edge_weights, Tag *th_part_weights,
78 int verb, const bool debug) {
79 return static_cast<MoFEM::Interface &>(cOre)
80 .getInterface<CommInterface>()
81 ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
82 th_edge_weights, th_part_weights, verb, debug);
83}
84
85MoFEMErrorCode ProblemsManager::buildProblem(const std::string name,
86 const bool square_matrix,
87 int verb) {
88
89 MoFEM::Interface &m_field = cOre;
91 if (!(cOre.getBuildMoFEM() & (1 << 0)))
92 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
93 if (!(cOre.getBuildMoFEM() & (1 << 1)))
94 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
95 if (!(cOre.getBuildMoFEM() & (1 << 2)))
96 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
97 const Problem *problem_ptr;
98 CHKERR m_field.get_problem(name, &problem_ptr);
99 CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
100 cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
101 // function knows what he is doing
103}
104
105MoFEMErrorCode ProblemsManager::buildProblem(Problem *problem_ptr,
106 const bool square_matrix,
107 int verb) {
108 MoFEM::Interface &m_field = cOre;
109 auto dofs_field_ptr = m_field.get_dofs();
110 auto ents_field_ptr = m_field.get_field_ents();
111 auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
113 PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
114
115 // Note: Only allowed changes on problem_ptr structure which not influence
116 // multi-index.
117
118 if (problem_ptr->getBitRefLevel().none()) {
119 SETERRQ(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
120 problem_ptr->getName().c_str());
121 }
122 CHKERR m_field.clear_problem(problem_ptr->getName());
123
124 // zero finite elements
125 problem_ptr->numeredFiniteElementsPtr->clear();
126
127 DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
128 auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
130 for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
131
132 const auto uid = (*it)->getLocalUniqueId();
133
134 auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
135 for (auto lo = r.first; lo != r.second; ++lo) {
136
137 if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
138 std::array<bool, 2> row_col = {false, false};
139
140 const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
141 const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
142 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
143
144 // if entity is not problem refinement level
145 if ((fe_bit & prb_mask) != fe_bit)
146 continue;
147 if ((fe_bit & prb_bit).none())
148 continue;
149
150 auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
151 auto dit = nb_dofs->lower_bound(uid);
152 decltype(dit) hi_dit;
153 if (dit != nb_dofs->end()) {
154 hi_dit = nb_dofs->upper_bound(
155 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
156 view.insert(dit, hi_dit);
157 row_col[rc] = true;
158 }
159 };
160
161 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
162 add_to_view(dofs_field_ptr, dofs_rows, ROW);
163
164 if (!square_matrix)
165 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
166 add_to_view(dofs_field_ptr, dofs_cols, COL);
167
168 if (square_matrix && row_col[ROW])
169 break;
170 else if (row_col[ROW] && row_col[COL])
171 break;
172 }
173 }
174 }
176 };
177
178 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
179
180 // Add row dofs to problem
181 {
182 // zero rows
183 problem_ptr->nbDofsRow = 0;
184 problem_ptr->nbLocDofsRow = 0;
185 problem_ptr->nbGhostDofsRow = 0;
186
187 // add dofs for rows
188 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
189 hi_miit;
190 hi_miit = dofs_rows.get<0>().end();
191
192 int count_dofs = dofs_rows.get<1>().count(true);
193 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
194 boost::shared_ptr<std::vector<NumeredDofEntity>>(
195 new std::vector<NumeredDofEntity>());
196 problem_ptr->getRowDofsSequence()->push_back(dofs_array);
197 dofs_array->reserve(count_dofs);
198 miit = dofs_rows.get<0>().begin();
199 for (; miit != hi_miit; miit++) {
200 if ((*miit)->getActive()) {
201 dofs_array->emplace_back(*miit);
202 dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
203 }
204 }
205 auto hint = problem_ptr->numeredRowDofsPtr->end();
206 for (auto &v : *dofs_array) {
207 hint = problem_ptr->numeredRowDofsPtr->emplace_hint(hint, dofs_array, &v);
208 }
209 }
210
211 // Add col dofs to problem
212 if (!square_matrix) {
213 // zero cols
214 problem_ptr->nbDofsCol = 0;
215 problem_ptr->nbLocDofsCol = 0;
216 problem_ptr->nbGhostDofsCol = 0;
217
218 // add dofs for cols
219 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
220 hi_miit;
221 hi_miit = dofs_cols.get<0>().end();
222
223 int count_dofs = 0;
224 miit = dofs_cols.get<0>().begin();
225 for (; miit != hi_miit; miit++) {
226 if (!(*miit)->getActive()) {
227 continue;
228 }
229 count_dofs++;
230 }
231
232 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
233 boost::shared_ptr<std::vector<NumeredDofEntity>>(
234 new std::vector<NumeredDofEntity>());
235 problem_ptr->getColDofsSequence()->push_back(dofs_array);
236 dofs_array->reserve(count_dofs);
237 miit = dofs_cols.get<0>().begin();
238 for (; miit != hi_miit; miit++) {
239 if (!(*miit)->getActive()) {
240 continue;
241 }
242 dofs_array->emplace_back(*miit);
243 dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
244 }
245 auto hint = problem_ptr->numeredColDofsPtr->end();
246 for (auto &v : *dofs_array) {
247 hint = problem_ptr->numeredColDofsPtr->emplace_hint(hint, dofs_array, &v);
248 }
249 } else {
250 problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
251 problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
252 problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
253 }
254
255 // job done, some debugging and postprocessing
256 if (verb >= VERBOSE) {
257 MOFEM_LOG("SYNC", Sev::verbose)
258 << problem_ptr->getName() << " Nb. local dofs "
259 << problem_ptr->numeredRowDofsPtr->size() << " by "
260 << problem_ptr->numeredColDofsPtr->size();
261 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
262 }
263
264 if (verb >= NOISY) {
265 MOFEM_LOG("SYNC", Sev::noisy)
266 << "FEs row dofs " << *problem_ptr << " Nb. row dof "
267 << problem_ptr->getNbDofsRow();
268 for (auto &miit : *problem_ptr->numeredRowDofsPtr)
269 MOFEM_LOG("SYNC", Sev::noisy) << *miit;
270
271 MOFEM_LOG("SYNC", Sev::noisy)
272 << "FEs col dofs " << *problem_ptr << " Nb. col dof "
273 << problem_ptr->getNbDofsCol();
274 for (auto &miit : *problem_ptr->numeredColDofsPtr)
275 MOFEM_LOG("SYNC", Sev::noisy) << *miit;
276 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
277 }
278
279 cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
280 // uses this function knows
281 // what he is doing
282
283 PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
284
286}
287
288MoFEMErrorCode ProblemsManager::buildProblemOnDistributedMesh(
289 const std::string name, const bool square_matrix, int verb) {
290 MoFEM::Interface &m_field = cOre;
292
293 if (!((cOre.getBuildMoFEM()) & Core::BUILD_FIELD))
294 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
295 if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
296 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
297 if (!((cOre.getBuildMoFEM()) & Core::BUILD_ADJ))
298 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
299
300 const Problem *problem_ptr;
301 CHKERR m_field.get_problem(name, &problem_ptr);
302 CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
303 square_matrix, verb);
304
305 cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM;
306 cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
307
309}
310
311MoFEMErrorCode ProblemsManager::buildProblemOnDistributedMesh(
312 Problem *problem_ptr, const bool square_matrix, int verb) {
313 MoFEM::Interface &m_field = cOre;
314 auto fields_ptr = m_field.get_fields();
315 auto fe_ptr = m_field.get_finite_elements();
316 auto dofs_field_ptr = m_field.get_dofs();
318 PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
319
320 // clear data structures
321 CHKERR m_field.clear_problem(problem_ptr->getName());
322
323 CHKERR getOptions();
324
325 if (problem_ptr->getBitRefLevel().none())
326 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
327 "problem <%s> refinement level not set",
328 problem_ptr->getName().c_str());
329
330 int loop_size = 2;
331 if (square_matrix) {
332 loop_size = 1;
333 problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
334 } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
335 problem_ptr->numeredColDofsPtr =
336 boost::shared_ptr<NumeredDofEntity_multiIndex>(
338 }
339
340 const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
341 const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
342
343 // // get rows and cols dofs view based on data on elements
344 DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
345
346 // Add DOFs to problem by visiting all elements and adding DOFs from
347 // elements to the problem
348 if (buildProblemFromFields == PETSC_FALSE) {
349
350 auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
351 auto ents_field_ptr = m_field.get_field_ents();
352 auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
354 for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
355 ++it) {
356
357 const auto uid = (*it)->getLocalUniqueId();
358
359 auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
360 for (auto lo = r.first; lo != r.second; ++lo) {
361
362 if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
363 std::array<bool, 2> row_col = {false, false};
364
365 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
366
367 // if entity is not problem refinement level
368 if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
369
370 auto add_to_view = [&](auto dofs, auto &view, auto rc) {
371 auto dit = dofs->lower_bound(uid);
372 decltype(dit) hi_dit;
373 if (dit != dofs->end()) {
374 hi_dit = dofs->upper_bound(
375 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
376 view.insert(dit, hi_dit);
377 row_col[rc] = true;
378 }
379 };
380
381 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
382 add_to_view(dofs_field_ptr, dofs_rows, ROW);
383
384 if (!square_matrix)
385 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
386 add_to_view(dofs_field_ptr, dofs_cols, COL);
387
388 if (square_matrix && row_col[ROW])
389 break;
390 else if (row_col[ROW] && row_col[COL])
391 break;
392 }
393 }
394 }
395 }
397 };
398
399 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
400 }
401
402 // Add DOFS to the problem by searching all the filedes, and adding to
403 // problem owned or shared DOFs
404 if (buildProblemFromFields == PETSC_TRUE) {
405 // Get fields IDs on elements
406 BitFieldId fields_ids_row, fields_ids_col;
407 for (auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
408 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
409 fields_ids_row |= fit->get()->getBitFieldIdRow();
410 fields_ids_col |= fit->get()->getBitFieldIdCol();
411 }
412 }
413 // Get fields DOFs
414 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
415 if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
416
417 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
418 FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
419 auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
420 FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
421
422 for (; dit != hi_dit; dit++) {
423
424 const int owner_proc = dit->get()->getOwnerProc();
425 if (owner_proc != m_field.get_comm_rank()) {
426 const unsigned char pstatus = dit->get()->getPStatus();
427 if (pstatus == 0) {
428 continue;
429 }
430 }
431
432 const auto &dof_bit = (*dit)->getBitRefLevel();
433 // if entity is not problem refinement level
434 if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
435
436 if ((fit->get()->getId() & fields_ids_row).any()) {
437 dofs_rows.insert(*dit);
438 }
439 if (!square_matrix) {
440 if ((fit->get()->getId() & fields_ids_col).any()) {
441 dofs_cols.insert(*dit);
442 }
443 }
444 }
445 }
446 }
447 }
448 }
449
450 if (synchroniseProblemEntities) {
451 // Get fields IDs on elements
452 BitFieldId fields_ids_row, fields_ids_col;
453 BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
454 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
455 fit != fe_ptr->end(); fit++) {
456 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
457 fields_ids_row |= fit->get()->getBitFieldIdRow();
458 fields_ids_col |= fit->get()->getBitFieldIdCol();
459 }
460 }
461
462 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
463 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
464 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
465 hi_miit;
466 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
467 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
468 Range ents_to_synchronise;
469 for (; miit != hi_miit; ++miit) {
470 if (miit->get()->getEntDofIdx() != 0)
471 continue;
472 ents_to_synchronise.insert(miit->get()->getEnt());
473 }
474 Range tmp_ents = ents_to_synchronise;
475 CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
476 ents_to_synchronise, nullptr, verb);
477 ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
478 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
479 if ((fit->get()->getId() & *fields_ids[ss]).any()) {
480 const auto bit_number = (*fit)->getBitNumber();
481 for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
482 pit != ents_to_synchronise.pair_end(); ++pit) {
483 const auto f = pit->first;
484 const auto s = pit->second;
485 const auto lo_uid =
486 FieldEntity::getLocalUniqueIdCalculate(bit_number, f);
487 const auto hi_uid =
488 FieldEntity::getLocalUniqueIdCalculate(bit_number, s);
489
490 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
491 auto hi_dit =
492 dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
493
494 dofs_ptr[ss]->insert(dit, hi_dit);
495 }
496 }
497 }
498 }
499 }
500
501 // add dofs for rows and cols and set ownership
502 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
503 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
504 problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
505 int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
506 int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
507 &problem_ptr->nbLocDofsCol};
508 int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
509 &problem_ptr->nbGhostDofsCol};
510 for (int ss = 0; ss < 2; ss++) {
511 *(nbdof_ptr[ss]) = 0;
512 *(local_nbdof_ptr[ss]) = 0;
513 *(ghost_nbdof_ptr[ss]) = 0;
514 }
515
516 // Loop over dofs on rows and columns and add to multi-indices in dofs
517 // problem structure, set partition for each dof
518 int nb_local_dofs[] = {0, 0};
519 for (int ss = 0; ss < loop_size; ss++) {
520 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
521 hi_miit;
522 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
523 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
524 for (; miit != hi_miit; miit++) {
525 int owner_proc = (*miit)->getOwnerProc();
526 if (owner_proc == m_field.get_comm_rank()) {
527 nb_local_dofs[ss]++;
528 }
529 }
530 }
531
532 // get layout
533 int start_ranges[2], end_ranges[2];
534 for (int ss = 0; ss != loop_size; ss++) {
535 PetscLayout layout;
536 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
537 CHKERR PetscLayoutSetBlockSize(layout, 1);
538 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
539 CHKERR PetscLayoutSetUp(layout);
540 CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
541 CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
542 CHKERR PetscLayoutDestroy(&layout);
543 }
544 if (square_matrix) {
545 nbdof_ptr[1] = nbdof_ptr[0];
546 nb_local_dofs[1] = nb_local_dofs[0];
547 start_ranges[1] = start_ranges[0];
548 end_ranges[1] = end_ranges[0];
549 }
550
551 // set local and global indices on own dofs
552 const size_t idx_data_type_size = sizeof(IdxDataType);
553 const size_t data_block_size = idx_data_type_size / sizeof(int);
554
555 if (sizeof(IdxDataType) % sizeof(int)) {
556 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
557 }
558 std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
559 m_field.get_comm_size()),
560 ids_data_packed_cols(m_field.get_comm_size());
561
562 // Loop over dofs on this processor and prepare those dofs to send on
563 // another proc
564 for (int ss = 0; ss != loop_size; ++ss) {
565
566 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
567 hi_miit;
568 hi_miit = dofs_ptr[ss]->get<0>().end();
569
570 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
571 boost::shared_ptr<std::vector<NumeredDofEntity>>(
572 new std::vector<NumeredDofEntity>());
573 int nb_dofs_to_add = 0;
574 miit = dofs_ptr[ss]->get<0>().begin();
575 for (; miit != hi_miit; ++miit) {
576 // Only set global idx for dofs on this processor part
577 if (!(miit->get()->getActive()))
578 continue;
579 ++nb_dofs_to_add;
580 }
581 dofs_array->reserve(nb_dofs_to_add);
582 if (ss == 0) {
583 problem_ptr->getRowDofsSequence()->push_back(dofs_array);
584 } else {
585 problem_ptr->getColDofsSequence()->push_back(dofs_array);
586 }
587
588 int &local_idx = *local_nbdof_ptr[ss];
589 miit = dofs_ptr[ss]->get<0>().begin();
590 for (; miit != hi_miit; ++miit) {
591
592 // Only set global idx for dofs on this processor part
593 if (!(miit->get()->getActive()))
594 continue;
595
596 dofs_array->emplace_back(*miit);
597
598 int owner_proc = dofs_array->back().getOwnerProc();
599 if (owner_proc < 0) {
600 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
601 "data inconsistency");
602 }
603
604 if (owner_proc != m_field.get_comm_rank()) {
605 dofs_array->back().pArt = owner_proc;
606 dofs_array->back().dofIdx = -1;
607 dofs_array->back().petscGloablDofIdx = -1;
608 dofs_array->back().petscLocalDofIdx = -1;
609 } else {
610
611 // Set part and indexes
612 int glob_idx = start_ranges[ss] + local_idx;
613 dofs_array->back().pArt = owner_proc;
614 dofs_array->back().dofIdx = glob_idx;
615 dofs_array->back().petscGloablDofIdx = glob_idx;
616 dofs_array->back().petscLocalDofIdx = local_idx;
617 local_idx++;
618
619 unsigned char pstatus = dofs_array->back().getPStatus();
620 // check id dof is shared, if that is a case global idx added to data
621 // structure and is sended to other processors
622 if (pstatus > 0) {
623
624 for (int proc = 0;
625 proc < MAX_SHARING_PROCS &&
626 -1 != dofs_array->back().getSharingProcsPtr()[proc];
627 proc++) {
628 // make it different for rows and columns when needed
629 if (ss == 0) {
630 ids_data_packed_rows[dofs_array->back()
631 .getSharingProcsPtr()[proc]]
632 .emplace_back(dofs_array->back().getGlobalUniqueId(),
633 glob_idx);
634 } else {
635 ids_data_packed_cols[dofs_array->back()
636 .getSharingProcsPtr()[proc]]
637 .emplace_back(dofs_array->back().getGlobalUniqueId(),
638 glob_idx);
639 }
640 if (!(pstatus & PSTATUS_MULTISHARED)) {
641 break;
642 }
643 }
644 }
645 }
646 }
647
648 auto hint = numered_dofs_ptr[ss]->end();
649 for (auto &v : *dofs_array)
650 hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
651 }
652 if (square_matrix) {
653 local_nbdof_ptr[1] = local_nbdof_ptr[0];
654 }
655
656 int nsends_rows = 0, nsends_cols = 0;
657 // Non zero lengths[i] represent a message to i of length lengths[i].
658 std::vector<int> lengths_rows(m_field.get_comm_size()),
659 lengths_cols(m_field.get_comm_size());
660 lengths_rows.clear();
661 lengths_cols.clear();
662 for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
663 lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
664 lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
665 if (!ids_data_packed_rows[proc].empty())
666 nsends_rows++;
667 if (!ids_data_packed_cols[proc].empty())
668 nsends_cols++;
669 }
670
671 MPI_Status *status;
672 CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
673
674 // Do rows
675 int nrecvs_rows; // number of messages received
676 int *onodes_rows; // list of node-ids from which messages are expected
677 int *olengths_rows; // corresponding message lengths
678 int **rbuf_row; // must bee freed by user
679
680 // make sure it is a PETSc comm
681 MPI_Comm comm;
682 CHKERR PetscCommDuplicate(m_field.get_comm(), &comm, NULL);
683
684 {
685
686 // rows
687
688 // Computes the number of messages a node expects to receive
689 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
690 &nrecvs_rows);
691 // std::cerr << nrecvs_rows << std::endl;
692
693 // Computes info about messages that a MPI-node will receive, including
694 // (from-id,length) pairs for each message.
695 CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
696 &lengths_rows[0], &onodes_rows,
697 &olengths_rows);
698
699 // Gets a unique new tag from a PETSc communicator. All processors that
700 // share the communicator MUST call this routine EXACTLY the same number
701 // of times. This tag should only be used with the current objects
702 // communicator; do NOT use it with any other MPI communicator.
703 int tag_row;
704 CHKERR PetscCommGetNewTag(comm, &tag_row);
705
706 // Allocate a buffer sufficient to hold messages of size specified in
707 // olengths. And post Irecvs on these buffers using node info from onodes
708 MPI_Request *r_waits_row; // must bee freed by user
709 // rbuf has a pointers to messeges. It has size of of nrecvs (number of
710 // messages) +1. In the first index a block is allocated,
711 // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
712
713 CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
714 olengths_rows, &rbuf_row, &r_waits_row);
715 CHKERR PetscFree(onodes_rows);
716
717 MPI_Request *s_waits_row; // status of sens messages
718 CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
719
720 // Send messeges
721 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
722 if (!lengths_rows[proc])
723 continue; // no message to send to this proc
724 CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
725 lengths_rows[proc], // message length
726 MPIU_INT, proc, // to proc
727 tag_row, comm, s_waits_row + kk);
728 kk++;
729 }
730
731 if (nrecvs_rows) {
732 CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
733 }
734 if (nsends_rows) {
735 CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
736 }
737
738 CHKERR PetscFree(r_waits_row);
739 CHKERR PetscFree(s_waits_row);
740 }
741
742 // cols
743 int nrecvs_cols = nrecvs_rows;
744 int *olengths_cols = olengths_rows;
745 PetscInt **rbuf_col = rbuf_row;
746 if (!square_matrix) {
747
748 // Computes the number of messages a node expects to receive
749 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
750 &nrecvs_cols);
751
752 // Computes info about messages that a MPI-node will receive, including
753 // (from-id,length) pairs for each message.
754 int *onodes_cols;
755 CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
756 &lengths_cols[0], &onodes_cols,
757 &olengths_cols);
758
759 // Gets a unique new tag from a PETSc communicator.
760 int tag_col;
761 CHKERR PetscCommGetNewTag(comm, &tag_col);
762
763 // Allocate a buffer sufficient to hold messages of size specified in
764 // olengths. And post Irecvs on these buffers using node info from onodes
765 MPI_Request *r_waits_col; // must bee freed by user
766 CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
767 olengths_cols, &rbuf_col, &r_waits_col);
768 CHKERR PetscFree(onodes_cols);
769
770 MPI_Request *s_waits_col; // status of sens messages
771 CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
772
773 // Send messages
774 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
775 if (!lengths_cols[proc])
776 continue; // no message to send to this proc
777 CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
778 lengths_cols[proc], // message length
779 MPIU_INT, proc, // to proc
780 tag_col, comm, s_waits_col + kk);
781 kk++;
782 }
783
784 if (nrecvs_cols) {
785 CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
786 }
787 if (nsends_cols) {
788 CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
789 }
790
791 CHKERR PetscFree(r_waits_col);
792 CHKERR PetscFree(s_waits_col);
793 }
794
795 CHKERR PetscCommDestroy(&comm);
796 CHKERR PetscFree(status);
797
798 DofEntity_multiIndex_global_uid_view dofs_glob_uid_view;
799 auto hint = dofs_glob_uid_view.begin();
800 for (auto dof : *m_field.get_dofs())
801 dofs_glob_uid_view.emplace_hint(hint, dof);
802
803 // set values received from other processors
804 for (int ss = 0; ss != loop_size; ++ss) {
805
806 int nrecvs;
807 int *olengths;
808 int **data_procs;
809 if (ss == 0) {
810 nrecvs = nrecvs_rows;
811 olengths = olengths_rows;
812 data_procs = rbuf_row;
813 } else {
814 nrecvs = nrecvs_cols;
815 olengths = olengths_cols;
816 data_procs = rbuf_col;
817 }
818
819 UId uid;
820 for (int kk = 0; kk != nrecvs; ++kk) {
821 int len = olengths[kk];
822 int *data_from_proc = data_procs[kk];
823 for (int dd = 0; dd < len; dd += data_block_size) {
824 uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
825 auto ddit = dofs_glob_uid_view.find(uid);
826 const auto owner_proc = FieldEntity::getOwnerFromUniqueId(uid);
827
828 if (owner_proc == m_field.get_comm_rank() &&
829 PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
830
831#ifndef NDEBUG
832 auto ents_field_ptr = m_field.get_field_ents();
833 MOFEM_LOG("SELF", Sev::error)
834 << "DOF is shared or multishared between processors. For example "
835 "if order of field on given entity is set in inconsistently, "
836 "has different value on two processor, error such as this is "
837 "triggered";
838
839 MOFEM_LOG("SELF", Sev::error) << "UId " << uid << " is not found";
840 MOFEM_LOG("SELF", Sev::error)
841 << "Problematic UId owner proc is " << owner_proc;
842 const auto uid_handle = FieldEntity::getHandleFromUniqueId(uid);
843 MOFEM_LOG("SELF", Sev::error)
844 << "Problematic UId entity owning processor handle is "
845 << uid_handle << " entity type "
846 << moab::CN::EntityTypeName(type_from_handle(uid_handle));
847 const auto uid_bit_number =
848 FieldEntity::getFieldBitNumberFromUniqueId(uid);
849 MOFEM_LOG("SELF", Sev::error)
850 << "Problematic UId field is "
851 << m_field.get_field_name(
852 field_bit_from_bit_number(uid_bit_number));
853
855 field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
856 auto fe_it = field_view.find(FieldEntity::getGlobalUniqueIdCalculate(
857 owner_proc, uid_bit_number, uid_handle));
858 if (fe_it == field_view.end()) {
859 MOFEM_LOG("SELF", Sev::error)
860 << "Also, no field entity in database for given global UId";
861 } else {
862 MOFEM_LOG("SELF", Sev::error) << "Field entity in databse exist "
863 "(but have no DOF wih give UId";
864 MOFEM_LOG("SELF", Sev::error) << **fe_it;
865
866 // Save file with missing entity
867 auto error_file_name =
868 "error_with_missing_entity_" +
869 boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
870 ".vtk";
871 MOFEM_LOG("SELF", Sev::error)
872 << "Look to file < " << error_file_name
873 << " > it contains entity with missing DOF.";
874
875 auto tmp_msh = get_temp_meshset_ptr(m_field.get_moab());
876 const auto local_fe_ent = (*fe_it)->getEnt();
877 CHKERR m_field.get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
878 CHKERR m_field.get_moab().write_file(error_file_name.c_str(), "VTK",
879 "", tmp_msh->get_ptr(), 1);
880 }
881#endif
882
883 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
884 "DOF with global UId not found (Compile code in Debug to "
885 "learn more about problem");
886 }
887
888 if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
889 continue;
890 }
891
892 auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
893
894 if (dit != numered_dofs_ptr[ss]->end()) {
895
896 int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
897#ifndef NDEBUG
898 if (PetscUnlikely(global_idx < 0))
899 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
900 "received negative dof");
901#endif
902 bool success;
903 success = numered_dofs_ptr[ss]->modify(
904 dit, NumeredDofEntity_mofem_index_change(global_idx));
905
906#ifndef NDEBUG
907 if (PetscUnlikely(!success))
908 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
909 "modification unsuccessful");
910#endif
911 success = numered_dofs_ptr[ss]->modify(
912 dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
913 global_idx));
914#ifndef NDEBUG
915 if (PetscUnlikely(!success))
916 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
917 "modification unsuccessful");
918#endif
919 };
920
921#ifndef NDEBUG
922 if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
923
924 // Dof is shared on this processor, however there is no element
925 // which have this dof. If DOF is not shared and received from other
926 // processor, but not marked as a shared on other that means that is
927 // data inconstancy and error should be thrown.
928
929 std::ostringstream zz;
930 zz << **ddit << std::endl;
931 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
932 "data inconsistency, dofs is not shared, but received "
933 "from other proc\n"
934 "%s",
935 zz.str().c_str());
936 }
937#endif
938 }
939 }
940 }
941
942 if (square_matrix) {
943 (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
944 (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
945 }
946
947 CHKERR PetscFree(olengths_rows);
948 CHKERR PetscFree(rbuf_row[0]);
949 CHKERR PetscFree(rbuf_row);
950 if (!square_matrix) {
951 CHKERR PetscFree(olengths_cols);
952 CHKERR PetscFree(rbuf_col[0]);
953 CHKERR PetscFree(rbuf_col);
954 }
955
956 if (square_matrix) {
957 if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
958 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
959 "data inconsistency for square_matrix %ld!=%ld",
960 numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
961 }
962 if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
963 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
964 "data inconsistency for square_matrix");
965 }
966 }
967
968 CHKERR printPartitionedProblem(problem_ptr, verb);
969 CHKERR debugPartitionedProblem(problem_ptr, verb);
970
971 PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
972
974}
975
976MoFEMErrorCode ProblemsManager::buildSubProblem(
977 const std::string out_name,
978
979 const std::vector<std::string> &fields_row,
980 const std::vector<std::string> &fields_col,
981
982 const std::string main_problem, const bool square_matrix,
983
984 const map<std::string, boost::shared_ptr<Range>> *entityMapRow,
985 const map<std::string, boost::shared_ptr<Range>> *entityMapCol,
986
987 int verb) {
988 MoFEM::Interface &m_field = cOre;
989 auto problems_ptr = m_field.get_problems();
991
992 CHKERR m_field.clear_problem(out_name);
993
994 // get reference to all problems
995 using ProblemByName = decltype(problems_ptr->get<Problem_mi_tag>());
996 auto &problems_by_name =
997 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
998
999 // get iterators to out problem, i.e. build problem
1000 auto out_problem_it = problems_by_name.find(out_name);
1001 if (out_problem_it == problems_by_name.end()) {
1002 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1003 "subproblem with name < %s > not defined (top tip check spelling)",
1004 out_name.c_str());
1005 }
1006 // get iterator to main problem, i.e. out problem is subproblem of main
1007 // problem
1008 auto main_problem_it = problems_by_name.find(main_problem);
1009 if (main_problem_it == problems_by_name.end()) {
1010 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1011 "problem of subproblem with name < %s > not defined (top tip "
1012 "check spelling)",
1013 main_problem.c_str());
1014 }
1015
1016 // get dofs for row & columns for out problem,
1017 boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1018 out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1019 // get dofs for row & columns for main problem
1020 boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1021 main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1022 // get local indices counter
1023 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1024 &out_problem_it->nbLocDofsCol};
1025 // get global indices counter
1026 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1027
1028 // set number of ghost nodes to zero
1029 {
1030 out_problem_it->nbGhostDofsRow = 0;
1031 out_problem_it->nbGhostDofsCol = 0;
1032 }
1033
1034 // put rows & columns field names in array
1035 std::vector<std::string> fields[] = {fields_row, fields_col};
1036 const map<std::string, boost::shared_ptr<Range>> *entityMap[] = {
1037 entityMapRow, entityMapCol};
1038
1039 // make data structure fos sub-problem data
1040 out_problem_it->subProblemData =
1041 boost::make_shared<Problem::SubProblemData>();
1042
1043 // Loop over rows and columns
1044 for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1045
1046 // reset dofs and columns counters
1047 (*nb_local_dofs[ss]) = 0;
1048 (*nb_dofs[ss]) = 0;
1049 // clear arrays
1050 out_problem_dofs[ss]->clear();
1051
1052 // If DOFs are cleared clear finite elements too.
1053 out_problem_it->numeredFiniteElementsPtr->clear();
1054
1055 // get dofs by field name and insert them in out problem multi-indices
1056 for (auto field : fields[ss]) {
1057
1058 // Following reserve memory in sequences, only two allocations are here,
1059 // once for array of objects, next for array of shared pointers
1060
1061 // aliased sequence of pointer is killed with element
1062 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1063 boost::make_shared<std::vector<NumeredDofEntity>>();
1064 // reserve memory for field dofs
1065 if (!ss)
1066 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1067 else
1068 out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1069
1070 // create elements objects
1071 auto bit_number = m_field.get_field_bit_number(field);
1072
1073 auto add_dit_to_dofs_array = [&](auto &dit) {
1074 if (dit->get()->getPetscGlobalDofIdx() >= 0)
1075 dofs_array->emplace_back(
1076 dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1077 dit->get()->getPetscGlobalDofIdx(),
1078 dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1079 };
1080
1081 auto get_dafult_dof_range = [&]() {
1082 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1083 FieldEntity::getLoBitNumberUId(bit_number));
1084 auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1085 FieldEntity::getHiBitNumberUId(bit_number));
1086 return std::make_pair(dit, hi_dit);
1087 };
1088
1089 auto check = [&](auto &field) -> boost::shared_ptr<Range> {
1090 if (entityMap[ss]) {
1091 auto mit = entityMap[ss]->find(field);
1092 if (mit != entityMap[ss]->end()) {
1093 return mit->second;
1094 } else {
1095 return nullptr;
1096 }
1097 } else {
1098 return nullptr;
1099 }
1100 };
1101
1102 if (auto r_ptr = check(field)) {
1103 for (auto p = r_ptr->pair_begin(); p != r_ptr->pair_end(); ++p) {
1104 const auto lo_ent = p->first;
1105 const auto hi_ent = p->second;
1106 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1107 DofEntity::getLoFieldEntityUId(bit_number, lo_ent));
1108 auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1109 DofEntity::getHiFieldEntityUId(bit_number, hi_ent));
1110 dofs_array->reserve(std::distance(dit, hi_dit));
1111 for (; dit != hi_dit; dit++) {
1112 add_dit_to_dofs_array(dit);
1113 }
1114 }
1115 } else {
1116 auto [dit, hi_dit] = get_dafult_dof_range();
1117 dofs_array->reserve(std::distance(dit, hi_dit));
1118 for (; dit != hi_dit; dit++)
1119 add_dit_to_dofs_array(dit);
1120 }
1121
1122 // fill multi-index
1123 auto hint = out_problem_dofs[ss]->end();
1124 for (auto &v : *dofs_array)
1125 hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1126 }
1127 // Set local indexes
1128 {
1129 auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1130 auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1131 for (; dit != hi_dit; dit++) {
1132 int idx = -1; // if dof is not part of partition, set local index to -1
1133 if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1134 idx = (*nb_local_dofs[ss])++;
1135 }
1136 bool success = out_problem_dofs[ss]->modify(
1137 out_problem_dofs[ss]->project<0>(dit),
1138 NumeredDofEntity_local_idx_change(idx));
1139 if (!success) {
1140 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1141 "operation unsuccessful");
1142 }
1143 };
1144 }
1145 // Set global indexes, compress global indices
1146 {
1147 auto dit =
1148 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1149 auto hi_dit =
1150 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1151 out_problem_dofs[ss]->size());
1152 const int nb = std::distance(dit, hi_dit);
1153 // get main problem global indices
1154 std::vector<int> main_indices(nb);
1155 for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1156 *it = dit->get()->getPetscGlobalDofIdx();
1157 }
1158 // create is with global dofs
1159 IS is;
1160 CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1161 PETSC_USE_POINTER, &is);
1162 // create map form main problem global indices to out problem global
1163 // indices
1164 AO ao;
1165 CHKERR AOCreateMappingIS(is, PETSC_NULLPTR, &ao);
1166 if (ss == 0) {
1167 IS is_dup;
1168 CHKERR ISDuplicate(is, &is_dup);
1169 out_problem_it->getSubData()->rowIs = SmartPetscObj<IS>(is_dup, false);
1170 out_problem_it->getSubData()->rowMap = SmartPetscObj<AO>(ao, true);
1171 } else {
1172 IS is_dup;
1173 CHKERR ISDuplicate(is, &is_dup);
1174 out_problem_it->getSubData()->colIs = SmartPetscObj<IS>(is_dup, false);
1175 out_problem_it->getSubData()->colMap = SmartPetscObj<AO>(ao, true);
1176 }
1177 CHKERR AOApplicationToPetscIS(ao, is);
1178 // set global number of DOFs
1179 CHKERR ISGetSize(is, nb_dofs[ss]);
1180 CHKERR ISDestroy(&is);
1181 // set out problem global indices after applying map
1182 dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1183 for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1184 dit++, it++) {
1185 bool success = out_problem_dofs[ss]->modify(
1186 out_problem_dofs[ss]->project<0>(dit),
1187 NumeredDofEntity_part_and_all_indices_change(
1188 dit->get()->getPart(), *it, *it,
1189 dit->get()->getPetscLocalDofIdx()));
1190 if (!success) {
1191 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1192 "operation unsuccessful");
1193 }
1194 }
1195 // set global indices to nodes not on this part
1196 {
1197 NumeredDofEntityByLocalIdx::iterator dit =
1198 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1199 NumeredDofEntityByLocalIdx::iterator hi_dit =
1200 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1201 const int nb = std::distance(dit, hi_dit);
1202 std::vector<int> main_indices_non_local(nb);
1203 for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1204 dit++, it++) {
1205 *it = dit->get()->getPetscGlobalDofIdx();
1206 }
1207 IS is;
1208 CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1209 &*main_indices_non_local.begin(),
1210 PETSC_USE_POINTER, &is);
1211 CHKERR AOApplicationToPetscIS(ao, is);
1212 CHKERR ISDestroy(&is);
1213 dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1214 for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1215 dit++, it++) {
1216 bool success = out_problem_dofs[ss]->modify(
1217 out_problem_dofs[ss]->project<0>(dit),
1218 NumeredDofEntity_part_and_all_indices_change(
1219 dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1220 dit->get()->getPetscLocalDofIdx()));
1221 if (!success) {
1222 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1223 "operation unsuccessful");
1224 }
1225 }
1226 }
1227 CHKERR AODestroy(&ao);
1228 }
1229 }
1230
1231 if (square_matrix) {
1232 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1233 out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1234 out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1235 out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1236 out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1237 }
1238
1239 CHKERR printPartitionedProblem(&*out_problem_it, verb);
1240 CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1241
1243}
1244
1245MoFEMErrorCode ProblemsManager::buildComposedProblem(
1246 const std::string out_name, const std::vector<std::string> add_row_problems,
1247 const std::vector<std::string> add_col_problems, const bool square_matrix,
1248 int verb) {
1249 if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
1250 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1251 if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1252 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1253 if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1254 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1255 MoFEM::Interface &m_field = cOre;
1256 auto problems_ptr = m_field.get_problems();
1258
1259 CHKERR m_field.clear_problem(out_name);
1260 // get reference to all problems
1261 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1262 ProblemByName &problems_by_name =
1263 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1264
1265 // Get iterators to out problem, i.e. build problem
1266 ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1267 if (out_problem_it == problems_by_name.end()) {
1268 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1269 "problem with name < %s > not defined (top tip check spelling)",
1270 out_name.c_str());
1271 }
1272 // Make data structure for composed-problem data
1273 out_problem_it->composedProblemsData =
1274 boost::make_shared<ComposedProblemsData>();
1275 boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1276 out_problem_it->getComposedProblemsData();
1277
1278 const std::vector<std::string> *add_prb[] = {&add_row_problems,
1279 &add_col_problems};
1280 std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1281 &cmp_prb_data->colProblemsAdd};
1282 std::vector<SmartPetscObj<IS>> *add_prb_is[] = {&cmp_prb_data->rowIs,
1283 &cmp_prb_data->colIs};
1284
1285 // Get local indices counter
1286 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1287 &out_problem_it->nbLocDofsCol};
1288 // Get global indices counter
1289 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1290
1291 // Set number of ghost nodes to zero
1292 {
1293 out_problem_it->nbDofsRow = 0;
1294 out_problem_it->nbDofsCol = 0;
1295 out_problem_it->nbLocDofsRow = 0;
1296 out_problem_it->nbLocDofsCol = 0;
1297 out_problem_it->nbGhostDofsRow = 0;
1298 out_problem_it->nbGhostDofsCol = 0;
1299 }
1300 int nb_dofs_reserve[] = {0, 0};
1301
1302 // Loop over rows and columns in the main problem and sub-problems
1303 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1304 add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1305 add_prb_is[ss]->reserve(add_prb[ss]->size());
1306 for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1307 vit != add_prb[ss]->end(); vit++) {
1308 ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1309 if (prb_it == problems_by_name.end()) {
1310 SETERRQ(
1311 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1312 "problem with name < %s > not defined (top tip check spelling)",
1313 vit->c_str());
1314 }
1315 add_prb_ptr[ss]->push_back(&*prb_it);
1316 // set number of dofs on rows and columns
1317 if (ss == 0) {
1318 // row
1319 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1320 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1321 nb_dofs_reserve[ss] +=
1322 add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1323 } else {
1324 // column
1325 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1326 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1327 nb_dofs_reserve[ss] +=
1328 add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1329 }
1330 }
1331 }
1332 // if squre problem, rows and columns are the same
1333 if (square_matrix) {
1334 add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1335 add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1336 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1337 *nb_dofs[1] = *nb_dofs[0];
1338 *nb_local_dofs[1] = *nb_local_dofs[0];
1339 }
1340
1341 // reserve memory for dofs
1342 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1343 // Reserve memory
1344 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1345 dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1346 dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1347 if (!ss)
1348 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1349 else
1350 out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1351 }
1352
1353 // Push back DOFs
1354 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1355 NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1356 dit,
1357 hi_dit;
1358 int shift_glob = 0;
1359 int shift_loc = 0;
1360 for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1361 PetscInt *dofs_out_idx_ptr;
1362 int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1363 CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1364 if (ss == 0) {
1365 dit = (*add_prb_ptr[ss])[pp]
1366 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1367 .begin();
1368 hi_dit = (*add_prb_ptr[ss])[pp]
1369 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1370 .end();
1371 } else {
1372 dit = (*add_prb_ptr[ss])[pp]
1373 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1374 .begin();
1375 hi_dit = (*add_prb_ptr[ss])[pp]
1376 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1377 .end();
1378 }
1379 int is_nb = 0;
1380 for (; dit != hi_dit; dit++) {
1381 const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
1382 const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
1383 const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
1384 if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1385 continue;
1386 const int rank = m_field.get_comm_rank();
1387 const int part = dit->get()->getPart();
1388 const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1389 const int loc_idx =
1390 (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1391 : -1;
1392 dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1393 glob_idx, loc_idx, part);
1394 if (part == rank) {
1395 dofs_out_idx_ptr[is_nb++] = glob_idx;
1396 }
1397 }
1398 if (is_nb > nb_local_dofs) {
1399 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1400 "Data inconsistency");
1401 }
1402 IS is;
1403 CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1404 PETSC_OWN_POINTER, &is);
1405 auto smart_is = SmartPetscObj<IS>(is);
1406 (*add_prb_is[ss]).push_back(smart_is);
1407 if (ss == 0) {
1408 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1409 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1410 } else {
1411 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1412 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1413 }
1414 if (square_matrix) {
1415 (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1416 (*add_prb_is[1]).push_back(smart_is);
1417 }
1418 }
1419 }
1420
1421 if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1422 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1423 }
1424
1425 // Insert DOFs to problem multi-index
1426 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1427 auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1428 : out_problem_it->numeredColDofsPtr->end();
1429 for (auto &v : *dofs_array[ss])
1430 hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1431 hint, dofs_array[ss], &v)
1432 : out_problem_it->numeredColDofsPtr->emplace_hint(
1433 hint, dofs_array[ss], &v);
1434 }
1435
1436 // Compress DOFs
1437 *nb_dofs[0] = 0;
1438 *nb_dofs[1] = 0;
1439 *nb_local_dofs[0] = 0;
1440 *nb_local_dofs[1] = 0;
1441 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1442
1443 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1444 if (ss == 0) {
1445 dofs_ptr = out_problem_it->numeredRowDofsPtr;
1446 } else {
1447 dofs_ptr = out_problem_it->numeredColDofsPtr;
1448 }
1449 NumeredDofEntityByUId::iterator dit, hi_dit;
1450 dit = dofs_ptr->get<Unique_mi_tag>().begin();
1451 hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1452 std::vector<int> idx;
1453 idx.reserve(std::distance(dit, hi_dit));
1454 // set dofs in order entity and dof number on entity
1455 for (; dit != hi_dit; dit++) {
1456 if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1457 bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1458 dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1459 if (!success) {
1460 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1461 "modification unsuccessful");
1462 }
1463 idx.push_back(dit->get()->getPetscGlobalDofIdx());
1464 } else {
1465 if (dit->get()->getPetscLocalDofIdx() != -1) {
1466 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1467 "local index should be negative");
1468 }
1469 }
1470 }
1471 if (square_matrix) {
1472 *nb_local_dofs[1] = *nb_local_dofs[0];
1473 }
1474
1475 // set new dofs mapping
1476 IS is;
1477 CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1478 PETSC_USE_POINTER, &is);
1479 CHKERR ISGetSize(is, nb_dofs[ss]);
1480 if (square_matrix) {
1481 *nb_dofs[1] = *nb_dofs[0];
1482 }
1483
1484 AO ao;
1485 CHKERR AOCreateMappingIS(is, PETSC_NULLPTR, &ao);
1486 for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1487 CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1488
1489 // Set DOFs numeration
1490 {
1491 std::vector<int> idx_new;
1492 idx_new.reserve(dofs_ptr->size());
1493 for (NumeredDofEntityByUId::iterator dit =
1494 dofs_ptr->get<Unique_mi_tag>().begin();
1495 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1496 idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1497 }
1498 // set new global dofs numeration
1499 IS is_new;
1500 CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1501 &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1502 CHKERR AOApplicationToPetscIS(ao, is_new);
1503 // set global indices to multi-index
1504 std::vector<int>::iterator vit = idx_new.begin();
1505 for (NumeredDofEntityByUId::iterator dit =
1506 dofs_ptr->get<Unique_mi_tag>().begin();
1507 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1508 bool success =
1509 dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1510 dit->get()->getPart(), *(vit++)));
1511 if (!success) {
1512 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1513 "modification unsuccessful");
1514 }
1515 }
1516 CHKERR ISDestroy(&is_new);
1517 }
1518 CHKERR ISDestroy(&is);
1519 CHKERR AODestroy(&ao);
1520 }
1521
1522 CHKERR printPartitionedProblem(&*out_problem_it, verb);
1523 CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1524
1525 // Inidcate that porble has been build
1526 cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM;
1527 cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
1528
1530}
1531
1532MoFEMErrorCode ProblemsManager::partitionSimpleProblem(const std::string name,
1533 int verb) {
1534
1535 MoFEM::Interface &m_field = cOre;
1536 auto problems_ptr = m_field.get_problems();
1538
1539 if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
1540 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1541 if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
1542 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1543 if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
1544 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1545 if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
1546 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1547 MOFEM_LOG("WORLD", Sev::verbose) << "Simple partition problem " << name;
1548
1549 // find p_miit
1550 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1551 ProblemByName &problems_set =
1552 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1553 ProblemByName::iterator p_miit = problems_set.find(name);
1554 if (p_miit == problems_set.end()) {
1555 SETERRQ(PETSC_COMM_SELF, 1,
1556 "problem < %s > is not found (top tip: check spelling)",
1557 name.c_str());
1558 }
1559 typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1560 Idx_mi_tag>::type NumeredDofEntitysByIdx;
1561 NumeredDofEntitysByIdx &dofs_row_by_idx =
1562 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>();
1563 NumeredDofEntitysByIdx &dofs_col_by_idx =
1564 p_miit->numeredColDofsPtr->get<Idx_mi_tag>();
1565 boost::multi_index::index<NumeredDofEntity_multiIndex,
1566 Idx_mi_tag>::type::iterator miit_row,
1567 hi_miit_row;
1568 boost::multi_index::index<NumeredDofEntity_multiIndex,
1569 Idx_mi_tag>::type::iterator miit_col,
1570 hi_miit_col;
1571 DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1572 DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1573 nb_row_local_dofs = 0;
1574 nb_row_ghost_dofs = 0;
1575 DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1576 DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1577 nb_col_local_dofs = 0;
1578 nb_col_ghost_dofs = 0;
1579
1580 bool square_matrix = false;
1581 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1582 square_matrix = true;
1583 }
1584
1585 // get row range of local indices
1586 PetscLayout layout_row;
1587 const int *ranges_row;
1588
1589 DofIdx nb_dofs_row = dofs_row_by_idx.size();
1590 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1591 CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1592 CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1593 CHKERR PetscLayoutSetUp(layout_row);
1594 CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1595 // get col range of local indices
1596 PetscLayout layout_col;
1597 const int *ranges_col;
1598 if (!square_matrix) {
1599 DofIdx nb_dofs_col = dofs_col_by_idx.size();
1600 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1601 CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1602 CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1603 CHKERR PetscLayoutSetUp(layout_col);
1604 CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1605 }
1606 for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1607 part++) {
1608 miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1609 hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1610 if (std::distance(miit_row, hi_miit_row) !=
1611 ranges_row[part + 1] - ranges_row[part]) {
1612 SETERRQ(
1613 PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1614 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1615 "rstart (%ld != %d - %d = %d) ",
1616 std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1617 ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1618 }
1619 // loop rows
1620 for (; miit_row != hi_miit_row; miit_row++) {
1621 bool success = dofs_row_by_idx.modify(
1622 miit_row,
1623 NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1624 if (!success)
1625 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1626 "modification unsuccessful");
1627 if (part == (unsigned int)m_field.get_comm_rank()) {
1628 success = dofs_row_by_idx.modify(
1629 miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1630 if (!success)
1631 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1632 "modification unsuccessful");
1633 }
1634 }
1635 if (!square_matrix) {
1636 miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1637 hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1638 if (std::distance(miit_col, hi_miit_col) !=
1639 ranges_col[part + 1] - ranges_col[part]) {
1640 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1641 "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1642 "rend - "
1643 "rstart (%ld != %d - %d = %d) ",
1644 std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1645 ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1646 }
1647 // loop cols
1648 for (; miit_col != hi_miit_col; miit_col++) {
1649 bool success = dofs_col_by_idx.modify(
1650 miit_col, NumeredDofEntity_part_and_glob_idx_change(
1651 part, (*miit_col)->dofIdx));
1652 if (!success)
1653 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1654 "modification unsuccessful");
1655 if (part == (unsigned int)m_field.get_comm_rank()) {
1656 success = dofs_col_by_idx.modify(
1657 miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1658 if (!success)
1659 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1660 "modification unsuccessful");
1661 }
1662 }
1663 }
1664 }
1665 CHKERR PetscLayoutDestroy(&layout_row);
1666 if (!square_matrix) {
1667 CHKERR PetscLayoutDestroy(&layout_col);
1668 }
1669 if (square_matrix) {
1670 nb_col_local_dofs = nb_row_local_dofs;
1671 nb_col_ghost_dofs = nb_row_ghost_dofs;
1672 }
1673 CHKERR printPartitionedProblem(&*p_miit, verb);
1674 cOre.getBuildMoFEM() |= Core::PARTITION_PROBLEM;
1676}
1677
1678MoFEMErrorCode ProblemsManager::partitionProblem(const std::string name,
1679 int verb) {
1680 MoFEM::Interface &m_field = cOre;
1681 auto problems_ptr = m_field.get_problems();
1683
1684 MOFEM_LOG("WORLD", Sev::noisy) << "Partition problem " << name;
1685
1686 using NumeredDofEntitysByIdx =
1687 NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type;
1688 using ProblemsByName = Problem_multiIndex::index<Problem_mi_tag>::type;
1689
1690 // Find problem pointer by name
1691 auto &problems_set =
1692 const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1693 auto p_miit = problems_set.find(name);
1694 if (p_miit == problems_set.end()) {
1695 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1696 "problem with name %s not defined (top tip check spelling)",
1697 name.c_str());
1698 }
1699 int nb_dofs_row = p_miit->getNbDofsRow();
1700
1701 if (m_field.get_comm_size() != 1) {
1702
1703 if (!(cOre.getBuildMoFEM() & (1 << 0)))
1704 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1705 if (!(cOre.getBuildMoFEM() & (1 << 1)))
1706 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1707 if (!(cOre.getBuildMoFEM() & (1 << 2)))
1708 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1709 "entFEAdjacencies not build");
1710 if (!(cOre.getBuildMoFEM() & (1 << 3)))
1711 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1712
1713 Mat Adj;
1714 CHKERR m_field.getInterface<MatrixManager>()
1715 ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1716
1717 int m, n;
1718 CHKERR MatGetSize(Adj, &m, &n);
1719 if (verb > VERY_VERBOSE)
1720 MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1721
1722 // partitioning
1723 MatPartitioning part;
1724 IS is;
1725 CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
1726 CHKERR MatPartitioningSetAdjacency(part, Adj);
1727 CHKERR MatPartitioningSetFromOptions(part);
1728 CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
1729 CHKERR MatPartitioningApply(part, &is);
1730 if (verb > VERY_VERBOSE)
1731 ISView(is, PETSC_VIEWER_STDOUT_WORLD);
1732
1733 // gather
1734 IS is_gather, is_num, is_gather_num;
1735 CHKERR ISAllGather(is, &is_gather);
1736 CHKERR ISPartitioningToNumbering(is, &is_num);
1737 CHKERR ISAllGather(is_num, &is_gather_num);
1738 const int *part_number, *petsc_idx;
1739 CHKERR ISGetIndices(is_gather, &part_number);
1740 CHKERR ISGetIndices(is_gather_num, &petsc_idx);
1741 int size_is_num, size_is_gather;
1742 CHKERR ISGetSize(is_gather, &size_is_gather);
1743 if (size_is_gather != (int)nb_dofs_row)
1744 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1745 "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
1746
1747 CHKERR ISGetSize(is_num, &size_is_num);
1748 if (size_is_num != (int)nb_dofs_row)
1749 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1750 "data inconsistency %d != %d", size_is_num, nb_dofs_row);
1751
1752 bool square_matrix = false;
1753 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1754 square_matrix = true;
1755
1756 // if (!square_matrix) {
1757 // // FIXME: This is for back compatibility, if deprecate interface
1758 // function
1759 // // build interfaces is removed, this part of the code will be obsolete
1760 // auto mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().begin();
1761 // auto hi_mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().end();
1762 // auto mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().begin();
1763 // auto hi_mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().end();
1764 // for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
1765 // if (mit_col == hi_mit_col) {
1766 // SETERRQ(
1767 // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1768 // "check finite element definition, nb. of rows is not equal to "
1769 // "number for columns");
1770 // }
1771 // if (mit_row->get()->getLocalUniqueId() !=
1772 // mit_col->get()->getLocalUniqueId()) {
1773 // SETERRQ(
1774 // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1775 // "check finite element definition, nb. of rows is not equal to "
1776 // "number for columns");
1777 // }
1778 // }
1779 // }
1780
1781 auto number_dofs = [&](auto &dofs_idx, auto &counter) {
1783 for (auto miit_dofs_row = dofs_idx.begin();
1784 miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
1785 const int part = part_number[(*miit_dofs_row)->dofIdx];
1786 if (part == (unsigned int)m_field.get_comm_rank()) {
1787 const bool success = dofs_idx.modify(
1788 miit_dofs_row,
1789 NumeredDofEntity_part_and_indices_change(
1790 part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
1791 if (!success) {
1792 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1793 "modification unsuccessful");
1794 }
1795 } else {
1796 const bool success = dofs_idx.modify(
1797 miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
1798 part, petsc_idx[(*miit_dofs_row)->dofIdx]));
1799 if (!success) {
1800 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1801 "modification unsuccessful");
1802 }
1803 }
1804 }
1806 };
1807
1808 // Set petsc global indices
1809 auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1810 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1811 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1812 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1813 nb_row_local_dofs = 0;
1814 nb_row_ghost_dofs = 0;
1815
1816 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1817
1818 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1819 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1820 if (square_matrix) {
1821 nb_col_local_dofs = nb_row_local_dofs;
1822 nb_col_ghost_dofs = nb_row_ghost_dofs;
1823 } else {
1824 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1825 const_cast<NumeredDofEntitysByIdx &>(
1826 p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1827 nb_col_local_dofs = 0;
1828 nb_col_ghost_dofs = 0;
1829 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1830 }
1831
1832 CHKERR ISRestoreIndices(is_gather, &part_number);
1833 CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
1834 CHKERR ISDestroy(&is_num);
1835 CHKERR ISDestroy(&is_gather_num);
1836 CHKERR ISDestroy(&is_gather);
1837 CHKERR ISDestroy(&is);
1838 CHKERR MatPartitioningDestroy(&part);
1839 CHKERR MatDestroy(&Adj);
1840 CHKERR printPartitionedProblem(&*p_miit, verb);
1841 } else {
1842
1843 auto number_dofs = [&](auto &dof_idx, auto &counter) {
1845 for (auto miit_dofs_row = dof_idx.begin(); miit_dofs_row != dof_idx.end();
1846 miit_dofs_row++) {
1847 const bool success = dof_idx.modify(
1848 miit_dofs_row,
1849 NumeredDofEntity_part_and_indices_change(0, counter, counter));
1850 ++counter;
1851 if (!success) {
1852 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1853 "modification unsuccessful");
1854 }
1855 }
1857 };
1858
1859 auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1860 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1861 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1862 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1863 nb_row_local_dofs = 0;
1864 nb_row_ghost_dofs = 0;
1865
1866 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1867
1868 bool square_matrix = false;
1869 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1870 square_matrix = true;
1871
1872 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1873 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1874 if (square_matrix) {
1875 nb_col_local_dofs = nb_row_local_dofs;
1876 nb_col_ghost_dofs = nb_row_ghost_dofs;
1877 } else {
1878 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1879 const_cast<NumeredDofEntitysByIdx &>(
1880 p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1881 nb_col_local_dofs = 0;
1882 nb_col_ghost_dofs = 0;
1883 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1884 }
1885 }
1886
1887 cOre.getBuildMoFEM() |= 1 << 4;
1889}
1890
1891MoFEMErrorCode ProblemsManager::inheritPartition(
1892 const std::string name, const std::string problem_for_rows, bool copy_rows,
1893 const std::string problem_for_cols, bool copy_cols, int verb) {
1894 MoFEM::Interface &m_field = cOre;
1895 auto problems_ptr = m_field.get_problems();
1897
1898 if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
1899 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
1900
1901 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1902
1903 // find p_miit
1904 ProblemByName &problems_by_name =
1905 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1906 ProblemByName::iterator p_miit = problems_by_name.find(name);
1907 if (p_miit == problems_by_name.end()) {
1908 SETERRQ(m_field.get_comm(), MOFEM_NOT_FOUND,
1909 "problem with name < %s > not defined (top tip check spelling)",
1910 name.c_str());
1911 }
1912 if (verb > QUIET)
1913 MOFEM_LOG("WORLD", Sev::inform)
1914 << p_miit->getName() << " from rows of " << problem_for_rows
1915 << " and columns of " << problem_for_cols;
1916
1917 // find p_miit_row
1918 ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
1919 if (p_miit_row == problems_by_name.end()) {
1920 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1921 "problem with name < %s > not defined (top tip check spelling)",
1922 problem_for_rows.c_str());
1923 }
1924 const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
1925 p_miit_row->numeredRowDofsPtr;
1926
1927 // find p_mit_col
1928 ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
1929 if (p_miit_col == problems_by_name.end()) {
1930 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1931 "problem with name < %s > not defined (top tip check spelling)",
1932 problem_for_cols.c_str());
1933 }
1934 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
1935 p_miit_col->numeredColDofsPtr;
1936
1937 bool copy[] = {copy_rows, copy_cols};
1938 boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
1939 p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
1940
1941 int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
1942 int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
1943 boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
1944 dofs_col};
1945
1946 for (int ss = 0; ss < 2; ss++) {
1947
1948 // build indices
1949 *nb_local_dofs[ss] = 0;
1950 if (!copy[ss]) {
1951
1952 // only copy indices which are belong to some elements if this problem
1953 std::vector<int> is_local, is_new;
1954
1955 NumeredDofEntityByUId &dofs_by_uid =
1956 copied_dofs[ss]->get<Unique_mi_tag>();
1957 for (NumeredDofEntity_multiIndex::iterator dit =
1958 composed_dofs[ss]->begin();
1959 dit != composed_dofs[ss]->end(); dit++) {
1960 NumeredDofEntityByUId::iterator diit =
1961 dofs_by_uid.find((*dit)->getLocalUniqueId());
1962 if (diit == dofs_by_uid.end()) {
1963 SETERRQ(
1965 "data inconsistency, could not find dof in composite problem");
1966 }
1967 int part_number = (*diit)->getPart(); // get part number
1968 int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
1969 bool success;
1970 success = composed_dofs[ss]->modify(
1971 dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
1972 petsc_global_dof));
1973 if (!success) {
1974 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1975 "modification unsuccessful");
1976 }
1977 if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
1978 success = composed_dofs[ss]->modify(
1979 dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1980 if (!success) {
1981 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1982 "modification unsuccessful");
1983 }
1984 is_local.push_back(petsc_global_dof);
1985 }
1986 }
1987
1988 AO ao;
1989 CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
1990 NULL, &ao);
1991
1992 // apply local to global mapping
1993 is_local.resize(0);
1994 for (NumeredDofEntity_multiIndex::iterator dit =
1995 composed_dofs[ss]->begin();
1996 dit != composed_dofs[ss]->end(); dit++) {
1997 is_local.push_back((*dit)->getPetscGlobalDofIdx());
1998 }
1999 CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2000 int idx2 = 0;
2001 for (NumeredDofEntity_multiIndex::iterator dit =
2002 composed_dofs[ss]->begin();
2003 dit != composed_dofs[ss]->end(); dit++) {
2004 int part_number = (*dit)->getPart(); // get part number
2005 int petsc_global_dof = is_local[idx2++];
2006 bool success;
2007 success = composed_dofs[ss]->modify(
2008 dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2009 petsc_global_dof));
2010 if (!success) {
2011 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2012 "modification unsuccessful");
2013 }
2014 }
2015
2016 CHKERR AODestroy(&ao);
2017
2018 } else {
2019
2020 for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2021 dit != copied_dofs[ss]->end(); dit++) {
2022 std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2023 p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2024 new NumeredDofEntity((*dit)->getDofEntityPtr())));
2025 if (p.second) {
2026 (*nb_dofs[ss])++;
2027 }
2028 int dof_idx = (*dit)->getDofIdx();
2029 int part_number = (*dit)->getPart(); // get part number
2030 int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2031 if (part_number == (unsigned int)m_field.get_comm_rank()) {
2032 const bool success = composed_dofs[ss]->modify(
2033 p.first, NumeredDofEntity_part_and_all_indices_change(
2034 part_number, dof_idx, petsc_global_dof,
2035 (*nb_local_dofs[ss])++));
2036 if (!success) {
2037 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2038 "modification unsuccessful");
2039 }
2040 } else {
2041 const bool success = composed_dofs[ss]->modify(
2042 p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
2043 part_number, dof_idx, petsc_global_dof));
2044 if (!success) {
2045 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2046 "modification unsuccessful");
2047 }
2048 }
2049 }
2050 }
2051 }
2052
2053 CHKERR printPartitionedProblem(&*p_miit, verb);
2054 CHKERR debugPartitionedProblem(&*p_miit, verb);
2055
2057}
2058
2060ProblemsManager::printPartitionedProblem(const Problem *problem_ptr, int verb) {
2061 MoFEM::Interface &m_field = cOre;
2063
2064 if (verb > QUIET) {
2065
2066 MOFEM_LOG("SYNC", Sev::inform)
2067 << problem_ptr->getName() << " Nb. local dof "
2068 << problem_ptr->getNbLocalDofsRow() << " by "
2069 << problem_ptr->getNbLocalDofsCol() << " nb global dofs "
2070 << problem_ptr->getNbDofsRow() << " by " << problem_ptr->getNbDofsCol();
2071
2073 }
2074
2076}
2077
2079ProblemsManager::debugPartitionedProblem(const Problem *problem_ptr, int verb) {
2080 MoFEM::Interface &m_field = cOre;
2082
2083 auto save_ent = [](moab::Interface &moab, const std::string name,
2084 const EntityHandle ent) {
2086 EntityHandle out_meshset;
2087 CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
2088 CHKERR moab.add_entities(out_meshset, &ent, 1);
2089 CHKERR moab.write_file(name.c_str(), "VTK", "", &out_meshset, 1);
2090 CHKERR moab.delete_entities(&out_meshset, 1);
2092 };
2093
2094 if (debug > 0) {
2095
2096 typedef NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type
2097 NumeredDofEntitysByIdx;
2098 NumeredDofEntitysByIdx::iterator dit, hi_dit;
2099 const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
2100 &(problem_ptr->numeredRowDofsPtr->get<Idx_mi_tag>()),
2101 &(problem_ptr->numeredColDofsPtr->get<Idx_mi_tag>())};
2102
2103 int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
2104 int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
2105 &problem_ptr->nbLocDofsCol};
2106
2107 for (int ss = 0; ss < 2; ss++) {
2108
2109 dit = numered_dofs_ptr[ss]->begin();
2110 hi_dit = numered_dofs_ptr[ss]->end();
2111 for (; dit != hi_dit; dit++) {
2112 if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
2113 if ((*dit)->getPetscLocalDofIdx() < 0) {
2114 std::ostringstream zz;
2115 zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2116 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2117 "local dof index for %d (0-row, 1-col) not set, i.e. has "
2118 "negative value\n %s",
2119 ss, zz.str().c_str());
2120 }
2121 if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
2122 std::ostringstream zz;
2123 zz << "rank " << m_field.get_comm_rank() << " " << **dit;
2124 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2125 "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
2126 zz.str().c_str());
2127 }
2128 } else {
2129 if ((*dit)->getPetscGlobalDofIdx() < 0) {
2130
2131 const EntityHandle ent = (*dit)->getEnt();
2132 CHKERR save_ent(
2133 m_field.get_moab(),
2134 "debug_part" +
2135 boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
2136 "_negative_global_index.vtk",
2137 ent);
2138
2139 std::ostringstream zz;
2140 zz << "rank " << m_field.get_comm_rank() << " "
2141 << dit->get()->getBitRefLevel() << " " << **dit;
2142 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2143 "global dof index for %d (0-row, 1-col) row not set, i.e. "
2144 "has negative value\n %s",
2145 ss, zz.str().c_str());
2146 }
2147 if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
2148 std::ostringstream zz;
2149 zz << "rank " << m_field.get_comm_rank() << " nb_dofs "
2150 << *nbdof_ptr[ss] << " " << **dit;
2151 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
2152 "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
2153 zz.str().c_str());
2154 }
2155 }
2156 }
2157 }
2158 }
2160}
2161
2162MoFEMErrorCode ProblemsManager::partitionFiniteElements(const std::string name,
2163 bool part_from_moab,
2164 int low_proc,
2165 int hi_proc, int verb) {
2166 MoFEM::Interface &m_field = cOre;
2167 auto problems_ptr = m_field.get_problems();
2168 auto fe_ent_ptr = m_field.get_ents_finite_elements();
2170
2171 if (!(cOre.getBuildMoFEM() & Core::BUILD_FIELD))
2172 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2173
2174 if (!(cOre.getBuildMoFEM() & Core::BUILD_FE))
2175 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2176
2177 if (!(cOre.getBuildMoFEM() & Core::BUILD_ADJ))
2178 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2179 "adjacencies not build");
2180
2181 if (!(cOre.getBuildMoFEM() & Core::BUILD_PROBLEM))
2182 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2183
2184 if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
2185 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2186 "problem not partitioned");
2187
2188 if (low_proc == -1)
2189 low_proc = m_field.get_comm_rank();
2190 if (hi_proc == -1)
2191 hi_proc = m_field.get_comm_rank();
2192
2193 // Find pointer to problem of given name
2194 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2195 auto &problems =
2196 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2197 ProblemByName::iterator p_miit = problems.find(name);
2198 if (p_miit == problems.end()) {
2199 SETERRQ(m_field.get_comm(), MOFEM_NOT_FOUND,
2200 "problem < %s > not found (top tip: check spelling)",
2201 name.c_str());
2202 }
2203
2204 // Get reference on finite elements multi-index on the problem
2205 NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2206 *p_miit->numeredFiniteElementsPtr;
2207
2208 // Clear all elements and data, build it again
2209 problem_finite_elements.clear();
2210
2211 // Check if dofs and columns are the same, i.e. structurally symmetric
2212 // problem
2213 bool do_cols_prob = true;
2214 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2215 do_cols_prob = false;
2216 }
2217
2218 auto get_good_elems = [&]() {
2219 auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
2220 good_elems.reserve(fe_ent_ptr->size());
2221
2222 const auto prb_bit = p_miit->getBitRefLevel();
2223 const auto prb_mask = p_miit->getBitRefLevelMask();
2224
2225 // Loop over all elements in database and if right element is there add it
2226 // to problem finite element multi-index
2227 for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2228
2229 // if element is not part of problem
2230 if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2231
2232 const auto &fe_bit = (*efit)->getBitRefLevel();
2233
2234 // if entity is not problem refinement level
2235 if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit).any())
2236 good_elems.emplace_back(efit);
2237 }
2238 }
2239
2240 return good_elems;
2241 };
2242
2243 auto good_elems = get_good_elems();
2244
2245 auto numbered_good_elems_ptr =
2246 boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2247 numbered_good_elems_ptr->reserve(good_elems.size());
2248 for (auto &efit : good_elems)
2249 numbered_good_elems_ptr->emplace_back(NumeredEntFiniteElement(*efit));
2250
2251 if (!do_cols_prob) {
2252 for (auto &fe : *numbered_good_elems_ptr) {
2253 if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2254 fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2255 }
2256 }
2257 }
2258
2259 if (part_from_moab) {
2260 for (auto &fe : *numbered_good_elems_ptr) {
2261 if (fe.getEntType() == MBENTITYSET) {
2262 fe.part = m_field.get_comm_rank();
2263 } else {
2264 // if partition is taken from moab partition
2265 int proc = fe.getPartProc();
2266 if (proc == -1 && fe.getEntType() == MBVERTEX)
2267 proc = fe.getOwnerProc();
2268 fe.part = proc;
2269 }
2270 }
2271 }
2272
2273 for (auto &fe : *numbered_good_elems_ptr) {
2274
2276 CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2277
2278 if (!part_from_moab) {
2279 if (fe.getEntType() == MBENTITYSET) {
2280 fe.part = m_field.get_comm_rank();
2281 } else {
2282 std::vector<int> parts(m_field.get_comm_size(), 0);
2283 for (auto &dof_ptr : rows_view)
2284 parts[dof_ptr->pArt]++;
2285 std::vector<int>::iterator pos =
2286 max_element(parts.begin(), parts.end());
2287 const auto max_part = std::distance(parts.begin(), pos);
2288 fe.part = max_part;
2289 }
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
2339 cOre.getBuildMoFEM() |= Core::PARTITION_FE;
2341}
2342
2343MoFEMErrorCode ProblemsManager::partitionGhostDofs(const std::string name,
2344 int verb) {
2345 MoFEM::Interface &m_field = cOre;
2346 auto problems_ptr = m_field.get_problems();
2348
2349 if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
2350 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2351 "partition of problem not build");
2352 if (!(cOre.getBuildMoFEM() & Core::PARTITION_FE))
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 SETERRQ(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
2472 cOre.getBuildMoFEM() |= Core::PARTITION_GHOST_DOFS;
2474}
2475
2477ProblemsManager::partitionGhostDofsOnDistributedMesh(const std::string name,
2478 int verb) {
2479 MoFEM::Interface &m_field = cOre;
2480 auto problems_ptr = m_field.get_problems();
2482
2483 if (!(cOre.getBuildMoFEM() & Core::PARTITION_PROBLEM))
2484 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2485 "partition of problem not build");
2486 if (!(cOre.getBuildMoFEM() & Core::PARTITION_FE))
2487 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2488 "partitions finite elements not build");
2489
2490 // get problem pointer
2491 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
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 // iterate 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 // iterate 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
2559 cOre.getBuildMoFEM() |= Core::PARTITION_GHOST_DOFS;
2561}
2562
2563MoFEMErrorCode ProblemsManager::getFEMeshset(const std::string prb_name,
2564 const std::string &fe_name,
2565 EntityHandle *meshset) const {
2566 MoFEM::Interface &m_field = cOre;
2567 const Problem *problem_ptr;
2568 const FiniteElement_multiIndex *fes_ptr;
2570
2571 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2572 CHKERR m_field.get_problem(prb_name, &problem_ptr);
2573 CHKERR m_field.get_finite_elements(&fes_ptr);
2574
2575 auto fe_miit = fes_ptr->get<FiniteElement_name_mi_tag>().find(fe_name);
2576 if (fe_miit != fes_ptr->get<FiniteElement_name_mi_tag>().end()) {
2577 auto fit =
2578 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
2579 EntFiniteElement::getLocalUniqueIdCalculate(
2580 0, (*fe_miit)->getFEUId()));
2581 auto hi_fe_it =
2582 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
2583 EntFiniteElement::getLocalUniqueIdCalculate(
2584 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
2585 std::vector<EntityHandle> fe_vec;
2586 fe_vec.reserve(std::distance(fit, hi_fe_it));
2587 for (; fit != hi_fe_it; fit++)
2588 fe_vec.push_back(fit->get()->getEnt());
2589 CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2590 fe_vec.size());
2591 }
2592
2594}
2595
2597ProblemsManager::getProblemElementsLayout(const std::string name,
2598 const std::string &fe_name,
2599 PetscLayout *layout) const {
2600 MoFEM::Interface &m_field = cOre;
2601 const Problem *problem_ptr;
2603 CHKERR m_field.get_problem(name, &problem_ptr);
2604 CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2605 fe_name, layout);
2607}
2608
2609MoFEMErrorCode ProblemsManager::getSideDofsOnBrokenSpaceEntities(
2610
2611 std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view,
2612
2613 const std::string problem_name, RowColData rc, const std::string field_name,
2614 const Range ents, int bridge_dim, const int lo_coeff, const int hi_coeff,
2615 const int lo_order, const int hi_order, int verb, const bool debug
2616
2617) const {
2618 MoFEM::Interface &m_field = cOre;
2620
2621 switch (rc) {
2622 case ROW:
2623 case COL:
2624 break;
2625 default:
2626 SETERRQ(m_field.get_comm(), MOFEM_INVALID_DATA, "invalid RowColData");
2627 break;
2628 }
2629
2630 const Problem *prb_ptr;
2631 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2632
2633 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2634 prb_ptr->numeredRowDofsPtr, nullptr};
2635 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2636 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2637
2638 Range bridge_ents;
2639 CHKERR m_field.get_moab().get_adjacencies(
2640 ents, bridge_dim, false, bridge_ents, moab::Interface::UNION);
2641
2642 auto &moab = m_field.get_moab();
2643 const auto bit_number = m_field.get_field_bit_number(field_name);
2644 const auto nb_coeffs =
2646 const auto &side_dof_map = m_field.get_field_structure(field_name)
2647 ->getDofSideMap();
2648
2649 using NumeredDofEntity_it_view_multiIndex = multi_index_container<
2650
2651 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2652
2653 >;
2654 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2655
2656 for (auto pit = bridge_ents.const_pair_begin();
2657 pit != bridge_ents.const_pair_end(); ++pit) {
2658 auto lo = numered_dofs[rc]->get<Unique_mi_tag>().lower_bound(
2659 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2660 auto hi = numered_dofs[rc]->get<Unique_mi_tag>().upper_bound(
2661 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2662
2663 auto bride_type = moab.type_from_handle(pit->first);
2664
2665 for (; lo != hi; ++lo) {
2666 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2667 (*lo)->getDofCoeffIdx() <= hi_coeff &&
2668 (*lo)->getDofOrder() >= lo_order &&
2669 (*lo)->getDofOrder() <= hi_order) {
2670 auto ent_dof_index = (*lo)->getEntDofIdx();
2671 auto side_it = side_dof_map.at(bride_type)
2672 .get<EntDofIdx_mi_tag>()
2673 .find(std::floor(static_cast<double>(ent_dof_index) /
2674 nb_coeffs));
2675 if (side_it !=
2676 side_dof_map.at(bride_type).get<EntDofIdx_mi_tag>().end()) {
2677 auto bridge_ent = (*lo)->getEnt();
2678 auto type = side_it->type;
2679 auto side = side_it->side;
2680 auto dim = CN::Dimension(type);
2681 EntityHandle side_ent = 0;
2682 CHKERR moab.side_element(bridge_ent, dim, side, side_ent);
2683 if (ents.find(side_ent) != ents.end()) {
2684 dofs_it_view.emplace_back(numered_dofs[rc]->project<0>(lo));
2685 }
2686 } else {
2687#ifndef NDEBUG
2688 MOFEM_LOG("SELF", Sev::error)
2689 << *m_field.get_field_structure(field_name);
2690 MOFEM_LOG("SELF", Sev::error)
2691 << "side not found for entity " << CN::EntityTypeName(bride_type);
2692 MOFEM_LOG("SELF", Sev::error)
2693 << "ent_dof_index / nb_coeffs "
2694 << std::floor(static_cast<double>(ent_dof_index) / nb_coeffs);
2695 MOFEM_LOG("SELF", Sev::error)
2696 << "side_dof_map.size() " << side_dof_map.size();
2697#endif // NDEBUG
2698 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2699 "side not found - you will get more information in debug");
2700 }
2701 }
2702 }
2703 }
2704
2705 if (verb > QUIET) {
2706 for (auto &dof : dofs_it_view)
2707 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2708 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
2709 }
2710
2711 vec_dof_view.reserve(vec_dof_view.size() + dofs_it_view.size());
2712 for (auto dit : dofs_it_view)
2713 vec_dof_view.push_back(*dit);
2714
2716}
2717
2718MoFEMErrorCode ProblemsManager::removeDofs(
2719 const std::string problem_name, RowColData rc,
2720 std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view, int verb,
2721 const bool debug) {
2722
2723 MoFEM::Interface &m_field = cOre;
2725
2726 switch (rc) {
2727 case ROW:
2728 case COL:
2729 break;
2730 default:
2731 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
2732 }
2733
2734 const Problem *prb_ptr;
2735 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2736
2737 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2738 prb_ptr->numeredRowDofsPtr, nullptr};
2739 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2740 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2741
2742 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2743 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2744 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2745
2746 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2747 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2748 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2749 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2750 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2751
2752 if (numered_dofs[rc]) {
2753
2754 if (verb >= NOISY)
2755 MOFEM_LOG_C("SYNC", Sev::noisy,
2756 "Number of DOFs in multi-index %d and to delete %d\n",
2757 numered_dofs[rc]->size(), vec_dof_view.size());
2758
2759 // erase dofs from problem
2760 for (auto weak_dit : vec_dof_view)
2761 if (auto dit = weak_dit.lock()) {
2762 numered_dofs[rc]->erase(dit->getLocalUniqueId());
2763 }
2764
2765 if (verb >= NOISY)
2766 MOFEM_LOG_C("SYNC", Sev::noisy,
2767 "Number of DOFs in multi-index after delete %d\n",
2768 numered_dofs[rc]->size());
2769
2770 // get current number of ghost dofs
2771 int nb_local_dofs = 0;
2772 int nb_ghost_dofs = 0;
2773 for (auto dit = numered_dofs[rc]->get<PetscLocalIdx_mi_tag>().begin();
2774 dit != numered_dofs[rc]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2775 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2776 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[rc]))
2777 ++nb_local_dofs;
2778 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[rc]))
2779 ++nb_ghost_dofs;
2780 else
2781 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2782 "Impossible case. You could run problem on no distributed "
2783 "mesh. That is not implemented. Dof local index is %d",
2784 (*dit)->getPetscLocalDofIdx());
2785 }
2786
2787 // get indices
2788 auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2789 const int nb_dofs = numered_dofs[rc]->size();
2790 indices.clear();
2791 indices.reserve(nb_dofs);
2792 for (auto dit = numered_dofs[rc]->get<decltype(tag)>().begin();
2793 dit != numered_dofs[rc]->get<decltype(tag)>().end(); ++dit) {
2794 bool add = true;
2795 if (only_local) {
2796 if ((*dit)->getPetscLocalDofIdx() < 0 ||
2797 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[rc])) {
2798 add = false;
2799 }
2800 }
2801 if (add)
2802 indices.push_back(decltype(tag)::get_index(dit));
2803 }
2804 };
2805
2806 auto get_indices_by_uid = [&](auto tag, auto &indices) {
2807 const int nb_dofs = numered_dofs[rc]->size();
2808 indices.clear();
2809 indices.reserve(nb_dofs);
2810 for (auto dit = numered_dofs[rc]->begin(); dit != numered_dofs[rc]->end();
2811 ++dit)
2812 indices.push_back(decltype(tag)::get_index(dit));
2813 };
2814
2815 auto get_sub_ao = [&](auto sub_data) {
2816 if (rc == 0) {
2817 return sub_data->getSmartRowMap();
2818 } else {
2819 return sub_data->getSmartColMap();
2820 }
2821 };
2822
2823 auto set_sub_is_and_ao = [&rc, &prb_ptr](auto sub_data, auto is, auto ao) {
2824 if (rc == ROW) {
2825 sub_data->rowIs = is;
2826 sub_data->rowMap = ao;
2827 sub_data->colIs = is; // if square problem col is equal to row
2828 sub_data->colMap = ao;
2829 } else {
2830 sub_data->colIs = is;
2831 sub_data->colMap = ao;
2832 }
2833 };
2834
2835 auto apply_symmetry = [&rc, &prb_ptr](auto sub_data) {
2836 if (rc == ROW) {
2837 if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
2838 sub_data->colIs = sub_data->getSmartRowIs();
2839 sub_data->colMap = sub_data->getSmartRowMap();
2840 }
2841 }
2842 };
2843
2844 auto concatenate_dofs = [&](auto tag, auto &indices,
2845 const auto local_only) {
2847 get_indices_by_tag(tag, indices, local_only);
2848
2849 SmartPetscObj<AO> ao;
2850 // Create AO from app indices (i.e. old), to pestc indices (new after
2851 // remove)
2852 if (local_only)
2853 ao = createAOMapping(m_field.get_comm(), indices.size(),
2854 &*indices.begin(), PETSC_NULLPTR);
2855 else
2856 ao = createAOMapping(PETSC_COMM_SELF, indices.size(), &*indices.begin(),
2857 PETSC_NULLPTR);
2858
2859 // Set mapping to sub dm data
2860 if (local_only) {
2861 if (auto sub_data = prb_ptr->getSubData()) {
2862 // create is and then map it to main problem of sub-problem
2863 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2864 &*indices.begin(), PETSC_COPY_VALUES);
2865 // get old app, i.e. oroginal befor sub indices, and ao, from app,
2866 // to petsc sub indices.
2867 auto sub_ao = get_sub_ao(sub_data);
2868 CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
2869 sub_ao = createAOMappingIS(sub_is, PETSC_NULLPTR);
2870 // set new sub ao
2871 set_sub_is_and_ao(sub_data, sub_is, sub_ao);
2872 apply_symmetry(sub_data);
2873 } else {
2874 // create sub data
2875 prb_ptr->getSubData() = boost::make_shared<Problem::SubProblemData>();
2876 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2877 &*indices.begin(), PETSC_COPY_VALUES);
2878 // set sub is ao
2879 set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
2880 apply_symmetry(prb_ptr->getSubData());
2881 }
2882 }
2883
2884 get_indices_by_uid(tag, indices);
2885 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
2886
2888 };
2889
2890 // set indices index
2891 auto set_concatenated_indices = [&]() {
2892 std::vector<int> global_indices;
2893 std::vector<int> local_indices;
2895 CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
2896 CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
2897 auto gi = global_indices.begin();
2898 auto li = local_indices.begin();
2899 for (auto dit = numered_dofs[rc]->begin(); dit != numered_dofs[rc]->end();
2900 ++dit) {
2901 auto mod = NumeredDofEntity_part_and_all_indices_change(
2902 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
2903 bool success = numered_dofs[rc]->modify(dit, mod);
2904 if (!success)
2905 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2906 "can not set negative indices");
2907 ++gi;
2908 ++li;
2909 }
2911 };
2912 CHKERR set_concatenated_indices();
2913
2914 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[rc], 1, MPI_INT, MPI_SUM,
2915 m_field.get_comm());
2916 *(local_nbdof_ptr[rc]) = nb_local_dofs;
2917 *(ghost_nbdof_ptr[rc]) = nb_ghost_dofs;
2918
2919 if (debug)
2920 for (auto dof : (*numered_dofs[rc])) {
2921 if (dof->getPetscGlobalDofIdx() < 0) {
2922 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2923 "Negative global idx");
2924 }
2925 if (dof->getPetscLocalDofIdx() < 0) {
2926 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2927 "Negative local idx");
2928 }
2929 }
2930
2931 } else {
2932
2933 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2934 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2935 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
2936 }
2937
2938 if (verb > QUIET) {
2940 "WORLD", Sev::inform,
2941 "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
2942 prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
2943 prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
2944 MOFEM_LOG_C("SYNC", Sev::verbose,
2945 "Removed DOFs from problem %s dofs [ %d / %d "
2946 "(before %d / %d) local, %d / %d (before %d / %d)]",
2947 prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
2948 prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
2949 nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
2950 prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
2951 nb_init_ghost_col_dofs);
2952 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2953 }
2954
2956}
2957
2958MoFEMErrorCode ProblemsManager::removeDofsOnEntities(
2959 const std::string problem_name, const std::string field_name,
2960 const Range ents, const int lo_coeff, const int hi_coeff,
2961 const int lo_order, const int hi_order, int verb, const bool debug) {
2962
2963 MoFEM::Interface &m_field = cOre;
2965
2966 const Problem *prb_ptr;
2967 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2968
2969 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2970 prb_ptr->numeredRowDofsPtr, nullptr};
2971 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2972 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2973
2974 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2975 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2976 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2977
2978 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2979 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2980 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2981 // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2982 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2983 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2984
2985 for (int s = 0; s != 2; ++s)
2986 if (numered_dofs[s]) {
2987
2988 using NumeredDofEntity_it_view_multiIndex = multi_index_container<
2989
2990 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2991
2992 >;
2993
2994 const auto bit_number = m_field.get_field_bit_number(field_name);
2995 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2996
2997 // Set -1 to global and local dofs indices
2998 for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2999 ++pit) {
3000 auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
3001 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
3002 auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
3003 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
3004
3005 for (; lo != hi; ++lo)
3006 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
3007 (*lo)->getDofCoeffIdx() <= hi_coeff &&
3008 (*lo)->getDofOrder() >= lo_order &&
3009 (*lo)->getDofOrder() <= hi_order)
3010 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
3011 }
3012
3013 if (verb > QUIET) {
3014 for (auto &dof : dofs_it_view)
3015 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
3016 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
3017 }
3018
3019 // create weak view
3020 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
3021 dofs_weak_view.reserve(dofs_it_view.size());
3022 for (auto dit : dofs_it_view)
3023 dofs_weak_view.push_back(*dit);
3024
3025 if (verb >= NOISY)
3026 MOFEM_LOG_C("SYNC", Sev::noisy,
3027 "Number of DOFs in multi-index %d and to delete %d\n",
3028 numered_dofs[s]->size(), dofs_it_view.size());
3029
3030 // erase dofs from problem
3031 for (auto weak_dit : dofs_weak_view)
3032 if (auto dit = weak_dit.lock()) {
3033 numered_dofs[s]->erase(dit->getLocalUniqueId());
3034 }
3035
3036 if (verb >= NOISY)
3037 MOFEM_LOG_C("SYNC", Sev::noisy,
3038 "Number of DOFs in multi-index after delete %d\n",
3039 numered_dofs[s]->size());
3040
3041 // get current number of ghost dofs
3042 int nb_local_dofs = 0;
3043 int nb_ghost_dofs = 0;
3044 for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
3045 dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
3046 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
3047 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
3048 ++nb_local_dofs;
3049 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
3050 ++nb_ghost_dofs;
3051 else
3052 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3053 "Impossible case. You could run problem on no distributed "
3054 "mesh. That is not implemented. Dof local index is %d",
3055 (*dit)->getPetscLocalDofIdx());
3056 }
3057
3058 // get indices
3059 auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
3060 const int nb_dofs = numered_dofs[s]->size();
3061 indices.clear();
3062 indices.reserve(nb_dofs);
3063 for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
3064 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3065 bool add = true;
3066 if (only_local) {
3067 if ((*dit)->getPetscLocalDofIdx() < 0 ||
3068 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
3069 add = false;
3070 }
3071 }
3072 if (add)
3073 indices.push_back(decltype(tag)::get_index(dit));
3074 }
3075 };
3076
3077 auto get_indices_by_uid = [&](auto tag, auto &indices) {
3078 const int nb_dofs = numered_dofs[s]->size();
3079 indices.clear();
3080 indices.reserve(nb_dofs);
3081 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3082 ++dit)
3083 indices.push_back(decltype(tag)::get_index(dit));
3084 };
3085
3086 auto get_sub_ao = [&](auto sub_data) {
3087 if (s == 0) {
3088 return sub_data->getSmartRowMap();
3089 } else {
3090 return sub_data->getSmartColMap();
3091 }
3092 };
3093
3094 auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
3095 if (s == 0) {
3096 sub_data->rowIs = is;
3097 sub_data->rowMap = ao;
3098 sub_data->colIs = is; // if square problem col is equal to row
3099 sub_data->colMap = ao;
3100 } else {
3101 sub_data->colIs = is;
3102 sub_data->colMap = ao;
3103 }
3104 };
3105
3106 auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
3107 if (s == 0) {
3108 if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
3109 sub_data->colIs = sub_data->getSmartRowIs();
3110 sub_data->colMap = sub_data->getSmartRowMap();
3111 }
3112 }
3113 };
3114
3115 auto concatenate_dofs = [&](auto tag, auto &indices,
3116 const auto local_only) {
3118 get_indices_by_tag(tag, indices, local_only);
3119
3120 SmartPetscObj<AO> ao;
3121 // Create AO from app indices (i.e. old), to pestc indices (new after
3122 // remove)
3123 if (local_only)
3124 ao = createAOMapping(m_field.get_comm(), indices.size(),
3125 &*indices.begin(), PETSC_NULLPTR);
3126 else
3127 ao = createAOMapping(PETSC_COMM_SELF, indices.size(),
3128 &*indices.begin(), PETSC_NULLPTR);
3129
3130 // Set mapping to sub dm data
3131 if (local_only) {
3132 if (auto sub_data = prb_ptr->getSubData()) {
3133 // create is and then map it to main problem of sub-problem
3134 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3135 &*indices.begin(), PETSC_COPY_VALUES);
3136 // get old app, i.e. oroginal befor sub indices, and ao, from app,
3137 // to petsc sub indices.
3138 auto sub_ao = get_sub_ao(sub_data);
3139 CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
3140 sub_ao = createAOMappingIS(sub_is, PETSC_NULLPTR);
3141 // set new sub ao
3142 set_sub_is_and_ao(sub_data, sub_is, sub_ao);
3143 apply_symmetry(sub_data);
3144 } else {
3145 // create sub data
3146 prb_ptr->getSubData() =
3147 boost::make_shared<Problem::SubProblemData>();
3148 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3149 &*indices.begin(), PETSC_COPY_VALUES);
3150 // set sub is ao
3151 set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
3152 apply_symmetry(prb_ptr->getSubData());
3153 }
3154 }
3155
3156 get_indices_by_uid(tag, indices);
3157 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3158
3160 };
3161
3162 // set indices index
3163 auto set_concatenated_indices = [&]() {
3164 std::vector<int> global_indices;
3165 std::vector<int> local_indices;
3167 CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3168 CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3169 auto gi = global_indices.begin();
3170 auto li = local_indices.begin();
3171 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3172 ++dit) {
3173 auto mod = NumeredDofEntity_part_and_all_indices_change(
3174 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3175 bool success = numered_dofs[s]->modify(dit, mod);
3176 if (!success)
3177 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3178 "can not set negative indices");
3179 ++gi;
3180 ++li;
3181 }
3183 };
3184 CHKERR set_concatenated_indices();
3185
3186 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3187 m_field.get_comm());
3188 *(local_nbdof_ptr[s]) = nb_local_dofs;
3189 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3190
3191 if (debug)
3192 for (auto dof : (*numered_dofs[s])) {
3193 if (dof->getPetscGlobalDofIdx() < 0) {
3194 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3195 "Negative global idx");
3196 }
3197 if (dof->getPetscLocalDofIdx() < 0) {
3198 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3199 "Negative local idx");
3200 }
3201 }
3202
3203 } else {
3204
3205 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3206 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3207 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3208 }
3209
3210 if (verb > QUIET) {
3212 "WORLD", Sev::inform,
3213 "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
3214 prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
3215 prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
3216 MOFEM_LOG_C("SYNC", Sev::verbose,
3217 "Removed DOFs from problem %s dofs [ %d / %d "
3218 "(before %d / %d) local, %d / %d (before %d / %d)]",
3219 prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
3220 prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
3221 nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
3222 prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
3223 nb_init_ghost_col_dofs);
3224 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
3225 }
3226
3228}
3229
3230MoFEMErrorCode ProblemsManager::removeDofsOnEntitiesNotDistributed(
3231 const std::string problem_name, const std::string field_name,
3232 const Range ents, const int lo_coeff, const int hi_coeff,
3233 const int lo_order, const int hi_order, int verb, const bool debug) {
3234
3235 MoFEM::Interface &m_field = cOre;
3237
3238 const Problem *prb_ptr;
3239 CHKERR m_field.get_problem(problem_name, &prb_ptr);
3240
3241 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
3242 prb_ptr->numeredRowDofsPtr, nullptr};
3243 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
3244 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
3245
3246 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
3247 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
3248 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
3249
3250 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
3251 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
3252 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
3253 // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
3254 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
3255 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
3256
3257 const std::array<int, 2> nb_init_dofs = {nb_init_row_dofs, nb_init_col_dofs};
3258
3259 for (int s = 0; s != 2; ++s)
3260 if (numered_dofs[s]) {
3261
3262 typedef multi_index_container<
3263
3264 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
3265
3266 >
3267 NumeredDofEntity_it_view_multiIndex;
3268
3269 const auto bit_number = m_field.get_field_bit_number(field_name);
3270 NumeredDofEntity_it_view_multiIndex dofs_it_view;
3271
3272 // Set -1 to global and local dofs indices
3273 for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
3274 ++pit) {
3275 auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
3276 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
3277 auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
3278 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
3279
3280 for (; lo != hi; ++lo)
3281 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
3282 (*lo)->getDofCoeffIdx() <= hi_coeff &&
3283 (*lo)->getDofOrder() >= lo_order &&
3284 (*lo)->getDofOrder() <= hi_order)
3285 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
3286 }
3287
3288 if (verb > QUIET) {
3289 for (auto &dof : dofs_it_view)
3290 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
3291 }
3292
3293 // set negative index
3294 auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
3295 for (auto dit : dofs_it_view) {
3296 bool success = numered_dofs[s]->modify(dit, mod);
3297 if (!success)
3298 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3299 "can not set negative indices");
3300 }
3301
3302 // create weak view
3303 std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
3304 dosf_weak_view.reserve(dofs_it_view.size());
3305 for (auto dit : dofs_it_view)
3306 dosf_weak_view.push_back(*dit);
3307
3308 if (verb >= NOISY)
3309 MOFEM_LOG_C("SYNC", Sev::noisy,
3310 "Number of DOFs in multi-index %d and to delete %d\n",
3311 numered_dofs[s]->size(), dofs_it_view.size());
3312
3313 // erase dofs from problem
3314 for (auto weak_dit : dosf_weak_view)
3315 if (auto dit = weak_dit.lock()) {
3316 numered_dofs[s]->erase(dit->getLocalUniqueId());
3317 }
3318
3319 if (verb >= NOISY)
3320 MOFEM_LOG_C("SYNC", Sev::noisy,
3321 "Number of DOFs in multi-index after delete %d\n",
3322 numered_dofs[s]->size());
3323
3324 // get current number of ghost dofs
3325 int nb_global_dof = 0;
3326 int nb_local_dofs = 0;
3327 int nb_ghost_dofs = 0;
3328
3329 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3330 ++dit) {
3331
3332 if ((*dit)->getDofIdx() >= 0) {
3333
3334 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
3335 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
3336 ++nb_local_dofs;
3337 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
3338 ++nb_ghost_dofs;
3339
3340 ++nb_global_dof;
3341 }
3342 }
3343
3344 if (debug) {
3345 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3346 m_field.get_comm());
3347 if (*(nbdof_ptr[s]) != nb_global_dof)
3348 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3349 "Number of local DOFs do not add up %d != %d",
3350 *(nbdof_ptr[s]), nb_global_dof);
3351 }
3352
3353 *(nbdof_ptr[s]) = nb_global_dof;
3354 *(local_nbdof_ptr[s]) = nb_local_dofs;
3355 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3356
3357 // get indices
3358 auto get_indices_by_tag = [&](auto tag) {
3359 std::vector<int> indices;
3360 indices.resize(nb_init_dofs[s], -1);
3361 for (auto dit = numered_dofs[s]->get<Idx_mi_tag>().lower_bound(0);
3362 dit != numered_dofs[s]->get<Idx_mi_tag>().end(); ++dit) {
3363 indices[(*dit)->getDofIdx()] = decltype(tag)::get_index(dit);
3364 }
3365 return indices;
3366 };
3367
3368 auto renumber = [&](auto tag, auto &indices) {
3370 int idx = 0;
3371 for (auto dit = numered_dofs[s]->get<decltype(tag)>().lower_bound(0);
3372 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3373 indices[(*dit)->getDofIdx()] = idx++;
3374 }
3376 };
3377
3378 auto get_sub_ao = [&](auto sub_data) {
3379 if (s == 0) {
3380 return sub_data->getSmartRowMap();
3381 } else {
3382 return sub_data->getSmartColMap();
3383 }
3384 };
3385
3386 auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
3387 if (s == 0) {
3388 sub_data->rowIs = is;
3389 sub_data->rowMap = ao;
3390 sub_data->colIs = is;
3391 sub_data->colMap = ao;
3392 } else {
3393 sub_data->colIs = is;
3394 sub_data->colMap = ao;
3395 }
3396 };
3397
3398 auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
3399 if (s == 0) {
3400 if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
3401 sub_data->colIs = sub_data->getSmartRowIs();
3402 sub_data->colMap = sub_data->getSmartRowMap();
3403 }
3404 }
3405 };
3406
3407 auto set_sub_data = [&](auto &indices) {
3409 if (auto sub_data = prb_ptr->getSubData()) {
3410 // create is and then map it to main problem of sub-problem
3411 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3412 &*indices.begin(), PETSC_COPY_VALUES);
3413 // get old app, i.e. oroginal befor sub indices, and ao, from
3414 // app, to petsc sub indices.
3415 auto sub_ao = get_sub_ao(sub_data);
3416 CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
3417 sub_ao = createAOMappingIS(sub_is, PETSC_NULLPTR);
3418 // set new sub ao
3419 set_sub_is_and_ao(sub_data, sub_is, sub_ao);
3420 apply_symmetry(sub_data);
3421 } else {
3422 prb_ptr->getSubData() = boost::make_shared<Problem::SubProblemData>();
3423 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3424 &*indices.begin(), PETSC_COPY_VALUES);
3425 auto sub_ao = createAOMappingIS(sub_is, PETSC_NULLPTR);
3426 // set sub is ao
3427 set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, sub_ao);
3428 apply_symmetry(prb_ptr->getSubData());
3429 }
3431 };
3432
3433 auto global_indices = get_indices_by_tag(PetscGlobalIdx_mi_tag());
3434 auto local_indices = get_indices_by_tag(PetscLocalIdx_mi_tag());
3435 CHKERR set_sub_data(global_indices);
3436 CHKERR renumber(PetscGlobalIdx_mi_tag(), global_indices);
3437 CHKERR renumber(PetscLocalIdx_mi_tag(), local_indices);
3438
3439 int i = 0;
3440 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3441 ++dit) {
3442 auto idx = (*dit)->getDofIdx();
3443 if (idx >= 0) {
3444 auto mod = NumeredDofEntity_part_and_all_indices_change(
3445 (*dit)->getPart(), i++, global_indices[idx], local_indices[idx]);
3446 bool success = numered_dofs[s]->modify(dit, mod);
3447 if (!success)
3448 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3449 "can not set negative indices");
3450 } else {
3451 auto mod = NumeredDofEntity_part_and_all_indices_change(
3452 (*dit)->getPart(), -1, -1, -1);
3453 bool success = numered_dofs[s]->modify(dit, mod);
3454 if (!success)
3455 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3456 "can not set negative indices");
3457 }
3458 };
3459
3460 if (debug) {
3461 for (auto dof : (*numered_dofs[s])) {
3462 if (dof->getDofIdx() >= 0 && dof->getPetscGlobalDofIdx() < 0) {
3463 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3464 "Negative global idx");
3465 }
3466 }
3467 }
3468
3469 } else {
3470
3471 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3472 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3473 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3474 }
3475
3476 if (verb >= NOISY)
3478
3479 if (verb > QUIET) {
3481 "WORLD", Sev::inform,
3482 "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
3483 prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
3484 prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
3485 MOFEM_LOG_C("SYNC", Sev::verbose,
3486 "Removed DOFs from problem %s dofs [ %d / %d "
3487 "(before %d / %d) local, %d / %d (before %d / %d)]",
3488 prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
3489 prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
3490 nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
3491 prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
3492 nb_init_ghost_col_dofs);
3494 }
3496}
3497
3498MoFEMErrorCode ProblemsManager::removeDofsOnEntities(
3499 const std::string problem_name, const std::string field_name,
3500 const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
3501 Range *ents_ptr, const int lo_coeff, const int hi_coeff, const int lo_order,
3502 const int hi_order, int verb, const bool debug) {
3503 MoFEM::Interface &m_field = cOre;
3505
3506 auto bit_manager = m_field.getInterface<BitRefManager>();
3507
3508 Range ents;
3509 if (ents_ptr) {
3510 ents = *ents_ptr;
3511 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3512 ents, verb);
3513 } else {
3514 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3515 verb);
3516 }
3517
3518 CHKERR removeDofsOnEntities(problem_name, field_name, ents, lo_coeff,
3519 hi_coeff, lo_order, hi_order, verb, debug);
3520
3522}
3523
3524MoFEMErrorCode ProblemsManager::removeDofsOnEntitiesNotDistributed(
3525 const std::string problem_name, const std::string field_name,
3526 const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask,
3527 Range *ents_ptr, const int lo_coeff, const int hi_coeff, const int lo_order,
3528 const int hi_order, int verb, const bool debug) {
3529 MoFEM::Interface &m_field = cOre;
3531
3532 auto bit_manager = m_field.getInterface<BitRefManager>();
3533
3534 Range ents;
3535 if (ents_ptr) {
3536 ents = *ents_ptr;
3537 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3538 ents, verb);
3539 } else {
3540 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3541 verb);
3542 }
3543
3544 CHKERR removeDofsOnEntitiesNotDistributed(problem_name, field_name, ents,
3545 lo_coeff, hi_coeff, lo_order,
3546 hi_order, verb, debug);
3547
3549}
3550
3552ProblemsManager::markDofs(const std::string problem_name, RowColData rc,
3553 const enum MarkOP op, const Range ents,
3554 std::vector<unsigned char> &marker) const {
3555
3556 Interface &m_field = cOre;
3557 const Problem *problem_ptr;
3559 CHKERR m_field.get_problem(problem_name, &problem_ptr);
3560 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3561 switch (rc) {
3562 case ROW:
3563 dofs = problem_ptr->getNumeredRowDofsPtr();
3564 break;
3565 case COL:
3566 dofs = problem_ptr->getNumeredColDofsPtr();
3567 default:
3568 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
3569 }
3570 marker.resize(dofs->size(), 0);
3571 std::vector<unsigned char> marker_tmp;
3572
3573 switch (op) {
3574 case MarkOP::OR:
3575 for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
3576 auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
3577 auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
3578 for (; lo != hi; ++lo) {
3579 if ((*lo)->getPetscLocalDofIdx() < 0)
3580 continue;
3581 marker[(*lo)->getPetscLocalDofIdx()] |= 1;
3582 }
3583 }
3584 break;
3585 case MarkOP::AND:
3586 marker_tmp.resize(dofs->size(), 0);
3587 for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
3588 auto lo = dofs->get<Ent_mi_tag>().lower_bound(p->first);
3589 auto hi = dofs->get<Ent_mi_tag>().upper_bound(p->second);
3590 for (; lo != hi; ++lo) {
3591 if ((*lo)->getPetscLocalDofIdx() < 0)
3592 continue;
3593 marker_tmp[(*lo)->getPetscLocalDofIdx()] = 1;
3594 }
3595 }
3596 for (int i = 0; i != marker.size(); ++i) {
3597 marker[i] &= marker_tmp[i];
3598 }
3599 break;
3600 }
3601
3603}
3604
3605MoFEMErrorCode ProblemsManager::modifyMarkDofs(
3606 const std::string problem_name, RowColData rc, const std::string field_name,
3607 const int lo, const int hi, const enum ProblemsManager::MarkOP op,
3608 const unsigned char c, std::vector<unsigned char> &marker) const {
3609
3610 Interface &m_field = cOre;
3611 const Problem *problem_ptr;
3613 CHKERR m_field.get_problem(problem_name, &problem_ptr);
3614 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3615 switch (rc) {
3616 case ROW:
3617 dofs = problem_ptr->getNumeredRowDofsPtr();
3618 break;
3619 case COL:
3620 dofs = problem_ptr->getNumeredColDofsPtr();
3621 default:
3622 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
3623 }
3624 marker.resize(dofs->size(), 0);
3625
3626 auto dof_lo = dofs->get<Unique_mi_tag>().lower_bound(
3627 FieldEntity::getLoBitNumberUId(m_field.get_field_bit_number(field_name)));
3628 auto dof_hi = dofs->get<Unique_mi_tag>().upper_bound(
3629 FieldEntity::getHiBitNumberUId(m_field.get_field_bit_number(field_name)));
3630
3631 switch (op) {
3632 case MarkOP::OR:
3633 for (; dof_lo != dof_hi; ++dof_lo) {
3634 if ((*dof_lo)->getPetscLocalDofIdx() < 0)
3635 continue;
3636 if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3637 (*dof_lo)->getDofCoeffIdx() <= hi)
3638 marker[(*dof_lo)->getPetscLocalDofIdx()] |= c;
3639 }
3640 break;
3641 case MarkOP::AND:
3642 for (; dof_lo != dof_hi; ++dof_lo) {
3643 if ((*dof_lo)->getPetscLocalDofIdx() < 0)
3644 continue;
3645 if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3646 (*dof_lo)->getDofCoeffIdx() <= hi)
3647 marker[(*dof_lo)->getPetscLocalDofIdx()] &= c;
3648 }
3649 break;
3650 }
3651
3653}
3654
3655MoFEMErrorCode ProblemsManager::markDofs(
3656
3657 const std::string problem_name, RowColData rc,
3658 std::vector<boost::weak_ptr<NumeredDofEntity>> &vec_dof_view,
3659 const enum MarkOP op, std::vector<unsigned char> &marker
3660
3661) const {
3662 Interface &m_field = cOre;
3664
3665 const Problem *problem_ptr;
3666 CHKERR m_field.get_problem(problem_name, &problem_ptr);
3667 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3668 switch (rc) {
3669 case ROW:
3670 dofs = problem_ptr->getNumeredRowDofsPtr();
3671 break;
3672 case COL:
3673 dofs = problem_ptr->getNumeredColDofsPtr();
3674 break;
3675 default:
3676 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Should be row or column");
3677 }
3678 marker.resize(dofs->size(), 0);
3679
3680 for (auto &dof : vec_dof_view) {
3681 if (auto dof_ptr = dof.lock()) {
3682 if (dof_ptr->getPetscLocalDofIdx() < 0)
3683 continue;
3684 if (op == MarkOP::OR)
3685 marker[dof_ptr->getPetscLocalDofIdx()] |= 1;
3686 else
3687 marker[dof_ptr->getPetscLocalDofIdx()] &= 1;
3688 }
3689 }
3690
3692}
3693
3695ProblemsManager::addFieldToEmptyFieldBlocks(const std::string problem_name,
3696 const std::string row_field,
3697 const std::string col_field) const {
3698
3699 Interface &m_field = cOre;
3701
3702 const auto problem_ptr = m_field.get_problem(problem_name);
3703 auto get_field_id = [&](const std::string field_name) {
3704 return m_field.get_field_structure(field_name)->getId();
3705 };
3706 const auto row_id = get_field_id(row_field);
3707 const auto col_id = get_field_id(col_field);
3708
3709 problem_ptr->addFieldToEmptyFieldBlocks(BlockFieldPair(row_id, col_id));
3710
3712}
3713
3714} // namespace MoFEM
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define ProblemManagerFunctionBegin
constexpr double a
@ VERY_VERBOSE
@ QUIET
@ NOISY
@ VERBOSE
RowColData
RowColData.
@ COL
@ ROW
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
static const bool debug
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< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
auto marker
set bit to marker
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
int DofIdx
Index of DOF.
Definition Types.hpp:18
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
uint128_t UId
Unique Id.
Definition Types.hpp:31
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
auto createISGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
Creates a data structure for an index set containing a list of integers.
auto createAOMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[])
Creates an application mapping using two integer arrays.
CoreTmp< 0 > Core
Definition Core.hpp:1148
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
MoFEM::LogManager::SeverityLevel Sev
Problem::BlockFieldPair BlockFieldPair
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
DeprecatedCoreInterface Interface
int r
Definition sdf.py:8
constexpr auto field_name
FTensor::Index< 'm', 3 > m
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual std::string get_field_name(const BitFieldId id) const =0
get field name from id
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
Deprecated interface functions.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients per DOF.
std::map< int, BaseFunction::DofsSideMap > & getDofSideMap() const
Get DOF side mapping for broken spaces.
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
BitProblemId getId() const
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.