24MoFEMErrorCode MatPartitioningApply_Parmetis_MoFEM(MatPartitioning part,
38 : cOre(const_cast<
MoFEM::
Core &>(core)), dEbug(false) {}
47 std::vector<std::vector<EntityHandle>> sbuffer(m_field.
get_comm_size());
49 Range::iterator eit = ents.begin();
50 for (; eit != ents.end(); eit++) {
52 auto meit = ref_ents_ptr->get<
Ent_mi_tag>().find(*eit);
53 if (meit == ref_ents_ptr->get<
Ent_mi_tag>().end()) {
56 "rank %d entity %lu not exist on database, local entity can not "
57 "be found for this owner",
61 unsigned char pstatus = (*meit)->getPStatus();
67 MOFEM_LOG(
"SYNC", Sev::noisy) <<
"pstatus " << std::bitset<8>(pstatus);
71 proc < MAX_SHARING_PROCS && -1 != (*meit)->getSharingProcsPtr()[proc];
73 if ((*meit)->getSharingProcsPtr()[proc] == -1)
75 "sharing processor not set");
77 if ((*meit)->getSharingProcsPtr()[proc] == m_field.
get_comm_rank())
80 EntityHandle handle_on_sharing_proc =
81 (*meit)->getSharingHandlersPtr()[proc];
82 sbuffer[(*meit)->getSharingProcsPtr()[proc]].push_back(
83 handle_on_sharing_proc);
85 MOFEM_LOG_C(
"SYNC", Sev::noisy,
"send %lu (%lu) to %d at %d\n",
86 (*meit)->getEnt(), handle_on_sharing_proc,
87 (*meit)->getSharingProcsPtr()[proc],
90 if (!(pstatus & PSTATUS_MULTISHARED))
96 std::vector<int> sbuffer_lengths(
98 const size_t block_size =
sizeof(EntityHandle) /
sizeof(
int);
101 if (!sbuffer[proc].empty()) {
103 sbuffer_lengths[proc] = sbuffer[proc].size() * block_size;
108 sbuffer_lengths[proc] = 0;
120 CHKERR PetscGatherNumberOfMessages(comm, NULL, &sbuffer_lengths[0], &nrecvs);
126 CHKERR PetscGatherMessageLengths(comm, nsends, nrecvs, &sbuffer_lengths[0],
134 CHKERR PetscCommGetNewTag(comm, &tag);
139 MPI_Request *r_waits;
143 CHKERR PetscPostIrecvInt(comm, tag, nrecvs, onodes, olengths, &rbuf,
146 MPI_Request *s_waits;
147 CHKERR PetscMalloc1(nsends, &s_waits);
150 for (
int proc = 0, kk = 0; proc < m_field.
get_comm_size(); proc++) {
151 if (!sbuffer_lengths[proc])
153 CHKERR MPI_Isend(&(sbuffer[proc])[0],
154 sbuffer_lengths[proc],
156 tag, comm, s_waits + kk);
162 CHKERR MPI_Waitall(nrecvs, r_waits, &status[0]);
166 CHKERR MPI_Waitall(nsends, s_waits, &status[0]);
169 MOFEM_LOG_C(
"SYNC", Sev::verbose,
"Rank %d nb. before ents %u\n",
174 for (
int kk = 0; kk < nrecvs; kk++) {
176 int len = olengths[kk];
177 int *data_from_proc = rbuf[kk];
179 for (
int ee = 0; ee < len; ee += block_size) {
182 bcopy(&data_from_proc[ee], &ent,
sizeof(EntityHandle));
183 auto meit = ref_ents_ptr->get<
Ent_mi_tag>().find(ent);
184 if (meit == ref_ents_ptr->get<
Ent_mi_tag>().end())
186 "rank %d entity %lu not exist on database, local entity can "
187 "not be found for this owner",
191 MOFEM_LOG_C(
"SYNC", Sev::verbose,
"received %ul (%ul) from %d at %d\n",
192 (*meit)->getEnt(), ent, onodes[kk],
195 ents.insert((*meit)->getEnt());
200 PetscSynchronizedPrintf(m_field.
get_comm(),
"Rank %d nb. after ents %u\n",
204 CHKERR PetscFree(s_waits);
205 CHKERR PetscFree(rbuf[0]);
207 CHKERR PetscFree(r_waits);
209 CHKERR PetscFree(olengths);
210 CHKERR PetscCommDestroy(&comm);
224 CHKERR m_field.
get_moab().get_entities_by_handle(idm, ents,
false);
231 const Problem *problem_ptr,
const std::string &fe_name,
int verb) {
234 ParallelComm *pcomm = ParallelComm::get_pcomm(
236 std::vector<int> shprocs(MAX_SHARING_PROCS, 0);
237 std::vector<EntityHandle> shhandles(MAX_SHARING_PROCS, 0);
239 Tag th_gid = m_field.
get_moab().globalId_tag();
244 CHKERR PetscLayoutGetRange(layout, &gid, &last_gid);
245 CHKERR PetscLayoutDestroy(&layout);
247 EntityHandle ent = (*fe_it)->getEnt();
249 unsigned int part = (*fe_it)->getPart();
250 CHKERR m_field.
get_moab().tag_set_data(pcomm->part_tag(), &ent, 1, &part);
251 if (part == pcomm->rank()) {
258 if (pcomm->size() > 1) {
260 unsigned char pstatus = 0;
261 if (pcomm->rank() != part) {
262 pstatus = PSTATUS_NOT_OWNED;
263 pstatus |= PSTATUS_GHOST;
266 if (pcomm->size() > 2) {
267 pstatus |= PSTATUS_SHARED;
268 pstatus |= PSTATUS_MULTISHARED;
270 pstatus |= PSTATUS_SHARED;
274 for (
size_t rr = 0; rr < pcomm->size(); ++rr) {
275 if (rr != pcomm->rank()) {
276 shhandles[rrr] = ent;
281 for (; rrr != pcomm->size(); ++rrr)
284 if (pstatus & PSTATUS_SHARED) {
285 CHKERR m_field.
get_moab().tag_set_data(pcomm->sharedp_tag(), &ent, 1,
287 CHKERR m_field.
get_moab().tag_set_data(pcomm->sharedh_tag(), &ent, 1,
291 if (pstatus & PSTATUS_MULTISHARED) {
292 CHKERR m_field.
get_moab().tag_set_data(pcomm->sharedps_tag(), &ent, 1,
294 CHKERR m_field.
get_moab().tag_set_data(pcomm->sharedhs_tag(), &ent, 1,
297 CHKERR m_field.
get_moab().tag_set_data(pcomm->pstatus_tag(), &ent, 1,
301 CHKERR pcomm->exchange_tags(th_gid, ents);
306 const std::string &name,
const std::string &fe_name,
int verb) {
317 const int num_entities,
318 const int owner_proc,
int verb) {
324 ParallelComm *pcomm = ParallelComm::get_pcomm(
327 Range all_ents_range;
328 all_ents_range.insert_list(entities, entities + num_entities);
330 auto get_tag = [&]() {
return m_field.
get_moab().globalId_tag(); };
332 auto delete_tag = [&](
auto &&th_gid) {
338 auto resolve_shared_ents = [&](
auto &&th_gid,
auto &all_ents_range) {
339 auto set_gid = [&](
auto &th_gid) {
340 std::vector<int> gids(num_entities);
341 for (
size_t g = 0;
g != all_ents_range.size(); ++
g)
349 auto get_skin_ents = [&](
auto &all_ents_range) {
350 std::array<Range, 4> proc_ents_skin;
351 proc_ents_skin[3] = all_ents_range.subset_by_dimension(3);
352 proc_ents_skin[2] = all_ents_range.subset_by_dimension(2);
353 proc_ents_skin[1] = all_ents_range.subset_by_dimension(1);
354 proc_ents_skin[0] = all_ents_range.subset_by_dimension(0);
355 return proc_ents_skin;
358 auto resolve_dim = [&](
auto &all_ents_range) {
359 for (
int resolve_dim = 3; resolve_dim >= 0; --resolve_dim) {
360 if (all_ents_range.num_of_dimension(resolve_dim))
366 auto get_proc_ent = [&](
auto &all_ents_range) {
369 proc_ent = all_ents_range;
373 auto resolve_shared_ents = [&](
auto &&proc_ents,
auto &&skin_ents) {
374 return pcomm->resolve_shared_ents(
375 0, proc_ents, resolve_dim(all_ents_range),
376 resolve_dim(all_ents_range), skin_ents.data(), set_gid(th_gid));
379 CHKERR resolve_shared_ents(get_proc_ent(all_ents_range),
380 get_skin_ents(all_ents_range));
385 CHKERR delete_tag(resolve_shared_ents(get_tag(), all_ents_range));
389 auto print_owner = [&](
const EntityHandle e) {
392 EntityHandle moab_owner_handle;
393 CHKERR pcomm->get_owner_handle(e, moab_owner_proc, moab_owner_handle);
395 unsigned char pstatus = 0;
397 CHKERR m_field.
get_moab().tag_get_data(pcomm->pstatus_tag(), &e, 1,
400 std::vector<int> shprocs(MAX_SHARING_PROCS, 0);
401 std::vector<EntityHandle> shhandles(MAX_SHARING_PROCS, 0);
403 CHKERR m_field.
get_moab().tag_get_data(pcomm->sharedp_tag(), &e, 1,
405 CHKERR m_field.
get_moab().tag_get_data(pcomm->sharedh_tag(), &e, 1,
407 if (pstatus & PSTATUS_MULTISHARED) {
408 CHKERR m_field.
get_moab().tag_get_data(pcomm->sharedps_tag(), &e, 1,
410 CHKERR m_field.
get_moab().tag_get_data(pcomm->sharedhs_tag(), &e, 1,
414 std::ostringstream ss;
417 if (!(pstatus & PSTATUS_NOT_OWNED))
419 if (pstatus & PSTATUS_SHARED)
420 ss <<
"PSTATUS_SHARED ";
421 if (pstatus & PSTATUS_MULTISHARED)
422 ss <<
"PSTATUS_MULTISHARED ";
424 ss <<
"owner " << moab_owner_proc <<
" (" << owner_proc <<
") ";
428 ss << shprocs[
r] <<
" ";
432 ss << shhandles[
r] <<
" ";
435 MOFEM_LOG(
"SYNC", Sev::noisy) << ss.str();
441 for (
auto e : all_ents_range)
450 const int owner_proc,
455 const int num_ents = entities.size();
456 std::vector<EntityHandle> vec_ents(num_ents);
457 std::copy(entities.begin(), entities.end(), vec_ents.begin());
466 const int owner_proc,
int verb) {
471 std::vector<EntityHandle> field_ents;
472 CHKERR m_field.
get_moab().get_entities_by_handle(field_meshset, field_ents,
486 Range exchange_ents_data_verts, exchange_ents_data;
495 for (
auto it = lo; it != hi; ++it)
498 ((*it)->getPStatus()) &&
500 (*it)->getNbDofsOnEnt()
503 if ((*it)->getEntType() == MBVERTEX)
504 exchange_ents_data_verts.insert((*it)->getEnt());
506 exchange_ents_data.insert((*it)->getEnt());
510 ParallelComm *pcomm = ParallelComm::get_pcomm(
513 auto exchange = [&](
const Range &ents, Tag
th) {
516 std::vector<Tag> tags;
518 CHKERR pcomm->exchange_tags(tags, tags, ents);
523 CHKERR exchange(exchange_ents_data_verts, field_ptr->th_FieldDataVerts);
524 CHKERR exchange(exchange_ents_data, field_ptr->th_FieldData);
531 const int adj_dim,
const int n_parts,
532 Tag *th_vertex_weights, Tag *th_edge_weights,
533 Tag *th_part_weights,
int verb,
const bool debug) {
538 int rstart, rend, nb_elems;
542 CHKERR PetscLayoutSetBlockSize(layout, 1);
543 CHKERR PetscLayoutSetSize(layout, ents.size());
544 CHKERR PetscLayoutSetUp(layout);
545 CHKERR PetscLayoutGetSize(layout, &nb_elems);
546 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
547 CHKERR PetscLayoutDestroy(&layout);
550 <<
"Finite elements in problem: row lower " << rstart <<
" row upper "
551 << rend <<
" nb. elems " << nb_elems <<
" ( " << ents.size() <<
" )";
556 std::vector<EntityHandle> weight_ents;
557 weight_ents.reserve(rend - rstart + 1);
561 std::vector<int> adj;
562 AdjBridge(
const EntityHandle ent, std::vector<int> &adj)
563 : ent(ent), adj(adj) {}
566 typedef multi_index_container<
570 hashed_unique<member<AdjBridge, EntityHandle, &AdjBridge::ent>>
575 auto get_it = [&](
auto i) {
576 auto it = ents.begin();
578 if (it == ents.end())
586 proc_ents.insert(get_it(rstart), get_it(rend));
587 if (proc_ents.size() != rend - rstart)
589 "Wrong number of elements in range %d != %d", proc_ents.size(),
594 proc_ents, adj_dim,
true, all_dim_ents, moab::Interface::UNION);
596 AdjBridgeMap adj_bridge_map;
597 auto hint = adj_bridge_map.begin();
598 std::vector<int> adj;
599 for (
auto ent : all_dim_ents) {
602 adj_ents = intersect(adj_ents, ents);
604 adj.reserve(adj_ents.size());
605 for (
auto a : adj_ents)
606 adj.emplace_back(ents.index(
a));
607 hint = adj_bridge_map.emplace_hint(hint, ent, adj);
613 const int nb_loc_elements = rend - rstart;
614 std::vector<int>
i(nb_loc_elements + 1, 0),
j;
616 std::vector<int> row_adj;
617 Range::iterator fe_it;
620 for (fe_it = proc_ents.begin(), ii = rstart, jj = 0, max_row_size = 0;
621 fe_it != proc_ents.end(); ++fe_it, ++ii) {
626 "not yet implemented, don't know what to do for meshset element");
630 CHKERR m_field.
get_moab().get_adjacencies(&*fe_it, 1, adj_dim,
false,
632 dim_ents = intersect(dim_ents, all_dim_ents);
635 for (
auto e : dim_ents) {
636 auto adj_it = adj_bridge_map.find(e);
637 if (adj_it != adj_bridge_map.end()) {
639 for (
const auto idx : adj_it->adj)
640 row_adj.push_back(idx);
647 std::sort(row_adj.begin(), row_adj.end());
648 auto end = std::unique(row_adj.begin(), row_adj.end());
650 size_t row_size = std::distance(row_adj.begin(), end);
651 max_row_size = std::max(max_row_size, row_size);
652 if (
j.capacity() < (nb_loc_elements - jj) * max_row_size)
653 j.reserve(nb_loc_elements * max_row_size);
656 auto diag = ents.index(*fe_it);
657 for (
auto it = row_adj.begin(); it != end; ++it)
664 if (th_vertex_weights != NULL)
665 weight_ents.push_back(*fe_it);
671 CHKERR PetscMalloc(
i.size() *
sizeof(
int), &_i);
672 CHKERR PetscMalloc(
j.size() *
sizeof(
int), &_j);
673 copy(
i.begin(),
i.end(), _i);
674 copy(
j.begin(),
j.end(), _j);
678 int *vertex_weights = NULL;
679 if (th_vertex_weights != NULL) {
680 CHKERR PetscMalloc(weight_ents.size() *
sizeof(
int), &vertex_weights);
682 &*weight_ents.begin(),
683 weight_ents.size(), vertex_weights);
689 CHKERR MatCreateMPIAdj(m_field.
get_comm(), rend - rstart, nb_elems, _i, _j,
691 CHKERR MatSetOption(Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
695 CHKERR MatConvert(Adj, MATMPIAIJ, MAT_INITIAL_MATRIX, &
A);
696 CHKERR MatView(
A, PETSC_VIEWER_DRAW_WORLD);
703 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Start";
705 MatPartitioning part;
709 CHKERR MatPartitioningSetAdjacency(part, Adj);
710 CHKERR MatPartitioningSetFromOptions(part);
711 CHKERR MatPartitioningSetNParts(part, n_parts);
712 if (th_vertex_weights != NULL) {
713 CHKERR MatPartitioningSetVertexWeights(part, vertex_weights);
716 PetscObjectTypeCompare((PetscObject)part, MATPARTITIONINGPARMETIS, &same);
719 CHKERR MatPartitioningApply_Parmetis_MoFEM(part, &is);
722 CHKERR MatPartitioningApply(part, &is);
725 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"End";
728 IS is_gather, is_num, is_gather_num;
729 CHKERR ISAllGather(is, &is_gather);
730 CHKERR ISPartitioningToNumbering(is, &is_num);
731 CHKERR ISAllGather(is_num, &is_gather_num);
733 const int *part_number, *gids;
734 CHKERR ISGetIndices(is_gather, &part_number);
735 CHKERR ISGetIndices(is_gather_num, &gids);
738 ParallelComm *pcomm = ParallelComm::get_pcomm(
740 Tag part_tag = pcomm->part_tag();
741 CHKERR m_field.
get_moab().tag_set_data(part_tag, ents, part_number);
742 Tag gid_tag = m_field.
get_moab().globalId_tag();
744 std::map<int, Range> parts_ents;
747 Range::iterator eit = ents.begin();
748 for (
int ii = 0; eit != ents.end(); eit++, ii++) {
749 parts_ents[part_number[ii]].insert(*eit);
753 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
754 moab::Interface::UNION);
755 if (!tagged_sets.empty())
758 if (n_parts > (
int)tagged_sets.size()) {
760 int num_new = n_parts - tagged_sets.size();
761 for (
int i = 0;
i < num_new;
i++) {
762 EntityHandle new_set;
764 tagged_sets.insert(new_set);
766 }
else if (n_parts < (
int)tagged_sets.size()) {
768 int num_del = tagged_sets.size() - n_parts;
769 for (
int i = 0;
i < num_del;
i++) {
770 EntityHandle old_set = tagged_sets.pop_back();
776 std::vector<int> dum_ids(n_parts);
777 for (
int i = 0;
i < n_parts;
i++)
784 for (
int pp = 0; pp != n_parts; pp++) {
785 Range dim_ents = parts_ents[pp].subset_by_dimension(
dim);
790 dim_ents,
dd,
false, adj_ents, moab::Interface::UNION);
795 parts_ents[pp].merge(adj_ents);
798 for (
int pp = 1; pp != n_parts; pp++) {
799 for (
int ppp = 0; ppp != pp; ppp++) {
800 parts_ents[pp] = subtract(parts_ents[pp], parts_ents[ppp]);
805 for (
int rr = 0; rr != n_parts; rr++) {
807 ss <<
"out_part_" << rr <<
".vtk";
809 <<
"Save debug part mesh " << ss.str();
810 EntityHandle meshset;
819 for (
int pp = 0; pp != n_parts; pp++) {
820 CHKERR m_field.
get_moab().add_entities(tagged_sets[pp], parts_ents[pp]);
830 for (
int pp = 0; pp != n_parts; pp++) {
831 Range type_ents = parts_ents[pp].subset_by_type(
t);
833 CHKERR m_field.
get_moab().tag_clear_data(part_tag, type_ents, &pp);
836 auto eit = type_ents.begin();
837 for (; eit != type_ents.end();) {
838 CHKERR m_field.
get_moab().tag_iterate(gid_tag, eit, type_ents.end(),
840 auto gid_tag_ptr =
static_cast<int *
>(ptr);
841 for (; count > 0; --count) {
852 CHKERR ISRestoreIndices(is_gather, &part_number);
853 CHKERR ISRestoreIndices(is_gather_num, &gids);
854 CHKERR ISDestroy(&is_num);
855 CHKERR ISDestroy(&is_gather_num);
856 CHKERR ISDestroy(&is_gather);
858 CHKERR MatPartitioningDestroy(&part);
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#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 ...
@ 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 ...
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 const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
MoFEMErrorCode synchroniseFieldEntities(const std::string name, int verb=DEFAULT_VERBOSITY)
virtual Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
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 resolveSharedFiniteElements(const Problem *problem_ptr, const std::string &fe_name, int verb=DEFAULT_VERBOSITY)
resolve shared entities for finite elements in the problem
#define _IT_NUMEREDFE_BY_NAME_FOR_LOOP_(PROBLEMPTR, NAME, IT)
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
const double r
rate factor
constexpr double t
plate stiffness
constexpr auto field_name
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
~CommInterface()
Destructor.
CommInterface(const MoFEM::Core &core)
MoFEMErrorCode makeEntitiesMultishared(const EntityHandle *entities, const int num_entities, const int owner_proc=0, int verb=DEFAULT_VERBOSITY)
make entities from proc 0 shared on all proc
MoFEMErrorCode synchroniseEntities(Range &ent, int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode makeFieldEntitiesMultishared(const std::string field_name, const int owner_proc=0, int verb=DEFAULT_VERBOSITY)
make field entities multi shared
MoFEMErrorCode exchangeFieldData(const std::string field_name, int verb=DEFAULT_VERBOSITY)
Exchange field data.
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 EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
Deprecated interface functions.
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
keeps basic data about problem
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.
base class for all interface classes