22#define ProblemManagerFunctionBegin \
24 MOFEM_LOG_CHANNEL("WORLD"); \
25 MOFEM_LOG_CHANNEL("SYNC"); \
26 MOFEM_LOG_FUNCTION(); \
27 MOFEM_LOG_TAG("SYNC", "ProblemsManager"); \
28 MOFEM_LOG_TAG("WORLD", "ProblemsManager")
32 bcopy(&uid,
dAta, 4 *
sizeof(
int));
43 int global_dof =
pTr[4];
47 unsigned int b0, b1, b2, b3;
48 bcopy(&
pTr[0], &b0,
sizeof(
int));
49 bcopy(&
pTr[1], &b1,
sizeof(
int));
50 bcopy(&
pTr[2], &b2,
sizeof(
int));
51 bcopy(&
pTr[3], &b3,
sizeof(
int));
52 UId uid =
static_cast<UId>(b0) |
static_cast<UId>(b1) << 8 *
sizeof(
int) |
53 static_cast<UId>(b2) << 16 *
sizeof(
int) |
54 static_cast<UId>(b3) << 24 *
sizeof(
int);
71 buildProblemFromFields(PETSC_FALSE),
72 synchroniseProblemEntities(PETSC_FALSE) {
79 CHKERR PetscOptionsBegin(m_field.
get_comm(),
"",
"Problem manager",
"none");
82 "-problem_build_from_fields",
83 "Add DOFs to problem directly from fields not through DOFs on elements",
86 ierr = PetscOptionsEnd();
92 const Range &ents,
const int dim,
const int adj_dim,
const int n_parts,
93 Tag *th_vertex_weights, Tag *th_edge_weights, Tag *th_part_weights,
94 int verb,
const bool debug) {
96 .getInterface<CommInterface>()
97 ->partitionMesh(ents,
dim, adj_dim, n_parts, th_vertex_weights,
98 th_edge_weights, th_part_weights, verb,
debug);
102 const bool square_matrix,
122 const bool square_matrix,
125 auto dofs_field_ptr = m_field.
get_dofs();
135 SETERRQ1(PETSC_COMM_SELF, 1,
"problem <%s> refinement level not set",
136 problem_ptr->
getName().c_str());
144 auto make_rows_and_cols_view = [&](
auto &dofs_rows,
auto &dofs_cols) {
146 for (
auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
148 const auto uid = (*it)->getLocalUniqueId();
151 for (
auto lo =
r.first; lo !=
r.second; ++lo) {
153 if ((lo->getBitFEId() & problem_ptr->
getBitFEId()).any()) {
154 std::array<bool, 2> row_col = {
false,
false};
158 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
161 if ((fe_bit & prb_mask) != fe_bit)
163 if ((fe_bit & prb_bit).none())
166 auto add_to_view = [&](
auto &nb_dofs,
auto &view,
auto rc) {
167 auto dit = nb_dofs->lower_bound(uid);
168 decltype(dit) hi_dit;
169 if (dit != nb_dofs->end()) {
170 hi_dit = nb_dofs->upper_bound(
172 view.insert(dit, hi_dit);
177 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
178 add_to_view(dofs_field_ptr, dofs_rows,
ROW);
181 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
182 add_to_view(dofs_field_ptr, dofs_cols,
COL);
184 if (square_matrix && row_col[
ROW])
186 else if (row_col[
ROW] && row_col[
COL])
194 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
204 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
206 hi_miit = dofs_rows.get<0>().end();
208 int count_dofs = dofs_rows.get<1>().count(
true);
209 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
210 boost::shared_ptr<std::vector<NumeredDofEntity>>(
211 new std::vector<NumeredDofEntity>());
213 dofs_array->reserve(count_dofs);
214 miit = dofs_rows.get<0>().begin();
215 for (; miit != hi_miit; miit++) {
216 if ((*miit)->getActive()) {
217 dofs_array->emplace_back(*miit);
218 dofs_array->back().dofIdx = (problem_ptr->
nbDofsRow)++;
222 for (
auto &
v : *dofs_array) {
228 if (!square_matrix) {
235 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
237 hi_miit = dofs_cols.get<0>().end();
240 miit = dofs_cols.get<0>().begin();
241 for (; miit != hi_miit; miit++) {
242 if (!(*miit)->getActive()) {
248 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
249 boost::shared_ptr<std::vector<NumeredDofEntity>>(
250 new std::vector<NumeredDofEntity>());
252 dofs_array->reserve(count_dofs);
253 miit = dofs_cols.get<0>().begin();
254 for (; miit != hi_miit; miit++) {
255 if (!(*miit)->getActive()) {
258 dofs_array->emplace_back(*miit);
259 dofs_array->back().dofIdx = problem_ptr->
nbDofsCol++;
262 for (
auto &
v : *dofs_array) {
274 << problem_ptr->
getName() <<
" Nb. local dofs "
282 <<
"FEs row dofs " << *problem_ptr <<
" Nb. row dof "
288 <<
"FEs col dofs " << *problem_ptr <<
" Nb. col dof "
305 const std::string name,
const bool square_matrix,
int verb) {
319 square_matrix, verb);
328 Problem *problem_ptr,
const bool square_matrix,
int verb) {
334 auto dofs_field_ptr = m_field.
get_dofs();
345 "problem <%s> refinement level not set",
346 problem_ptr->
getName().c_str());
354 boost::shared_ptr<NumeredDofEntity_multiIndex>(
368 auto make_rows_and_cols_view = [&](
auto &dofs_rows,
auto &dofs_cols) {
372 for (
auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
375 const auto uid = (*it)->getLocalUniqueId();
378 for (
auto lo =
r.first; lo !=
r.second; ++lo) {
380 if ((lo->getBitFEId() & problem_ptr->
getBitFEId()).any()) {
381 std::array<bool, 2> row_col = {
false,
false};
383 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
386 if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
388 auto add_to_view = [&](
auto dofs,
auto &view,
auto rc) {
389 auto dit = dofs->lower_bound(uid);
390 decltype(dit) hi_dit;
391 if (dit != dofs->end()) {
392 hi_dit = dofs->upper_bound(
394 view.insert(dit, hi_dit);
399 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
400 add_to_view(dofs_field_ptr, dofs_rows,
ROW);
403 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
404 add_to_view(dofs_field_ptr, dofs_cols,
COL);
406 if (square_matrix && row_col[
ROW])
408 else if (row_col[
ROW] && row_col[
COL])
417 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
425 for (
auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
426 if ((fit->get()->getId() & problem_ptr->
getBitFEId()).any()) {
427 fields_ids_row |= fit->get()->getBitFieldIdRow();
428 fields_ids_col |= fit->get()->getBitFieldIdCol();
432 for (
auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
433 if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
437 auto hi_dit = dofs_field_ptr->get<
Unique_mi_tag>().upper_bound(
440 for (; dit != hi_dit; dit++) {
442 const int owner_proc = dit->get()->getOwnerProc();
444 const unsigned char pstatus = dit->get()->getPStatus();
450 const auto &dof_bit = (*dit)->getBitRefLevel();
452 if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
454 if ((fit->get()->getId() & fields_ids_row).any()) {
455 dofs_rows.insert(*dit);
457 if (!square_matrix) {
458 if ((fit->get()->getId() & fields_ids_col).any()) {
459 dofs_cols.insert(*dit);
472 BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
473 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
474 fit != fe_ptr->end(); fit++) {
475 if ((fit->get()->getId() & problem_ptr->
getBitFEId()).any()) {
476 fields_ids_row |= fit->get()->getBitFieldIdRow();
477 fields_ids_col |= fit->get()->getBitFieldIdCol();
482 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
483 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
485 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
486 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
487 Range ents_to_synchronise;
488 for (; miit != hi_miit; ++miit) {
489 if (miit->get()->getEntDofIdx() != 0)
491 ents_to_synchronise.insert(miit->get()->getEnt());
493 Range tmp_ents = ents_to_synchronise;
495 ents_to_synchronise, verb);
496 ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
497 for (
auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
498 if ((fit->get()->getId() & *fields_ids[ss]).any()) {
499 const auto bit_number = (*fit)->getBitNumber();
500 for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
501 pit != ents_to_synchronise.pair_end(); ++pit) {
502 const auto f = pit->first;
503 const auto s = pit->second;
509 auto dit = dofs_field_ptr->get<
Unique_mi_tag>().lower_bound(lo_uid);
513 dofs_ptr[ss]->insert(dit, hi_dit);
522 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
529 for (
int ss = 0; ss < 2; ss++) {
530 *(nbdof_ptr[ss]) = 0;
531 *(local_nbdof_ptr[ss]) = 0;
532 *(ghost_nbdof_ptr[ss]) = 0;
537 int nb_local_dofs[] = {0, 0};
538 for (
int ss = 0; ss < loop_size; ss++) {
539 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
541 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
542 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
543 for (; miit != hi_miit; miit++) {
544 int owner_proc = (*miit)->getOwnerProc();
552 int start_ranges[2], end_ranges[2];
553 for (
int ss = 0; ss != loop_size; ss++) {
556 CHKERR PetscLayoutSetBlockSize(layout, 1);
557 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
558 CHKERR PetscLayoutSetUp(layout);
559 CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
560 CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
561 CHKERR PetscLayoutDestroy(&layout);
564 nbdof_ptr[1] = nbdof_ptr[0];
565 nb_local_dofs[1] = nb_local_dofs[0];
566 start_ranges[1] = start_ranges[0];
567 end_ranges[1] = end_ranges[0];
579 const size_t idx_data_type_size =
sizeof(
IdxDataType);
580 const size_t data_block_size = idx_data_type_size /
sizeof(
int);
585 std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
591 for (
int ss = 0; ss != loop_size; ++ss) {
593 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
595 hi_miit = dofs_ptr[ss]->get<0>().end();
597 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
598 boost::shared_ptr<std::vector<NumeredDofEntity>>(
599 new std::vector<NumeredDofEntity>());
600 int nb_dofs_to_add = 0;
601 miit = dofs_ptr[ss]->get<0>().begin();
602 for (; miit != hi_miit; ++miit) {
604 if (!(miit->get()->getActive()))
608 dofs_array->reserve(nb_dofs_to_add);
615 int &local_idx = *local_nbdof_ptr[ss];
616 miit = dofs_ptr[ss]->get<0>().begin();
617 for (; miit != hi_miit; ++miit) {
620 if (!(miit->get()->getActive()))
623 dofs_array->emplace_back(*miit);
625 int owner_proc = dofs_array->back().getOwnerProc();
626 if (owner_proc < 0) {
628 "data inconsistency");
632 dofs_array->back().pArt = owner_proc;
633 dofs_array->back().dofIdx = -1;
634 dofs_array->back().petscGloablDofIdx = -1;
635 dofs_array->back().petscLocalDofIdx = -1;
639 int glob_idx = start_ranges[ss] + local_idx;
640 dofs_array->back().pArt = owner_proc;
641 dofs_array->back().dofIdx = glob_idx;
642 dofs_array->back().petscGloablDofIdx = glob_idx;
643 dofs_array->back().petscLocalDofIdx = local_idx;
646 unsigned char pstatus = dofs_array->back().getPStatus();
652 proc < MAX_SHARING_PROCS &&
653 -1 != dofs_array->back().getSharingProcsPtr()[proc];
657 ids_data_packed_rows[dofs_array->back()
658 .getSharingProcsPtr()[proc]]
659 .emplace_back(dofs_array->back().getGlobalUniqueId(),
662 ids_data_packed_cols[dofs_array->back()
663 .getSharingProcsPtr()[proc]]
664 .emplace_back(dofs_array->back().getGlobalUniqueId(),
667 if (!(pstatus & PSTATUS_MULTISHARED)) {
675 auto hint = numered_dofs_ptr[ss]->end();
676 for (
auto &
v : *dofs_array)
677 hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &
v);
680 local_nbdof_ptr[1] = local_nbdof_ptr[0];
683 int nsends_rows = 0, nsends_cols = 0;
687 lengths_rows.clear();
688 lengths_cols.clear();
690 lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
691 lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
692 if (!ids_data_packed_rows[proc].empty())
694 if (!ids_data_packed_cols[proc].empty())
716 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
722 CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
723 &lengths_rows[0], &onodes_rows,
731 CHKERR PetscCommGetNewTag(comm, &tag_row);
735 MPI_Request *r_waits_row;
740 CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
741 olengths_rows, &rbuf_row, &r_waits_row);
742 CHKERR PetscFree(onodes_rows);
744 MPI_Request *s_waits_row;
745 CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
748 for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
749 if (!lengths_rows[proc])
751 CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0],
754 tag_row, comm, s_waits_row + kk);
759 CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
762 CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
765 CHKERR PetscFree(r_waits_row);
766 CHKERR PetscFree(s_waits_row);
770 int nrecvs_cols = nrecvs_rows;
771 int *olengths_cols = olengths_rows;
772 PetscInt **rbuf_col = rbuf_row;
773 if (!square_matrix) {
776 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
782 CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
783 &lengths_cols[0], &onodes_cols,
788 CHKERR PetscCommGetNewTag(comm, &tag_col);
792 MPI_Request *r_waits_col;
793 CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
794 olengths_cols, &rbuf_col, &r_waits_col);
795 CHKERR PetscFree(onodes_cols);
797 MPI_Request *s_waits_col;
798 CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
801 for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
802 if (!lengths_cols[proc])
804 CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0],
807 tag_col, comm, s_waits_col + kk);
812 CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
815 CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
818 CHKERR PetscFree(r_waits_col);
819 CHKERR PetscFree(s_waits_col);
822 CHKERR PetscCommDestroy(&comm);
826 auto hint = dofs_glob_uid_view.begin();
827 for (
auto dof : *m_field.
get_dofs())
828 dofs_glob_uid_view.emplace_hint(hint, dof);
831 for (
int ss = 0; ss != loop_size; ++ss) {
837 nrecvs = nrecvs_rows;
838 olengths = olengths_rows;
839 data_procs = rbuf_row;
841 nrecvs = nrecvs_cols;
842 olengths = olengths_cols;
843 data_procs = rbuf_col;
847 for (
int kk = 0; kk != nrecvs; ++kk) {
848 int len = olengths[kk];
849 int *data_from_proc = data_procs[kk];
850 for (
int dd = 0;
dd < len;
dd += data_block_size) {
852 auto ddit = dofs_glob_uid_view.find(uid);
854 if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
858 <<
"DOF is shared or multishared between processors. For example "
859 "if order of field on given entity is set in inconsistently, "
860 "has different value on two processor, error such as this is "
863 MOFEM_LOG(
"SELF", Sev::error) <<
"UId " << uid <<
" is not found";
866 <<
"Problematic UId owner proc is " << owner_proc;
869 <<
"Problematic UId entity owning processor handle is "
870 << uid_handle <<
" entity type "
872 const auto uid_bit_number =
875 <<
"Problematic UId field is "
880 field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
882 owner_proc, uid_bit_number, uid_handle));
883 if (fe_it == field_view.end()) {
885 <<
"Also, no field entity in database for given global UId";
887 MOFEM_LOG(
"SELF", Sev::error) <<
"Field entity in databse exist "
888 "(but have no DOF wih give UId";
889 MOFEM_LOG(
"SELF", Sev::error) << **fe_it;
892 auto error_file_name =
893 "error_with_missing_entity_" +
897 <<
"Look to file < " << error_file_name
898 <<
" > it contains entity with missing DOF.";
901 const auto local_fe_ent = (*fe_it)->getEnt();
902 CHKERR m_field.
get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
903 CHKERR m_field.
get_moab().write_file(error_file_name.c_str(),
"VTK",
904 "", tmp_msh->get_ptr(), 1);
909 "DOF with global UId not found (Compile code in Debug to "
910 "learn more about problem");
913 auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
915 if (dit != numered_dofs_ptr[ss]->end()) {
919 if (PetscUnlikely(global_idx < 0))
921 "received negative dof");
924 success = numered_dofs_ptr[ss]->modify(
928 if (PetscUnlikely(!success))
930 "modification unsuccessful");
932 success = numered_dofs_ptr[ss]->modify(
936 if (PetscUnlikely(!success))
938 "modification unsuccessful");
943 if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
950 std::ostringstream zz;
951 zz << **ddit << std::endl;
953 "data inconsistency, dofs is not shared, but received "
968 CHKERR PetscFree(olengths_rows);
969 CHKERR PetscFree(rbuf_row[0]);
970 CHKERR PetscFree(rbuf_row);
971 if (!square_matrix) {
972 CHKERR PetscFree(olengths_cols);
973 CHKERR PetscFree(rbuf_col[0]);
974 CHKERR PetscFree(rbuf_col);
978 if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
980 "data inconsistency for square_matrix %d!=%d",
981 numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
985 "data inconsistency for square_matrix");
998 const std::string out_name,
const std::vector<std::string> &fields_row,
999 const std::vector<std::string> &fields_col,
const std::string main_problem,
1000 const bool square_matrix,
1001 const map<std::string, std::pair<EntityType, EntityType>> *entityMapRow,
1002 const map<std::string, std::pair<EntityType, EntityType>> *entityMapCol,
1011 using ProblemByName =
decltype(problems_ptr->get<
Problem_mi_tag>());
1012 auto &problems_by_name =
1013 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1016 auto out_problem_it = problems_by_name.find(out_name);
1017 if (out_problem_it == problems_by_name.end()) {
1019 "subproblem with name < %s > not defined (top tip check spelling)",
1024 auto main_problem_it = problems_by_name.find(main_problem);
1025 if (main_problem_it == problems_by_name.end()) {
1027 "problem of subproblem with name < %s > not defined (top tip "
1029 main_problem.c_str());
1033 boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1034 out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1036 boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1037 main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1039 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1040 &out_problem_it->nbLocDofsCol};
1042 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1046 out_problem_it->nbGhostDofsRow = 0;
1047 out_problem_it->nbGhostDofsCol = 0;
1051 std::vector<std::string> fields[] = {fields_row, fields_col};
1052 const map<std::string, std::pair<EntityType, EntityType>> *entityMap[] = {
1053 entityMapRow, entityMapCol};
1056 out_problem_it->subProblemData =
1057 boost::make_shared<Problem::SubProblemData>();
1060 for (
int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1063 (*nb_local_dofs[ss]) = 0;
1066 out_problem_dofs[ss]->clear();
1069 out_problem_it->numeredFiniteElementsPtr->clear();
1072 for (
auto field : fields[ss]) {
1078 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1079 boost::make_shared<std::vector<NumeredDofEntity>>();
1082 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1084 out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1088 auto dit = main_problem_dofs[ss]->get<
Unique_mi_tag>().lower_bound(
1090 auto hi_dit = main_problem_dofs[ss]->get<
Unique_mi_tag>().upper_bound(
1093 auto add_dit_to_dofs_array = [&](
auto &dit) {
1094 if (dit->get()->getPetscGlobalDofIdx() >= 0)
1095 dofs_array->emplace_back(
1096 dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1097 dit->get()->getPetscGlobalDofIdx(),
1098 dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1101 if (entityMap[ss]) {
1102 auto mit = entityMap[ss]->find(field);
1103 if (mit != entityMap[ss]->end()) {
1107 for (
auto diit = dit; diit != hi_dit; ++diit) {
1109 if (ent_type >= lo_type && ent_type <= hi_type)
1112 dofs_array->reserve(count);
1113 for (; dit != hi_dit; ++dit) {
1115 if (ent_type >= lo_type && ent_type <= hi_type)
1116 add_dit_to_dofs_array(dit);
1119 dofs_array->reserve(std::distance(dit, hi_dit));
1120 for (; dit != hi_dit; dit++)
1121 add_dit_to_dofs_array(dit);
1124 dofs_array->reserve(std::distance(dit, hi_dit));
1125 for (; dit != hi_dit; dit++)
1126 add_dit_to_dofs_array(dit);
1130 auto hint = out_problem_dofs[ss]->end();
1131 for (
auto &
v : *dofs_array)
1132 hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &
v);
1136 auto dit = out_problem_dofs[ss]->get<
Idx_mi_tag>().begin();
1137 auto hi_dit = out_problem_dofs[ss]->get<
Idx_mi_tag>().end();
1138 for (; dit != hi_dit; dit++) {
1140 if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
1141 idx = (*nb_local_dofs[ss])++;
1143 bool success = out_problem_dofs[ss]->modify(
1144 out_problem_dofs[ss]->project<0>(dit),
1148 "operation unsuccessful");
1158 out_problem_dofs[ss]->size());
1159 const int nb = std::distance(dit, hi_dit);
1161 std::vector<int> main_indices(nb);
1162 for (
auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1163 *it = dit->get()->getPetscGlobalDofIdx();
1167 CHKERR ISCreateGeneral(m_field.
get_comm(), nb, &*main_indices.begin(),
1168 PETSC_USE_POINTER, &is);
1172 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1175 CHKERR ISDuplicate(is, &is_dup);
1180 CHKERR ISDuplicate(is, &is_dup);
1184 CHKERR AOApplicationToPetscIS(ao, is);
1186 CHKERR ISGetSize(is, nb_dofs[ss]);
1190 for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1192 bool success = out_problem_dofs[ss]->modify(
1193 out_problem_dofs[ss]->project<0>(dit),
1195 dit->get()->getPart(), *it, *it,
1196 dit->get()->getPetscLocalDofIdx()));
1199 "operation unsuccessful");
1204 NumeredDofEntityByLocalIdx::iterator dit =
1206 NumeredDofEntityByLocalIdx::iterator hi_dit =
1208 const int nb = std::distance(dit, hi_dit);
1209 std::vector<int> main_indices_non_local(nb);
1210 for (
auto it = main_indices_non_local.begin(); dit != hi_dit;
1212 *it = dit->get()->getPetscGlobalDofIdx();
1216 &*main_indices_non_local.begin(),
1217 PETSC_USE_POINTER, &is);
1218 CHKERR AOApplicationToPetscIS(ao, is);
1221 for (
auto it = main_indices_non_local.begin(); dit != hi_dit;
1223 bool success = out_problem_dofs[ss]->modify(
1224 out_problem_dofs[ss]->project<0>(dit),
1226 dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1227 dit->get()->getPetscLocalDofIdx()));
1230 "operation unsuccessful");
1238 if (square_matrix) {
1239 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1240 out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1241 out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1242 out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1243 out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1253 const std::string out_name,
const std::vector<std::string> add_row_problems,
1254 const std::vector<std::string> add_col_problems,
const bool square_matrix,
1269 ProblemByName &problems_by_name =
1270 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1273 ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1274 if (out_problem_it == problems_by_name.end()) {
1276 "problem with name < %s > not defined (top tip check spelling)",
1280 out_problem_it->composedProblemsData =
1281 boost::make_shared<ComposedProblemsData>();
1282 boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1283 out_problem_it->getComposedProblemsData();
1285 const std::vector<std::string> *add_prb[] = {&add_row_problems,
1287 std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1288 &cmp_prb_data->colProblemsAdd};
1289 std::vector<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1292 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1293 &out_problem_it->nbLocDofsCol};
1295 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1299 out_problem_it->nbDofsRow = 0;
1300 out_problem_it->nbDofsCol = 0;
1301 out_problem_it->nbLocDofsRow = 0;
1302 out_problem_it->nbLocDofsCol = 0;
1303 out_problem_it->nbGhostDofsRow = 0;
1304 out_problem_it->nbGhostDofsCol = 0;
1306 int nb_dofs_reserve[] = {0, 0};
1309 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1313 add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1314 add_prb_is[ss]->reserve(add_prb[ss]->size());
1315 for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1316 vit != add_prb[ss]->end(); vit++) {
1317 ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1318 if (prb_it == problems_by_name.end()) {
1321 "problem with name < %s > not defined (top tip check spelling)",
1324 add_prb_ptr[ss]->push_back(&*prb_it);
1328 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1329 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1330 nb_dofs_reserve[ss] +=
1331 add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1334 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1335 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1336 nb_dofs_reserve[ss] +=
1337 add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1342 if (square_matrix) {
1343 add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1344 add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1345 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1346 *nb_dofs[1] = *nb_dofs[0];
1347 *nb_local_dofs[1] = *nb_local_dofs[0];
1351 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1353 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1354 dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1355 dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1357 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1359 out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1363 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1364 NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1369 for (
unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1370 PetscInt *dofs_out_idx_ptr;
1371 int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1372 CHKERR PetscMalloc(nb_local_dofs *
sizeof(
int), &dofs_out_idx_ptr);
1374 dit = (*add_prb_ptr[ss])[pp]
1377 hi_dit = (*add_prb_ptr[ss])[pp]
1381 dit = (*add_prb_ptr[ss])[pp]
1384 hi_dit = (*add_prb_ptr[ss])[pp]
1389 for (; dit != hi_dit; dit++) {
1390 const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
1391 const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
1392 const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
1393 if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1396 const int part = dit->get()->getPart();
1397 const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1399 (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1401 dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1402 glob_idx, loc_idx, part);
1404 dofs_out_idx_ptr[is_nb++] = glob_idx;
1407 if (is_nb > nb_local_dofs) {
1409 "Data inconsistency");
1412 CHKERR ISCreateGeneral(m_field.
get_comm(), is_nb, dofs_out_idx_ptr,
1413 PETSC_OWN_POINTER, &is);
1414 (*add_prb_is[ss]).push_back(is);
1416 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1417 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1419 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1420 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1422 if (square_matrix) {
1423 (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1424 (*add_prb_is[1]).push_back(is);
1425 CHKERR PetscObjectReference((PetscObject)is);
1430 if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1435 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1436 auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1437 : out_problem_it->numeredColDofsPtr->end();
1438 for (
auto &
v : *dofs_array[ss])
1439 hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1440 hint, dofs_array[ss], &
v)
1441 : out_problem_it->numeredColDofsPtr->emplace_hint(
1442 hint, dofs_array[ss], &
v);
1448 *nb_local_dofs[0] = 0;
1449 *nb_local_dofs[1] = 0;
1450 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1452 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1454 dofs_ptr = out_problem_it->numeredRowDofsPtr;
1456 dofs_ptr = out_problem_it->numeredColDofsPtr;
1458 NumeredDofEntityByUId::iterator dit, hi_dit;
1461 std::vector<int> idx;
1462 idx.reserve(std::distance(dit, hi_dit));
1464 for (; dit != hi_dit; dit++) {
1465 if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
1470 "modification unsuccessful");
1472 idx.push_back(dit->get()->getPetscGlobalDofIdx());
1474 if (dit->get()->getPetscLocalDofIdx() != -1) {
1476 "local index should be negative");
1480 if (square_matrix) {
1481 *nb_local_dofs[1] = *nb_local_dofs[0];
1486 CHKERR ISCreateGeneral(m_field.
get_comm(), idx.size(), &*idx.begin(),
1487 PETSC_USE_POINTER, &is);
1488 CHKERR ISGetSize(is, nb_dofs[ss]);
1489 if (square_matrix) {
1490 *nb_dofs[1] = *nb_dofs[0];
1494 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1495 for (
unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1496 CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1500 std::vector<int> idx_new;
1501 idx_new.reserve(dofs_ptr->size());
1502 for (NumeredDofEntityByUId::iterator dit =
1505 idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1510 &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1511 CHKERR AOApplicationToPetscIS(ao, is_new);
1513 std::vector<int>::iterator vit = idx_new.begin();
1514 for (NumeredDofEntityByUId::iterator dit =
1519 dit->get()->getPart(), *(vit++)));
1522 "modification unsuccessful");
1525 CHKERR ISDestroy(&is_new);
1556 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Simple partition problem " << name;
1560 ProblemByName &problems_set =
1561 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1562 ProblemByName::iterator p_miit = problems_set.find(name);
1563 if (p_miit == problems_set.end()) {
1564 SETERRQ1(PETSC_COMM_SELF, 1,
1565 "problem < %s > is not found (top tip: check spelling)",
1570 NumeredDofEntitysByIdx &dofs_row_by_idx =
1571 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>();
1572 NumeredDofEntitysByIdx &dofs_col_by_idx =
1573 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>();
1580 DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1581 DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1582 nb_row_local_dofs = 0;
1583 nb_row_ghost_dofs = 0;
1584 DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1585 DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1586 nb_col_local_dofs = 0;
1587 nb_col_ghost_dofs = 0;
1589 bool square_matrix =
false;
1590 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1591 square_matrix =
true;
1595 PetscLayout layout_row;
1596 const int *ranges_row;
1598 DofIdx nb_dofs_row = dofs_row_by_idx.size();
1600 CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1601 CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1602 CHKERR PetscLayoutSetUp(layout_row);
1603 CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1605 PetscLayout layout_col;
1606 const int *ranges_col;
1607 if (!square_matrix) {
1608 DofIdx nb_dofs_col = dofs_col_by_idx.size();
1610 CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1611 CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1612 CHKERR PetscLayoutSetUp(layout_col);
1613 CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1617 miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1618 hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1619 if (std::distance(miit_row, hi_miit_row) !=
1620 ranges_row[part + 1] - ranges_row[part]) {
1622 PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1623 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1624 "rstart (%d != %d - %d = %d) ",
1625 std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1626 ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1629 for (; miit_row != hi_miit_row; miit_row++) {
1630 bool success = dofs_row_by_idx.modify(
1635 "modification unsuccessful");
1637 success = dofs_row_by_idx.modify(
1641 "modification unsuccessful");
1644 if (!square_matrix) {
1645 miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1646 hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1647 if (std::distance(miit_col, hi_miit_col) !=
1648 ranges_col[part + 1] - ranges_col[part]) {
1649 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1650 "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1652 "rstart (%d != %d - %d = %d) ",
1653 std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1654 ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1657 for (; miit_col != hi_miit_col; miit_col++) {
1658 bool success = dofs_col_by_idx.modify(
1660 part, (*miit_col)->dofIdx));
1663 "modification unsuccessful");
1665 success = dofs_col_by_idx.modify(
1669 "modification unsuccessful");
1674 CHKERR PetscLayoutDestroy(&layout_row);
1675 if (!square_matrix) {
1676 CHKERR PetscLayoutDestroy(&layout_col);
1678 if (square_matrix) {
1679 nb_col_local_dofs = nb_row_local_dofs;
1680 nb_col_ghost_dofs = nb_row_ghost_dofs;
1693 MOFEM_LOG(
"WORLD", Sev::noisy) <<
"Partition problem " << name;
1695 using NumeredDofEntitysByIdx =
1700 auto &problems_set =
1701 const_cast<ProblemsByName &
>(problems_ptr->get<
Problem_mi_tag>());
1702 auto p_miit = problems_set.find(name);
1703 if (p_miit == problems_set.end()) {
1705 "problem with name %s not defined (top tip check spelling)",
1708 int nb_dofs_row = p_miit->getNbDofsRow();
1718 "entFEAdjacencies not build");
1724 ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1729 MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1732 MatPartitioning part;
1735 CHKERR MatPartitioningSetAdjacency(part, Adj);
1736 CHKERR MatPartitioningSetFromOptions(part);
1738 CHKERR MatPartitioningApply(part, &is);
1740 ISView(is, PETSC_VIEWER_STDOUT_WORLD);
1743 IS is_gather, is_num, is_gather_num;
1744 CHKERR ISAllGather(is, &is_gather);
1745 CHKERR ISPartitioningToNumbering(is, &is_num);
1746 CHKERR ISAllGather(is_num, &is_gather_num);
1747 const int *part_number, *petsc_idx;
1748 CHKERR ISGetIndices(is_gather, &part_number);
1749 CHKERR ISGetIndices(is_gather_num, &petsc_idx);
1750 int size_is_num, size_is_gather;
1751 CHKERR ISGetSize(is_gather, &size_is_gather);
1752 if (size_is_gather != (
int)nb_dofs_row)
1753 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1754 "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
1756 CHKERR ISGetSize(is_num, &size_is_num);
1757 if (size_is_num != (
int)nb_dofs_row)
1758 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1759 "data inconsistency %d != %d", size_is_num, nb_dofs_row);
1761 bool square_matrix =
false;
1762 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1763 square_matrix =
true;
1790 auto number_dofs = [&](
auto &dofs_idx,
auto &counter) {
1792 for (
auto miit_dofs_row = dofs_idx.begin();
1793 miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
1794 const int part = part_number[(*miit_dofs_row)->dofIdx];
1796 const bool success = dofs_idx.modify(
1799 part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
1802 "modification unsuccessful");
1805 const bool success = dofs_idx.modify(
1807 part, petsc_idx[(*miit_dofs_row)->dofIdx]));
1810 "modification unsuccessful");
1818 auto &dofs_row_by_idx_no_const =
const_cast<NumeredDofEntitysByIdx &
>(
1819 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>());
1820 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1821 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1822 nb_row_local_dofs = 0;
1823 nb_row_ghost_dofs = 0;
1825 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1827 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1828 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1829 if (square_matrix) {
1830 nb_col_local_dofs = nb_row_local_dofs;
1831 nb_col_ghost_dofs = nb_row_ghost_dofs;
1833 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1834 const_cast<NumeredDofEntitysByIdx &
>(
1835 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>());
1836 nb_col_local_dofs = 0;
1837 nb_col_ghost_dofs = 0;
1838 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1841 CHKERR ISRestoreIndices(is_gather, &part_number);
1842 CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
1843 CHKERR ISDestroy(&is_num);
1844 CHKERR ISDestroy(&is_gather_num);
1845 CHKERR ISDestroy(&is_gather);
1847 CHKERR MatPartitioningDestroy(&part);
1852 auto number_dofs = [&](
auto &dof_idx,
auto &counter) {
1854 for (
auto miit_dofs_row = dof_idx.begin(); miit_dofs_row != dof_idx.end();
1856 const bool success = dof_idx.modify(
1862 "modification unsuccessful");
1868 auto &dofs_row_by_idx_no_const =
const_cast<NumeredDofEntitysByIdx &
>(
1869 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>());
1870 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1871 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1872 nb_row_local_dofs = 0;
1873 nb_row_ghost_dofs = 0;
1875 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1877 bool square_matrix =
false;
1878 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1879 square_matrix =
true;
1881 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1882 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1883 if (square_matrix) {
1884 nb_col_local_dofs = nb_row_local_dofs;
1885 nb_col_ghost_dofs = nb_row_ghost_dofs;
1887 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1888 const_cast<NumeredDofEntitysByIdx &
>(
1889 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>());
1890 nb_col_local_dofs = 0;
1891 nb_col_ghost_dofs = 0;
1892 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1901 const std::string name,
const std::string problem_for_rows,
bool copy_rows,
1902 const std::string problem_for_cols,
bool copy_cols,
int verb) {
1913 ProblemByName &problems_by_name =
1914 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1915 ProblemByName::iterator p_miit = problems_by_name.find(name);
1916 if (p_miit == problems_by_name.end()) {
1918 "problem with name < %s > not defined (top tip check spelling)",
1923 << p_miit->getName() <<
" from rows of " << problem_for_rows
1924 <<
" and columns of " << problem_for_cols;
1927 ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
1928 if (p_miit_row == problems_by_name.end()) {
1930 "problem with name < %s > not defined (top tip check spelling)",
1931 problem_for_rows.c_str());
1933 const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
1934 p_miit_row->numeredRowDofsPtr;
1937 ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
1938 if (p_miit_col == problems_by_name.end()) {
1940 "problem with name < %s > not defined (top tip check spelling)",
1941 problem_for_cols.c_str());
1943 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
1944 p_miit_col->numeredColDofsPtr;
1946 bool copy[] = {copy_rows, copy_cols};
1947 boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
1948 p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
1950 int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
1951 int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
1952 boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
1955 for (
int ss = 0; ss < 2; ss++) {
1958 *nb_local_dofs[ss] = 0;
1962 std::vector<int> is_local, is_new;
1966 for (NumeredDofEntity_multiIndex::iterator dit =
1967 composed_dofs[ss]->begin();
1968 dit != composed_dofs[ss]->end(); dit++) {
1969 NumeredDofEntityByUId::iterator diit =
1970 dofs_by_uid.find((*dit)->getLocalUniqueId());
1971 if (diit == dofs_by_uid.end()) {
1974 "data inconsistency, could not find dof in composite problem");
1976 int part_number = (*diit)->getPart();
1977 int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
1979 success = composed_dofs[ss]->modify(
1984 "modification unsuccessful");
1986 if ((*dit)->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
1987 success = composed_dofs[ss]->modify(
1991 "modification unsuccessful");
1993 is_local.push_back(petsc_global_dof);
1998 CHKERR AOCreateMapping(m_field.
get_comm(), is_local.size(), &is_local[0],
2003 for (NumeredDofEntity_multiIndex::iterator dit =
2004 composed_dofs[ss]->begin();
2005 dit != composed_dofs[ss]->end(); dit++) {
2006 is_local.push_back((*dit)->getPetscGlobalDofIdx());
2008 CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2010 for (NumeredDofEntity_multiIndex::iterator dit =
2011 composed_dofs[ss]->begin();
2012 dit != composed_dofs[ss]->end(); dit++) {
2013 int part_number = (*dit)->getPart();
2014 int petsc_global_dof = is_local[idx2++];
2016 success = composed_dofs[ss]->modify(
2021 "modification unsuccessful");
2029 for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2030 dit != copied_dofs[ss]->end(); dit++) {
2031 std::pair<NumeredDofEntity_multiIndex::iterator, bool>
p;
2032 p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2037 int dof_idx = (*dit)->getDofIdx();
2038 int part_number = (*dit)->getPart();
2039 int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2041 const bool success = composed_dofs[ss]->modify(
2043 part_number, dof_idx, petsc_global_dof,
2044 (*nb_local_dofs[ss])++));
2047 "modification unsuccessful");
2050 const bool success = composed_dofs[ss]->modify(
2052 part_number, dof_idx, petsc_global_dof));
2055 "modification unsuccessful");
2076 << problem_ptr->
getName() <<
" Nb. local dof "
2093 const EntityHandle ent) {
2095 EntityHandle out_meshset;
2096 CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
2097 CHKERR moab.add_entities(out_meshset, &ent, 1);
2098 CHKERR moab.write_file(name.c_str(),
"VTK",
"", &out_meshset, 1);
2099 CHKERR moab.delete_entities(&out_meshset, 1);
2106 NumeredDofEntitysByIdx;
2107 NumeredDofEntitysByIdx::iterator dit, hi_dit;
2108 const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
2116 for (
int ss = 0; ss < 2; ss++) {
2118 dit = numered_dofs_ptr[ss]->begin();
2119 hi_dit = numered_dofs_ptr[ss]->end();
2120 for (; dit != hi_dit; dit++) {
2122 if ((*dit)->getPetscLocalDofIdx() < 0) {
2123 std::ostringstream zz;
2126 "local dof index for %d (0-row, 1-col) not set, i.e. has "
2127 "negative value\n %s",
2128 ss, zz.str().c_str());
2130 if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
2131 std::ostringstream zz;
2134 "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
2138 if ((*dit)->getPetscGlobalDofIdx() < 0) {
2140 const EntityHandle ent = (*dit)->getEnt();
2145 "_negative_global_index.vtk",
2148 std::ostringstream zz;
2150 << dit->get()->getBitRefLevel() <<
" " << **dit;
2152 "global dof index for %d (0-row, 1-col) row not set, i.e. "
2153 "has negative value\n %s",
2154 ss, zz.str().c_str());
2156 if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
2157 std::ostringstream zz;
2159 << *nbdof_ptr[ss] <<
" " << **dit;
2161 "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
2172 bool part_from_moab,
2174 int hi_proc,
int verb) {
2188 "adjacencies not build");
2195 "problem not partitioned");
2205 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
2206 ProblemByName::iterator p_miit = problems.find(name);
2207 if (p_miit == problems.end()) {
2209 "problem < %s > not found (top tip: check spelling)",
2215 *p_miit->numeredFiniteElementsPtr;
2218 problem_finite_elements.clear();
2222 bool do_cols_prob =
true;
2223 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2224 do_cols_prob =
false;
2227 auto get_good_elems = [&]() {
2228 auto good_elems = std::vector<
decltype(fe_ent_ptr->begin())>();
2229 good_elems.reserve(fe_ent_ptr->size());
2231 const auto prb_bit = p_miit->getBitRefLevel();
2232 const auto prb_mask = p_miit->getBitRefLevelMask();
2236 for (
auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2239 if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2241 const auto &fe_bit = (*efit)->getBitRefLevel();
2244 if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit).any())
2245 good_elems.emplace_back(efit);
2252 auto good_elems = get_good_elems();
2254 auto numbered_good_elems_ptr =
2255 boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2256 numbered_good_elems_ptr->reserve(good_elems.size());
2257 for (
auto &efit : good_elems)
2260 if (!do_cols_prob) {
2261 for (
auto &fe : *numbered_good_elems_ptr) {
2262 if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2263 fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2268 if (part_from_moab) {
2269 for (
auto &fe : *numbered_good_elems_ptr) {
2271 int proc = fe.getPartProc();
2272 if (proc == -1 && fe.getEntType() == MBVERTEX)
2273 proc = fe.getOwnerProc();
2278 for (
auto &fe : *numbered_good_elems_ptr) {
2281 CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2283 if (!part_from_moab) {
2285 for (
auto &dof_ptr : rows_view)
2286 parts[dof_ptr->pArt]++;
2287 std::vector<int>::iterator pos = max_element(parts.begin(), parts.end());
2288 const auto max_part = std::distance(parts.begin(), pos);
2293 for (
auto &fe : *numbered_good_elems_ptr) {
2295 auto check_fields_and_dofs = [&]() {
2296 if (!part_from_moab) {
2297 if (fe.getBitFieldIdRow().none() && m_field.
get_comm_size() == 0) {
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());
2308 if (check_fields_and_dofs()) {
2310 auto p = problem_finite_elements.insert(
2311 boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
2319 auto elements_on_rank =
2320 problem_finite_elements.get<
Part_mi_tag>().equal_range(
2323 << p_miit->getName() <<
" nb. elems "
2324 << std::distance(elements_on_rank.first, elements_on_rank.second);
2326 for (
auto &fe : *fe_ptr) {
2332 <<
"Element " << fe->getName() <<
" nb. elems "
2333 << std::distance(e_range.first, e_range.second);
2351 "partition of problem not build");
2354 "partitions finite elements not build");
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;
2376 p_miit->numeredFiniteElementsPtr->get<
Part_mi_tag>().equal_range(
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();
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) {
2396 fe_vec_view.clear();
2398 *(p_miit->getNumeredRowDofsPtr()),
2399 fe_vec_view, Inserter());
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);
2409 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2411 auto hint_c = ghost_idx_col_view.begin();
2412 for (
auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2414 fe_vec_view.clear();
2416 *(p_miit->getNumeredColDofsPtr()),
2417 fe_vec_view, Inserter());
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);
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};
2431 &ghost_idx_row_view, &ghost_idx_col_view};
2437 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
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))
2448 "inconsistent data, ghost dof already set");
2449 bool success = dof_by_uid_no_const[ss]->modify(
2451 if (PetscUnlikely(!success))
2453 "modification unsuccessful");
2454 (*nb_ghost_dofs[ss])++;
2457 if (loop_size == 1) {
2458 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
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();
2485 "partition of problem not build");
2488 "partitions finite elements not build");
2493 ProblemsByName &problems_set =
2494 const_cast<ProblemsByName &
>(problems_ptr->get<
Problem_mi_tag>());
2495 ProblemsByName::iterator p_miit = problems_set.find(name);
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;
2510 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2514 typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2515 NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2516 p_miit->numeredColDofsPtr};
2519 for (
int ss = 0; ss != loop_size; ++ss) {
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));
2529 auto cmp = [](
auto a,
auto b) {
2530 return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2532 sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2535 for (
auto gid_it : ghost_idx_view) {
2536 bool success = numered_dofs[ss]->modify(
2540 "modification unsuccessful");
2541 ++(*nb_ghost_dofs[ss]);
2544 if (loop_size == 1) {
2545 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
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();
2564 const std::string fe_name,
2565 EntityHandle *meshset)
const {
2574 .lower_bound(fe_name);
2577 .upper_bound(fe_name);
2578 std::vector<EntityHandle> fe_vec;
2579 fe_vec.reserve(std::distance(fit, hi_fe_it));
2580 for (; fit != hi_fe_it; fit++)
2581 fe_vec.push_back(fit->get()->getEnt());
2582 CHKERR m_field.
get_moab().add_entities(*meshset, &*fe_vec.begin(),
2589 const std::string fe_name,
2590 PetscLayout *layout)
const {
2601 const std::string problem_name,
const std::string
field_name,
2602 const Range ents,
const int lo_coeff,
const int hi_coeff,
2603 const int lo_order,
const int hi_order,
int verb,
const bool debug) {
2627 for (
int s = 0; s != 2; ++s)
2628 if (numered_dofs[s]) {
2630 typedef multi_index_container<
2632 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2635 NumeredDofEntity_it_view_multiIndex;
2638 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2641 for (
auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2643 auto lo = numered_dofs[s]->get<
Unique_mi_tag>().lower_bound(
2645 auto hi = numered_dofs[s]->get<
Unique_mi_tag>().upper_bound(
2648 for (; lo != hi; ++lo)
2649 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2650 (*lo)->getDofCoeffIdx() <= hi_coeff &&
2651 (*lo)->getDofOrder() >= lo_order &&
2652 (*lo)->getDofOrder() <= hi_order)
2653 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2657 for (
auto &dof : dofs_it_view)
2665 for (
auto dit : dofs_it_view) {
2666 bool success = numered_dofs[s]->modify(dit, mod);
2669 "can not set negative indices");
2673 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
2674 dofs_weak_view.reserve(dofs_it_view.size());
2675 for (
auto dit : dofs_it_view)
2676 dofs_weak_view.push_back(*dit);
2680 "Number of DOFs in multi-index %d and to delete %d\n",
2681 numered_dofs[s]->size(), dofs_it_view.size());
2684 for (
auto weak_dit : dofs_weak_view)
2685 if (
auto dit = weak_dit.lock()) {
2686 numered_dofs[s]->erase(dit->getLocalUniqueId());
2691 "Number of DOFs in multi-index after delete %d\n",
2692 numered_dofs[s]->size());
2695 int nb_local_dofs = 0;
2696 int nb_ghost_dofs = 0;
2697 for (
auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2699 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2700 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2702 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2706 "Impossible case. You could run problem on no distributed "
2707 "mesh. That is not implemented. Dof local index is %d",
2708 (*dit)->getPetscLocalDofIdx());
2712 auto get_indices_by_tag = [&](
auto tag,
auto &indices,
bool only_local) {
2713 const int nb_dofs = numered_dofs[s]->size();
2715 indices.reserve(nb_dofs);
2716 for (
auto dit = numered_dofs[s]->get<
decltype(tag)>().begin();
2717 dit != numered_dofs[s]->get<
decltype(tag)>().end(); ++dit) {
2720 if ((*dit)->getPetscLocalDofIdx() < 0 ||
2721 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2726 indices.push_back(
decltype(tag)::get_index(dit));
2730 auto get_indices_by_uid = [&](
auto tag,
auto &indices) {
2731 const int nb_dofs = numered_dofs[s]->size();
2733 indices.reserve(nb_dofs);
2734 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2736 indices.push_back(
decltype(tag)::get_index(dit));
2739 auto concatenate_dofs = [&](
auto tag,
auto &indices,
2740 const auto local_only) {
2742 get_indices_by_tag(tag, indices, local_only);
2747 &*indices.begin(), PETSC_NULL, &ao);
2749 CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
2750 &*indices.begin(), PETSC_NULL, &ao);
2752 get_indices_by_uid(tag, indices);
2753 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
2759 auto set_concatinated_indices = [&]() {
2760 std::vector<int> global_indices;
2761 std::vector<int> local_indices;
2765 auto gi = global_indices.begin();
2766 auto li = local_indices.begin();
2767 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2770 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
2771 bool success = numered_dofs[s]->modify(dit, mod);
2774 "can not set negative indices");
2780 CHKERR set_concatinated_indices();
2782 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2784 *(local_nbdof_ptr[s]) = nb_local_dofs;
2785 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
2788 for (
auto dof : (*numered_dofs[s])) {
2789 if (dof->getPetscGlobalDofIdx() < 0) {
2791 "Negative global idx");
2793 if (dof->getPetscLocalDofIdx() < 0) {
2795 "Negative local idx");
2801 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2802 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2803 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
2808 "removed ents on rank %d from problem %s dofs [ %d / %d "
2810 "%d) local, %d / %d (before %d / %d) "
2811 "ghost, %d / %d (before %d / %d) global]",
2814 nb_init_row_dofs, nb_init_col_dofs,
2816 nb_init_ghost_row_dofs, nb_init_ghost_col_dofs,
2818 nb_init_loc_row_dofs, nb_init_loc_col_dofs);
2826 const std::string problem_name,
const std::string
field_name,
2827 const Range ents,
const int lo_coeff,
const int hi_coeff,
2828 const int lo_order,
const int hi_order,
int verb,
const bool debug) {
2852 for (
int s = 0; s != 2; ++s)
2853 if (numered_dofs[s]) {
2855 typedef multi_index_container<
2857 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2860 NumeredDofEntity_it_view_multiIndex;
2863 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2866 for (
auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2868 auto lo = numered_dofs[s]->get<
Unique_mi_tag>().lower_bound(
2870 auto hi = numered_dofs[s]->get<
Unique_mi_tag>().upper_bound(
2873 for (; lo != hi; ++lo)
2874 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2875 (*lo)->getDofCoeffIdx() <= hi_coeff &&
2876 (*lo)->getDofOrder() >= lo_order &&
2877 (*lo)->getDofOrder() <= hi_order)
2878 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2882 for (
auto &dof : dofs_it_view)
2888 for (
auto dit : dofs_it_view) {
2889 bool success = numered_dofs[s]->modify(dit, mod);
2892 "can not set negative indices");
2896 std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2897 dosf_weak_view.reserve(dofs_it_view.size());
2898 for (
auto dit : dofs_it_view)
2899 dosf_weak_view.push_back(*dit);
2903 "Number of DOFs in multi-index %d and to delete %d\n",
2904 numered_dofs[s]->size(), dofs_it_view.size());
2907 for (
auto weak_dit : dosf_weak_view)
2908 if (
auto dit = weak_dit.lock()) {
2909 numered_dofs[s]->erase(dit->getLocalUniqueId());
2914 "Number of DOFs in multi-index after delete %d\n",
2915 numered_dofs[s]->size());
2918 int nb_global_dof = 0;
2919 int nb_local_dofs = 0;
2920 int nb_ghost_dofs = 0;
2922 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2925 if ((*dit)->getDofIdx() >= 0) {
2927 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2928 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2930 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2938 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2940 if (*(nbdof_ptr[s]) != nb_global_dof)
2942 "Number of local DOFs do not add up %d != %d",
2943 *(nbdof_ptr[s]), nb_global_dof);
2946 *(nbdof_ptr[s]) = nb_global_dof;
2947 *(local_nbdof_ptr[s]) = nb_local_dofs;
2948 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
2951 auto get_indices_by_tag = [&](
auto tag) {
2952 std::vector<int> indices;
2953 const int nb_dofs = numered_dofs[s]->size();
2954 indices.resize(nb_dofs, -1);
2956 for (
auto dit = numered_dofs[s]->get<
decltype(tag)>().begin();
2957 dit != numered_dofs[s]->get<
decltype(tag)>().end(); ++dit) {
2958 int idx = std::distance(numered_dofs[s]->begin(),
2959 numered_dofs[s]->project<0>(dit));
2960 int current_idx =
decltype(tag)::get_index(dit);
2961 if (current_idx >= 0)
2971 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2975 (*dit)->getPart(), (*dit)->getDofIdx(), global_indices[idx],
2976 local_indices[idx]);
2977 bool success = numered_dofs[s]->modify(dit, mod);
2980 "can not set negative indices");
2984 for (
auto dof : (*numered_dofs[s])) {
2985 if (dof->getDofIdx() >= 0 && dof->getPetscGlobalDofIdx() < 0) {
2987 "Negative global idx");
2992 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2993 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2994 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3001 auto log = [&](
auto str,
auto level) {
3003 "removed ents on rank %d from problem %s dofs [ %d / %d "
3005 "%d) local, %d / %d (before %d / %d) "
3006 "ghost, %d / %d (before %d / %d) global]",
3009 nb_init_row_dofs, nb_init_col_dofs,
3011 nb_init_ghost_row_dofs, nb_init_ghost_col_dofs,
3013 nb_init_loc_row_dofs, nb_init_loc_col_dofs);
3016 log(
"WORLD", Sev::verbose);
3018 log(
"SYNC", Sev::noisy);
3026 const std::string problem_name,
const std::string
field_name,
3028 Range *ents_ptr,
const int lo_coeff,
const int hi_coeff,
const int lo_order,
3029 const int hi_order,
int verb,
const bool debug) {
3038 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3041 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3046 hi_coeff, lo_order, hi_order, verb,
debug);
3052 const std::string problem_name,
const std::string
field_name,
3054 Range *ents_ptr,
const int lo_coeff,
const int hi_coeff,
const int lo_order,
3055 const int hi_order,
int verb,
const bool debug) {
3064 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3067 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3072 lo_coeff, hi_coeff, lo_order,
3073 hi_order, verb,
debug);
3080 const enum MarkOP op,
const Range ents,
3081 std::vector<unsigned char> &
marker)
const {
3087 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3097 marker.resize(dofs->size(), 0);
3098 std::vector<unsigned char> marker_tmp;
3102 for (
auto p = ents.pair_begin();
p != ents.pair_end(); ++
p) {
3103 auto lo = dofs->get<
Ent_mi_tag>().lower_bound(
p->first);
3104 auto hi = dofs->get<
Ent_mi_tag>().upper_bound(
p->second);
3105 for (; lo != hi; ++lo)
3106 marker[(*lo)->getPetscLocalDofIdx()] |= 1;
3110 marker_tmp.resize(dofs->size(), 0);
3111 for (
auto p = ents.pair_begin();
p != ents.pair_end(); ++
p) {
3112 auto lo = dofs->get<
Ent_mi_tag>().lower_bound(
p->first);
3113 auto hi = dofs->get<
Ent_mi_tag>().upper_bound(
p->second);
3114 for (; lo != hi; ++lo)
3115 marker_tmp[(*lo)->getPetscLocalDofIdx()] = 1;
3117 for (
int i = 0;
i !=
marker.size(); ++
i) {
3129 const unsigned char c, std::vector<unsigned char> &
marker)
const {
3135 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3145 marker.resize(dofs->size(), 0);
3152 auto marker_ref = [
marker](
auto &it) ->
unsigned int & {
3153 return marker[(*it)->getPetscLocalDofIdx()];
3158 for (; dof_lo != dof_hi; ++dof_lo)
3159 if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3160 (*dof_lo)->getDofCoeffIdx() <= hi)
3161 marker[(*dof_lo)->getPetscLocalDofIdx()] |=
c;
3164 for (; dof_lo != dof_hi; ++dof_lo)
3165 if ((*dof_lo)->getDofCoeffIdx() >= lo &&
3166 (*dof_lo)->getDofCoeffIdx() <= hi)
3167 marker[(*dof_lo)->getPetscLocalDofIdx()] &=
c;
3176 const std::string row_field,
3177 const std::string col_field)
const {
3182 const auto problem_ptr = m_field.
get_problem(problem_name);
3183 auto get_field_id = [&](
const std::string
field_name) {
3186 const auto row_id = get_field_id(row_field);
3187 const auto col_id = get_field_id(col_field);
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define ProblemManagerFunctionBegin
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#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 ...
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FieldEntity, UId, &FieldEntity::getGlobalUniqueId > > > > FieldEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getLocalUniqueId > >, ordered_non_unique< const_mem_fun< DofEntity, char, &DofEntity::getActive > > > > DofEntity_multiIndex_active_view
multi-index view on DofEntity activity
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getGlobalUniqueId > > > > DofEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
MoFEMErrorCode getProblemElementsLayout(const std::string name, const std::string fe_name, PetscLayout *layout) const
Get layout of elements in the problem.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildComposedProblem(const std::string out_name, const std::vector< std::string > add_row_problems, const std::vector< std::string > add_col_problems, const bool square_matrix=true, int verb=1)
build composite problem
MoFEMErrorCode getFEMeshset(const std::string prb_name, const std::string fe_name, EntityHandle *meshset) const
create add entities of finite element in the problem
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, std::pair< EntityType, EntityType > > *entityMapRow=nullptr, const map< std::string, std::pair< EntityType, EntityType > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
DEPRECATED MoFEMErrorCode partitionMesh(const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
Set partition tag to each finite element in the problem.
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionGhostDofsOnDistributedMesh(const std::string name, int verb=VERBOSE)
determine ghost nodes on distributed meshes
MoFEMErrorCode inheritPartition(const std::string name, const std::string problem_for_rows, bool copy_rows, const std::string problem_for_cols, bool copy_cols, int verb=VERBOSE)
build indexing and partition problem inheriting indexing and partitioning from two other problems
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
auto marker
set bit to marker
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
Problem::EmptyFieldBlocks EmptyFieldBlocks
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary 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
const double r
rate factor
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
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static MoFEMErrorCode getDofView(const FE_ENTS &fe_ents_view, const MOFEM_DOFS &mofem_dofs, MOFEM_DOFS_VIEW &dofs_view, INSERTER &&inserter)
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static auto getHandleFromUniqueId(const UId uid)
static auto getFieldBitNumberFromUniqueId(const UId uid)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
static auto getOwnerFromUniqueId(const UId uid)
UId getGlobalUniqueIdCalculate() const
Calculate global UId.
const BitFieldId & getId() const
Get unique field id.
IdxDataType(const UId uid, const int dof)
IdxDataTypePtr(const int *ptr)
Matrix manager is used to build and partition problems.
keeps information about indexed dofs for the problem
Partitioned (Indexed) Finite Element in Problem.
keeps basic data about problem
DofIdx getNbLocalDofsRow() const
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
DofIdx getNbDofsRow() const
BitRefLevel getBitRefLevel() const
auto & getColDofsSequence() const
Get reference to sequence data numbered dof container.
BitRefLevel getBitRefLevelMask() const
DofIdx nbGhostDofsCol
Number of ghost DOFs in col.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
DofIdx getNbGhostDofsCol() const
DofIdx nbLocDofsRow
Local number of DOFs in row.
BitFEId getBitFEId() const
DofIdx nbDofsCol
Global number of DOFs in col.
DofIdx getNbDofsCol() const
DofIdx getNbLocalDofsCol() const
DofIdx nbGhostDofsRow
Number of ghost DOFs in row.
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
DofIdx nbLocDofsCol
Local number of DOFs in colIs.
DofIdx nbDofsRow
Global number of DOFs in row.
DofIdx getNbGhostDofsRow() const
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
auto & getRowDofsSequence() const
Get reference to sequence data numbered dof container.
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
ProblemsManager(const MoFEM::Core &core)
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string problem_name, const std::string row_field, const std::string col_field) const
add empty block to problem
MoFEMErrorCode modifyMarkDofs(const std::string problem_name, RowColData rc, const std::string field_name, const int lo, const int hi, const enum MarkOP op, const unsigned char c, std::vector< unsigned char > &marker) const
Mark DOFs.
PetscLogEvent MOFEM_EVENT_ProblemsManager
PetscBool buildProblemFromFields
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode removeDofsOnEntitiesNotDistributed(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
MoFEMErrorCode getOptions()
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.