24 MoFEMErrorCode MatPartitioningApply_Parmetis_MoFEM(MatPartitioning part,
29 #define ProblemManagerFunctionBegin \
31 MOFEM_LOG_CHANNEL("WORLD"); \
32 MOFEM_LOG_CHANNEL("SYNC"); \
33 MOFEM_LOG_FUNCTION(); \
34 MOFEM_LOG_TAG("SYNC", "ProblemsManager"); \
35 MOFEM_LOG_TAG("WORLD", "ProblemsManager")
39 bcopy(&uid,
dAta, 4 *
sizeof(
int));
50 int global_dof =
pTr[4];
54 unsigned int b0, b1, b2, b3;
55 bcopy(&
pTr[0], &b0,
sizeof(
int));
56 bcopy(&
pTr[1], &b1,
sizeof(
int));
57 bcopy(&
pTr[2], &b2,
sizeof(
int));
58 bcopy(&
pTr[3], &b3,
sizeof(
int));
59 UId uid =
static_cast<UId>(b0) |
static_cast<UId>(b1) << 8 *
sizeof(
int) |
60 static_cast<UId>(b2) << 16 *
sizeof(
int) |
61 static_cast<UId>(b3) << 24 *
sizeof(
int);
84 buildProblemFromFields(PETSC_FALSE),
85 synchroniseProblemEntities(PETSC_FALSE) {
92 CHKERR PetscOptionsBegin(m_field.
get_comm(),
"",
"Problem manager",
"none");
95 "-problem_build_from_fields",
96 "Add DOFs to problem directly from fields not through DOFs on elements",
99 ierr = PetscOptionsEnd();
105 const Range &ents,
const int dim,
const int adj_dim,
const int n_parts,
106 Tag *th_vertex_weights, Tag *th_edge_weights, Tag *th_part_weights,
107 int verb,
const bool debug) {
112 int rstart, rend, nb_elems;
116 CHKERR PetscLayoutSetBlockSize(layout, 1);
117 CHKERR PetscLayoutSetSize(layout, ents.size());
118 CHKERR PetscLayoutSetUp(layout);
119 CHKERR PetscLayoutGetSize(layout, &nb_elems);
120 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
121 CHKERR PetscLayoutDestroy(&layout);
124 <<
"Finite elements in problem: row lower " << rstart <<
" row upper "
125 << rend <<
" nb. elems " << nb_elems <<
" ( " << ents.size() <<
" )";
130 std::vector<EntityHandle> weight_ents;
131 weight_ents.reserve(rend - rstart + 1);
135 std::vector<int> adj;
136 AdjBridge(
const EntityHandle ent, std::vector<int> &adj)
137 : ent(ent), adj(adj) {}
140 typedef multi_index_container<
144 hashed_unique<member<AdjBridge, EntityHandle, &AdjBridge::ent>>
150 CHKERR m_field.
get_moab().get_adjacencies(ents, adj_dim,
true, all_dim_ents,
151 moab::Interface::UNION);
153 AdjBridgeMap adj_bridge_map;
154 auto hint = adj_bridge_map.begin();
155 std::vector<int> adj;
156 for (
auto ent : all_dim_ents) {
159 adj_ents = intersect(adj_ents, ents);
161 adj.reserve(adj_ents.size());
162 for (
auto a : adj_ents)
163 adj.emplace_back(ents.index(a));
164 hint = adj_bridge_map.emplace_hint(hint, ent, adj);
170 const int nb_loc_elements = rend - rstart;
171 std::vector<int>
i(nb_loc_elements + 1, 0),
j;
173 std::vector<int> row_adj;
174 Range::iterator fe_it;
179 fe_it = ents.begin(), ii = 0, jj = 0, max_row_size = 0;
181 fe_it != ents.end(); ++fe_it, ++ii) {
188 if (m_field.
get_moab().type_from_handle(*fe_it) == MBENTITYSET) {
191 "not yet implemented, don't know what to do for meshset element");
195 CHKERR m_field.
get_moab().get_adjacencies(&*fe_it, 1, adj_dim,
false,
197 dim_ents = intersect(dim_ents, all_dim_ents);
200 for (
auto e : dim_ents) {
201 auto adj_it = adj_bridge_map.find(e);
202 if (adj_it != adj_bridge_map.end()) {
204 for (
const auto idx : adj_it->adj)
205 row_adj.push_back(idx);
212 std::sort(row_adj.begin(), row_adj.end());
213 auto end = std::unique(row_adj.begin(), row_adj.end());
215 size_t row_size = std::distance(row_adj.begin(), end);
216 max_row_size = std::max(max_row_size, row_size);
217 if (
j.capacity() < (nb_loc_elements - jj) * max_row_size)
218 j.reserve(nb_loc_elements * max_row_size);
221 auto diag = ents.index(*fe_it);
222 for (
auto it = row_adj.begin(); it != end; ++it)
229 if (th_vertex_weights != NULL)
230 weight_ents.push_back(*fe_it);
236 CHKERR PetscMalloc(
i.size() *
sizeof(
int), &_i);
237 CHKERR PetscMalloc(
j.size() *
sizeof(
int), &_j);
238 copy(
i.begin(),
i.end(), _i);
239 copy(
j.begin(),
j.end(), _j);
243 int *vertex_weights = NULL;
244 if (th_vertex_weights != NULL) {
245 CHKERR PetscMalloc(weight_ents.size() *
sizeof(
int), &vertex_weights);
247 &*weight_ents.begin(),
248 weight_ents.size(), vertex_weights);
254 CHKERR MatCreateMPIAdj(m_field.
get_comm(), rend - rstart, nb_elems, _i, _j,
256 CHKERR MatSetOption(Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
260 CHKERR MatConvert(Adj, MATMPIAIJ, MAT_INITIAL_MATRIX, &A);
261 CHKERR MatView(A, PETSC_VIEWER_DRAW_WORLD);
268 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Start";
270 MatPartitioning part;
274 CHKERR MatPartitioningSetAdjacency(part, Adj);
275 CHKERR MatPartitioningSetFromOptions(part);
276 CHKERR MatPartitioningSetNParts(part, n_parts);
277 if (th_vertex_weights != NULL) {
278 CHKERR MatPartitioningSetVertexWeights(part, vertex_weights);
281 PetscObjectTypeCompare((PetscObject)part, MATPARTITIONINGPARMETIS, &same);
284 CHKERR MatPartitioningApply_Parmetis_MoFEM(part, &is);
287 CHKERR MatPartitioningApply(part, &is);
290 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"End";
293 IS is_gather, is_num, is_gather_num;
294 CHKERR ISAllGather(is, &is_gather);
295 CHKERR ISPartitioningToNumbering(is, &is_num);
296 CHKERR ISAllGather(is_num, &is_gather_num);
298 const int *part_number, *gids;
299 CHKERR ISGetIndices(is_gather, &part_number);
300 CHKERR ISGetIndices(is_gather_num, &gids);
303 ParallelComm *pcomm = ParallelComm::get_pcomm(
305 Tag part_tag = pcomm->part_tag();
306 CHKERR m_field.
get_moab().tag_set_data(part_tag, ents, part_number);
307 auto get_gid_tag = [&]() {
309 rval = m_field.
get_moab().tag_get_handle(GLOBAL_ID_TAG_NAME, gid_tag);
310 if (
rval != MB_SUCCESS) {
313 GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag,
314 MB_TAG_DENSE | MB_TAG_CREAT, &zero);
318 Tag gid_tag = get_gid_tag();
320 std::map<int, Range> parts_ents;
323 Range::iterator eit = ents.begin();
324 for (
int ii = 0; eit != ents.end(); eit++, ii++) {
325 parts_ents[part_number[ii]].insert(*eit);
329 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
330 moab::Interface::UNION);
331 if (!tagged_sets.empty())
334 if (n_parts > (
int)tagged_sets.size()) {
336 int num_new = n_parts - tagged_sets.size();
337 for (
int i = 0;
i < num_new;
i++) {
338 EntityHandle new_set;
340 tagged_sets.insert(new_set);
342 }
else if (n_parts < (
int)tagged_sets.size()) {
344 int num_del = tagged_sets.size() - n_parts;
345 for (
int i = 0;
i < num_del;
i++) {
346 EntityHandle old_set = tagged_sets.pop_back();
352 std::vector<int> dum_ids(n_parts);
353 for (
int i = 0;
i < n_parts;
i++)
360 for (
int pp = 0; pp != n_parts; pp++) {
361 Range dim_ents = parts_ents[pp].subset_by_dimension(
dim);
366 dim_ents,
dd,
true, adj_ents, moab::Interface::UNION);
371 parts_ents[pp].merge(adj_ents);
374 for (
int pp = 1; pp != n_parts; pp++) {
375 for (
int ppp = 0; ppp != pp; ppp++) {
376 parts_ents[pp] = subtract(parts_ents[pp], parts_ents[ppp]);
382 ss <<
"out_part_" << rr <<
".vtk";
383 EntityHandle meshset;
391 for (
int pp = 0; pp != n_parts; pp++) {
392 CHKERR m_field.
get_moab().add_entities(tagged_sets[pp], parts_ents[pp]);
398 for (
int pp = 0; pp != n_parts; pp++) {
399 Range dim_ents = parts_ents[pp].subset_by_dimension(
dd);
401 for (Range::iterator eit = dim_ents.begin(); eit != dim_ents.end();
413 CHKERR ISRestoreIndices(is_gather, &part_number);
414 CHKERR ISRestoreIndices(is_gather_num, &gids);
415 CHKERR ISDestroy(&is_num);
416 CHKERR ISDestroy(&is_gather_num);
417 CHKERR ISDestroy(&is_gather);
419 CHKERR MatPartitioningDestroy(&part);
429 const bool square_matrix,
449 const bool square_matrix,
452 auto dofs_field_ptr = m_field.
get_dofs();
462 SETERRQ1(PETSC_COMM_SELF, 1,
"problem <%s> refinement level not set",
463 problem_ptr->
getName().c_str());
471 auto make_rows_and_cols_view = [&](
auto &dofs_rows,
auto &dofs_cols) {
473 for (
auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
475 const auto uid = (*it)->getLocalUniqueId();
478 for (
auto lo =
r.first; lo !=
r.second; ++lo) {
480 if ((lo->getBitFEId() & problem_ptr->
getBitFEId()).any()) {
481 std::array<bool, 2> row_col = {
false,
false};
485 const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
488 if ((fe_bit & prb_mask) != fe_bit)
490 if ((fe_bit & prb_bit) != prb_bit)
493 auto add_to_view = [&](
auto &nb_dofs,
auto &view,
auto rc) {
494 auto dit = nb_dofs->lower_bound(uid);
495 decltype(dit) hi_dit;
496 if (dit != nb_dofs->end()) {
497 hi_dit = nb_dofs->upper_bound(
499 view.insert(dit, hi_dit);
504 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
505 add_to_view(dofs_field_ptr, dofs_rows,
ROW);
508 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
509 add_to_view(dofs_field_ptr, dofs_cols,
COL);
511 if (square_matrix && row_col[
ROW])
513 else if (row_col[
ROW] && row_col[
COL])
521 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
531 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
533 hi_miit = dofs_rows.get<0>().end();
535 int count_dofs = dofs_rows.get<1>().count(
true);
536 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
537 boost::shared_ptr<std::vector<NumeredDofEntity>>(
538 new std::vector<NumeredDofEntity>());
540 dofs_array->reserve(count_dofs);
541 miit = dofs_rows.get<0>().begin();
542 for (; miit != hi_miit; miit++) {
543 if ((*miit)->getActive()) {
544 dofs_array->emplace_back(*miit);
545 dofs_array->back().dofIdx = (problem_ptr->
nbDofsRow)++;
549 for (
auto &v : *dofs_array) {
555 if (!square_matrix) {
562 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
564 hi_miit = dofs_cols.get<0>().end();
567 miit = dofs_cols.get<0>().begin();
568 for (; miit != hi_miit; miit++) {
569 if (!(*miit)->getActive()) {
575 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
576 boost::shared_ptr<std::vector<NumeredDofEntity>>(
577 new std::vector<NumeredDofEntity>());
579 dofs_array->reserve(count_dofs);
580 miit = dofs_cols.get<0>().begin();
581 for (; miit != hi_miit; miit++) {
582 if (!(*miit)->getActive()) {
585 dofs_array->emplace_back(*miit);
586 dofs_array->back().dofIdx = problem_ptr->
nbDofsCol++;
589 for (
auto &v : *dofs_array) {
601 << problem_ptr->
getName() <<
" Nb. local dofs "
609 <<
"FEs row dofs " << *problem_ptr <<
" Nb. row dof "
615 <<
"FEs col dofs " << *problem_ptr <<
" Nb. col dof "
632 const std::string name,
const bool square_matrix,
int verb) {
646 square_matrix, verb);
655 Problem *problem_ptr,
const bool square_matrix,
int verb) {
660 auto dofs_field_ptr = m_field.
get_dofs();
671 "problem <%s> refinement level not set",
672 problem_ptr->
getName().c_str());
680 boost::shared_ptr<NumeredDofEntity_multiIndex>(
694 auto make_rows_and_cols_view = [&](
auto &dofs_rows,
auto &dofs_cols) {
698 for (
auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
701 const auto uid = (*it)->getLocalUniqueId();
704 for (
auto lo =
r.first; lo !=
r.second; ++lo) {
706 if ((lo->getBitFEId() & problem_ptr->
getBitFEId()).any()) {
707 std::array<bool, 2> row_col = {
false,
false};
711 const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
714 if ((fe_bit & prb_mask) != fe_bit)
716 if ((fe_bit & prb_bit) != prb_bit)
719 auto add_to_view = [&](
auto &nb_dofs,
auto &view,
auto rc) {
720 auto dit = nb_dofs->lower_bound(uid);
721 decltype(dit) hi_dit;
722 if (dit != nb_dofs->end()) {
723 hi_dit = nb_dofs->upper_bound(
725 view.insert(dit, hi_dit);
730 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
731 add_to_view(dofs_field_ptr, dofs_rows,
ROW);
734 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
735 add_to_view(dofs_field_ptr, dofs_cols,
COL);
737 if (square_matrix && row_col[
ROW])
739 else if (row_col[
ROW] && row_col[
COL])
747 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
755 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
756 fit != fe_ptr->end(); fit++) {
757 if ((fit->get()->getId() & problem_ptr->
getBitFEId()).any()) {
758 fields_ids_row |= fit->get()->getBitFieldIdRow();
759 fields_ids_col |= fit->get()->getBitFieldIdCol();
763 for (Field_multiIndex::iterator fit = fields_ptr->begin();
764 fit != fields_ptr->end(); fit++) {
765 if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
769 auto hi_dit = dofs_field_ptr->get<
Unique_mi_tag>().upper_bound(
772 for (; dit != hi_dit; dit++) {
774 const int owner_proc = dit->get()->getOwnerProc();
776 const unsigned char pstatus = dit->get()->getPStatus();
781 const BitRefLevel dof_bit = (*dit)->getBitRefLevel();
783 if ((dof_bit & prb_mask) != dof_bit)
785 if ((dof_bit & prb_bit) != prb_bit)
787 if ((fit->get()->getId() & fields_ids_row).any()) {
788 dofs_rows.insert(*dit);
790 if (!square_matrix) {
791 if ((fit->get()->getId() & fields_ids_col).any()) {
792 dofs_cols.insert(*dit);
803 BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
804 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
805 fit != fe_ptr->end(); fit++) {
806 if ((fit->get()->getId() & problem_ptr->
getBitFEId()).any()) {
807 fields_ids_row |= fit->get()->getBitFieldIdRow();
808 fields_ids_col |= fit->get()->getBitFieldIdCol();
813 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
814 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
816 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
817 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
818 Range ents_to_synchronise;
819 for (; miit != hi_miit; ++miit) {
820 if (miit->get()->getEntDofIdx() != 0)
822 ents_to_synchronise.insert(miit->get()->getEnt());
824 Range tmp_ents = ents_to_synchronise;
826 ents_to_synchronise, verb);
827 ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
828 for (
auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
829 if ((fit->get()->getId() & *fields_ids[ss]).any()) {
830 const auto bit_number = (*fit)->getBitNumber();
831 for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
832 pit != ents_to_synchronise.pair_end(); ++pit) {
833 const auto f = pit->first;
834 const auto s = pit->second;
840 auto dit = dofs_field_ptr->get<
Unique_mi_tag>().lower_bound(lo_uid);
844 dofs_ptr[ss]->insert(dit, hi_dit);
853 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
860 for (
int ss = 0; ss < 2; ss++) {
861 *(nbdof_ptr[ss]) = 0;
862 *(local_nbdof_ptr[ss]) = 0;
863 *(ghost_nbdof_ptr[ss]) = 0;
868 int nb_local_dofs[] = {0, 0};
869 for (
int ss = 0; ss < loop_size; ss++) {
870 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
872 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
873 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
874 for (; miit != hi_miit; miit++) {
875 int owner_proc = (*miit)->getOwnerProc();
883 int start_ranges[2], end_ranges[2];
884 for (
int ss = 0; ss != loop_size; ss++) {
887 CHKERR PetscLayoutSetBlockSize(layout, 1);
888 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
889 CHKERR PetscLayoutSetUp(layout);
890 CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
891 CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
892 CHKERR PetscLayoutDestroy(&layout);
895 nbdof_ptr[1] = nbdof_ptr[0];
896 nb_local_dofs[1] = nb_local_dofs[0];
897 start_ranges[1] = start_ranges[0];
898 end_ranges[1] = end_ranges[0];
910 const size_t idx_data_type_size =
sizeof(
IdxDataType);
911 const size_t data_block_size = idx_data_type_size /
sizeof(
int);
916 std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
922 for (
int ss = 0; ss != loop_size; ++ss) {
924 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
926 hi_miit = dofs_ptr[ss]->get<0>().end();
928 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
929 boost::shared_ptr<std::vector<NumeredDofEntity>>(
930 new std::vector<NumeredDofEntity>());
931 int nb_dofs_to_add = 0;
932 miit = dofs_ptr[ss]->get<0>().begin();
933 for (; miit != hi_miit; ++miit) {
935 if (!(miit->get()->getActive()))
939 dofs_array->reserve(nb_dofs_to_add);
946 int &local_idx = *local_nbdof_ptr[ss];
947 miit = dofs_ptr[ss]->get<0>().begin();
948 for (; miit != hi_miit; ++miit) {
951 if (!(miit->get()->getActive()))
954 dofs_array->emplace_back(*miit);
956 int owner_proc = dofs_array->back().getOwnerProc();
957 if (owner_proc < 0) {
959 "data inconsistency");
963 dofs_array->back().pArt = owner_proc;
964 dofs_array->back().dofIdx = -1;
965 dofs_array->back().petscGloablDofIdx = -1;
966 dofs_array->back().petscLocalDofIdx = -1;
970 int glob_idx = start_ranges[ss] + local_idx;
971 dofs_array->back().pArt = owner_proc;
972 dofs_array->back().dofIdx = glob_idx;
973 dofs_array->back().petscGloablDofIdx = glob_idx;
974 dofs_array->back().petscLocalDofIdx = local_idx;
977 unsigned char pstatus = dofs_array->back().getPStatus();
983 proc < MAX_SHARING_PROCS &&
984 -1 != dofs_array->back().getSharingProcsPtr()[proc];
988 ids_data_packed_rows[dofs_array->back()
989 .getSharingProcsPtr()[proc]]
990 .emplace_back(dofs_array->back().getGlobalUniqueId(),
993 ids_data_packed_cols[dofs_array->back()
994 .getSharingProcsPtr()[proc]]
995 .emplace_back(dofs_array->back().getGlobalUniqueId(),
998 if (!(pstatus & PSTATUS_MULTISHARED)) {
1006 auto hint = numered_dofs_ptr[ss]->end();
1007 for (
auto &v : *dofs_array)
1008 hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
1010 if (square_matrix) {
1011 local_nbdof_ptr[1] = local_nbdof_ptr[0];
1014 int nsends_rows = 0, nsends_cols = 0;
1018 lengths_rows.clear();
1019 lengths_cols.clear();
1020 for (
int proc = 0; proc < m_field.
get_comm_size(); proc++) {
1021 lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
1022 lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
1023 if (!ids_data_packed_rows[proc].empty())
1025 if (!ids_data_packed_cols[proc].empty())
1047 &lengths_rows[0], &nrecvs_rows);
1052 CHKERR PetscGatherMessageLengths(m_field.
get_comm(), nsends_rows,
1053 nrecvs_rows, &lengths_rows[0],
1054 &onodes_rows, &olengths_rows);
1065 MPI_Request *r_waits_row;
1070 CHKERR PetscPostIrecvInt(m_field.
get_comm(), tag_row, nrecvs_rows,
1071 onodes_rows, olengths_rows, &rbuf_row,
1073 CHKERR PetscFree(onodes_rows);
1075 MPI_Request *s_waits_row;
1076 CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
1079 for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
1080 if (!lengths_rows[proc])
1082 CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0],
1085 tag_row, m_field.
get_comm(), s_waits_row + kk);
1090 CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
1093 CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
1096 CHKERR PetscFree(r_waits_row);
1097 CHKERR PetscFree(s_waits_row);
1101 int nrecvs_cols = nrecvs_rows;
1102 int *olengths_cols = olengths_rows;
1103 PetscInt **rbuf_col = rbuf_row;
1104 if (!square_matrix) {
1108 &lengths_cols[0], &nrecvs_cols);
1113 CHKERR PetscGatherMessageLengths(m_field.
get_comm(), nsends_cols,
1114 nrecvs_cols, &lengths_cols[0],
1115 &onodes_cols, &olengths_cols);
1123 MPI_Request *r_waits_col;
1124 CHKERR PetscPostIrecvInt(m_field.
get_comm(), tag_col, nrecvs_cols,
1125 onodes_cols, olengths_cols, &rbuf_col,
1127 CHKERR PetscFree(onodes_cols);
1129 MPI_Request *s_waits_col;
1130 CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
1133 for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
1134 if (!lengths_cols[proc])
1136 CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0],
1139 tag_col, m_field.
get_comm(), s_waits_col + kk);
1144 CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
1147 CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
1150 CHKERR PetscFree(r_waits_col);
1151 CHKERR PetscFree(s_waits_col);
1154 CHKERR PetscFree(status);
1157 auto hint = dofs_glob_uid_view.begin();
1158 for (
auto dof : *m_field.
get_dofs())
1159 dofs_glob_uid_view.emplace_hint(hint, dof);
1162 for (
int ss = 0; ss != loop_size; ++ss) {
1168 nrecvs = nrecvs_rows;
1169 olengths = olengths_rows;
1170 data_procs = rbuf_row;
1172 nrecvs = nrecvs_cols;
1173 olengths = olengths_cols;
1174 data_procs = rbuf_col;
1178 for (
int kk = 0; kk != nrecvs; ++kk) {
1179 int len = olengths[kk];
1180 int *data_from_proc = data_procs[kk];
1181 for (
int dd = 0;
dd < len;
dd += data_block_size) {
1183 auto ddit = dofs_glob_uid_view.find(uid);
1184 if (ddit == dofs_glob_uid_view.end()) {
1185 std::ostringstream zz;
1186 zz << uid << std::endl;
1188 "no such dof %s in mofem database", zz.str().c_str());
1191 auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
1193 if (dit != numered_dofs_ptr[ss]->end()) {
1198 "received negative dof");
1200 success = numered_dofs_ptr[ss]->modify(
1204 "modification unsuccessful");
1205 success = numered_dofs_ptr[ss]->modify(
1210 "modification unsuccessful");
1212 }
else if (ddit->get()->getPStatus() == 0) {
1219 std::ostringstream zz;
1220 zz << **ddit << std::endl;
1222 "data inconsistency, dofs is not shared, but received "
1231 if (square_matrix) {
1236 CHKERR PetscFree(olengths_rows);
1237 CHKERR PetscFree(rbuf_row[0]);
1238 CHKERR PetscFree(rbuf_row);
1239 if (!square_matrix) {
1240 CHKERR PetscFree(olengths_cols);
1241 CHKERR PetscFree(rbuf_col[0]);
1242 CHKERR PetscFree(rbuf_col);
1245 if (square_matrix) {
1246 if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
1248 "data inconsistency for square_matrix %d!=%d",
1249 numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
1253 "data inconsistency for square_matrix");
1266 const std::string out_name,
const std::vector<std::string> &fields_row,
1267 const std::vector<std::string> &fields_col,
const std::string main_problem,
1268 const bool square_matrix,
1269 const map<std::string, std::pair<EntityType, EntityType>> *entityMapRow,
1270 const map<std::string, std::pair<EntityType, EntityType>> *entityMapCol,
1280 ProblemByName &problems_by_name =
1281 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1284 ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1285 if (out_problem_it == problems_by_name.end()) {
1287 "subproblem with name < %s > not defined (top tip check spelling)",
1292 ProblemByName::iterator main_problem_it = problems_by_name.find(main_problem);
1293 if (main_problem_it == problems_by_name.end()) {
1295 "problem of subproblem with name < %s > not defined (top tip "
1297 main_problem.c_str());
1301 boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1302 out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1304 boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1305 main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1307 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1308 &out_problem_it->nbLocDofsCol};
1310 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1314 out_problem_it->nbGhostDofsRow = 0;
1315 out_problem_it->nbGhostDofsCol = 0;
1319 std::vector<std::string> fields[] = {fields_row, fields_col};
1320 const map<std::string, std::pair<EntityType, EntityType>> *entityMap[] = {
1321 entityMapRow, entityMapCol};
1324 out_problem_it->subProblemData =
1325 boost::make_shared<Problem::SubProblemData>();
1328 for (
int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1331 (*nb_local_dofs[ss]) = 0;
1334 out_problem_dofs[ss]->clear();
1337 out_problem_it->numeredFiniteElementsPtr->clear();
1340 for (
auto field : fields[ss]) {
1346 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1347 boost::make_shared<std::vector<NumeredDofEntity>>();
1350 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1352 out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1356 auto dit = main_problem_dofs[ss]->get<
Unique_mi_tag>().lower_bound(
1358 auto hi_dit = main_problem_dofs[ss]->get<
Unique_mi_tag>().upper_bound(
1361 auto add_dit_to_dofs_array = [&](
auto &dit) {
1362 dofs_array->emplace_back(
1363 dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1364 dit->get()->getPetscGlobalDofIdx(),
1365 dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1368 if (entityMap[ss]) {
1369 auto mit = entityMap[ss]->find(field);
1370 if (mit != entityMap[ss]->end()) {
1371 EntityType lo_type = mit->second.first;
1372 EntityType hi_type = mit->second.second;
1374 for (
auto diit = dit; diit != hi_dit; ++diit) {
1375 EntityType ent_type = (*diit)->getEntType();
1376 if (ent_type >= lo_type && ent_type <= hi_type)
1379 dofs_array->reserve(count);
1380 for (; dit != hi_dit; ++dit) {
1381 EntityType ent_type = (*dit)->getEntType();
1382 if (ent_type >= lo_type && ent_type <= hi_type)
1383 add_dit_to_dofs_array(dit);
1386 dofs_array->reserve(std::distance(dit, hi_dit));
1387 for (; dit != hi_dit; dit++)
1388 add_dit_to_dofs_array(dit);
1391 dofs_array->reserve(std::distance(dit, hi_dit));
1392 for (; dit != hi_dit; dit++)
1393 add_dit_to_dofs_array(dit);
1397 auto hint = out_problem_dofs[ss]->end();
1398 for (
auto &v : *dofs_array)
1399 hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1403 auto dit = out_problem_dofs[ss]->get<
Idx_mi_tag>().begin();
1404 auto hi_dit = out_problem_dofs[ss]->get<
Idx_mi_tag>().end();
1405 for (; dit != hi_dit; dit++) {
1407 if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
1408 idx = (*nb_local_dofs[ss])++;
1410 bool success = out_problem_dofs[ss]->modify(
1411 out_problem_dofs[ss]->project<0>(dit),
1415 "operation unsuccessful");
1425 out_problem_dofs[ss]->size());
1426 const int nb = std::distance(dit, hi_dit);
1428 std::vector<int> main_indices(nb);
1429 for (
auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1430 *it = dit->get()->getPetscGlobalDofIdx();
1434 CHKERR ISCreateGeneral(m_field.
get_comm(), nb, &*main_indices.begin(),
1435 PETSC_USE_POINTER, &is);
1439 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1442 CHKERR ISDuplicate(is, &is_dup);
1447 CHKERR ISDuplicate(is, &is_dup);
1451 CHKERR AOApplicationToPetscIS(ao, is);
1453 CHKERR ISGetSize(is, nb_dofs[ss]);
1457 for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1459 bool success = out_problem_dofs[ss]->modify(
1460 out_problem_dofs[ss]->project<0>(dit),
1462 dit->get()->getPart(), *it, *it,
1463 dit->get()->getPetscLocalDofIdx()));
1466 "operation unsuccessful");
1471 NumeredDofEntityByLocalIdx::iterator dit =
1473 NumeredDofEntityByLocalIdx::iterator hi_dit =
1475 const int nb = std::distance(dit, hi_dit);
1476 std::vector<int> main_indices_non_local(nb);
1477 for (
auto it = main_indices_non_local.begin(); dit != hi_dit;
1479 *it = dit->get()->getPetscGlobalDofIdx();
1483 &*main_indices_non_local.begin(),
1484 PETSC_USE_POINTER, &is);
1485 CHKERR AOApplicationToPetscIS(ao, is);
1488 for (
auto it = main_indices_non_local.begin(); dit != hi_dit;
1490 bool success = out_problem_dofs[ss]->modify(
1491 out_problem_dofs[ss]->project<0>(dit),
1493 dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1494 dit->get()->getPetscLocalDofIdx()));
1497 "operation unsuccessful");
1505 if (square_matrix) {
1506 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1507 out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1508 out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1509 out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1510 out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1511 CHKERR PetscObjectReference(
1512 (PetscObject)out_problem_it->getSubData()->rowIs);
1522 const std::string out_name,
const std::vector<std::string> add_row_problems,
1523 const std::vector<std::string> add_col_problems,
const bool square_matrix,
1538 ProblemByName &problems_by_name =
1539 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1542 ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1543 if (out_problem_it == problems_by_name.end()) {
1545 "problem with name < %s > not defined (top tip check spelling)",
1549 out_problem_it->composedProblemsData =
1550 boost::make_shared<ComposedProblemsData>();
1551 boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1552 out_problem_it->getComposedProblemsData();
1554 const std::vector<std::string> *add_prb[] = {&add_row_problems,
1556 std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1557 &cmp_prb_data->colProblemsAdd};
1558 std::vector<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1561 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1562 &out_problem_it->nbLocDofsCol};
1564 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1568 out_problem_it->nbDofsRow = 0;
1569 out_problem_it->nbDofsCol = 0;
1570 out_problem_it->nbLocDofsRow = 0;
1571 out_problem_it->nbLocDofsCol = 0;
1572 out_problem_it->nbGhostDofsRow = 0;
1573 out_problem_it->nbGhostDofsCol = 0;
1575 int nb_dofs_reserve[] = {0, 0};
1578 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1582 add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1583 add_prb_is[ss]->reserve(add_prb[ss]->size());
1584 for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1585 vit != add_prb[ss]->end(); vit++) {
1586 ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1587 if (prb_it == problems_by_name.end()) {
1590 "problem with name < %s > not defined (top tip check spelling)",
1593 add_prb_ptr[ss]->push_back(&*prb_it);
1597 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1598 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1599 nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1602 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1603 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1604 nb_dofs_reserve[ss] += add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1609 if (square_matrix) {
1610 add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1611 add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1612 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1613 *nb_dofs[1] = *nb_dofs[0];
1614 *nb_local_dofs[1] = *nb_local_dofs[0];
1618 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1620 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1621 dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1622 dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1624 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1626 out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1630 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1631 NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1636 for (
unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1637 PetscInt *dofs_out_idx_ptr;
1638 int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1639 CHKERR PetscMalloc(nb_local_dofs *
sizeof(
int), &dofs_out_idx_ptr);
1641 dit = (*add_prb_ptr[ss])[pp]
1644 hi_dit = (*add_prb_ptr[ss])[pp]
1648 dit = (*add_prb_ptr[ss])[pp]
1651 hi_dit = (*add_prb_ptr[ss])[pp]
1656 for (; dit != hi_dit; dit++) {
1657 BitRefLevel prb_bit = out_problem_it->getBitRefLevel();
1658 BitRefLevel prb_mask = out_problem_it->getMaskBitRefLevel();
1659 BitRefLevel dof_bit = dit->get()->getBitRefLevel();
1660 if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1663 const int part = dit->get()->getPart();
1664 const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1666 (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1668 dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1669 glob_idx, loc_idx, part);
1671 dofs_out_idx_ptr[is_nb++] = glob_idx;
1674 if (is_nb > nb_local_dofs) {
1676 "Data inconsistency");
1679 CHKERR ISCreateGeneral(m_field.
get_comm(), is_nb, dofs_out_idx_ptr,
1680 PETSC_OWN_POINTER, &is);
1681 (*add_prb_is[ss]).push_back(is);
1683 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1684 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1686 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1687 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1689 if (square_matrix) {
1690 (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1691 (*add_prb_is[1]).push_back(is);
1692 CHKERR PetscObjectReference((PetscObject)is);
1697 if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1702 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1703 auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1704 : out_problem_it->numeredColDofsPtr->end();
1705 for (
auto &v : *dofs_array[ss])
1706 hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1707 hint, dofs_array[ss], &v)
1708 : out_problem_it->numeredColDofsPtr->emplace_hint(
1709 hint, dofs_array[ss], &v);
1715 *nb_local_dofs[0] = 0;
1716 *nb_local_dofs[1] = 0;
1717 for (
int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1719 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1721 dofs_ptr = out_problem_it->numeredRowDofsPtr;
1723 dofs_ptr = out_problem_it->numeredColDofsPtr;
1725 NumeredDofEntityByUId::iterator dit, hi_dit;
1728 std::vector<int> idx;
1729 idx.reserve(std::distance(dit, hi_dit));
1731 for (; dit != hi_dit; dit++) {
1732 if (dit->get()->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
1737 "modification unsuccessful");
1739 idx.push_back(dit->get()->getPetscGlobalDofIdx());
1741 if (dit->get()->getPetscLocalDofIdx() != -1) {
1743 "local index should be negative");
1747 if (square_matrix) {
1748 *nb_local_dofs[1] = *nb_local_dofs[0];
1753 CHKERR ISCreateGeneral(m_field.
get_comm(), idx.size(), &*idx.begin(),
1754 PETSC_USE_POINTER, &is);
1755 CHKERR ISGetSize(is, nb_dofs[ss]);
1756 if (square_matrix) {
1757 *nb_dofs[1] = *nb_dofs[0];
1761 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1762 for (
unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1763 CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1767 std::vector<int> idx_new;
1768 idx_new.reserve(dofs_ptr->size());
1769 for (NumeredDofEntityByUId::iterator dit =
1772 idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1777 &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1778 CHKERR AOApplicationToPetscIS(ao, is_new);
1780 std::vector<int>::iterator vit = idx_new.begin();
1781 for (NumeredDofEntityByUId::iterator dit =
1786 dit->get()->getPart(), *(vit++)));
1789 "modification unsuccessful");
1792 CHKERR ISDestroy(&is_new);
1823 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Simple partition problem " << name;
1827 ProblemByName &problems_set =
1828 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
1829 ProblemByName::iterator p_miit = problems_set.find(name);
1830 if (p_miit == problems_set.end()) {
1831 SETERRQ1(PETSC_COMM_SELF, 1,
1832 "problem < %s > is not found (top tip: check spelling)",
1837 NumeredDofEntitysByIdx &dofs_row_by_idx =
1838 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>();
1839 NumeredDofEntitysByIdx &dofs_col_by_idx =
1840 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>();
1847 DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1848 DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1849 nb_row_local_dofs = 0;
1850 nb_row_ghost_dofs = 0;
1851 DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1852 DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1853 nb_col_local_dofs = 0;
1854 nb_col_ghost_dofs = 0;
1856 bool square_matrix =
false;
1857 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1858 square_matrix =
true;
1862 PetscLayout layout_row;
1863 const int *ranges_row;
1865 DofIdx nb_dofs_row = dofs_row_by_idx.size();
1867 CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1868 CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1869 CHKERR PetscLayoutSetUp(layout_row);
1870 CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1872 PetscLayout layout_col;
1873 const int *ranges_col;
1874 if (!square_matrix) {
1875 DofIdx nb_dofs_col = dofs_col_by_idx.size();
1877 CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1878 CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1879 CHKERR PetscLayoutSetUp(layout_col);
1880 CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1884 miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1885 hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1886 if (std::distance(miit_row, hi_miit_row) !=
1887 ranges_row[part + 1] - ranges_row[part]) {
1889 PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1890 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1891 "rstart (%d != %d - %d = %d) ",
1892 std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1893 ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1896 for (; miit_row != hi_miit_row; miit_row++) {
1897 bool success = dofs_row_by_idx.modify(
1902 "modification unsuccessful");
1904 success = dofs_row_by_idx.modify(
1908 "modification unsuccessful");
1911 if (!square_matrix) {
1912 miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1913 hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1914 if (std::distance(miit_col, hi_miit_col) !=
1915 ranges_col[part + 1] - ranges_col[part]) {
1916 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1917 "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1919 "rstart (%d != %d - %d = %d) ",
1920 std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1921 ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1924 for (; miit_col != hi_miit_col; miit_col++) {
1925 bool success = dofs_col_by_idx.modify(
1927 part, (*miit_col)->dofIdx));
1930 "modification unsuccessful");
1932 success = dofs_col_by_idx.modify(
1936 "modification unsuccessful");
1941 CHKERR PetscLayoutDestroy(&layout_row);
1942 if (!square_matrix) {
1943 CHKERR PetscLayoutDestroy(&layout_col);
1945 if (square_matrix) {
1946 nb_col_local_dofs = nb_row_local_dofs;
1947 nb_col_ghost_dofs = nb_row_ghost_dofs;
1960 MOFEM_LOG(
"WORLD", Sev::noisy) <<
"Partition problem " << name;
1962 using NumeredDofEntitysByIdx =
1967 auto &problems_set =
1968 const_cast<ProblemsByName &
>(problems_ptr->get<
Problem_mi_tag>());
1969 auto p_miit = problems_set.find(name);
1970 if (p_miit == problems_set.end()) {
1972 "problem with name %s not defined (top tip check spelling)",
1975 int nb_dofs_row = p_miit->getNbDofsRow();
1985 "entFEAdjacencies not build");
1991 ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1996 MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1999 MatPartitioning part;
2002 CHKERR MatPartitioningSetAdjacency(part, Adj);
2003 CHKERR MatPartitioningSetFromOptions(part);
2005 CHKERR MatPartitioningApply(part, &is);
2007 ISView(is, PETSC_VIEWER_STDOUT_WORLD);
2010 IS is_gather, is_num, is_gather_num;
2011 CHKERR ISAllGather(is, &is_gather);
2012 CHKERR ISPartitioningToNumbering(is, &is_num);
2013 CHKERR ISAllGather(is_num, &is_gather_num);
2014 const int *part_number, *petsc_idx;
2015 CHKERR ISGetIndices(is_gather, &part_number);
2016 CHKERR ISGetIndices(is_gather_num, &petsc_idx);
2017 int size_is_num, size_is_gather;
2018 CHKERR ISGetSize(is_gather, &size_is_gather);
2019 if (size_is_gather != (
int)nb_dofs_row)
2020 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
2021 "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
2023 CHKERR ISGetSize(is_num, &size_is_num);
2024 if (size_is_num != (
int)nb_dofs_row)
2025 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
2026 "data inconsistency %d != %d", size_is_num, nb_dofs_row);
2028 bool square_matrix =
false;
2029 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
2030 square_matrix =
true;
2056 auto number_dofs = [&](
auto &dofs_idx,
auto &counter) {
2058 for (
auto miit_dofs_row = dofs_idx.begin();
2059 miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
2060 const int part = part_number[(*miit_dofs_row)->dofIdx];
2062 const bool success = dofs_idx.modify(
2065 part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
2068 "modification unsuccessful");
2071 const bool success = dofs_idx.modify(
2073 part, petsc_idx[(*miit_dofs_row)->dofIdx]));
2076 "modification unsuccessful");
2084 auto &dofs_row_by_idx_no_const =
const_cast<NumeredDofEntitysByIdx &
>(
2085 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>());
2086 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2087 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2088 nb_row_local_dofs = 0;
2089 nb_row_ghost_dofs = 0;
2091 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
2093 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2094 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2095 if (square_matrix) {
2096 nb_col_local_dofs = nb_row_local_dofs;
2097 nb_col_ghost_dofs = nb_row_ghost_dofs;
2099 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2100 const_cast<NumeredDofEntitysByIdx &
>(
2101 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>());
2102 nb_col_local_dofs = 0;
2103 nb_col_ghost_dofs = 0;
2104 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
2107 CHKERR ISRestoreIndices(is_gather, &part_number);
2108 CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
2109 CHKERR ISDestroy(&is_num);
2110 CHKERR ISDestroy(&is_gather_num);
2111 CHKERR ISDestroy(&is_gather);
2113 CHKERR MatPartitioningDestroy(&part);
2118 auto number_dofs = [&](
auto &dof_idx,
auto &counter) {
2120 for (
auto miit_dofs_row = dof_idx.begin();
2121 miit_dofs_row != dof_idx.end(); miit_dofs_row++) {
2122 const bool success = dof_idx.modify(
2128 "modification unsuccessful");
2134 auto &dofs_row_by_idx_no_const =
const_cast<NumeredDofEntitysByIdx &
>(
2135 p_miit->numeredRowDofsPtr->get<
Idx_mi_tag>());
2136 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
2137 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2138 nb_row_local_dofs = 0;
2139 nb_row_ghost_dofs = 0;
2141 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
2143 bool square_matrix =
false;
2144 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
2145 square_matrix =
true;
2147 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
2148 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2149 if (square_matrix) {
2150 nb_col_local_dofs = nb_row_local_dofs;
2151 nb_col_ghost_dofs = nb_row_ghost_dofs;
2153 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
2154 const_cast<NumeredDofEntitysByIdx &
>(
2155 p_miit->numeredColDofsPtr->get<
Idx_mi_tag>());
2156 nb_col_local_dofs = 0;
2157 nb_col_ghost_dofs = 0;
2158 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
2167 const std::string name,
const std::string problem_for_rows,
bool copy_rows,
2168 const std::string problem_for_cols,
bool copy_cols,
int verb) {
2179 ProblemByName &problems_by_name =
2180 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
2181 ProblemByName::iterator p_miit = problems_by_name.find(name);
2182 if (p_miit == problems_by_name.end()) {
2184 "problem with name < %s > not defined (top tip check spelling)",
2189 << p_miit->getName() <<
" from rows of " << problem_for_rows
2190 <<
" and columns of " << problem_for_cols;
2193 ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
2194 if (p_miit_row == problems_by_name.end()) {
2196 "problem with name < %s > not defined (top tip check spelling)",
2197 problem_for_rows.c_str());
2199 const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
2200 p_miit_row->numeredRowDofsPtr;
2203 ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
2204 if (p_miit_col == problems_by_name.end()) {
2206 "problem with name < %s > not defined (top tip check spelling)",
2207 problem_for_cols.c_str());
2209 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
2210 p_miit_col->numeredColDofsPtr;
2212 bool copy[] = {copy_rows, copy_cols};
2213 boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
2214 p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
2216 int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
2217 int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
2218 boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
2221 for (
int ss = 0; ss < 2; ss++) {
2224 *nb_local_dofs[ss] = 0;
2228 std::vector<int> is_local, is_new;
2232 for (NumeredDofEntity_multiIndex::iterator dit =
2233 composed_dofs[ss]->begin();
2234 dit != composed_dofs[ss]->end(); dit++) {
2235 NumeredDofEntityByUId::iterator diit =
2236 dofs_by_uid.find((*dit)->getLocalUniqueId());
2237 if (diit == dofs_by_uid.end()) {
2240 "data inconsistency, could not find dof in composite problem");
2242 int part_number = (*diit)->getPart();
2243 int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
2245 success = composed_dofs[ss]->modify(
2250 "modification unsuccessful");
2252 if ((*dit)->getPart() == (
unsigned int)m_field.
get_comm_rank()) {
2253 success = composed_dofs[ss]->modify(
2257 "modification unsuccessful");
2259 is_local.push_back(petsc_global_dof);
2264 CHKERR AOCreateMapping(m_field.
get_comm(), is_local.size(), &is_local[0],
2269 for (NumeredDofEntity_multiIndex::iterator dit =
2270 composed_dofs[ss]->begin();
2271 dit != composed_dofs[ss]->end(); dit++) {
2272 is_local.push_back((*dit)->getPetscGlobalDofIdx());
2274 CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2276 for (NumeredDofEntity_multiIndex::iterator dit =
2277 composed_dofs[ss]->begin();
2278 dit != composed_dofs[ss]->end(); dit++) {
2279 int part_number = (*dit)->getPart();
2280 int petsc_global_dof = is_local[idx2++];
2282 success = composed_dofs[ss]->modify(
2287 "modification unsuccessful");
2295 for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2296 dit != copied_dofs[ss]->end(); dit++) {
2297 std::pair<NumeredDofEntity_multiIndex::iterator, bool>
p;
2298 p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2303 int dof_idx = (*dit)->getDofIdx();
2304 int part_number = (*dit)->getPart();
2305 int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2307 const bool success = composed_dofs[ss]->modify(
2309 part_number, dof_idx, petsc_global_dof,
2310 (*nb_local_dofs[ss])++));
2313 "modification unsuccessful");
2316 const bool success = composed_dofs[ss]->modify(
2318 part_number, dof_idx, petsc_global_dof));
2321 "modification unsuccessful");
2342 << problem_ptr->
getName() <<
" Nb. local dof "
2360 NumeredDofEntitysByIdx;
2361 NumeredDofEntitysByIdx::iterator dit, hi_dit;
2362 const NumeredDofEntitysByIdx *numered_dofs_ptr[] = {
2370 for (
int ss = 0; ss < 2; ss++) {
2372 dit = numered_dofs_ptr[ss]->begin();
2373 hi_dit = numered_dofs_ptr[ss]->end();
2374 for (; dit != hi_dit; dit++) {
2376 if ((*dit)->getPetscLocalDofIdx() < 0) {
2377 std::ostringstream zz;
2380 "local dof index for %d (0-row, 1-col) not set, i.e. has "
2381 "negative value\n %s",
2382 ss, zz.str().c_str());
2384 if ((*dit)->getPetscLocalDofIdx() >= *local_nbdof_ptr[ss]) {
2385 std::ostringstream zz;
2388 "local dofs for %d (0-row, 1-col) out of range\n %s", ss,
2392 if ((*dit)->getPetscGlobalDofIdx() < 0) {
2393 std::ostringstream zz;
2395 << dit->get()->getBitRefLevel() <<
" " << **dit;
2397 "global dof index for %d (0-row, 1-col) row not set, i.e. "
2398 "has negative value\n %s",
2399 ss, zz.str().c_str());
2401 if ((*dit)->getPetscGlobalDofIdx() >= *nbdof_ptr[ss]) {
2402 std::ostringstream zz;
2404 << *nbdof_ptr[ss] <<
" " << **dit;
2406 "global dofs for %d (0-row, 1-col) out of range\n %s", ss,
2417 bool part_from_moab,
2419 int hi_proc,
int verb) {
2433 "adjacencies not build");
2440 "problem not partitioned");
2450 const_cast<ProblemByName &
>(problems_ptr->get<
Problem_mi_tag>());
2451 ProblemByName::iterator p_miit = problems.find(name);
2452 if (p_miit == problems.end()) {
2454 "problem < %s > not found (top tip: check spelling)",
2460 *p_miit->numeredFiniteElementsPtr;
2463 problem_finite_elements.clear();
2467 bool do_cols_prob =
true;
2468 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2469 do_cols_prob =
false;
2472 auto get_good_elems = [&]() {
2473 auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
2474 good_elems.reserve(fe_ent_ptr->size());
2476 const auto prb_bit = p_miit->getBitRefLevel();
2477 const auto prb_mask = p_miit->getMaskBitRefLevel();
2481 for (
auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2484 if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2486 const auto fe_bit = (*efit)->getBitRefLevel();
2489 if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit) == prb_bit)
2490 good_elems.emplace_back(efit);
2497 auto good_elems = get_good_elems();
2499 auto numbered_good_elems_ptr =
2500 boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2501 numbered_good_elems_ptr->reserve(good_elems.size());
2502 for (
auto &efit : good_elems)
2505 if (!do_cols_prob) {
2506 for (
auto &fe : *numbered_good_elems_ptr) {
2507 if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2508 fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2513 if (part_from_moab) {
2514 for (
auto &fe : *numbered_good_elems_ptr) {
2516 int proc = fe.getPartProc();
2517 if (proc == -1 && fe.getEntType() == MBVERTEX)
2518 proc = fe.getOwnerProc();
2523 for (
auto &fe : *numbered_good_elems_ptr) {
2526 CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2528 if (!part_from_moab) {
2530 for (
auto &dof_ptr : rows_view)
2531 parts[dof_ptr->pArt]++;
2532 std::vector<int>::iterator pos = max_element(parts.begin(), parts.end());
2533 const auto max_part = std::distance(parts.begin(), pos);
2538 for (
auto &fe : *numbered_good_elems_ptr) {
2540 auto check_fields_and_dofs = [&]() {
2541 if (!part_from_moab) {
2543 if (fe.getBitFieldIdRow().none() && m_field.
get_comm_size() == 0) {
2545 <<
"At least one field has to be added to element row to "
2546 "determine partition of finite element. Check element " +
2547 boost::lexical_cast<std::string>(fe.getName());
2556 return (!fe.sPtr->getRowFieldEnts().empty() ||
2557 !fe.sPtr->getColFieldEnts().empty())
2561 (fe.getBitFieldIdRow().none() || fe.getBitFieldIdCol().none());
2564 if (check_fields_and_dofs()) {
2566 auto p = problem_finite_elements.insert(
2567 boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
2575 auto elements_on_rank =
2576 problem_finite_elements.get<
Part_mi_tag>().equal_range(
2579 << p_miit->getName() <<
" nb. elems "
2580 << std::distance(elements_on_rank.first, elements_on_rank.second);
2582 for (
auto &fe : *fe_ptr) {
2588 <<
"Element " << fe->getName() <<
" nb. elems "
2589 << std::distance(e_range.first, e_range.second);
2607 "partition of problem not build");
2610 "partitions finite elements not build");
2619 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2620 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2621 nb_row_ghost_dofs = 0;
2622 nb_col_ghost_dofs = 0;
2632 p_miit->numeredFiniteElementsPtr->get<
Part_mi_tag>().equal_range(
2638 using Vec = std::vector<boost::shared_ptr<NumeredDofEntity>>;
2639 using It = Vec::iterator;
2640 It operator()(
Vec &dofs_view, It &hint,
2641 boost::shared_ptr<NumeredDofEntity> &&dof) {
2642 dofs_view.emplace_back(dof);
2643 return dofs_view.end();
2648 std::vector<boost::shared_ptr<NumeredDofEntity>> fe_vec_view;
2649 auto hint_r = ghost_idx_row_view.begin();
2650 for (
auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2652 fe_vec_view.clear();
2654 *(p_miit->getNumeredRowDofsPtr()),
2655 fe_vec_view, Inserter());
2657 for (
auto &dof_ptr : fe_vec_view) {
2658 if (dof_ptr->getPart() != (
unsigned int)m_field.
get_comm_rank()) {
2659 hint_r = ghost_idx_row_view.emplace_hint(hint_r, dof_ptr);
2665 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2667 auto hint_c = ghost_idx_col_view.begin();
2668 for (
auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2670 fe_vec_view.clear();
2672 (*fe_ptr)->getColFieldEnts(), *(p_miit->getNumeredColDofsPtr()),
2673 fe_vec_view, Inserter());
2675 for (
auto &dof_ptr : fe_vec_view) {
2676 if (dof_ptr->getPart() != (
unsigned int)m_field.
get_comm_rank()) {
2677 hint_c = ghost_idx_col_view.emplace_hint(hint_c, dof_ptr);
2683 int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2684 int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2687 &ghost_idx_row_view, &ghost_idx_col_view};
2693 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2698 for (
int ss = 0; ss != loop_size; ++ss) {
2699 for (
auto &gid : *ghost_idx_view[ss]) {
2700 NumeredDofEntityByUId::iterator dof =
2701 dof_by_uid_no_const[ss]->find(gid->getLocalUniqueId());
2702 if (PetscUnlikely((*dof)->petscLocalDofIdx != (
DofIdx)-1))
2704 "inconsistent data, ghost dof already set");
2705 bool success = dof_by_uid_no_const[ss]->modify(
2707 if (PetscUnlikely(!success))
2709 "modification unsuccessful");
2710 (*nb_ghost_dofs[ss])++;
2713 if (loop_size == 1) {
2714 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2720 <<
" FEs ghost dofs on problem " << p_miit->getName()
2721 <<
" Nb. ghost dof " << p_miit->getNbGhostDofsRow() <<
" by "
2722 << p_miit->getNbGhostDofsCol() <<
" Nb. local dof "
2723 << p_miit->getNbLocalDofsCol() <<
" by " << p_miit->getNbLocalDofsCol();
2741 "partition of problem not build");
2744 "partitions finite elements not build");
2749 ProblemsByName &problems_set =
2750 const_cast<ProblemsByName &
>(problems_ptr->get<
Problem_mi_tag>());
2751 ProblemsByName::iterator p_miit = problems_set.find(name);
2755 DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2756 &(p_miit->nbGhostDofsCol)};
2757 DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2758 for (
int ss = 0; ss != 2; ++ss) {
2759 (*nb_ghost_dofs[ss]) = 0;
2766 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2770 typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2771 NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2772 p_miit->numeredColDofsPtr};
2775 for (
int ss = 0; ss != loop_size; ++ss) {
2780 std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2781 ghost_idx_view.reserve(std::distance(
r.first,
r.second));
2782 for (;
r.first !=
r.second; ++
r.first)
2783 ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(
r.first));
2785 auto cmp = [](
auto a,
auto b) {
2786 return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2788 sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2791 for (
auto gid_it : ghost_idx_view) {
2792 bool success = numered_dofs[ss]->modify(
2796 "modification unsuccessful");
2797 ++(*nb_ghost_dofs[ss]);
2800 if (loop_size == 1) {
2801 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2807 <<
" FEs ghost dofs on problem " << p_miit->getName()
2808 <<
" Nb. ghost dof " << p_miit->getNbGhostDofsRow() <<
" by "
2809 << p_miit->getNbGhostDofsCol() <<
" Nb. local dof "
2810 << p_miit->getNbLocalDofsCol() <<
" by " << p_miit->getNbLocalDofsCol();
2820 const std::string fe_name,
2821 EntityHandle *meshset)
const {
2830 .lower_bound(fe_name);
2833 .upper_bound(fe_name);
2834 std::vector<EntityHandle> fe_vec;
2835 fe_vec.reserve(std::distance(fit, hi_fe_it));
2836 for (; fit != hi_fe_it; fit++)
2837 fe_vec.push_back(fit->get()->getEnt());
2838 CHKERR m_field.
get_moab().add_entities(*meshset, &*fe_vec.begin(),
2845 const std::string fe_name,
2846 PetscLayout *layout)
const {
2857 const std::string problem_name,
const std::string field_name,
2858 const Range ents,
const int lo_coeff,
const int hi_coeff,
int verb,
2883 for (
int s = 0; s != 2; ++s)
2884 if (numered_dofs[s]) {
2886 typedef multi_index_container<
2888 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2891 NumeredDofEntity_it_view_multiIndex;
2894 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2897 for (
auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2899 auto lo = numered_dofs[s]->get<
Unique_mi_tag>().lower_bound(
2901 auto hi = numered_dofs[s]->get<
Unique_mi_tag>().upper_bound(
2904 for (; lo != hi; ++lo)
2905 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2906 (*lo)->getDofCoeffIdx() <= hi_coeff)
2907 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2911 for (
auto &dof : dofs_it_view)
2919 for (
auto dit : dofs_it_view) {
2920 bool success = numered_dofs[s]->modify(dit, mod);
2923 "can not set negative indices");
2927 std::vector<boost::weak_ptr<NumeredDofEntity>> dosf_weak_view;
2928 dosf_weak_view.reserve(dofs_it_view.size());
2929 for (
auto dit : dofs_it_view)
2930 dosf_weak_view.push_back(*dit);
2934 "Number of DOFs in multi-index %d and to delete %d\n",
2935 numered_dofs[s]->size(), dofs_it_view.size());
2938 for (
auto weak_dit : dosf_weak_view)
2939 if (
auto dit = weak_dit.lock()) {
2940 numered_dofs[s]->erase(dit->getLocalUniqueId());
2945 "Number of DOFs in multi-index after delete %d\n",
2946 numered_dofs[s]->size());
2949 int nb_local_dofs = 0;
2950 int nb_ghost_dofs = 0;
2951 for (
auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2953 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2954 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2956 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2960 "Imposible case. You could run problem on no distributed "
2961 "mesh. That is not implemented.");
2965 auto get_indices_by_tag = [&](
auto tag,
auto &indices,
bool only_local) {
2966 const int nb_dofs = numered_dofs[s]->size();
2968 indices.reserve(nb_dofs);
2969 for (
auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2970 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2973 if ((*dit)->getPetscLocalDofIdx() < 0 ||
2974 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2979 indices.push_back(decltype(tag)::get_index(dit));
2983 auto get_indices_by_uid = [&](
auto tag,
auto &indices) {
2984 const int nb_dofs = numered_dofs[s]->size();
2986 indices.reserve(nb_dofs);
2987 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2989 indices.push_back(decltype(tag)::get_index(dit));
2992 auto concatenate_dofs = [&](
auto tag,
auto &indices,
2993 const auto local_only) {
2995 get_indices_by_tag(tag, indices, local_only);
3000 &*indices.begin(), PETSC_NULL, &ao);
3002 CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
3003 &*indices.begin(), PETSC_NULL, &ao);
3005 get_indices_by_uid(tag, indices);
3006 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3012 auto set_concatinated_indices = [&]() {
3013 std::vector<int> global_indices;
3014 std::vector<int> local_indices;
3018 auto gi = global_indices.begin();
3019 auto li = local_indices.begin();
3020 for (
auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3023 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3024 bool success = numered_dofs[s]->modify(dit, mod);
3027 "can not set negative indices");
3033 CHKERR set_concatinated_indices();
3035 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3037 *(local_nbdof_ptr[s]) = nb_local_dofs;
3038 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3041 for (
auto dof : (*numered_dofs[s])) {
3042 if (dof->getPetscGlobalDofIdx() < 0) {
3044 "Negative global idx");
3046 if (dof->getPetscLocalDofIdx() < 0) {
3048 "Negative local idx");
3054 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3055 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3056 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3061 "removed ents on rank %d from problem %s dofs [ %d / %d "
3063 "%d) local, %d / %d (before %d / %d) "
3064 "ghost, %d / %d (before %d / %d) global]",
3067 nb_init_row_dofs, nb_init_col_dofs,
3069 nb_init_ghost_row_dofs, nb_init_ghost_col_dofs,
3071 nb_init_loc_row_dofs, nb_init_loc_col_dofs);
3080 std::vector<bool> &marker) {
3086 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
3096 marker.resize(dofs->size());
3097 std::fill(marker.begin(), marker.end(),
false);
3098 for (
auto p = ents.pair_begin();
p != ents.pair_end(); ++
p) {
3099 auto lo = dofs->get<
Ent_mi_tag>().lower_bound(
p->first);
3100 auto hi = dofs->get<
Ent_mi_tag>().upper_bound(
p->second);
3101 for (; lo != hi; ++lo)
3102 marker[(*lo)->getPetscLocalDofIdx()] =
true;
#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< 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 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 Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
#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 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 partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildCompsedProblem(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 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, int verb=VERBOSE, const bool debug=false)
Remove DOFs from 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 buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
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
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
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
static const MOFEMuuid IDD_MOFEMProblemsManager
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
const double r
rate factor
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 boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
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)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
IdxDataType(const UId uid, const int dof)
IdxDataTypePtr(const int *ptr)
MoFEM interface unique ID.
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
BitRefLevel getMaskBitRefLevel() const
std::string getName() const
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
DofIdx nbGhostDofsCol
Number of ghost DOFs in col.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
boost::shared_ptr< SequenceDofContainer > & getRowDofsSequence() const
Get reference to sequence data numbered dof container.
boost::shared_ptr< SequenceDofContainer > & getColDofsSequence() const
Get reference to sequence data numbered dof container.
DofIdx getNbGhostDofsCol() const
boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
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.
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
Problem manager is used to build and partition problems.
ProblemsManager(const MoFEM::Core &core)
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
PetscLogEvent MOFEM_EVENT_ProblemsManager
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
PetscBool buildProblemFromFields
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const Range ents, std::vector< bool > &marker)
Create vector with marked indices.
MoFEMErrorCode getOptions()
base class for all interface classes
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.