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