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