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