1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "../config/config.hpp"
13
14 #ifdef MFEM_USE_MPI
15
16 #include "mesh_headers.hpp"
17 #include "pncmesh.hpp"
18 #include "../general/binaryio.hpp"
19
20 #include <map>
21 #include <climits> // INT_MIN, INT_MAX
22
23 namespace mfem
24 {
25
26 using namespace bin_io;
27
ParNCMesh(MPI_Comm comm,const NCMesh & ncmesh,int * part)28 ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh, int *part)
29 : NCMesh(ncmesh)
30 {
31 MyComm = comm;
32 MPI_Comm_size(MyComm, &NRanks);
33 MPI_Comm_rank(MyComm, &MyRank);
34
35 // assign leaf elements to the processors by simply splitting the
36 // sequence of leaf elements into 'NRanks' parts
37 for (int i = 0; i < leaf_elements.Size(); i++)
38 {
39 elements[leaf_elements[i]].rank = part ? part[i] : InitialPartition(i);
40 }
41
42 Update();
43
44 // note that at this point all processors still have all the leaf elements;
45 // we however may now start pruning the refinement tree to get rid of
46 // branches that only contain someone else's leaves (see Prune())
47 }
48
ParNCMesh(MPI_Comm comm,std::istream & input,int version,int & curved,int & is_nc)49 ParNCMesh::ParNCMesh(MPI_Comm comm, std::istream &input, int version,
50 int &curved, int &is_nc)
51 : NCMesh(input, version, curved, is_nc)
52 {
53 MyComm = comm;
54 MPI_Comm_size(MyComm, &NRanks);
55
56 int my_rank;
57 MPI_Comm_rank(MyComm, &my_rank);
58
59 int max_rank = 0;
60 for (int i = 0; i < leaf_elements.Size(); i++)
61 {
62 max_rank = std::max(elements[leaf_elements[i]].rank, max_rank);
63 }
64
65 MFEM_VERIFY((my_rank == MyRank) && (max_rank < NRanks),
66 "Parallel mesh file doesn't seem to match current MPI setup. "
67 "Loading a parallel NC mesh with a non-matching communicator "
68 "size is not supported.");
69
70 bool iso = Iso;
71 MPI_Allreduce(&iso, &Iso, 1, MPI_C_BOOL, MPI_LAND, MyComm);
72
73 Update();
74 }
75
ParNCMesh(const ParNCMesh & other)76 ParNCMesh::ParNCMesh(const ParNCMesh &other)
77 // copy primary data only
78 : NCMesh(other)
79 , MyComm(other.MyComm)
80 , NRanks(other.NRanks)
81 {
82 Update(); // mark all secondary stuff for recalculation
83 }
84
~ParNCMesh()85 ParNCMesh::~ParNCMesh()
86 {
87 ClearAuxPM();
88 }
89
Update()90 void ParNCMesh::Update()
91 {
92 NCMesh::Update();
93
94 groups.clear();
95 group_id.clear();
96
97 CommGroup self;
98 self.push_back(MyRank);
99 groups.push_back(self);
100 group_id[self] = 0;
101
102 for (int i = 0; i < 3; i++)
103 {
104 entity_owner[i].DeleteAll();
105 entity_pmat_group[i].DeleteAll();
106 entity_index_rank[i].DeleteAll();
107 }
108
109 shared_vertices.Clear();
110 shared_edges.Clear();
111 shared_faces.Clear();
112
113 element_type.SetSize(0);
114 ghost_layer.SetSize(0);
115 boundary_layer.SetSize(0);
116 }
117
ElementSharesFace(int elem,int local,int face)118 void ParNCMesh::ElementSharesFace(int elem, int local, int face)
119 {
120 // Analogous to ElementSharesEdge.
121
122 Element &el = elements[elem];
123 int f_index = faces[face].index;
124
125 int &owner = tmp_owner[f_index];
126 owner = std::min(owner, el.rank);
127
128 char &flag = tmp_shared_flag[f_index];
129 flag |= (el.rank == MyRank) ? 0x1 : 0x2;
130
131 entity_index_rank[2].Append(Connection(f_index, el.rank));
132
133 // derive globally consistent face ID from the global element sequence
134 int &el_loc = entity_elem_local[2][f_index];
135 if (el_loc < 0 || leaf_sfc_index[el.index] < leaf_sfc_index[(el_loc >> 4)])
136 {
137 el_loc = (el.index << 4) | local;
138 }
139 }
140
BuildFaceList()141 void ParNCMesh::BuildFaceList()
142 {
143 if (HaveTets()) { GetEdgeList(); } // needed by TraverseTetEdge()
144
145 // This is an extension of NCMesh::BuildFaceList() which also determines
146 // face ownership and prepares face processor groups.
147
148 face_list.Clear();
149 shared_faces.Clear();
150 boundary_faces.SetSize(0);
151
152 if (Dim < 3 || !leaf_elements.Size()) { return; }
153
154 int nfaces = NFaces + NGhostFaces;
155
156 tmp_owner.SetSize(nfaces);
157 tmp_owner = INT_MAX;
158
159 tmp_shared_flag.SetSize(nfaces);
160 tmp_shared_flag = 0;
161
162 entity_index_rank[2].SetSize(6*leaf_elements.Size() * 3/2);
163 entity_index_rank[2].SetSize(0);
164
165 entity_elem_local[2].SetSize(nfaces);
166 entity_elem_local[2] = -1;
167
168 NCMesh::BuildFaceList();
169
170 InitOwners(nfaces, entity_owner[2]);
171 MakeSharedList(face_list, shared_faces);
172
173 tmp_owner.DeleteAll();
174 tmp_shared_flag.DeleteAll();
175
176 // create simple conforming (cut-mesh) groups now
177 CreateGroups(NFaces, entity_index_rank[2], entity_conf_group[2]);
178 // NOTE: entity_index_rank[2] is not deleted until CalculatePMatrixGroups
179
180 CalcFaceOrientations();
181 }
182
ElementSharesEdge(int elem,int local,int enode)183 void ParNCMesh::ElementSharesEdge(int elem, int local, int enode)
184 {
185 // Called by NCMesh::BuildEdgeList when an edge is visited in a leaf element.
186 // This allows us to determine edge ownership and whether it is shared
187 // without duplicating all the HashTable lookups in NCMesh::BuildEdgeList().
188
189 Element &el= elements[elem];
190 int e_index = nodes[enode].edge_index;
191
192 int &owner = tmp_owner[e_index];
193 owner = std::min(owner, el.rank);
194
195 char &flag = tmp_shared_flag[e_index];
196 flag |= (el.rank == MyRank) ? 0x1 : 0x2;
197
198 entity_index_rank[1].Append(Connection(e_index, el.rank));
199
200 // derive globally consistent edge ID from the global element sequence
201 int &el_loc = entity_elem_local[1][e_index];
202 if (el_loc < 0 || leaf_sfc_index[el.index] < leaf_sfc_index[(el_loc >> 4)])
203 {
204 el_loc = (el.index << 4) | local;
205 }
206 }
207
BuildEdgeList()208 void ParNCMesh::BuildEdgeList()
209 {
210 // This is an extension of NCMesh::BuildEdgeList() which also determines
211 // edge ownership and prepares edge processor groups.
212
213 edge_list.Clear();
214 shared_edges.Clear();
215 if (Dim < 3) { boundary_faces.SetSize(0); }
216
217 if (Dim < 2 || !leaf_elements.Size()) { return; }
218
219 int nedges = NEdges + NGhostEdges;
220
221 tmp_owner.SetSize(nedges);
222 tmp_owner = INT_MAX;
223
224 tmp_shared_flag.SetSize(nedges);
225 tmp_shared_flag = 0;
226
227 entity_index_rank[1].SetSize(12*leaf_elements.Size() * 3/2);
228 entity_index_rank[1].SetSize(0);
229
230 entity_elem_local[1].SetSize(nedges);
231 entity_elem_local[1] = -1;
232
233 NCMesh::BuildEdgeList();
234
235 InitOwners(nedges, entity_owner[1]);
236 MakeSharedList(edge_list, shared_edges);
237
238 tmp_owner.DeleteAll();
239 tmp_shared_flag.DeleteAll();
240
241 // create simple conforming (cut-mesh) groups now
242 CreateGroups(NEdges, entity_index_rank[1], entity_conf_group[1]);
243 // NOTE: entity_index_rank[1] is not deleted until CalculatePMatrixGroups
244 }
245
ElementSharesVertex(int elem,int local,int vnode)246 void ParNCMesh::ElementSharesVertex(int elem, int local, int vnode)
247 {
248 // Analogous to ElementSharesEdge.
249
250 Element &el = elements[elem];
251 int v_index = nodes[vnode].vert_index;
252
253 int &owner = tmp_owner[v_index];
254 owner = std::min(owner, el.rank);
255
256 char &flag = tmp_shared_flag[v_index];
257 flag |= (el.rank == MyRank) ? 0x1 : 0x2;
258
259 entity_index_rank[0].Append(Connection(v_index, el.rank));
260
261 // derive globally consistent vertex ID from the global element sequence
262 int &el_loc = entity_elem_local[0][v_index];
263 if (el_loc < 0 || leaf_sfc_index[el.index] < leaf_sfc_index[(el_loc >> 4)])
264 {
265 el_loc = (el.index << 4) | local;
266 }
267 }
268
BuildVertexList()269 void ParNCMesh::BuildVertexList()
270 {
271 // This is an extension of NCMesh::BuildVertexList() which also determines
272 // vertex ownership and creates vertex processor groups.
273
274 int nvertices = NVertices + NGhostVertices;
275
276 tmp_owner.SetSize(nvertices);
277 tmp_owner = INT_MAX;
278
279 tmp_shared_flag.SetSize(nvertices);
280 tmp_shared_flag = 0;
281
282 entity_index_rank[0].SetSize(8*leaf_elements.Size());
283 entity_index_rank[0].SetSize(0);
284
285 entity_elem_local[0].SetSize(nvertices);
286 entity_elem_local[0] = -1;
287
288 NCMesh::BuildVertexList();
289
290 InitOwners(nvertices, entity_owner[0]);
291 MakeSharedList(vertex_list, shared_vertices);
292
293 tmp_owner.DeleteAll();
294 tmp_shared_flag.DeleteAll();
295
296 // create simple conforming (cut-mesh) groups now
297 CreateGroups(NVertices, entity_index_rank[0], entity_conf_group[0]);
298 // NOTE: entity_index_rank[0] is not deleted until CalculatePMatrixGroups
299 }
300
InitOwners(int num,Array<GroupId> & entity_owner)301 void ParNCMesh::InitOwners(int num, Array<GroupId> &entity_owner)
302 {
303 entity_owner.SetSize(num);
304 for (int i = 0; i < num; i++)
305 {
306 entity_owner[i] =
307 (tmp_owner[i] != INT_MAX) ? GetSingletonGroup(tmp_owner[i]) : 0;
308 }
309 }
310
MakeSharedList(const NCList & list,NCList & shared)311 void ParNCMesh::MakeSharedList(const NCList &list, NCList &shared)
312 {
313 MFEM_VERIFY(tmp_shared_flag.Size(), "wrong code path");
314
315 // combine flags of masters and slaves
316 for (int i = 0; i < list.masters.Size(); i++)
317 {
318 const Master &master = list.masters[i];
319 char &master_flag = tmp_shared_flag[master.index];
320 char master_old_flag = master_flag;
321
322 for (int j = master.slaves_begin; j < master.slaves_end; j++)
323 {
324 int si = list.slaves[j].index;
325 if (si >= 0)
326 {
327 char &slave_flag = tmp_shared_flag[si];
328 master_flag |= slave_flag;
329 slave_flag |= master_old_flag;
330 }
331 else // special case: prism edge-face constraint
332 {
333 if (entity_owner[1][-1-si] != MyRank)
334 {
335 master_flag |= 0x2;
336 }
337 }
338 }
339 }
340
341 shared.Clear();
342
343 for (int i = 0; i < list.conforming.Size(); i++)
344 {
345 if (tmp_shared_flag[list.conforming[i].index] == 0x3)
346 {
347 shared.conforming.Append(list.conforming[i]);
348 }
349 }
350 for (int i = 0; i < list.masters.Size(); i++)
351 {
352 if (tmp_shared_flag[list.masters[i].index] == 0x3)
353 {
354 shared.masters.Append(list.masters[i]);
355 }
356 }
357 for (int i = 0; i < list.slaves.Size(); i++)
358 {
359 int si = list.slaves[i].index;
360 if (si >= 0 && tmp_shared_flag[si] == 0x3)
361 {
362 shared.slaves.Append(list.slaves[i]);
363 }
364 }
365 }
366
operator <(const ParNCMesh::CommGroup & lhs,const ParNCMesh::CommGroup & rhs)367 bool operator<(const ParNCMesh::CommGroup &lhs, const ParNCMesh::CommGroup &rhs)
368 {
369 if (lhs.size() == rhs.size())
370 {
371 for (unsigned i = 0; i < lhs.size(); i++)
372 {
373 if (lhs[i] < rhs[i]) { return true; }
374 }
375 return false;
376 }
377 return lhs.size() < rhs.size();
378 }
379
380 #ifdef MFEM_DEBUG
group_sorted(const ParNCMesh::CommGroup & group)381 static bool group_sorted(const ParNCMesh::CommGroup &group)
382 {
383 for (unsigned i = 1; i < group.size(); i++)
384 {
385 if (group[i] <= group[i-1]) { return false; }
386 }
387 return true;
388 }
389 #endif
390
GetGroupId(const CommGroup & group)391 ParNCMesh::GroupId ParNCMesh::GetGroupId(const CommGroup &group)
392 {
393 if (group.size() == 1 && group[0] == MyRank)
394 {
395 return 0;
396 }
397 MFEM_ASSERT(group_sorted(group), "invalid group");
398 GroupId &id = group_id[group];
399 if (!id)
400 {
401 id = groups.size();
402 groups.push_back(group);
403 }
404 return id;
405 }
406
GetSingletonGroup(int rank)407 ParNCMesh::GroupId ParNCMesh::GetSingletonGroup(int rank)
408 {
409 MFEM_ASSERT(rank != INT_MAX, "invalid rank");
410 static std::vector<int> group;
411 group.resize(1);
412 group[0] = rank;
413 return GetGroupId(group);
414 }
415
GroupContains(GroupId id,int rank) const416 bool ParNCMesh::GroupContains(GroupId id, int rank) const
417 {
418 // TODO: would std::lower_bound() pay off here? Groups are usually small.
419 const CommGroup &group = groups[id];
420 for (unsigned i = 0; i < group.size(); i++)
421 {
422 if (group[i] == rank) { return true; }
423 }
424 return false;
425 }
426
CreateGroups(int nentities,Array<Connection> & index_rank,Array<GroupId> & entity_group)427 void ParNCMesh::CreateGroups(int nentities, Array<Connection> &index_rank,
428 Array<GroupId> &entity_group)
429 {
430 index_rank.Sort();
431 index_rank.Unique();
432
433 entity_group.SetSize(nentities);
434 entity_group = 0;
435
436 CommGroup group;
437 group.reserve(128);
438
439 int begin = 0, end = 0;
440 while (begin < index_rank.Size())
441 {
442 int index = index_rank[begin].from;
443 if (index >= nentities)
444 {
445 break; // probably creating entity_conf_group (no ghosts)
446 }
447 while (end < index_rank.Size() && index_rank[end].from == index)
448 {
449 end++;
450 }
451 group.resize(end - begin);
452 for (int i = begin; i < end; i++)
453 {
454 group[i - begin] = index_rank[i].to;
455 }
456 entity_group[index] = GetGroupId(group);
457 begin = end;
458 }
459 }
460
AddConnections(int entity,int index,const Array<int> & ranks)461 void ParNCMesh::AddConnections(int entity, int index, const Array<int> &ranks)
462 {
463 for (int i = 0; i < ranks.Size(); i++)
464 {
465 entity_index_rank[entity].Append(Connection(index, ranks[i]));
466 }
467 }
468
CalculatePMatrixGroups()469 void ParNCMesh::CalculatePMatrixGroups()
470 {
471 // make sure all entity_index_rank[i] arrays are filled
472 GetSharedVertices();
473 GetSharedEdges();
474 GetSharedFaces();
475
476 int v[4], e[4], eo[4];
477
478 Array<int> ranks;
479 ranks.Reserve(256);
480
481 // connect slave edges to master edges and their vertices
482 for (int i = 0; i < shared_edges.masters.Size(); i++)
483 {
484 const Master &master_edge = shared_edges.masters[i];
485 ranks.SetSize(0);
486 for (int j = master_edge.slaves_begin; j < master_edge.slaves_end; j++)
487 {
488 int owner = entity_owner[1][edge_list.slaves[j].index];
489 ranks.Append(groups[owner][0]);
490 }
491 ranks.Sort();
492 ranks.Unique();
493
494 AddConnections(1, master_edge.index, ranks);
495
496 GetEdgeVertices(master_edge, v);
497 for (int j = 0; j < 2; j++)
498 {
499 AddConnections(0, v[j], ranks);
500 }
501 }
502
503 // connect slave faces to master faces and their edges and vertices
504 for (int i = 0; i < shared_faces.masters.Size(); i++)
505 {
506 const Master &master_face = shared_faces.masters[i];
507 ranks.SetSize(0);
508 for (int j = master_face.slaves_begin; j < master_face.slaves_end; j++)
509 {
510 int si = face_list.slaves[j].index;
511 int owner = (si >= 0) ? entity_owner[2][si] // standard face dependency
512 /* */ : entity_owner[1][-1 - si]; // prism edge-face dep
513 ranks.Append(groups[owner][0]);
514 }
515 ranks.Sort();
516 ranks.Unique();
517
518 AddConnections(2, master_face.index, ranks);
519
520 int nfv = GetFaceVerticesEdges(master_face, v, e, eo);
521 for (int j = 0; j < nfv; j++)
522 {
523 AddConnections(0, v[j], ranks);
524 AddConnections(1, e[j], ranks);
525 }
526 }
527
528 int nentities[3] =
529 {
530 NVertices + NGhostVertices,
531 NEdges + NGhostEdges,
532 NFaces + NGhostFaces
533 };
534
535 // compress the index-rank arrays into group representation
536 for (int i = 0; i < 3; i++)
537 {
538 CreateGroups(nentities[i], entity_index_rank[i], entity_pmat_group[i]);
539 entity_index_rank[i].DeleteAll();
540 }
541 }
542
get_face_orientation(Face & face,Element & e1,Element & e2,int local[2])543 int ParNCMesh::get_face_orientation(Face &face, Element &e1, Element &e2,
544 int local[2])
545 {
546 // Return face orientation in e2, assuming the face has orientation 0 in e1.
547 int ids[2][4];
548 Element* e[2] = { &e1, &e2 };
549 for (int i = 0; i < 2; i++)
550 {
551 // get local face number (remember that p1, p2, p3 are not in order, and
552 // p4 is not stored)
553 int lf = find_local_face(e[i]->Geom(),
554 find_node(*e[i], face.p1),
555 find_node(*e[i], face.p2),
556 find_node(*e[i], face.p3));
557 // optional output
558 if (local) { local[i] = lf; }
559
560 // get node IDs for the face as seen from e[i]
561 const int* fv = GI[e[i]->Geom()].faces[lf];
562 for (int j = 0; j < 4; j++)
563 {
564 ids[i][j] = e[i]->node[fv[j]];
565 }
566 }
567
568 return (ids[0][3] >= 0) ? Mesh::GetQuadOrientation(ids[0], ids[1])
569 /* */ : Mesh::GetTriOrientation(ids[0], ids[1]);
570 }
571
CalcFaceOrientations()572 void ParNCMesh::CalcFaceOrientations()
573 {
574 if (Dim < 3) { return; }
575
576 // Calculate orientation of shared conforming faces.
577 // NOTE: face orientation is calculated relative to its lower rank element.
578 // Thanks to the ghost layer this can be done locally, without communication.
579
580 face_orient.SetSize(NFaces);
581 face_orient = 0;
582
583 for (auto face = faces.begin(); face != faces.end(); ++face)
584 {
585 if (face->elem[0] >= 0 && face->elem[1] >= 0 && face->index < NFaces)
586 {
587 Element *e1 = &elements[face->elem[0]];
588 Element *e2 = &elements[face->elem[1]];
589
590 if (e1->rank == e2->rank) { continue; }
591 if (e1->rank > e2->rank) { std::swap(e1, e2); }
592
593 face_orient[face->index] = get_face_orientation(*face, *e1, *e2);
594 }
595 }
596 }
597
GetBoundaryClosure(const Array<int> & bdr_attr_is_ess,Array<int> & bdr_vertices,Array<int> & bdr_edges)598 void ParNCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
599 Array<int> &bdr_vertices,
600 Array<int> &bdr_edges)
601 {
602 NCMesh::GetBoundaryClosure(bdr_attr_is_ess, bdr_vertices, bdr_edges);
603
604 int i, j;
605 // filter out ghost vertices
606 for (i = j = 0; i < bdr_vertices.Size(); i++)
607 {
608 if (bdr_vertices[i] < NVertices) { bdr_vertices[j++] = bdr_vertices[i]; }
609 }
610 bdr_vertices.SetSize(j);
611
612 // filter out ghost edges
613 for (i = j = 0; i < bdr_edges.Size(); i++)
614 {
615 if (bdr_edges[i] < NEdges) { bdr_edges[j++] = bdr_edges[i]; }
616 }
617 bdr_edges.SetSize(j);
618 }
619
620
621 //// Neighbors /////////////////////////////////////////////////////////////////
622
UpdateLayers()623 void ParNCMesh::UpdateLayers()
624 {
625 if (element_type.Size()) { return; }
626
627 int nleaves = leaf_elements.Size();
628
629 element_type.SetSize(nleaves);
630 for (int i = 0; i < nleaves; i++)
631 {
632 element_type[i] = (elements[leaf_elements[i]].rank == MyRank) ? 1 : 0;
633 }
634
635 // determine the ghost layer
636 Array<char> ghost_set;
637 FindSetNeighbors(element_type, NULL, &ghost_set);
638
639 // find the neighbors of the ghost layer
640 Array<char> boundary_set;
641 FindSetNeighbors(ghost_set, NULL, &boundary_set);
642
643 ghost_layer.SetSize(0);
644 boundary_layer.SetSize(0);
645 for (int i = 0; i < nleaves; i++)
646 {
647 char &etype = element_type[i];
648 if (ghost_set[i])
649 {
650 etype = 2;
651 ghost_layer.Append(leaf_elements[i]);
652 }
653 else if (boundary_set[i] && etype)
654 {
655 etype = 3;
656 boundary_layer.Append(leaf_elements[i]);
657 }
658 }
659 }
660
CheckElementType(int elem,int type)661 bool ParNCMesh::CheckElementType(int elem, int type)
662 {
663 Element &el = elements[elem];
664 if (!el.ref_type)
665 {
666 return (element_type[el.index] == type);
667 }
668 else
669 {
670 for (int i = 0; i < 8 && el.child[i] >= 0; i++)
671 {
672 if (!CheckElementType(el.child[i], type)) { return false; }
673 }
674 return true;
675 }
676 }
677
ElementNeighborProcessors(int elem,Array<int> & ranks)678 void ParNCMesh::ElementNeighborProcessors(int elem, Array<int> &ranks)
679 {
680 ranks.SetSize(0); // preserve capacity
681
682 // big shortcut: there are no neighbors if element_type == 1
683 if (CheckElementType(elem, 1)) { return; }
684
685 // ok, we do need to look for neighbors;
686 // at least we can only search in the ghost layer
687 tmp_neighbors.SetSize(0);
688 FindNeighbors(elem, tmp_neighbors, &ghost_layer);
689
690 // return a list of processors
691 for (int i = 0; i < tmp_neighbors.Size(); i++)
692 {
693 ranks.Append(elements[tmp_neighbors[i]].rank);
694 }
695 ranks.Sort();
696 ranks.Unique();
697 }
698
699 template<class T>
set_to_array(const std::set<T> & set,Array<T> & array)700 static void set_to_array(const std::set<T> &set, Array<T> &array)
701 {
702 array.Reserve(set.size());
703 array.SetSize(0);
704 for (std::set<int>::iterator it = set.begin(); it != set.end(); ++it)
705 {
706 array.Append(*it);
707 }
708 }
709
NeighborProcessors(Array<int> & neighbors)710 void ParNCMesh::NeighborProcessors(Array<int> &neighbors)
711 {
712 UpdateLayers();
713
714 // TODO: look at groups instead?
715
716 std::set<int> ranks;
717 for (int i = 0; i < ghost_layer.Size(); i++)
718 {
719 ranks.insert(elements[ghost_layer[i]].rank);
720 }
721 set_to_array(ranks, neighbors);
722 }
723
724
725 //// ParMesh compatibility /////////////////////////////////////////////////////
726
MakeSharedTable(int ngroups,int ent,Array<int> & shared_local,Table & group_shared,Array<char> * entity_geom,char geom)727 void ParNCMesh::MakeSharedTable(int ngroups, int ent, Array<int> &shared_local,
728 Table &group_shared, Array<char> *entity_geom,
729 char geom)
730 {
731 const Array<GroupId> &conf_group = entity_conf_group[ent];
732
733 group_shared.MakeI(ngroups-1);
734
735 // count shared entities
736 int num_shared = 0;
737 for (int i = 0; i < conf_group.Size(); i++)
738 {
739 if (conf_group[i])
740 {
741 if (entity_geom && (*entity_geom)[i] != geom) { continue; }
742
743 num_shared++;
744 group_shared.AddAColumnInRow(conf_group[i]-1);
745 }
746 }
747
748 shared_local.SetSize(num_shared);
749 group_shared.MakeJ();
750
751 // fill shared_local and group_shared
752 for (int i = 0, j = 0; i < conf_group.Size(); i++)
753 {
754 if (conf_group[i])
755 {
756 if (entity_geom && (*entity_geom)[i] != geom) { continue; }
757
758 shared_local[j] = i;
759 group_shared.AddConnection(conf_group[i]-1, j);
760 j++;
761 }
762 }
763 group_shared.ShiftUpI();
764
765 // sort the groups consistently across processors
766 for (int i = 0; i < group_shared.Size(); i++)
767 {
768 int size = group_shared.RowSize(i);
769 int *row = group_shared.GetRow(i);
770
771 Array<int> ref_row(row, size);
772 ref_row.Sort([&](const int a, const int b)
773 {
774 int el_loc_a = entity_elem_local[ent][shared_local[a]];
775 int el_loc_b = entity_elem_local[ent][shared_local[b]];
776
777 int lsi_a = leaf_sfc_index[el_loc_a >> 4];
778 int lsi_b = leaf_sfc_index[el_loc_b >> 4];
779
780 if (lsi_a != lsi_b) { return lsi_a < lsi_b; }
781
782 return (el_loc_a & 0xf) < (el_loc_b & 0xf);
783 });
784 }
785 }
786
GetConformingSharedStructures(ParMesh & pmesh)787 void ParNCMesh::GetConformingSharedStructures(ParMesh &pmesh)
788 {
789 if (leaf_elements.Size())
790 {
791 // make sure we have entity_conf_group[x] and the ordering arrays
792 for (int ent = 0; ent < Dim; ent++)
793 {
794 GetSharedList(ent);
795 MFEM_VERIFY(entity_conf_group[ent].Size(), "internal error");
796 MFEM_VERIFY(entity_elem_local[ent].Size(), "internal error");
797 }
798 }
799
800 // create ParMesh groups, and the map (ncmesh_group -> pmesh_group)
801 Array<int> group_map(groups.size());
802 {
803 group_map = 0;
804 IntegerSet iset;
805 ListOfIntegerSets int_groups;
806 for (unsigned i = 0; i < groups.size(); i++)
807 {
808 if (groups[i].size() > 1 || !i) // skip singleton groups
809 {
810 iset.Recreate(groups[i].size(), groups[i].data());
811 group_map[i] = int_groups.Insert(iset);
812 }
813 }
814 pmesh.gtopo.Create(int_groups, 822);
815 }
816
817 // renumber groups in entity_conf_group[] (due to missing singletons)
818 for (int ent = 0; ent < 3; ent++)
819 {
820 for (int i = 0; i < entity_conf_group[ent].Size(); i++)
821 {
822 GroupId &ecg = entity_conf_group[ent][i];
823 ecg = group_map[ecg];
824 }
825 }
826
827 // create shared to local index mappings and group tables
828 int ng = pmesh.gtopo.NGroups();
829 MakeSharedTable(ng, 0, pmesh.svert_lvert, pmesh.group_svert);
830 MakeSharedTable(ng, 1, pmesh.sedge_ledge, pmesh.group_sedge);
831
832 Array<int> slt, slq;
833 MakeSharedTable(ng, 2, slt, pmesh.group_stria, &face_geom, Geometry::TRIANGLE);
834 MakeSharedTable(ng, 2, slq, pmesh.group_squad, &face_geom, Geometry::SQUARE);
835
836 pmesh.sface_lface = slt;
837 pmesh.sface_lface.Append(slq);
838
839 // create shared_edges
840 for (int i = 0; i < pmesh.shared_edges.Size(); i++)
841 {
842 delete pmesh.shared_edges[i];
843 }
844 pmesh.shared_edges.SetSize(pmesh.sedge_ledge.Size());
845 for (int i = 0; i < pmesh.shared_edges.Size(); i++)
846 {
847 int el_loc = entity_elem_local[1][pmesh.sedge_ledge[i]];
848 MeshId edge_id(-1, leaf_elements[(el_loc >> 4)], (el_loc & 0xf));
849
850 int v[2];
851 GetEdgeVertices(edge_id, v, false);
852 pmesh.shared_edges[i] = new Segment(v, 1);
853 }
854
855 // create shared_trias
856 pmesh.shared_trias.SetSize(slt.Size());
857 for (int i = 0; i < slt.Size(); i++)
858 {
859 int el_loc = entity_elem_local[2][slt[i]];
860 MeshId face_id(-1, leaf_elements[(el_loc >> 4)], (el_loc & 0xf));
861
862 int v[4], e[4], eo[4];
863 GetFaceVerticesEdges(face_id, v, e, eo);
864 pmesh.shared_trias[i].Set(v);
865 }
866
867 // create shared_quads
868 pmesh.shared_quads.SetSize(slq.Size());
869 for (int i = 0; i < slq.Size(); i++)
870 {
871 int el_loc = entity_elem_local[2][slq[i]];
872 MeshId face_id(-1, leaf_elements[(el_loc >> 4)], (el_loc & 0xf));
873
874 int e[4], eo[4];
875 GetFaceVerticesEdges(face_id, pmesh.shared_quads[i].v, e, eo);
876 }
877
878 // free the arrays, they're not needed anymore (until next mesh update)
879 for (int ent = 0; ent < Dim; ent++)
880 {
881 entity_conf_group[ent].DeleteAll();
882 entity_elem_local[ent].DeleteAll();
883 }
884 }
885
GetFaceNeighbors(ParMesh & pmesh)886 void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
887 {
888 ClearAuxPM();
889
890 const NCList &shared = (Dim == 3) ? GetSharedFaces() : GetSharedEdges();
891 const NCList &full_list = (Dim == 3) ? GetFaceList() : GetEdgeList();
892
893 Array<Element*> fnbr;
894 Array<Connection> send_elems;
895
896 int bound = shared.conforming.Size() + shared.slaves.Size();
897
898 fnbr.Reserve(bound);
899 send_elems.Reserve(bound);
900
901 // go over all shared faces and collect face neighbor elements
902 for (int i = 0; i < shared.conforming.Size(); i++)
903 {
904 const MeshId &cf = shared.conforming[i];
905 Face* face = GetFace(elements[cf.element], cf.local);
906 MFEM_ASSERT(face != NULL, "");
907
908 MFEM_ASSERT(face->elem[0] >= 0 && face->elem[1] >= 0, "");
909 Element* e[2] = { &elements[face->elem[0]], &elements[face->elem[1]] };
910
911 if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
912 MFEM_ASSERT(e[0]->rank != MyRank && e[1]->rank == MyRank, "");
913
914 fnbr.Append(e[0]);
915 send_elems.Append(Connection(e[0]->rank, e[1]->index));
916 }
917
918 for (int i = 0; i < shared.masters.Size(); i++)
919 {
920 const Master &mf = shared.masters[i];
921 for (int j = mf.slaves_begin; j < mf.slaves_end; j++)
922 {
923 const Slave &sf = full_list.slaves[j];
924 if (sf.element < 0) { continue; }
925
926 MFEM_ASSERT(mf.element >= 0, "");
927 Element* e[2] = { &elements[mf.element], &elements[sf.element] };
928
929 bool loc0 = (e[0]->rank == MyRank);
930 bool loc1 = (e[1]->rank == MyRank);
931 if (loc0 == loc1) { continue; }
932 if (loc0) { std::swap(e[0], e[1]); }
933
934 fnbr.Append(e[0]);
935 send_elems.Append(Connection(e[0]->rank, e[1]->index));
936 }
937 }
938
939 MFEM_ASSERT(fnbr.Size() <= bound, "oops, bad upper bound");
940
941 // remove duplicate face neighbor elements and sort them by rank & index
942 // (note that the send table is sorted the same way and the order is also the
943 // same on different processors, this is important for ExchangeFaceNbrData)
944 fnbr.Sort();
945 fnbr.Unique();
946 fnbr.Sort([](const Element* a, const Element* b)
947 {
948 return (a->rank != b->rank) ? a->rank < b->rank
949 /* */ : a->index < b->index;
950 });
951
952 // put the ranks into 'face_nbr_group'
953 for (int i = 0; i < fnbr.Size(); i++)
954 {
955 if (!i || fnbr[i]->rank != pmesh.face_nbr_group.Last())
956 {
957 pmesh.face_nbr_group.Append(fnbr[i]->rank);
958 }
959 }
960 int nranks = pmesh.face_nbr_group.Size();
961
962 // create a new mfem::Element for each face neighbor element
963 pmesh.face_nbr_elements.SetSize(0);
964 pmesh.face_nbr_elements.Reserve(fnbr.Size());
965
966 pmesh.face_nbr_elements_offset.SetSize(0);
967 pmesh.face_nbr_elements_offset.Reserve(pmesh.face_nbr_group.Size()+1);
968
969 Array<int> fnbr_index(NGhostElements);
970 fnbr_index = -1;
971
972 std::map<int, int> vert_map;
973 for (int i = 0; i < fnbr.Size(); i++)
974 {
975 Element* elem = fnbr[i];
976 mfem::Element* fne = NewMeshElement(elem->geom);
977 fne->SetAttribute(elem->attribute);
978 pmesh.face_nbr_elements.Append(fne);
979
980 GeomInfo& gi = GI[(int) elem->geom];
981 for (int k = 0; k < gi.nv; k++)
982 {
983 int &v = vert_map[elem->node[k]];
984 if (!v) { v = vert_map.size(); }
985 fne->GetVertices()[k] = v-1;
986 }
987
988 if (!i || elem->rank != fnbr[i-1]->rank)
989 {
990 pmesh.face_nbr_elements_offset.Append(i);
991 }
992
993 MFEM_ASSERT(elem->index >= NElements, "not a ghost element");
994 fnbr_index[elem->index - NElements] = i;
995 }
996 pmesh.face_nbr_elements_offset.Append(fnbr.Size());
997
998 // create vertices in 'face_nbr_vertices'
999 {
1000 pmesh.face_nbr_vertices.SetSize(vert_map.size());
1001 if (coordinates.Size())
1002 {
1003 tmp_vertex = new TmpVertex[nodes.NumIds()]; // TODO: something cheaper?
1004 for (auto it = vert_map.begin(); it != vert_map.end(); ++it)
1005 {
1006 pmesh.face_nbr_vertices[it->second-1].SetCoords(
1007 spaceDim, CalcVertexPos(it->first));
1008 }
1009 delete [] tmp_vertex;
1010 }
1011 }
1012
1013 // make the 'send_face_nbr_elements' table
1014 send_elems.Sort();
1015 send_elems.Unique();
1016
1017 for (int i = 0, last_rank = -1; i < send_elems.Size(); i++)
1018 {
1019 Connection &c = send_elems[i];
1020 if (c.from != last_rank)
1021 {
1022 // renumber rank to position in 'face_nbr_group'
1023 last_rank = c.from;
1024 c.from = pmesh.face_nbr_group.Find(c.from);
1025 }
1026 else
1027 {
1028 c.from = send_elems[i-1].from; // avoid search
1029 }
1030 }
1031 pmesh.send_face_nbr_elements.MakeFromList(nranks, send_elems);
1032
1033 // go over the shared faces again and modify their Mesh::FaceInfo
1034 for (int i = 0; i < shared.conforming.Size(); i++)
1035 {
1036 const MeshId &cf = shared.conforming[i];
1037 Face* face = GetFace(elements[cf.element], cf.local);
1038
1039 Element* e[2] = { &elements[face->elem[0]], &elements[face->elem[1]] };
1040 if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
1041
1042 Mesh::FaceInfo &fi = pmesh.faces_info[cf.index];
1043 fi.Elem2No = -1 - fnbr_index[e[0]->index - NElements];
1044
1045 if (Dim == 3)
1046 {
1047 int local[2];
1048 int o = get_face_orientation(*face, *e[1], *e[0], local);
1049 fi.Elem2Inf = 64*local[1] + o;
1050 }
1051 else
1052 {
1053 fi.Elem2Inf = 64*find_element_edge(*e[0], face->p1, face->p3) + 1;
1054 }
1055 }
1056
1057 if (shared.slaves.Size())
1058 {
1059 int nfaces = NFaces, nghosts = NGhostFaces;
1060 if (Dim <= 2) { nfaces = NEdges, nghosts = NGhostEdges; }
1061
1062 // enlarge Mesh::faces_info for ghost slaves
1063 MFEM_ASSERT(pmesh.faces_info.Size() == nfaces, "");
1064 MFEM_ASSERT(pmesh.GetNumFaces() == nfaces, "");
1065 pmesh.faces_info.SetSize(nfaces + nghosts);
1066 for (int i = nfaces; i < pmesh.faces_info.Size(); i++)
1067 {
1068 Mesh::FaceInfo &fi = pmesh.faces_info[i];
1069 fi.Elem1No = fi.Elem2No = -1;
1070 fi.Elem1Inf = fi.Elem2Inf = -1;
1071 fi.NCFace = -1;
1072 }
1073 // Note that some of the indices i >= nfaces in pmesh.faces_info will
1074 // remain untouched below and they will have Elem1No == -1, in particular.
1075
1076 // fill in FaceInfo for shared slave faces
1077 for (int i = 0; i < shared.masters.Size(); i++)
1078 {
1079 const Master &mf = shared.masters[i];
1080 for (int j = mf.slaves_begin; j < mf.slaves_end; j++)
1081 {
1082 const Slave &sf = full_list.slaves[j];
1083 if (sf.element < 0) { continue; }
1084
1085 MFEM_ASSERT(mf.element >= 0, "");
1086 Element &sfe = elements[sf.element];
1087 Element &mfe = elements[mf.element];
1088
1089 bool sloc = (sfe.rank == MyRank);
1090 bool mloc = (mfe.rank == MyRank);
1091 if (sloc == mloc) { continue; }
1092
1093 Mesh::FaceInfo &fi = pmesh.faces_info[sf.index];
1094 fi.Elem1No = sfe.index;
1095 fi.Elem2No = mfe.index;
1096 fi.Elem1Inf = 64 * sf.local;
1097 fi.Elem2Inf = 64 * mf.local;
1098
1099 if (!sloc)
1100 {
1101 // 'fi' is the info for a ghost slave face with index:
1102 // sf.index >= nfaces
1103 std::swap(fi.Elem1No, fi.Elem2No);
1104 std::swap(fi.Elem1Inf, fi.Elem2Inf);
1105 // After the above swap, Elem1No refers to the local, master-side
1106 // element. In other words, side 1 IS NOT the side that generated
1107 // the face.
1108 }
1109 else
1110 {
1111 // 'fi' is the info for a local slave face with index:
1112 // sf.index < nfaces
1113 // Here, Elem1No refers to the local, slave-side element.
1114 // In other words, side 1 IS the side that generated the face.
1115 }
1116 MFEM_ASSERT(fi.Elem2No >= NElements, "");
1117 fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
1118
1119 const DenseMatrix* pm = full_list.point_matrices[sf.geom][sf.matrix];
1120 if (!sloc && Dim == 3)
1121 {
1122 // TODO: does this handle triangle faces correctly?
1123
1124 // ghost slave in 3D needs flipping orientation
1125 DenseMatrix* pm2 = new DenseMatrix(*pm);
1126 std::swap((*pm2)(0,1), (*pm2)(0,3));
1127 std::swap((*pm2)(1,1), (*pm2)(1,3));
1128 aux_pm_store.Append(pm2);
1129
1130 fi.Elem2Inf ^= 1;
1131 pm = pm2;
1132
1133 // The problem is that sf.point_matrix is designed for P matrix
1134 // construction and always has orientation relative to the slave
1135 // face. In ParMesh::GetSharedFaceTransformations the result
1136 // would therefore be the same on both processors, which is not
1137 // how that function works for conforming faces. The orientation
1138 // of Loc1, Loc2 and Face needs to always be relative to Element
1139 // 1, which is the element containing the slave face on one
1140 // processor, but on the other it is the element containing the
1141 // master face. In the latter case we need to flip the pm.
1142 }
1143 else if (!sloc && Dim == 2)
1144 {
1145 fi.Elem2Inf ^= 1; // set orientation to 1
1146 // The point matrix (used to define "side 1" which is the same as
1147 // "parent side" in this case) does not require a flip since it
1148 // is aligned with the parent side, so NO flip is performed in
1149 // Mesh::ApplyLocalSlaveTransformation.
1150 }
1151
1152 MFEM_ASSERT(fi.NCFace < 0, "");
1153 fi.NCFace = pmesh.nc_faces_info.Size();
1154 pmesh.nc_faces_info.Append(Mesh::NCFaceInfo(true, sf.master, pm));
1155 }
1156 }
1157 }
1158
1159 // NOTE: this function skips ParMesh::send_face_nbr_vertices and
1160 // ParMesh::face_nbr_vertices_offset, these are not used outside of ParMesh
1161 }
1162
ClearAuxPM()1163 void ParNCMesh::ClearAuxPM()
1164 {
1165 for (int i = 0; i < aux_pm_store.Size(); i++)
1166 {
1167 delete aux_pm_store[i];
1168 }
1169 aux_pm_store.DeleteAll();
1170 }
1171
1172
1173 //// Prune, Refine, Derefine ///////////////////////////////////////////////////
1174
PruneTree(int elem)1175 bool ParNCMesh::PruneTree(int elem)
1176 {
1177 Element &el = elements[elem];
1178 if (el.ref_type)
1179 {
1180 bool remove[8];
1181 bool removeAll = true;
1182
1183 // determine which subtrees can be removed (and whether it's all of them)
1184 for (int i = 0; i < 8; i++)
1185 {
1186 remove[i] = false;
1187 if (el.child[i] >= 0)
1188 {
1189 remove[i] = PruneTree(el.child[i]);
1190 if (!remove[i]) { removeAll = false; }
1191 }
1192 }
1193
1194 // all children can be removed, let the (maybe indirect) parent do it
1195 if (removeAll) { return true; }
1196
1197 // not all children can be removed, but remove those that can be
1198 for (int i = 0; i < 8; i++)
1199 {
1200 if (remove[i]) { DerefineElement(el.child[i]); }
1201 }
1202
1203 return false; // need to keep this element and up
1204 }
1205 else
1206 {
1207 // return true if this leaf can be removed
1208 return el.rank < 0;
1209 }
1210 }
1211
Prune()1212 void ParNCMesh::Prune()
1213 {
1214 if (!Iso && Dim == 3)
1215 {
1216 if (MyRank == 0)
1217 {
1218 MFEM_WARNING("Can't prune 3D aniso meshes yet.");
1219 }
1220 return;
1221 }
1222
1223 UpdateLayers();
1224
1225 for (int i = 0; i < leaf_elements.Size(); i++)
1226 {
1227 // rank of elements beyond the ghost layer is unknown / not updated
1228 if (element_type[i] == 0)
1229 {
1230 elements[leaf_elements[i]].rank = -1;
1231 // NOTE: rank == -1 will make the element disappear from leaf_elements
1232 // on next Update, see NCMesh::CollectLeafElements
1233 }
1234 }
1235
1236 // derefine subtrees whose leaves are all unneeded
1237 for (int i = 0; i < root_state.Size(); i++)
1238 {
1239 if (PruneTree(i)) { DerefineElement(i); }
1240 }
1241
1242 Update();
1243 }
1244
1245
Refine(const Array<Refinement> & refinements)1246 void ParNCMesh::Refine(const Array<Refinement> &refinements)
1247 {
1248 if (NRanks == 1)
1249 {
1250 NCMesh::Refine(refinements);
1251 return;
1252 }
1253
1254 for (int i = 0; i < refinements.Size(); i++)
1255 {
1256 const Refinement &ref = refinements[i];
1257 MFEM_VERIFY(ref.ref_type == 7 || Dim < 3,
1258 "anisotropic parallel refinement not supported yet in 3D.");
1259 }
1260 MFEM_VERIFY(Iso || Dim < 3,
1261 "parallel refinement of 3D aniso meshes not supported yet.");
1262
1263 NeighborRefinementMessage::Map send_ref;
1264
1265 // create refinement messages to all neighbors (NOTE: some may be empty)
1266 Array<int> neighbors;
1267 NeighborProcessors(neighbors);
1268 for (int i = 0; i < neighbors.Size(); i++)
1269 {
1270 send_ref[neighbors[i]].SetNCMesh(this);
1271 }
1272
1273 // populate messages: all refinements that occur next to the processor
1274 // boundary need to be sent to the adjoining neighbors so they can keep
1275 // their ghost layer up to date
1276 Array<int> ranks;
1277 ranks.Reserve(64);
1278 for (int i = 0; i < refinements.Size(); i++)
1279 {
1280 const Refinement &ref = refinements[i];
1281 MFEM_ASSERT(ref.index < NElements, "");
1282 int elem = leaf_elements[ref.index];
1283 ElementNeighborProcessors(elem, ranks);
1284 for (int j = 0; j < ranks.Size(); j++)
1285 {
1286 send_ref[ranks[j]].AddRefinement(elem, ref.ref_type);
1287 }
1288 }
1289
1290 // send the messages (overlap with local refinements)
1291 NeighborRefinementMessage::IsendAll(send_ref, MyComm);
1292
1293 // do local refinements
1294 for (int i = 0; i < refinements.Size(); i++)
1295 {
1296 const Refinement &ref = refinements[i];
1297 NCMesh::RefineElement(leaf_elements[ref.index], ref.ref_type);
1298 }
1299
1300 // receive (ghost layer) refinements from all neighbors
1301 for (int j = 0; j < neighbors.Size(); j++)
1302 {
1303 int rank, size;
1304 NeighborRefinementMessage::Probe(rank, size, MyComm);
1305
1306 NeighborRefinementMessage msg;
1307 msg.SetNCMesh(this);
1308 msg.Recv(rank, size, MyComm);
1309
1310 // do the ghost refinements
1311 for (int i = 0; i < msg.Size(); i++)
1312 {
1313 NCMesh::RefineElement(msg.elements[i], msg.values[i]);
1314 }
1315 }
1316
1317 Update();
1318
1319 // make sure we can delete the send buffers
1320 NeighborRefinementMessage::WaitAllSent(send_ref);
1321 }
1322
1323
LimitNCLevel(int max_nc_level)1324 void ParNCMesh::LimitNCLevel(int max_nc_level)
1325 {
1326 MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
1327
1328 while (1)
1329 {
1330 Array<Refinement> refinements;
1331 GetLimitRefinements(refinements, max_nc_level);
1332
1333 long size = refinements.Size(), glob_size;
1334 MPI_Allreduce(&size, &glob_size, 1, MPI_LONG, MPI_SUM, MyComm);
1335
1336 if (!glob_size) { break; }
1337
1338 Refine(refinements);
1339 }
1340 }
1341
Derefine(const Array<int> & derefs)1342 void ParNCMesh::Derefine(const Array<int> &derefs)
1343 {
1344 MFEM_VERIFY(Dim < 3 || Iso,
1345 "derefinement of 3D anisotropic meshes not implemented yet.");
1346
1347 InitDerefTransforms();
1348
1349 // store fine element ranks
1350 old_index_or_rank.SetSize(leaf_elements.Size());
1351 for (int i = 0; i < leaf_elements.Size(); i++)
1352 {
1353 old_index_or_rank[i] = elements[leaf_elements[i]].rank;
1354 }
1355
1356 // back up the leaf_elements array
1357 Array<int> old_elements;
1358 leaf_elements.Copy(old_elements);
1359
1360 // *** STEP 1: redistribute elements to avoid complex derefinements ***
1361
1362 Array<int> new_ranks(leaf_elements.Size());
1363 for (int i = 0; i < leaf_elements.Size(); i++)
1364 {
1365 new_ranks[i] = elements[leaf_elements[i]].rank;
1366 }
1367
1368 // make the lowest rank get all the fine elements for each derefinement
1369 for (int i = 0; i < derefs.Size(); i++)
1370 {
1371 int row = derefs[i];
1372 MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
1373 "invalid derefinement number.");
1374
1375 const int* fine = derefinements.GetRow(row);
1376 int size = derefinements.RowSize(row);
1377
1378 int coarse_rank = INT_MAX;
1379 for (int j = 0; j < size; j++)
1380 {
1381 int fine_rank = elements[leaf_elements[fine[j]]].rank;
1382 coarse_rank = std::min(coarse_rank, fine_rank);
1383 }
1384 for (int j = 0; j < size; j++)
1385 {
1386 new_ranks[fine[j]] = coarse_rank;
1387 }
1388 }
1389
1390 int target_elements = 0;
1391 for (int i = 0; i < new_ranks.Size(); i++)
1392 {
1393 if (new_ranks[i] == MyRank) { target_elements++; }
1394 }
1395
1396 // redistribute elements slightly to get rid of complex derefinements
1397 // straddling processor boundaries *and* update the ghost layer
1398 RedistributeElements(new_ranks, target_elements, false);
1399
1400 // *** STEP 2: derefine now, communication similar to Refine() ***
1401
1402 NeighborDerefinementMessage::Map send_deref;
1403
1404 // create derefinement messages to all neighbors (NOTE: some may be empty)
1405 Array<int> neighbors;
1406 NeighborProcessors(neighbors);
1407 for (int i = 0; i < neighbors.Size(); i++)
1408 {
1409 send_deref[neighbors[i]].SetNCMesh(this);
1410 }
1411
1412 // derefinements that occur next to the processor boundary need to be sent
1413 // to the adjoining neighbors to keep their ghost layers in sync
1414 Array<int> ranks;
1415 ranks.Reserve(64);
1416 for (int i = 0; i < derefs.Size(); i++)
1417 {
1418 const int* fine = derefinements.GetRow(derefs[i]);
1419 int parent = elements[old_elements[fine[0]]].parent;
1420
1421 // send derefinement to neighbors
1422 ElementNeighborProcessors(parent, ranks);
1423 for (int j = 0; j < ranks.Size(); j++)
1424 {
1425 send_deref[ranks[j]].AddDerefinement(parent, new_ranks[fine[0]]);
1426 }
1427 }
1428 NeighborDerefinementMessage::IsendAll(send_deref, MyComm);
1429
1430 // restore old (pre-redistribution) element indices, for SetDerefMatrixCodes
1431 for (int i = 0; i < leaf_elements.Size(); i++)
1432 {
1433 elements[leaf_elements[i]].index = -1;
1434 }
1435 for (int i = 0; i < old_elements.Size(); i++)
1436 {
1437 elements[old_elements[i]].index = i;
1438 }
1439
1440 // do local derefinements
1441 Array<int> coarse;
1442 old_elements.Copy(coarse);
1443 for (int i = 0; i < derefs.Size(); i++)
1444 {
1445 const int* fine = derefinements.GetRow(derefs[i]);
1446 int parent = elements[old_elements[fine[0]]].parent;
1447
1448 // record the relation of the fine elements to their parent
1449 SetDerefMatrixCodes(parent, coarse);
1450
1451 NCMesh::DerefineElement(parent);
1452 }
1453
1454 // receive ghost layer derefinements from all neighbors
1455 for (int j = 0; j < neighbors.Size(); j++)
1456 {
1457 int rank, size;
1458 NeighborDerefinementMessage::Probe(rank, size, MyComm);
1459
1460 NeighborDerefinementMessage msg;
1461 msg.SetNCMesh(this);
1462 msg.Recv(rank, size, MyComm);
1463
1464 // do the ghost derefinements
1465 for (int i = 0; i < msg.Size(); i++)
1466 {
1467 int elem = msg.elements[i];
1468 if (elements[elem].ref_type)
1469 {
1470 SetDerefMatrixCodes(elem, coarse);
1471 NCMesh::DerefineElement(elem);
1472 }
1473 elements[elem].rank = msg.values[i];
1474 }
1475 }
1476
1477 // update leaf_elements, Element::index etc.
1478 Update();
1479
1480 UpdateLayers();
1481
1482 // link old fine elements to the new coarse elements
1483 for (int i = 0; i < coarse.Size(); i++)
1484 {
1485 int index = elements[coarse[i]].index;
1486 if (element_type[index] == 0)
1487 {
1488 // this coarse element will get pruned, encode who owns it now
1489 index = -1 - elements[coarse[i]].rank;
1490 }
1491 transforms.embeddings[i].parent = index;
1492 }
1493
1494 leaf_elements.Copy(old_elements);
1495
1496 Prune();
1497
1498 // renumber coarse element indices after pruning
1499 for (int i = 0; i < coarse.Size(); i++)
1500 {
1501 int &index = transforms.embeddings[i].parent;
1502 if (index >= 0)
1503 {
1504 index = elements[old_elements[index]].index;
1505 }
1506 }
1507
1508 // make sure we can delete all send buffers
1509 NeighborDerefinementMessage::WaitAllSent(send_deref);
1510 }
1511
1512
1513 template<typename Type>
SynchronizeDerefinementData(Array<Type> & elem_data,const Table & deref_table)1514 void ParNCMesh::SynchronizeDerefinementData(Array<Type> &elem_data,
1515 const Table &deref_table)
1516 {
1517 const MPI_Datatype datatype = MPITypeMap<Type>::mpi_type;
1518
1519 Array<MPI_Request*> requests;
1520 Array<int> neigh;
1521
1522 requests.Reserve(64);
1523 neigh.Reserve(8);
1524
1525 // make room for ghost values (indices beyond NumElements)
1526 elem_data.SetSize(leaf_elements.Size(), 0);
1527
1528 for (int i = 0; i < deref_table.Size(); i++)
1529 {
1530 const int* fine = deref_table.GetRow(i);
1531 int size = deref_table.RowSize(i);
1532 MFEM_ASSERT(size <= 8, "");
1533
1534 int ranks[8], min_rank = INT_MAX, max_rank = INT_MIN;
1535 for (int j = 0; j < size; j++)
1536 {
1537 ranks[j] = elements[leaf_elements[fine[j]]].rank;
1538 min_rank = std::min(min_rank, ranks[j]);
1539 max_rank = std::max(max_rank, ranks[j]);
1540 }
1541
1542 // exchange values for derefinements that straddle processor boundaries
1543 if (min_rank != max_rank)
1544 {
1545 neigh.SetSize(0);
1546 for (int j = 0; j < size; j++)
1547 {
1548 if (ranks[j] != MyRank) { neigh.Append(ranks[j]); }
1549 }
1550 neigh.Sort();
1551 neigh.Unique();
1552
1553 for (int j = 0; j < size; j++/*pass*/)
1554 {
1555 Type *data = &elem_data[fine[j]];
1556
1557 int rnk = ranks[j], len = 1; /*j;
1558 do { j++; } while (j < size && ranks[j] == rnk);
1559 len = j - len;*/
1560
1561 if (rnk == MyRank)
1562 {
1563 for (int k = 0; k < neigh.Size(); k++)
1564 {
1565 MPI_Request* req = new MPI_Request;
1566 MPI_Isend(data, len, datatype, neigh[k], 292, MyComm, req);
1567 requests.Append(req);
1568 }
1569 }
1570 else
1571 {
1572 MPI_Request* req = new MPI_Request;
1573 MPI_Irecv(data, len, datatype, rnk, 292, MyComm, req);
1574 requests.Append(req);
1575 }
1576 }
1577 }
1578 }
1579
1580 for (int i = 0; i < requests.Size(); i++)
1581 {
1582 MPI_Wait(requests[i], MPI_STATUS_IGNORE);
1583 delete requests[i];
1584 }
1585 }
1586
1587 // instantiate SynchronizeDerefinementData for int and double
1588 template void
1589 ParNCMesh::SynchronizeDerefinementData<int>(Array<int> &, const Table &);
1590 template void
1591 ParNCMesh::SynchronizeDerefinementData<double>(Array<double> &, const Table &);
1592
1593
CheckDerefinementNCLevel(const Table & deref_table,Array<int> & level_ok,int max_nc_level)1594 void ParNCMesh::CheckDerefinementNCLevel(const Table &deref_table,
1595 Array<int> &level_ok, int max_nc_level)
1596 {
1597 Array<int> leaf_ok(leaf_elements.Size());
1598 leaf_ok = 1;
1599
1600 // check elements that we own
1601 for (int i = 0; i < deref_table.Size(); i++)
1602 {
1603 const int *fine = deref_table.GetRow(i),
1604 size = deref_table.RowSize(i);
1605
1606 int parent = elements[leaf_elements[fine[0]]].parent;
1607 Element &pa = elements[parent];
1608
1609 for (int j = 0; j < size; j++)
1610 {
1611 int child = leaf_elements[fine[j]];
1612 if (elements[child].rank == MyRank)
1613 {
1614 int splits[3];
1615 CountSplits(child, splits);
1616
1617 for (int k = 0; k < Dim; k++)
1618 {
1619 if ((pa.ref_type & (1 << k)) &&
1620 splits[k] >= max_nc_level)
1621 {
1622 leaf_ok[fine[j]] = 0; break;
1623 }
1624 }
1625 }
1626 }
1627 }
1628
1629 SynchronizeDerefinementData(leaf_ok, deref_table);
1630
1631 level_ok.SetSize(deref_table.Size());
1632 level_ok = 1;
1633
1634 for (int i = 0; i < deref_table.Size(); i++)
1635 {
1636 const int* fine = deref_table.GetRow(i),
1637 size = deref_table.RowSize(i);
1638
1639 for (int j = 0; j < size; j++)
1640 {
1641 if (!leaf_ok[fine[j]])
1642 {
1643 level_ok[i] = 0; break;
1644 }
1645 }
1646 }
1647 }
1648
1649
1650 //// Rebalance /////////////////////////////////////////////////////////////////
1651
Rebalance(const Array<int> * custom_partition)1652 void ParNCMesh::Rebalance(const Array<int> *custom_partition)
1653 {
1654 send_rebalance_dofs.clear();
1655 recv_rebalance_dofs.clear();
1656
1657 Array<int> old_elements;
1658 leaf_elements.GetSubArray(0, NElements, old_elements);
1659
1660 if (!custom_partition) // SFC based partitioning
1661 {
1662 Array<int> new_ranks(leaf_elements.Size());
1663 new_ranks = -1;
1664
1665 // figure out new assignments for Element::rank
1666 long local_elems = NElements, total_elems = 0;
1667 MPI_Allreduce(&local_elems, &total_elems, 1, MPI_LONG, MPI_SUM, MyComm);
1668
1669 long first_elem_global = 0;
1670 MPI_Scan(&local_elems, &first_elem_global, 1, MPI_LONG, MPI_SUM, MyComm);
1671 first_elem_global -= local_elems;
1672
1673 for (int i = 0, j = 0; i < leaf_elements.Size(); i++)
1674 {
1675 if (elements[leaf_elements[i]].rank == MyRank)
1676 {
1677 new_ranks[i] = Partition(first_elem_global + (j++), total_elems);
1678 }
1679 }
1680
1681 int target_elements = PartitionFirstIndex(MyRank+1, total_elems)
1682 - PartitionFirstIndex(MyRank, total_elems);
1683
1684 // assign the new ranks and send elements (plus ghosts) to new owners
1685 RedistributeElements(new_ranks, target_elements, true);
1686 }
1687 else // whatever partitioning the user has passed
1688 {
1689 MFEM_VERIFY(custom_partition->Size() == NElements,
1690 "Size of the partition array must match the number "
1691 "of local mesh elements (ParMesh::GetNE()).");
1692
1693 Array<int> new_ranks;
1694 custom_partition->Copy(new_ranks);
1695 new_ranks.SetSize(leaf_elements.Size(), -1); // make room for ghosts
1696
1697 RedistributeElements(new_ranks, -1, true);
1698 }
1699
1700 // set up the old index array
1701 old_index_or_rank.SetSize(NElements);
1702 old_index_or_rank = -1;
1703 for (int i = 0; i < old_elements.Size(); i++)
1704 {
1705 Element &el = elements[old_elements[i]];
1706 if (el.rank == MyRank) { old_index_or_rank[el.index] = i; }
1707 }
1708
1709 // get rid of elements beyond the new ghost layer
1710 Prune();
1711 }
1712
RedistributeElements(Array<int> & new_ranks,int target_elements,bool record_comm)1713 void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
1714 bool record_comm)
1715 {
1716 bool sfc = (target_elements >= 0);
1717
1718 UpdateLayers();
1719
1720 // *** STEP 1: communicate new rank assignments for the ghost layer ***
1721
1722 NeighborElementRankMessage::Map send_ghost_ranks, recv_ghost_ranks;
1723
1724 ghost_layer.Sort([&](const int a, const int b)
1725 {
1726 return elements[a].rank < elements[b].rank;
1727 });
1728
1729 {
1730 Array<int> rank_neighbors;
1731
1732 // loop over neighbor ranks and their elements
1733 int begin = 0, end = 0;
1734 while (end < ghost_layer.Size())
1735 {
1736 // find range of elements belonging to one rank
1737 int rank = elements[ghost_layer[begin]].rank;
1738 while (end < ghost_layer.Size() &&
1739 elements[ghost_layer[end]].rank == rank) { end++; }
1740
1741 Array<int> rank_elems;
1742 rank_elems.MakeRef(&ghost_layer[begin], end - begin);
1743
1744 // find elements within boundary_layer that are neighbors to 'rank'
1745 rank_neighbors.SetSize(0);
1746 NeighborExpand(rank_elems, rank_neighbors, &boundary_layer);
1747
1748 // send a message with new rank assignments within 'rank_neighbors'
1749 NeighborElementRankMessage& msg = send_ghost_ranks[rank];
1750 msg.SetNCMesh(this);
1751
1752 msg.Reserve(rank_neighbors.Size());
1753 for (int i = 0; i < rank_neighbors.Size(); i++)
1754 {
1755 int elem = rank_neighbors[i];
1756 msg.AddElementRank(elem, new_ranks[elements[elem].index]);
1757 }
1758
1759 msg.Isend(rank, MyComm);
1760
1761 // prepare to receive a message from the neighbor too, these will
1762 // be new the new rank assignments for our ghost layer
1763 recv_ghost_ranks[rank].SetNCMesh(this);
1764
1765 begin = end;
1766 }
1767 }
1768
1769 NeighborElementRankMessage::RecvAll(recv_ghost_ranks, MyComm);
1770
1771 // read new ranks for the ghost layer from messages received
1772 NeighborElementRankMessage::Map::iterator it;
1773 for (it = recv_ghost_ranks.begin(); it != recv_ghost_ranks.end(); ++it)
1774 {
1775 NeighborElementRankMessage &msg = it->second;
1776 for (int i = 0; i < msg.Size(); i++)
1777 {
1778 int ghost_index = elements[msg.elements[i]].index;
1779 MFEM_ASSERT(element_type[ghost_index] == 2, "");
1780 new_ranks[ghost_index] = msg.values[i];
1781 }
1782 }
1783
1784 recv_ghost_ranks.clear();
1785
1786 // *** STEP 2: send elements that no longer belong to us to new assignees ***
1787
1788 /* The result thus far is just the array 'new_ranks' containing new owners
1789 for elements that we currently own plus new owners for the ghost layer.
1790 Next we keep elements that still belong to us and send ElementSets with
1791 the remaining elements to their new owners. Each batch of elements needs
1792 to be sent together with their neighbors so the receiver also gets a
1793 ghost layer that is up to date (this is why we needed Step 1). */
1794
1795 int received_elements = 0;
1796 for (int i = 0; i < leaf_elements.Size(); i++)
1797 {
1798 Element &el = elements[leaf_elements[i]];
1799 if (el.rank == MyRank && new_ranks[i] == MyRank)
1800 {
1801 received_elements++; // initialize to number of elements we're keeping
1802 }
1803 el.rank = new_ranks[i];
1804 }
1805
1806 int nsent = 0, nrecv = 0; // for debug check
1807
1808 RebalanceMessage::Map send_elems;
1809 {
1810 // sort elements we own by the new rank
1811 Array<int> owned_elements;
1812 owned_elements.MakeRef(leaf_elements.GetData(), NElements);
1813 owned_elements.Sort([&](const int a, const int b)
1814 {
1815 return elements[a].rank < elements[b].rank;
1816 });
1817
1818 Array<int> batch;
1819 batch.Reserve(1024);
1820
1821 // send elements to new owners
1822 int begin = 0, end = 0;
1823 while (end < NElements)
1824 {
1825 // find range of elements belonging to one rank
1826 int rank = elements[owned_elements[begin]].rank;
1827 while (end < owned_elements.Size() &&
1828 elements[owned_elements[end]].rank == rank) { end++; }
1829
1830 if (rank != MyRank)
1831 {
1832 Array<int> rank_elems;
1833 rank_elems.MakeRef(&owned_elements[begin], end - begin);
1834
1835 // expand the 'rank_elems' set by its neighbor elements (ghosts)
1836 batch.SetSize(0);
1837 NeighborExpand(rank_elems, batch);
1838
1839 // send the batch
1840 RebalanceMessage &msg = send_elems[rank];
1841 msg.SetNCMesh(this);
1842
1843 msg.Reserve(batch.Size());
1844 for (int i = 0; i < batch.Size(); i++)
1845 {
1846 int elem = batch[i];
1847 Element &el = elements[elem];
1848
1849 if ((element_type[el.index] & 1) || el.rank != rank)
1850 {
1851 msg.AddElementRank(elem, el.rank);
1852 }
1853 // NOTE: we skip 'ghosts' that are of the receiver's rank because
1854 // they are not really ghosts and would get sent multiple times,
1855 // disrupting the termination mechanism in Step 4.
1856 }
1857
1858 if (sfc)
1859 {
1860 msg.Isend(rank, MyComm);
1861 }
1862 else
1863 {
1864 // custom partitioning needs synchronous sends
1865 msg.Issend(rank, MyComm);
1866 }
1867 nsent++;
1868
1869 // also: record what elements we sent (excluding the ghosts)
1870 // so that SendRebalanceDofs can later send data for them
1871 if (record_comm)
1872 {
1873 send_rebalance_dofs[rank].SetElements(rank_elems, this);
1874 }
1875 }
1876
1877 begin = end;
1878 }
1879 }
1880
1881 // *** STEP 3: receive elements from others ***
1882
1883 RebalanceMessage msg;
1884 msg.SetNCMesh(this);
1885
1886 if (sfc)
1887 {
1888 /* We don't know from whom we're going to receive, so we need to probe.
1889 However, for the default SFC partitioning, we do know how many elements
1890 we're going to own eventually, so the termination condition is easy. */
1891
1892 while (received_elements < target_elements)
1893 {
1894 int rank, size;
1895 RebalanceMessage::Probe(rank, size, MyComm);
1896
1897 // receive message; note: elements are created as the message is decoded
1898 msg.Recv(rank, size, MyComm);
1899 nrecv++;
1900
1901 for (int i = 0; i < msg.Size(); i++)
1902 {
1903 int elem_rank = msg.values[i];
1904 elements[msg.elements[i]].rank = elem_rank;
1905
1906 if (elem_rank == MyRank) { received_elements++; }
1907 }
1908
1909 // save the ranks we received from, for later use in RecvRebalanceDofs
1910 if (record_comm)
1911 {
1912 recv_rebalance_dofs[rank].SetNCMesh(this);
1913 }
1914 }
1915
1916 Update();
1917
1918 RebalanceMessage::WaitAllSent(send_elems);
1919 }
1920 else
1921 {
1922 /* The case (target_elements < 0) is used for custom partitioning.
1923 Here we need to employ the "non-blocking consensus" algorithm
1924 (https://scorec.rpi.edu/REPORTS/2015-9.pdf) to determine when the
1925 element exchange is finished. The algorithm uses a non-blocking
1926 barrier. */
1927
1928 MPI_Request barrier = MPI_REQUEST_NULL;
1929 int done = 0;
1930
1931 while (!done)
1932 {
1933 int rank, size;
1934 while (RebalanceMessage::IProbe(rank, size, MyComm))
1935 {
1936 // receive message; note: elements are created as the msg is decoded
1937 msg.Recv(rank, size, MyComm);
1938 nrecv++;
1939
1940 for (int i = 0; i < msg.Size(); i++)
1941 {
1942 elements[msg.elements[i]].rank = msg.values[i];
1943 }
1944
1945 // save the ranks we received from, for later use in RecvRebalanceDofs
1946 if (record_comm)
1947 {
1948 recv_rebalance_dofs[rank].SetNCMesh(this);
1949 }
1950 }
1951
1952 if (barrier != MPI_REQUEST_NULL)
1953 {
1954 MPI_Test(&barrier, &done, MPI_STATUS_IGNORE);
1955 }
1956 else
1957 {
1958 if (RebalanceMessage::TestAllSent(send_elems))
1959 {
1960 int err = MPI_Ibarrier(MyComm, &barrier);
1961
1962 MFEM_VERIFY(err == MPI_SUCCESS, "");
1963 MFEM_VERIFY(barrier != MPI_REQUEST_NULL, "");
1964 }
1965 }
1966 }
1967
1968 Update();
1969 }
1970
1971 NeighborElementRankMessage::WaitAllSent(send_ghost_ranks);
1972
1973 #ifdef MFEM_DEBUG
1974 int glob_sent, glob_recv;
1975 MPI_Reduce(&nsent, &glob_sent, 1, MPI_INT, MPI_SUM, 0, MyComm);
1976 MPI_Reduce(&nrecv, &glob_recv, 1, MPI_INT, MPI_SUM, 0, MyComm);
1977
1978 if (MyRank == 0)
1979 {
1980 MFEM_ASSERT(glob_sent == glob_recv,
1981 "(glob_sent, glob_recv) = ("
1982 << glob_sent << ", " << glob_recv << ")");
1983 }
1984 #endif
1985 }
1986
1987
SendRebalanceDofs(int old_ndofs,const Table & old_element_dofs,long old_global_offset,FiniteElementSpace * space)1988 void ParNCMesh::SendRebalanceDofs(int old_ndofs,
1989 const Table &old_element_dofs,
1990 long old_global_offset,
1991 FiniteElementSpace *space)
1992 {
1993 Array<int> dofs;
1994 int vdim = space->GetVDim();
1995
1996 // fill messages (prepared by Rebalance) with element DOFs
1997 RebalanceDofMessage::Map::iterator it;
1998 for (it = send_rebalance_dofs.begin(); it != send_rebalance_dofs.end(); ++it)
1999 {
2000 RebalanceDofMessage &msg = it->second;
2001 msg.dofs.clear();
2002 int ne = msg.elem_ids.size();
2003 if (ne)
2004 {
2005 msg.dofs.reserve(old_element_dofs.RowSize(msg.elem_ids[0]) * ne * vdim);
2006 }
2007 for (int i = 0; i < ne; i++)
2008 {
2009 old_element_dofs.GetRow(msg.elem_ids[i], dofs);
2010 space->DofsToVDofs(dofs, old_ndofs);
2011 msg.dofs.insert(msg.dofs.end(), dofs.begin(), dofs.end());
2012 }
2013 msg.dof_offset = old_global_offset;
2014 }
2015
2016 // send the DOFs to element recipients from last Rebalance()
2017 RebalanceDofMessage::IsendAll(send_rebalance_dofs, MyComm);
2018 }
2019
2020
RecvRebalanceDofs(Array<int> & elements,Array<long> & dofs)2021 void ParNCMesh::RecvRebalanceDofs(Array<int> &elements, Array<long> &dofs)
2022 {
2023 // receive from the same ranks as in last Rebalance()
2024 RebalanceDofMessage::RecvAll(recv_rebalance_dofs, MyComm);
2025
2026 // count the size of the result
2027 int ne = 0, nd = 0;
2028 RebalanceDofMessage::Map::iterator it;
2029 for (it = recv_rebalance_dofs.begin(); it != recv_rebalance_dofs.end(); ++it)
2030 {
2031 RebalanceDofMessage &msg = it->second;
2032 ne += msg.elem_ids.size();
2033 nd += msg.dofs.size();
2034 }
2035
2036 elements.SetSize(ne);
2037 dofs.SetSize(nd);
2038
2039 // copy element indices and their DOFs
2040 ne = nd = 0;
2041 for (it = recv_rebalance_dofs.begin(); it != recv_rebalance_dofs.end(); ++it)
2042 {
2043 RebalanceDofMessage &msg = it->second;
2044 for (unsigned i = 0; i < msg.elem_ids.size(); i++)
2045 {
2046 elements[ne++] = msg.elem_ids[i];
2047 }
2048 for (unsigned i = 0; i < msg.dofs.size(); i++)
2049 {
2050 dofs[nd++] = msg.dof_offset + msg.dofs[i];
2051 }
2052 }
2053
2054 RebalanceDofMessage::WaitAllSent(send_rebalance_dofs);
2055 }
2056
2057
2058 //// ElementSet ////////////////////////////////////////////////////////////////
2059
ElementSet(const ElementSet & other)2060 ParNCMesh::ElementSet::ElementSet(const ElementSet &other)
2061 : ncmesh(other.ncmesh), include_ref_types(other.include_ref_types)
2062 {
2063 other.data.Copy(data);
2064 }
2065
WriteInt(int value)2066 void ParNCMesh::ElementSet::WriteInt(int value)
2067 {
2068 // helper to put an int to the data array
2069 data.Append(value & 0xff);
2070 data.Append((value >> 8) & 0xff);
2071 data.Append((value >> 16) & 0xff);
2072 data.Append((value >> 24) & 0xff);
2073 }
2074
GetInt(int pos) const2075 int ParNCMesh::ElementSet::GetInt(int pos) const
2076 {
2077 // helper to get an int from the data array
2078 return (int) data[pos] +
2079 ((int) data[pos+1] << 8) +
2080 ((int) data[pos+2] << 16) +
2081 ((int) data[pos+3] << 24);
2082 }
2083
FlagElements(const Array<int> & elements,char flag)2084 void ParNCMesh::ElementSet::FlagElements(const Array<int> &elements, char flag)
2085 {
2086 for (int i = 0; i < elements.Size(); i++)
2087 {
2088 int elem = elements[i];
2089 while (elem >= 0)
2090 {
2091 Element &el = ncmesh->elements[elem];
2092 if (el.flag == flag) { break; }
2093 el.flag = flag;
2094 elem = el.parent;
2095 }
2096 }
2097 }
2098
EncodeTree(int elem)2099 void ParNCMesh::ElementSet::EncodeTree(int elem)
2100 {
2101 Element &el = ncmesh->elements[elem];
2102 if (!el.ref_type)
2103 {
2104 // we reached a leaf, mark this as zero child mask
2105 data.Append(0);
2106 }
2107 else
2108 {
2109 // check which subtrees contain marked elements
2110 int mask = 0;
2111 for (int i = 0; i < 8; i++)
2112 {
2113 if (el.child[i] >= 0 && ncmesh->elements[el.child[i]].flag)
2114 {
2115 mask |= 1 << i;
2116 }
2117 }
2118
2119 // write the bit mask and visit the subtrees
2120 data.Append(mask);
2121 if (include_ref_types)
2122 {
2123 data.Append(el.ref_type);
2124 }
2125
2126 for (int i = 0; i < 8; i++)
2127 {
2128 if (mask & (1 << i))
2129 {
2130 EncodeTree(el.child[i]);
2131 }
2132 }
2133 }
2134 }
2135
Encode(const Array<int> & elements)2136 void ParNCMesh::ElementSet::Encode(const Array<int> &elements)
2137 {
2138 FlagElements(elements, 1);
2139
2140 // Each refinement tree that contains at least one element from the set
2141 // is encoded as HEADER + TREE, where HEADER is the root element number and
2142 // TREE is the output of EncodeTree().
2143 for (int i = 0; i < ncmesh->root_state.Size(); i++)
2144 {
2145 if (ncmesh->elements[i].flag)
2146 {
2147 WriteInt(i);
2148 EncodeTree(i);
2149 }
2150 }
2151 WriteInt(-1); // mark end of data
2152
2153 FlagElements(elements, 0);
2154 }
2155
2156 #ifdef MFEM_DEBUG
RefPath() const2157 std::string ParNCMesh::ElementSet::RefPath() const
2158 {
2159 std::ostringstream oss;
2160 for (int i = 0; i < ref_path.Size(); i++)
2161 {
2162 oss << " elem " << ref_path[i] << " (";
2163 const Element &el = ncmesh->elements[ref_path[i]];
2164 for (int j = 0; j < GI[el.Geom()].nv; j++)
2165 {
2166 if (j) { oss << ", "; }
2167 oss << ncmesh->RetrieveNode(el, j);
2168 }
2169 oss << ")\n";
2170 }
2171 return oss.str();
2172 }
2173 #endif
2174
DecodeTree(int elem,int & pos,Array<int> & elements) const2175 void ParNCMesh::ElementSet::DecodeTree(int elem, int &pos,
2176 Array<int> &elements) const
2177 {
2178 #ifdef MFEM_DEBUG
2179 ref_path.Append(elem);
2180 #endif
2181 int mask = data[pos++];
2182 if (!mask)
2183 {
2184 elements.Append(elem);
2185 }
2186 else
2187 {
2188 Element &el = ncmesh->elements[elem];
2189 if (include_ref_types)
2190 {
2191 int ref_type = data[pos++];
2192 if (!el.ref_type)
2193 {
2194 ncmesh->RefineElement(elem, ref_type);
2195 }
2196 else { MFEM_ASSERT(ref_type == el.ref_type, "") }
2197 }
2198 else
2199 {
2200 MFEM_ASSERT(el.ref_type != 0, "Path not found:\n"
2201 << RefPath() << " mask = " << mask);
2202 }
2203
2204 for (int i = 0; i < 8; i++)
2205 {
2206 if (mask & (1 << i))
2207 {
2208 DecodeTree(el.child[i], pos, elements);
2209 }
2210 }
2211 }
2212 #ifdef MFEM_DEBUG
2213 ref_path.DeleteLast();
2214 #endif
2215 }
2216
Decode(Array<int> & elements) const2217 void ParNCMesh::ElementSet::Decode(Array<int> &elements) const
2218 {
2219 int root, pos = 0;
2220 while ((root = GetInt(pos)) >= 0)
2221 {
2222 pos += 4;
2223 DecodeTree(root, pos, elements);
2224 }
2225 }
2226
Dump(std::ostream & os) const2227 void ParNCMesh::ElementSet::Dump(std::ostream &os) const
2228 {
2229 write<int>(os, data.Size());
2230 os.write((const char*) data.GetData(), data.Size());
2231 }
2232
Load(std::istream & is)2233 void ParNCMesh::ElementSet::Load(std::istream &is)
2234 {
2235 data.SetSize(read<int>(is));
2236 is.read((char*) data.GetData(), data.Size());
2237 }
2238
2239
2240 //// EncodeMeshIds/DecodeMeshIds ///////////////////////////////////////////////
2241
AdjustMeshIds(Array<MeshId> ids[],int rank)2242 void ParNCMesh::AdjustMeshIds(Array<MeshId> ids[], int rank)
2243 {
2244 GetSharedVertices();
2245 GetSharedEdges();
2246 GetSharedFaces();
2247
2248 if (!shared_edges.masters.Size() &&
2249 !shared_faces.masters.Size()) { return; }
2250
2251 Array<bool> contains_rank(groups.size());
2252 for (unsigned i = 0; i < groups.size(); i++)
2253 {
2254 contains_rank[i] = GroupContains(i, rank);
2255 }
2256
2257 Array<Pair<int, int> > find_v(ids[0].Size());
2258 for (int i = 0; i < ids[0].Size(); i++)
2259 {
2260 find_v[i].one = ids[0][i].index;
2261 find_v[i].two = i;
2262 }
2263 find_v.Sort();
2264
2265 // find vertices of master edges shared with 'rank', and modify their
2266 // MeshIds so their element/local matches the element of the master edge
2267 for (int i = 0; i < shared_edges.masters.Size(); i++)
2268 {
2269 const MeshId &edge_id = shared_edges.masters[i];
2270 if (contains_rank[entity_pmat_group[1][edge_id.index]])
2271 {
2272 int v[2], pos, k;
2273 GetEdgeVertices(edge_id, v);
2274 for (int j = 0; j < 2; j++)
2275 {
2276 if ((pos = find_v.FindSorted(Pair<int, int>(v[j], 0))) != -1)
2277 {
2278 // switch to an element/local that is safe for 'rank'
2279 k = find_v[pos].two;
2280 ChangeVertexMeshIdElement(ids[0][k], edge_id.element);
2281 ChangeRemainingMeshIds(ids[0], pos, find_v);
2282 }
2283 }
2284 }
2285 }
2286
2287 if (!shared_faces.masters.Size()) { return; }
2288
2289 Array<Pair<int, int> > find_e(ids[1].Size());
2290 for (int i = 0; i < ids[1].Size(); i++)
2291 {
2292 find_e[i].one = ids[1][i].index;
2293 find_e[i].two = i;
2294 }
2295 find_e.Sort();
2296
2297 // find vertices/edges of master faces shared with 'rank', and modify their
2298 // MeshIds so their element/local matches the element of the master face
2299 for (int i = 0; i < shared_faces.masters.Size(); i++)
2300 {
2301 const MeshId &face_id = shared_faces.masters[i];
2302 if (contains_rank[entity_pmat_group[2][face_id.index]])
2303 {
2304 int v[4], e[4], eo[4], pos, k;
2305 int nfv = GetFaceVerticesEdges(face_id, v, e, eo);
2306 for (int j = 0; j < nfv; j++)
2307 {
2308 if ((pos = find_v.FindSorted(Pair<int, int>(v[j], 0))) != -1)
2309 {
2310 k = find_v[pos].two;
2311 ChangeVertexMeshIdElement(ids[0][k], face_id.element);
2312 ChangeRemainingMeshIds(ids[0], pos, find_v);
2313 }
2314 if ((pos = find_e.FindSorted(Pair<int, int>(e[j], 0))) != -1)
2315 {
2316 k = find_e[pos].two;
2317 ChangeEdgeMeshIdElement(ids[1][k], face_id.element);
2318 ChangeRemainingMeshIds(ids[1], pos, find_e);
2319 }
2320 }
2321 }
2322 }
2323 }
2324
ChangeVertexMeshIdElement(NCMesh::MeshId & id,int elem)2325 void ParNCMesh::ChangeVertexMeshIdElement(NCMesh::MeshId &id, int elem)
2326 {
2327 Element &el = elements[elem];
2328 MFEM_ASSERT(el.ref_type == 0, "");
2329
2330 GeomInfo& gi = GI[el.Geom()];
2331 for (int i = 0; i < gi.nv; i++)
2332 {
2333 if (nodes[el.node[i]].vert_index == id.index)
2334 {
2335 id.local = i;
2336 id.element = elem;
2337 return;
2338 }
2339 }
2340 MFEM_ABORT("Vertex not found.");
2341 }
2342
ChangeEdgeMeshIdElement(NCMesh::MeshId & id,int elem)2343 void ParNCMesh::ChangeEdgeMeshIdElement(NCMesh::MeshId &id, int elem)
2344 {
2345 Element &old = elements[id.element];
2346 const int *ev = GI[old.Geom()].edges[(int) id.local];
2347 Node* node = nodes.Find(old.node[ev[0]], old.node[ev[1]]);
2348 MFEM_ASSERT(node != NULL, "Edge not found.");
2349
2350 Element &el = elements[elem];
2351 MFEM_ASSERT(el.ref_type == 0, "");
2352
2353 GeomInfo& gi = GI[el.Geom()];
2354 for (int i = 0; i < gi.ne; i++)
2355 {
2356 const int* ev = gi.edges[i];
2357 if ((el.node[ev[0]] == node->p1 && el.node[ev[1]] == node->p2) ||
2358 (el.node[ev[1]] == node->p1 && el.node[ev[0]] == node->p2))
2359 {
2360 id.local = i;
2361 id.element = elem;
2362 return;
2363 }
2364
2365 }
2366 MFEM_ABORT("Edge not found.");
2367 }
2368
ChangeRemainingMeshIds(Array<MeshId> & ids,int pos,const Array<Pair<int,int>> & find)2369 void ParNCMesh::ChangeRemainingMeshIds(Array<MeshId> &ids, int pos,
2370 const Array<Pair<int, int> > &find)
2371 {
2372 const MeshId &first = ids[find[pos].two];
2373 while (++pos < find.Size() && ids[find[pos].two].index == first.index)
2374 {
2375 MeshId &other = ids[find[pos].two];
2376 other.element = first.element;
2377 other.local = first.local;
2378 }
2379 }
2380
EncodeMeshIds(std::ostream & os,Array<MeshId> ids[])2381 void ParNCMesh::EncodeMeshIds(std::ostream &os, Array<MeshId> ids[])
2382 {
2383 std::map<int, int> stream_id;
2384
2385 // get a list of elements involved, dump them to 'os' and create the mapping
2386 // element_id: (Element index -> stream ID)
2387 {
2388 Array<int> elements;
2389 for (int type = 0; type < 3; type++)
2390 {
2391 for (int i = 0; i < ids[type].Size(); i++)
2392 {
2393 elements.Append(ids[type][i].element);
2394 }
2395 }
2396
2397 ElementSet eset(this);
2398 eset.Encode(elements);
2399 eset.Dump(os);
2400
2401 Array<int> decoded;
2402 decoded.Reserve(elements.Size());
2403 eset.Decode(decoded);
2404
2405 for (int i = 0; i < decoded.Size(); i++)
2406 {
2407 stream_id[decoded[i]] = i;
2408 }
2409 }
2410
2411 // write the IDs as element/local pairs
2412 for (int type = 0; type < 3; type++)
2413 {
2414 write<int>(os, ids[type].Size());
2415 for (int i = 0; i < ids[type].Size(); i++)
2416 {
2417 const MeshId& id = ids[type][i];
2418 write<int>(os, stream_id[id.element]); // TODO: variable 1-4 bytes
2419 write<char>(os, id.local);
2420 }
2421 }
2422 }
2423
DecodeMeshIds(std::istream & is,Array<MeshId> ids[])2424 void ParNCMesh::DecodeMeshIds(std::istream &is, Array<MeshId> ids[])
2425 {
2426 // read the list of elements
2427 ElementSet eset(this);
2428 eset.Load(is);
2429
2430 Array<int> elems;
2431 eset.Decode(elems);
2432
2433 // read vertex/edge/face IDs
2434 for (int type = 0; type < 3; type++)
2435 {
2436 int ne = read<int>(is);
2437 ids[type].SetSize(ne);
2438
2439 for (int i = 0; i < ne; i++)
2440 {
2441 int el_num = read<int>(is);
2442 int elem = elems[el_num];
2443 Element &el = elements[elem];
2444
2445 MFEM_VERIFY(!el.ref_type, "not a leaf element: " << el_num);
2446
2447 MeshId &id = ids[type][i];
2448 id.element = elem;
2449 id.local = read<char>(is);
2450
2451 // find vertex/edge/face index
2452 GeomInfo &gi = GI[el.Geom()];
2453 switch (type)
2454 {
2455 case 0:
2456 {
2457 id.index = nodes[el.node[(int) id.local]].vert_index;
2458 break;
2459 }
2460 case 1:
2461 {
2462 const int* ev = gi.edges[(int) id.local];
2463 Node* node = nodes.Find(el.node[ev[0]], el.node[ev[1]]);
2464 MFEM_ASSERT(node && node->HasEdge(), "edge not found.");
2465 id.index = node->edge_index;
2466 break;
2467 }
2468 default:
2469 {
2470 const int* fv = gi.faces[(int) id.local];
2471 Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
2472 el.node[fv[2]], el.node[fv[3]]);
2473 MFEM_ASSERT(face, "face not found.");
2474 id.index = face->index;
2475 }
2476 }
2477 }
2478 }
2479 }
2480
EncodeGroups(std::ostream & os,const Array<GroupId> & ids)2481 void ParNCMesh::EncodeGroups(std::ostream &os, const Array<GroupId> &ids)
2482 {
2483 // get a list of unique GroupIds
2484 std::map<GroupId, GroupId> stream_id;
2485 for (int i = 0; i < ids.Size(); i++)
2486 {
2487 if (i && ids[i] == ids[i-1]) { continue; }
2488 unsigned size = stream_id.size();
2489 GroupId &sid = stream_id[ids[i]];
2490 if (size != stream_id.size()) { sid = size; }
2491 }
2492
2493 // write the unique groups
2494 write<short>(os, stream_id.size());
2495 for (std::map<GroupId, GroupId>::iterator
2496 it = stream_id.begin(); it != stream_id.end(); ++it)
2497 {
2498 write<GroupId>(os, it->second);
2499 if (it->first >= 0)
2500 {
2501 const CommGroup &group = groups[it->first];
2502 write<short>(os, group.size());
2503 for (unsigned i = 0; i < group.size(); i++)
2504 {
2505 write<int>(os, group[i]);
2506 }
2507 }
2508 else
2509 {
2510 // special "invalid" group, marks forwarded rows
2511 write<short>(os, -1);
2512 }
2513 }
2514
2515 // write the list of all GroupIds
2516 write<int>(os, ids.Size());
2517 for (int i = 0; i < ids.Size(); i++)
2518 {
2519 write<GroupId>(os, stream_id[ids[i]]);
2520 }
2521 }
2522
DecodeGroups(std::istream & is,Array<GroupId> & ids)2523 void ParNCMesh::DecodeGroups(std::istream &is, Array<GroupId> &ids)
2524 {
2525 int ngroups = read<short>(is);
2526 Array<GroupId> groups(ngroups);
2527
2528 // read stream groups, convert to our groups
2529 CommGroup ranks;
2530 ranks.reserve(128);
2531 for (int i = 0; i < ngroups; i++)
2532 {
2533 int id = read<GroupId>(is);
2534 int size = read<short>(is);
2535 if (size >= 0)
2536 {
2537 ranks.resize(size);
2538 for (int i = 0; i < size; i++)
2539 {
2540 ranks[i] = read<int>(is);
2541 }
2542 groups[id] = GetGroupId(ranks);
2543 }
2544 else
2545 {
2546 groups[id] = -1; // forwarded
2547 }
2548 }
2549
2550 // read the list of IDs
2551 ids.SetSize(read<int>(is));
2552 for (int i = 0; i < ids.Size(); i++)
2553 {
2554 ids[i] = groups[read<GroupId>(is)];
2555 }
2556 }
2557
2558
2559 //// Messages //////////////////////////////////////////////////////////////////
2560
2561 template<class ValueType, bool RefTypes, int Tag>
Encode(int)2562 void ParNCMesh::ElementValueMessage<ValueType, RefTypes, Tag>::Encode(int)
2563 {
2564 std::ostringstream ostream;
2565
2566 Array<int> tmp_elements;
2567 tmp_elements.MakeRef(elements.data(), elements.size());
2568
2569 ElementSet eset(pncmesh, RefTypes);
2570 eset.Encode(tmp_elements);
2571 eset.Dump(ostream);
2572
2573 // decode the element set to obtain a local numbering of elements
2574 Array<int> decoded;
2575 decoded.Reserve(tmp_elements.Size());
2576 eset.Decode(decoded);
2577
2578 std::map<int, int> element_index;
2579 for (int i = 0; i < decoded.Size(); i++)
2580 {
2581 element_index[decoded[i]] = i;
2582 }
2583
2584 write<int>(ostream, values.size());
2585 MFEM_ASSERT(elements.size() == values.size(), "");
2586
2587 for (unsigned i = 0; i < values.size(); i++)
2588 {
2589 write<int>(ostream, element_index[elements[i]]); // element number
2590 write<ValueType>(ostream, values[i]);
2591 }
2592
2593 ostream.str().swap(data);
2594 }
2595
2596 template<class ValueType, bool RefTypes, int Tag>
Decode(int)2597 void ParNCMesh::ElementValueMessage<ValueType, RefTypes, Tag>::Decode(int)
2598 {
2599 std::istringstream istream(data);
2600
2601 ElementSet eset(pncmesh, RefTypes);
2602 eset.Load(istream);
2603
2604 Array<int> tmp_elements;
2605 eset.Decode(tmp_elements);
2606
2607 int* el = tmp_elements.GetData();
2608 elements.assign(el, el + tmp_elements.Size());
2609 values.resize(elements.size());
2610
2611 int count = read<int>(istream);
2612 for (int i = 0; i < count; i++)
2613 {
2614 int index = read<int>(istream);
2615 MFEM_ASSERT(index >= 0 && (size_t) index < values.size(), "");
2616 values[index] = read<ValueType>(istream);
2617 }
2618
2619 // no longer need the raw data
2620 data.clear();
2621 }
2622
SetElements(const Array<int> & elems,NCMesh * ncmesh)2623 void ParNCMesh::RebalanceDofMessage::SetElements(const Array<int> &elems,
2624 NCMesh *ncmesh)
2625 {
2626 eset.SetNCMesh(ncmesh);
2627 eset.Encode(elems);
2628
2629 Array<int> decoded;
2630 decoded.Reserve(elems.Size());
2631 eset.Decode(decoded);
2632
2633 elem_ids.resize(decoded.Size());
2634 for (int i = 0; i < decoded.Size(); i++)
2635 {
2636 elem_ids[i] = eset.GetNCMesh()->elements[decoded[i]].index;
2637 }
2638 }
2639
write_dofs(std::ostream & os,const std::vector<int> & dofs)2640 static void write_dofs(std::ostream &os, const std::vector<int> &dofs)
2641 {
2642 write<int>(os, dofs.size());
2643 // TODO: we should compress the ints, mostly they are contiguous ranges
2644 os.write((const char*) dofs.data(), dofs.size() * sizeof(int));
2645 }
2646
read_dofs(std::istream & is,std::vector<int> & dofs)2647 static void read_dofs(std::istream &is, std::vector<int> &dofs)
2648 {
2649 dofs.resize(read<int>(is));
2650 is.read((char*) dofs.data(), dofs.size() * sizeof(int));
2651 }
2652
Encode(int)2653 void ParNCMesh::RebalanceDofMessage::Encode(int)
2654 {
2655 std::ostringstream stream;
2656
2657 eset.Dump(stream);
2658 write<long>(stream, dof_offset);
2659 write_dofs(stream, dofs);
2660
2661 stream.str().swap(data);
2662 }
2663
Decode(int)2664 void ParNCMesh::RebalanceDofMessage::Decode(int)
2665 {
2666 std::istringstream stream(data);
2667
2668 eset.Load(stream);
2669 dof_offset = read<long>(stream);
2670 read_dofs(stream, dofs);
2671
2672 data.clear();
2673
2674 Array<int> elems;
2675 eset.Decode(elems);
2676
2677 elem_ids.resize(elems.Size());
2678 for (int i = 0; i < elems.Size(); i++)
2679 {
2680 elem_ids[i] = eset.GetNCMesh()->elements[elems[i]].index;
2681 }
2682 }
2683
2684
2685 //// Utility ///////////////////////////////////////////////////////////////////
2686
GetDebugMesh(Mesh & debug_mesh) const2687 void ParNCMesh::GetDebugMesh(Mesh &debug_mesh) const
2688 {
2689 // create a serial NCMesh containing all our elements (ghosts and all)
2690 NCMesh* copy = new NCMesh(*this);
2691
2692 Array<int> &cle = copy->leaf_elements;
2693 for (int i = 0; i < cle.Size(); i++)
2694 {
2695 Element &el = copy->elements[cle[i]];
2696 el.attribute = el.rank + 1;
2697 }
2698
2699 debug_mesh.InitFromNCMesh(*copy);
2700 debug_mesh.SetAttributes();
2701 debug_mesh.ncmesh = copy;
2702 }
2703
Trim()2704 void ParNCMesh::Trim()
2705 {
2706 NCMesh::Trim();
2707
2708 shared_vertices.Clear();
2709 shared_edges.Clear();
2710 shared_faces.Clear();
2711
2712 for (int i = 0; i < 3; i++)
2713 {
2714 entity_owner[i].DeleteAll();
2715 entity_pmat_group[i].DeleteAll();
2716 entity_index_rank[i].DeleteAll();
2717 }
2718
2719 send_rebalance_dofs.clear();
2720 recv_rebalance_dofs.clear();
2721
2722 old_index_or_rank.DeleteAll();
2723
2724 ClearAuxPM();
2725 }
2726
MemoryUsage() const2727 long ParNCMesh::RebalanceDofMessage::MemoryUsage() const
2728 {
2729 return (elem_ids.capacity() + dofs.capacity()) * sizeof(int);
2730 }
2731
2732 template<typename K, typename V>
map_memory_usage(const std::map<K,V> & map)2733 static long map_memory_usage(const std::map<K, V> &map)
2734 {
2735 long result = 0;
2736 for (typename std::map<K, V>::const_iterator
2737 it = map.begin(); it != map.end(); ++it)
2738 {
2739 result += it->second.MemoryUsage();
2740 result += sizeof(std::pair<K, V>) + 3*sizeof(void*) + sizeof(bool);
2741 }
2742 return result;
2743 }
2744
GroupsMemoryUsage() const2745 long ParNCMesh::GroupsMemoryUsage() const
2746 {
2747 long groups_size = groups.capacity() * sizeof(CommGroup);
2748 for (unsigned i = 0; i < groups.size(); i++)
2749 {
2750 groups_size += groups[i].capacity() * sizeof(int);
2751 }
2752 const int approx_node_size =
2753 sizeof(std::pair<CommGroup, GroupId>) + 3*sizeof(void*) + sizeof(bool);
2754 return groups_size + group_id.size() * approx_node_size;
2755 }
2756
2757 template<typename Type, int Size>
arrays_memory_usage(const Array<Type> (& arrays)[Size])2758 static long arrays_memory_usage(const Array<Type> (&arrays)[Size])
2759 {
2760 long total = 0;
2761 for (int i = 0; i < Size; i++)
2762 {
2763 total += arrays[i].MemoryUsage();
2764 }
2765 return total;
2766 }
2767
MemoryUsage(bool with_base) const2768 long ParNCMesh::MemoryUsage(bool with_base) const
2769 {
2770 long total_groups_owners = 0;
2771 for (int i = 0; i < 3; i++)
2772 {
2773 total_groups_owners += entity_owner[i].MemoryUsage() +
2774 entity_pmat_group[i].MemoryUsage() +
2775 entity_index_rank[i].MemoryUsage();
2776 }
2777
2778 return (with_base ? NCMesh::MemoryUsage() : 0) +
2779 GroupsMemoryUsage() +
2780 arrays_memory_usage(entity_owner) +
2781 arrays_memory_usage(entity_pmat_group) +
2782 arrays_memory_usage(entity_conf_group) +
2783 arrays_memory_usage(entity_elem_local) +
2784 shared_vertices.MemoryUsage() +
2785 shared_edges.MemoryUsage() +
2786 shared_faces.MemoryUsage() +
2787 face_orient.MemoryUsage() +
2788 element_type.MemoryUsage() +
2789 ghost_layer.MemoryUsage() +
2790 boundary_layer.MemoryUsage() +
2791 tmp_owner.MemoryUsage() +
2792 tmp_shared_flag.MemoryUsage() +
2793 arrays_memory_usage(entity_index_rank) +
2794 tmp_neighbors.MemoryUsage() +
2795 map_memory_usage(send_rebalance_dofs) +
2796 map_memory_usage(recv_rebalance_dofs) +
2797 old_index_or_rank.MemoryUsage() +
2798 aux_pm_store.MemoryUsage() +
2799 sizeof(ParNCMesh) - sizeof(NCMesh);
2800 }
2801
PrintMemoryDetail(bool with_base) const2802 int ParNCMesh::PrintMemoryDetail(bool with_base) const
2803 {
2804 if (with_base) { NCMesh::PrintMemoryDetail(); }
2805
2806 mfem::out << GroupsMemoryUsage() << " groups\n"
2807 << arrays_memory_usage(entity_owner) << " entity_owner\n"
2808 << arrays_memory_usage(entity_pmat_group) << " entity_pmat_group\n"
2809 << arrays_memory_usage(entity_conf_group) << " entity_conf_group\n"
2810 << arrays_memory_usage(entity_elem_local) << " entity_elem_local\n"
2811 << shared_vertices.MemoryUsage() << " shared_vertices\n"
2812 << shared_edges.MemoryUsage() << " shared_edges\n"
2813 << shared_faces.MemoryUsage() << " shared_faces\n"
2814 << face_orient.MemoryUsage() << " face_orient\n"
2815 << element_type.MemoryUsage() << " element_type\n"
2816 << ghost_layer.MemoryUsage() << " ghost_layer\n"
2817 << boundary_layer.MemoryUsage() << " boundary_layer\n"
2818 << tmp_owner.MemoryUsage() << " tmp_owner\n"
2819 << tmp_shared_flag.MemoryUsage() << " tmp_shared_flag\n"
2820 << arrays_memory_usage(entity_index_rank) << " entity_index_rank\n"
2821 << tmp_neighbors.MemoryUsage() << " tmp_neighbors\n"
2822 << map_memory_usage(send_rebalance_dofs) << " send_rebalance_dofs\n"
2823 << map_memory_usage(recv_rebalance_dofs) << " recv_rebalance_dofs\n"
2824 << old_index_or_rank.MemoryUsage() << " old_index_or_rank\n"
2825 << aux_pm_store.MemoryUsage() << " aux_pm_store\n"
2826 << sizeof(ParNCMesh) - sizeof(NCMesh) << " ParNCMesh" << std::endl;
2827
2828 return leaf_elements.Size();
2829 }
2830
2831 } // namespace mfem
2832
2833 #endif // MFEM_USE_MPI
2834