MFEM  v3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pncmesh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "mesh_headers.hpp"
17 #include "pncmesh.hpp"
18 #include "../fem/fe_coll.hpp"
19 #include "../fem/fespace.hpp"
20 
21 #include <map>
22 #include <climits> // INT_MIN, INT_MAX
23 
24 namespace mfem
25 {
26 
27 ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh)
28  : NCMesh(ncmesh)
29 {
30  MyComm = comm;
31  MPI_Comm_size(MyComm, &NRanks);
32  MPI_Comm_rank(MyComm, &MyRank);
33 
34  // assign leaf elements to the processors by simply splitting the
35  // sequence of leaf elements into 'NRanks' parts
36  for (int i = 0; i < leaf_elements.Size(); i++)
37  {
38  leaf_elements[i]->rank = InitialPartition(i);
39  }
40 
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 {
51  ClearAuxPM();
52 }
53 
55 {
57 
61 
63  ghost_layer.SetSize(0);
64  boundary_layer.SetSize(0);
65 }
66 
68 {
69  // This is an override of NCMesh::AssignLeafIndices(). The difference is
70  // that we shift all elements we own to the beginning of the array
71  // 'leaf_elements' and assign all ghost elements indices >= NElements. This
72  // will make the ghosts skipped in NCMesh::GetMeshComponents.
73 
74  // Also note that the ordering of ghosts and non-ghosts is preserved here,
75  // which is important for ParNCMesh::GetFaceNeighbors.
76 
77  Array<Element*> ghosts;
78  ghosts.Reserve(leaf_elements.Size());
79 
80  NElements = 0;
81  for (int i = 0; i < leaf_elements.Size(); i++)
82  {
83  Element* elem = leaf_elements[i];
84  if (elem->rank == MyRank)
85  {
86  leaf_elements[NElements++] = elem;
87  }
88  else
89  {
90  ghosts.Append(elem);
91  }
92  }
93  NGhostElements = ghosts.Size();
94 
95  leaf_elements.SetSize(NElements);
96  leaf_elements.Append(ghosts);
97 
99 }
100 
102 {
103  // This is an override of NCMesh::UpdateVertices. This version first
104  // assigns Vertex::index to vertices of elements of our rank. Only these
105  // vertices then make it to the Mesh in NCMesh::GetMeshComponents.
106  // The remaining (ghost) vertices are assigned indices greater or equal to
107  // Mesh::GetNV().
108 
109  for (HashTable<Node>::Iterator it(nodes); it; ++it)
110  {
111  if (it->vertex) { it->vertex->index = -1; }
112  }
113 
114  NVertices = 0;
115  for (int i = 0; i < leaf_elements.Size(); i++)
116  {
117  Element* elem = leaf_elements[i];
118  if (elem->rank == MyRank)
119  {
120  for (int j = 0; j < GI[(int) elem->geom].nv; j++)
121  {
122  int &vindex = elem->node[j]->vertex->index;
123  if (vindex < 0) { vindex = NVertices++; }
124  }
125  }
126  }
127 
129  for (HashTable<Node>::Iterator it(nodes); it; ++it)
130  {
131  if (it->vertex && it->vertex->index >= 0)
132  {
133  vertex_nodeId[it->vertex->index] = it->id;
134  }
135  }
136 
137  NGhostVertices = 0;
138  for (HashTable<Node>::Iterator it(nodes); it; ++it)
139  {
140  if (it->vertex && it->vertex->index < 0)
141  {
142  it->vertex->index = NVertices + (NGhostVertices++);
143  }
144  }
145 }
146 
148 {
149  // This is an override (or extension of) NCMesh::OnMeshUpdated().
150  // In addition to getting edge/face indices from 'mesh', we also
151  // assign indices to ghost edges/faces that don't exist in the 'mesh'.
152 
153  // clear Edge:: and Face::index
154  for (HashTable<Node>::Iterator it(nodes); it; ++it)
155  {
156  if (it->edge) { it->edge->index = -1; }
157  }
158  for (HashTable<Face>::Iterator it(faces); it; ++it)
159  {
160  it->index = -1;
161  }
162 
163  // go assign existing edge/face indices
164  NCMesh::OnMeshUpdated(mesh);
165 
166  // assign ghost edge indices
167  NEdges = mesh->GetNEdges();
168  NGhostEdges = 0;
169  for (HashTable<Node>::Iterator it(nodes); it; ++it)
170  {
171  if (it->edge && it->edge->index < 0)
172  {
173  it->edge->index = NEdges + (NGhostEdges++);
174  }
175  }
176 
177  // assign ghost face indices
178  NFaces = mesh->GetNumFaces();
179  NGhostFaces = 0;
180  for (HashTable<Face>::Iterator it(faces); it; ++it)
181  {
182  if (it->index < 0) { it->index = NFaces + (NGhostFaces++); }
183  }
184 
185  if (Dim == 2)
186  {
187  MFEM_ASSERT(NFaces == NEdges, "");
188  MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
189  }
190 }
191 
193 {
194  // Called by NCMesh::BuildEdgeList when an edge is visited in a leaf element.
195  // This allows us to determine edge ownership and processors that share it
196  // without duplicating all the HashTable lookups in NCMesh::BuildEdgeList().
197 
198  int &owner = edge_owner[edge->index];
199  owner = std::min(owner, elem->rank);
200 
201  index_rank.Append(Connection(edge->index, elem->rank));
202 }
203 
205 {
206  // Analogous to ElementHasEdge.
207 
208  int &owner = face_owner[face->index];
209  owner = std::min(owner, elem->rank);
210 
211  index_rank.Append(Connection(face->index, elem->rank));
212 }
213 
215 {
216  // This is an extension of NCMesh::BuildEdgeList() which also determines
217  // edge ownership, creates edge processor groups and lists shared edges.
218 
219  int nedges = NEdges + NGhostEdges;
220  edge_owner.SetSize(nedges);
221  edge_owner = INT_MAX;
222 
223  index_rank.SetSize(12*leaf_elements.Size() * 3/2);
224  index_rank.SetSize(0);
225 
227 
229 
230  index_rank.Sort();
231  index_rank.Unique();
233  index_rank.DeleteAll();
234 
236 }
237 
239 {
240  // This is an extension of NCMesh::BuildFaceList() which also determines
241  // face ownership, creates face processor groups and lists shared faces.
242 
243  int nfaces = NFaces + NGhostFaces;
244  face_owner.SetSize(nfaces);
245  face_owner = INT_MAX;
246 
247  index_rank.SetSize(6*leaf_elements.Size() * 3/2);
248  index_rank.SetSize(0);
249 
251 
253 
254  index_rank.Sort();
255  index_rank.Unique();
257  index_rank.DeleteAll();
258 
260 
262 }
263 
264 struct MasterSlaveInfo
265 {
266  int master; // master index if this is a slave
267  int slaves_begin, slaves_end; // slave list if this is a master
268  MasterSlaveInfo() : master(-1), slaves_begin(0), slaves_end(0) {}
269 };
270 
271 void ParNCMesh::AddMasterSlaveRanks(int nitems, const NCList& list)
272 {
273  // create an auxiliary structure for each edge/face
274  std::vector<MasterSlaveInfo> info(nitems);
275 
276  for (unsigned i = 0; i < list.masters.size(); i++)
277  {
278  const Master &mf = list.masters[i];
279  info[mf.index].slaves_begin = mf.slaves_begin;
280  info[mf.index].slaves_end = mf.slaves_end;
281  }
282  for (unsigned i = 0; i < list.slaves.size(); i++)
283  {
284  const Slave& sf = list.slaves[i];
285  info[sf.index].master = sf.master;
286  }
287 
288  // We need the processor groups of master edges/faces to contain the ranks of
289  // their slaves (so that master DOFs get sent to those who share the slaves).
290  // Conversely, we need the groups of slave edges/faces to contain the ranks
291  // of their masters. Both can be done by appending more items to the
292  // 'index_rank' array, before it is sorted and converted to the group table.
293  // (Note that a master/slave edge can be shared by more than one processor.)
294 
295  int size = index_rank.Size();
296  for (int i = 0; i < size; i++)
297  {
298  int index = index_rank[i].from;
299  int rank = index_rank[i].to;
300 
301  const MasterSlaveInfo &msi = info[index];
302  if (msi.master >= 0)
303  {
304  // 'index' is a slave, add its rank to the master's group
305  index_rank.Append(Connection(msi.master, rank));
306  }
307  else
308  {
309  for (int j = msi.slaves_begin; j < msi.slaves_end; j++)
310  {
311  // 'index' is a master, add its rank to the groups of the slaves
312  index_rank.Append(Connection(list.slaves[j].index, rank));
313  }
314  }
315  }
316 }
317 
318 static bool is_shared(const Table& groups, int index, int MyRank)
319 {
320  // A vertex/edge/face is shared if its group contains more than one processor
321  // and at the same time one of them is ourselves.
322 
323  int size = groups.RowSize(index);
324  if (size <= 1)
325  {
326  return false;
327  }
328 
329  const int* group = groups.GetRow(index);
330  for (int i = 0; i < size; i++)
331  {
332  if (group[i] == MyRank) { return true; }
333  }
334 
335  return false;
336 }
337 
338 void ParNCMesh::MakeShared(const Table &groups, const NCList &list,
339  NCList &shared)
340 {
341  shared.Clear();
342 
343  for (unsigned i = 0; i < list.conforming.size(); i++)
344  {
345  if (is_shared(groups, list.conforming[i].index, MyRank))
346  {
347  shared.conforming.push_back(list.conforming[i]);
348  }
349  }
350  for (unsigned i = 0; i < list.masters.size(); i++)
351  {
352  if (is_shared(groups, list.masters[i].index, MyRank))
353  {
354  shared.masters.push_back(list.masters[i]);
355  }
356  }
357  for (unsigned i = 0; i < list.slaves.size(); i++)
358  {
359  if (is_shared(groups, list.slaves[i].index, MyRank))
360  {
361  shared.slaves.push_back(list.slaves[i]);
362  }
363  }
364 }
365 
367 {
368  int nvertices = NVertices + NGhostVertices;
369  vertex_owner.SetSize(nvertices);
370  vertex_owner = INT_MAX;
371 
372  index_rank.SetSize(8*leaf_elements.Size());
373  index_rank.SetSize(0);
374 
375  Array<MeshId> vertex_id(nvertices);
376 
377  // similarly to edges/faces, we loop over the vertices of all leaf elements
378  // to determine which processors share each vertex
379  for (int i = 0; i < leaf_elements.Size(); i++)
380  {
381  Element* elem = leaf_elements[i];
382  for (int j = 0; j < GI[(int) elem->geom].nv; j++)
383  {
384  Node* node = elem->node[j];
385  int index = node->vertex->index;
386 
387  int &owner = vertex_owner[index];
388  owner = std::min(owner, elem->rank);
389 
390  index_rank.Append(Connection(index, elem->rank));
391 
392  MeshId &id = vertex_id[index];
393  id.index = (node->edge ? -1 : index);
394  id.element = elem;
395  id.local = j;
396  }
397  }
398 
399  index_rank.Sort();
400  index_rank.Unique();
402  index_rank.DeleteAll();
403 
404  // create a list of shared vertices, skip obviously slave vertices
405  // (for simplicity, we don't guarantee to skip all slave vertices)
407  for (int i = 0; i < nvertices; i++)
408  {
409  if (is_shared(vertex_group, i, MyRank) && vertex_id[i].index >= 0)
410  {
411  shared_vertices.conforming.push_back(vertex_id[i]);
412  }
413  }
414 }
415 
417  int local[2])
418 {
419  // Return face orientation in e2, assuming the face has orientation 0 in e1.
420  int ids[2][4];
421  Element* e[2] = { e1, e2 };
422  for (int i = 0; i < 2; i++)
423  {
424  int lf = find_hex_face(find_node(e[i], face->p1),
425  find_node(e[i], face->p2),
426  find_node(e[i], face->p3));
427  if (local) { local[i] = lf; }
428 
429  // get node IDs for the face as seen from e[i]
430  const int* fv = GI[Geometry::CUBE].faces[lf];
431  for (int j = 0; j < 4; j++)
432  {
433  ids[i][j] = e[i]->node[fv[j]]->id;
434  }
435  }
436  return Mesh::GetQuadOrientation(ids[0], ids[1]);
437 }
438 
440 {
441  if (Dim < 3) { return; }
442 
443  // Calculate orientation of shared conforming faces.
444  // NOTE: face orientation is calculated relative to its lower rank element.
445  // Thanks to the ghost layer this can be done locally, without communication.
446 
448  face_orient = 0;
449 
450  for (HashTable<Face>::Iterator it(faces); it; ++it)
451  {
452  if (it->ref_count == 2 && it->index < NFaces)
453  {
454  Element* e[2] = { it->elem[0], it->elem[1] };
455  if (e[0]->rank == e[1]->rank) { continue; }
456  if (e[0]->rank > e[1]->rank) { std::swap(e[0], e[1]); }
457 
458  face_orient[it->index] = get_face_orientation(it, e[0], e[1]);
459  }
460  }
461 }
462 
463 void ParNCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
464  Array<int> &bdr_vertices,
465  Array<int> &bdr_edges)
466 {
467  NCMesh::GetBoundaryClosure(bdr_attr_is_ess, bdr_vertices, bdr_edges);
468 
469  int i, j;
470  // filter out ghost vertices
471  for (i = j = 0; i < bdr_vertices.Size(); i++)
472  {
473  if (bdr_vertices[i] < NVertices) { bdr_vertices[j++] = bdr_vertices[i]; }
474  }
475  bdr_vertices.SetSize(j);
476 
477  // filter out ghost edges
478  for (i = j = 0; i < bdr_edges.Size(); i++)
479  {
480  if (bdr_edges[i] < NEdges) { bdr_edges[j++] = bdr_edges[i]; }
481  }
482  bdr_edges.SetSize(j);
483 }
484 
485 
487 
489 {
490  if (element_type.Size()) { return; }
491 
492  int nleaves = leaf_elements.Size();
493 
494  element_type.SetSize(nleaves);
495  for (int i = 0; i < nleaves; i++)
496  {
497  element_type[i] = (leaf_elements[i]->rank == MyRank) ? 1 : 0;
498  }
499 
500  // determine the ghost layer
501  Array<char> ghost_set;
502  FindSetNeighbors(element_type, NULL, &ghost_set);
503 
504  // find the neighbors of the ghost layer
505  Array<char> boundary_set;
506  FindSetNeighbors(ghost_set, NULL, &boundary_set);
507 
508  ghost_layer.SetSize(0);
509  boundary_layer.SetSize(0);
510  for (int i = 0; i < nleaves; i++)
511  {
512  char &etype = element_type[i];
513  if (ghost_set[i])
514  {
515  etype = 2;
516  ghost_layer.Append(leaf_elements[i]);
517  }
518  else if (boundary_set[i] && etype)
519  {
520  etype = 3;
521  boundary_layer.Append(leaf_elements[i]);
522  }
523  }
524 }
525 
527 {
528  if (!elem->ref_type)
529  {
530  return (element_type[elem->index] == type);
531  }
532  else
533  {
534  for (int i = 0; i < 8; i++)
535  {
536  if (elem->child[i] &&
537  !CheckElementType(elem->child[i], type)) { return false; }
538  }
539  return true;
540  }
541 }
542 
544  Array<int> &ranks)
545 {
546  ranks.SetSize(0); // preserve capacity
547 
548  // big shortcut: there are no neighbors if element_type == 1
549  if (CheckElementType(elem, 1)) { return; }
550 
551  // ok, we do need to look for neighbors;
552  // at least we can only search in the ghost layer
553  tmp_neighbors.SetSize(0);
555 
556  // return a list of processors
557  for (int i = 0; i < tmp_neighbors.Size(); i++)
558  {
559  ranks.Append(tmp_neighbors[i]->rank);
560  }
561  ranks.Sort();
562  ranks.Unique();
563 }
564 
565 template<class T>
566 static void set_to_array(const std::set<T> &set, Array<T> &array)
567 {
568  array.Reserve(set.size());
569  array.SetSize(0);
570  for (std::set<int>::iterator it = set.begin(); it != set.end(); ++it)
571  {
572  array.Append(*it);
573  }
574 }
575 
577 {
578  UpdateLayers();
579 
580  std::set<int> ranks;
581  for (int i = 0; i < ghost_layer.Size(); i++)
582  {
583  ranks.insert(ghost_layer[i]->rank);
584  }
585  set_to_array(ranks, neighbors);
586 }
587 
589 {
590  return (a->rank != b->rank) ? a->rank < b->rank
591  /* */ : a->index < b->index;
592 }
593 
595 {
596  ClearAuxPM();
597 
598  const NCList &shared = (Dim == 3) ? GetSharedFaces() : GetSharedEdges();
599  const NCList &full_list = (Dim == 3) ? GetFaceList() : GetEdgeList();
600 
601  Array<Element*> fnbr;
602  Array<Connection> send_elems;
603 
604  int bound = shared.conforming.size() + shared.slaves.size();
605 
606  fnbr.Reserve(bound);
607  send_elems.Reserve(bound);
608 
609  // go over all shared faces and collect face neighbor elements
610  for (unsigned i = 0; i < shared.conforming.size(); i++)
611  {
612  const MeshId &cf = shared.conforming[i];
613  Face* face = GetFace(cf.element, cf.local);
614  MFEM_ASSERT(face != NULL, "");
615 
616  Element* e[2] = { face->elem[0], face->elem[1] };
617  MFEM_ASSERT(e[0] != NULL && e[1] != NULL, "");
618 
619  if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
620  MFEM_ASSERT(e[0]->rank != MyRank && e[1]->rank == MyRank, "");
621 
622  fnbr.Append(e[0]);
623  send_elems.Append(Connection(e[0]->rank, e[1]->index));
624  }
625 
626  for (unsigned i = 0; i < shared.masters.size(); i++)
627  {
628  const Master &mf = shared.masters[i];
629  for (int j = mf.slaves_begin; j < mf.slaves_end; j++)
630  {
631  const Slave &sf = full_list.slaves[j];
632 
633  Element* e[2] = { mf.element, sf.element };
634  MFEM_ASSERT(e[0] != NULL && e[1] != NULL, "");
635 
636  bool loc0 = (e[0]->rank == MyRank);
637  bool loc1 = (e[1]->rank == MyRank);
638  if (loc0 == loc1) { continue; }
639  if (loc0) { std::swap(e[0], e[1]); }
640 
641  fnbr.Append(e[0]);
642  send_elems.Append(Connection(e[0]->rank, e[1]->index));
643  }
644  }
645 
646  MFEM_ASSERT(fnbr.Size() <= bound, "oops, bad upper bound");
647 
648  // remove duplicate face neighbor elements and sort them by rank & index
649  // (note that the send table is sorted the same way and the order is also the
650  // same on different processors, this is important for ExchangeFaceNbrData)
651  fnbr.Sort();
652  fnbr.Unique();
654 
655  // put the ranks into 'face_nbr_group'
656  for (int i = 0; i < fnbr.Size(); i++)
657  {
658  if (!i || fnbr[i]->rank != pmesh.face_nbr_group.Last())
659  {
660  pmesh.face_nbr_group.Append(fnbr[i]->rank);
661  }
662  }
663  int nranks = pmesh.face_nbr_group.Size();
664 
665  // create a new mfem::Element for each face neighbor element
666  pmesh.face_nbr_elements.SetSize(0);
667  pmesh.face_nbr_elements.Reserve(fnbr.Size());
668 
671 
672  Array<int> fnbr_index(NGhostElements);
673  fnbr_index = -1;
674 
675  std::map<Vertex*, int> vert_map;
676  for (int i = 0; i < fnbr.Size(); i++)
677  {
678  Element* elem = fnbr[i];
679  mfem::Element* fne = NewMeshElement(elem->geom);
680  fne->SetAttribute(elem->attribute);
681  pmesh.face_nbr_elements.Append(fne);
682 
683  GeomInfo& gi = GI[(int) elem->geom];
684  for (int k = 0; k < gi.nv; k++)
685  {
686  int &v = vert_map[elem->node[k]->vertex];
687  if (!v) { v = vert_map.size(); }
688  fne->GetVertices()[k] = v-1;
689  }
690 
691  if (!i || elem->rank != fnbr[i-1]->rank)
692  {
694  }
695 
696  MFEM_ASSERT(elem->index >= NElements, "not a ghost element");
697  fnbr_index[elem->index - NElements] = i;
698  }
699  pmesh.face_nbr_elements_offset.Append(fnbr.Size());
700 
701  // create vertices in 'face_nbr_vertices'
702  pmesh.face_nbr_vertices.SetSize(vert_map.size());
703  std::map<Vertex*, int>::iterator it;
704  for (it = vert_map.begin(); it != vert_map.end(); ++it)
705  {
706  pmesh.face_nbr_vertices[it->second-1].SetCoords(it->first->pos);
707  }
708 
709  // make the 'send_face_nbr_elements' table
710  send_elems.Sort();
711  send_elems.Unique();
712 
713  for (int i = 0, last_rank = -1; i < send_elems.Size(); i++)
714  {
715  Connection &c = send_elems[i];
716  if (c.from != last_rank)
717  {
718  // renumber rank to position in 'face_nbr_group'
719  last_rank = c.from;
720  c.from = pmesh.face_nbr_group.Find(c.from);
721  }
722  else
723  {
724  c.from = send_elems[i-1].from; // avoid search
725  }
726  }
727  pmesh.send_face_nbr_elements.MakeFromList(nranks, send_elems);
728 
729  // go over the shared faces again and modify their Mesh::FaceInfo
730  for (unsigned i = 0; i < shared.conforming.size(); i++)
731  {
732  const MeshId &cf = shared.conforming[i];
733  Face* face = GetFace(cf.element, cf.local);
734 
735  Element* e[2] = { face->elem[0], face->elem[1] };
736  if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
737 
738  Mesh::FaceInfo &fi = pmesh.faces_info[cf.index];
739  fi.Elem2No = -1 - fnbr_index[e[0]->index - NElements];
740 
741  if (Dim == 3)
742  {
743  int local[2];
744  int o = get_face_orientation(face, e[1], e[0], local);
745  fi.Elem2Inf = 64*local[1] + o;
746  }
747  else
748  {
749  fi.Elem2Inf = 64*find_element_edge(e[0], face->p1, face->p3) + 1;
750  }
751  }
752 
753  if (shared.slaves.size())
754  {
755  int nfaces = NFaces, nghosts = NGhostFaces;
756  if (Dim <= 2) { nfaces = NEdges, nghosts = NGhostEdges; }
757 
758  // enlarge Mesh::faces_info for ghost slaves
759  pmesh.faces_info.SetSize(nfaces + nghosts);
760  for (int i = nfaces; i < pmesh.faces_info.Size(); i++)
761  {
762  Mesh::FaceInfo &fi = pmesh.faces_info[i];
763  fi.Elem1No = fi.Elem2No = -1;
764  fi.Elem1Inf = fi.Elem2Inf = 0;
765  fi.NCFace = -1;
766  }
767 
768  // fill in FaceInfo for shared slave faces
769  for (unsigned i = 0; i < shared.masters.size(); i++)
770  {
771  const Master &mf = shared.masters[i];
772  for (int j = mf.slaves_begin; j < mf.slaves_end; j++)
773  {
774  const Slave &sf = full_list.slaves[j];
775 
776  MFEM_ASSERT(sf.element && mf.element, "");
777  bool sloc = (sf.element->rank == MyRank);
778  bool mloc = (mf.element->rank == MyRank);
779  if (sloc == mloc) { continue; }
780 
781  Mesh::FaceInfo &fi = pmesh.faces_info[sf.index];
782  fi.Elem1No = sf.element->index;
783  fi.Elem2No = mf.element->index;
784  fi.Elem1Inf = 64 * sf.local;
785  fi.Elem2Inf = 64 * mf.local;
786 
787  if (!sloc)
788  {
789  std::swap(fi.Elem1No, fi.Elem2No);
790  std::swap(fi.Elem1Inf, fi.Elem2Inf);
791  }
792  MFEM_ASSERT(fi.Elem2No >= NElements, "");
793  fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
794 
795  const DenseMatrix* pm = &sf.point_matrix;
796  if (!sloc && Dim == 3)
797  {
798  // ghost slave in 3D needs flipping orientation
799  DenseMatrix* pm2 = new DenseMatrix(*pm);
800  std::swap((*pm2)(0,1), (*pm2)(0,3));
801  std::swap((*pm2)(1,1), (*pm2)(1,3));
802  aux_pm_store.Append(pm2);
803 
804  fi.Elem2Inf ^= 1;
805  pm = pm2;
806 
807  // The problem is that sf.point_matrix is designed for P matrix
808  // construction and always has orientation relative to the slave
809  // face. In ParMesh::GetSharedFaceTransformations the result
810  // would therefore be the same on both processors, which is not
811  // how that function works for conforming faces. The orientation
812  // of Loc1, Loc2 and Face needs to always be relative to Element
813  // 1, which is the element containing the slave face on one
814  // processor, but on the other it is the element containing the
815  // master face. In the latter case we need to flip the pm.
816  }
817 
818  MFEM_ASSERT(fi.NCFace < 0, "");
819  fi.NCFace = pmesh.nc_faces_info.Size();
820  pmesh.nc_faces_info.Append(Mesh::NCFaceInfo(true, sf.master, pm));
821  }
822  }
823  }
824 
825  // NOTE: this function skips ParMesh::send_face_nbr_vertices and
826  // ParMesh::face_nbr_vertices_offset, these are not used outside of ParMesh
827 }
828 
830 {
831  for (int i = 0; i < aux_pm_store.Size(); i++)
832  {
833  delete aux_pm_store[i];
834  }
835  aux_pm_store.DeleteAll();
836 }
837 
838 
840 
842 {
843  if (elem->ref_type)
844  {
845  bool remove[8];
846  bool removeAll = true;
847 
848  // determine which subtrees can be removed (and whether it's all of them)
849  for (int i = 0; i < 8; i++)
850  {
851  remove[i] = false;
852  if (elem->child[i])
853  {
854  remove[i] = PruneTree(elem->child[i]);
855  if (!remove[i]) { removeAll = false; }
856  }
857  }
858 
859  // all children can be removed, let the (maybe indirect) parent do it
860  if (removeAll) { return true; }
861 
862  // not all children can be removed, but remove those that can be
863  for (int i = 0; i < 8; i++)
864  {
865  if (remove[i]) { DerefineElement(elem->child[i]); }
866  }
867 
868  return false; // need to keep this element and up
869  }
870  else
871  {
872  // return true if this leaf can be removed
873  return elem->rank < 0;
874  }
875 }
876 
878 {
879  if (!Iso && Dim == 3)
880  {
881  MFEM_WARNING("Can't prune 3D aniso meshes yet.");
882  return;
883  }
884 
885  UpdateLayers();
886 
887  for (int i = 0; i < leaf_elements.Size(); i++)
888  {
889  // rank of elements beyond the ghost layer is unknown / not updated
890  if (element_type[i] == 0)
891  {
892  leaf_elements[i]->rank = -1;
893  // NOTE: rank == -1 will make the element disappear from leaf_elements
894  // on next Update, see NCMesh::CollectLeafElements
895  }
896  }
897 
898  // derefine subtrees whose leaves are all unneeded
899  for (int i = 0; i < root_elements.Size(); i++)
900  {
901  if (PruneTree(root_elements[i]))
902  {
904  }
905  }
906 
907  Update();
908 }
909 
910 
911 void ParNCMesh::Refine(const Array<Refinement> &refinements)
912 {
913  for (int i = 0; i < refinements.Size(); i++)
914  {
915  const Refinement &ref = refinements[i];
916  MFEM_VERIFY(ref.ref_type == 7 || Dim < 3,
917  "anisotropic parallel refinement not supported yet in 3D.");
918  }
919  MFEM_VERIFY(Iso || Dim < 3,
920  "parallel refinement of 3D aniso meshes not supported yet.");
921 
923 
924  // create refinement messages to all neighbors (NOTE: some may be empty)
925  Array<int> neighbors;
926  NeighborProcessors(neighbors);
927  for (int i = 0; i < neighbors.Size(); i++)
928  {
929  send_ref[neighbors[i]].SetNCMesh(this);
930  }
931 
932  // populate messages: all refinements that occur next to the processor
933  // boundary need to be sent to the adjoining neighbors so they can keep
934  // their ghost layer up to date
935  Array<int> ranks;
936  ranks.Reserve(64);
937  for (int i = 0; i < refinements.Size(); i++)
938  {
939  const Refinement &ref = refinements[i];
940  MFEM_ASSERT(ref.index < NElements, "");
941  Element* elem = leaf_elements[ref.index];
942  ElementNeighborProcessors(elem, ranks);
943  for (int j = 0; j < ranks.Size(); j++)
944  {
945  send_ref[ranks[j]].AddRefinement(elem, ref.ref_type);
946  }
947  }
948 
949  // send the messages (overlap with local refinements)
951 
952  // do local refinements
953  for (int i = 0; i < refinements.Size(); i++)
954  {
955  const Refinement &ref = refinements[i];
957  }
958 
959  // receive (ghost layer) refinements from all neighbors
960  for (int j = 0; j < neighbors.Size(); j++)
961  {
962  int rank, size;
964 
966  msg.SetNCMesh(this);
967  msg.Recv(rank, size, MyComm);
968 
969  // do the ghost refinements
970  for (int i = 0; i < msg.Size(); i++)
971  {
972  NCMesh::RefineElement(msg.elements[i], msg.values[i]);
973  }
974  }
975 
976  Update();
977 
978  // make sure we can delete the send buffers
980 }
981 
982 
983 void ParNCMesh::LimitNCLevel(int max_nc_level)
984 {
985  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
986 
987  while (1)
988  {
989  Array<Refinement> refinements;
990  GetLimitRefinements(refinements, max_nc_level);
991 
992  long size = refinements.Size(), glob_size;
993  MPI_Allreduce(&size, &glob_size, 1, MPI_LONG, MPI_SUM, MyComm);
994 
995  if (!glob_size) { break; }
996 
997  Refine(refinements);
998  }
999 }
1000 
1001 void ParNCMesh::Derefine(const Array<int> &derefs)
1002 {
1003  MFEM_VERIFY(Dim < 3 || Iso,
1004  "derefinement of 3D anisotropic meshes not implemented yet.");
1005 
1007 
1008  // store fine element ranks
1010  for (int i = 0; i < leaf_elements.Size(); i++)
1011  {
1012  old_index_or_rank[i] = leaf_elements[i]->rank;
1013  }
1014 
1015  // back up the leaf_elements array
1016  Array<Element*> old_elements;
1017  leaf_elements.Copy(old_elements);
1018 
1019  // *** STEP 1: redistribute elements to avoid complex derefinements ***
1020 
1021  Array<int> new_ranks(leaf_elements.Size());
1022  for (int i = 0; i < leaf_elements.Size(); i++)
1023  {
1024  new_ranks[i] = leaf_elements[i]->rank;
1025  }
1026 
1027  // make the lowest rank get all the fine elements for each derefinement
1028  for (int i = 0; i < derefs.Size(); i++)
1029  {
1030  int row = derefs[i];
1031  MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
1032  "invalid derefinement number.");
1033 
1034  const int* fine = derefinements.GetRow(row);
1035  int size = derefinements.RowSize(row);
1036 
1037  int coarse_rank = INT_MAX;
1038  for (int j = 0; j < size; j++)
1039  {
1040  int fine_rank = leaf_elements[fine[j]]->rank;
1041  coarse_rank = std::min(coarse_rank, fine_rank);
1042  }
1043  for (int j = 0; j < size; j++)
1044  {
1045  new_ranks[fine[j]] = coarse_rank;
1046  }
1047  }
1048 
1049  int target_elements = 0;
1050  for (int i = 0; i < new_ranks.Size(); i++)
1051  {
1052  if (new_ranks[i] == MyRank) { target_elements++; }
1053  }
1054 
1055  // redistribute elements slightly to get rid of complex derefinements
1056  // straddling processor boundaries *and* update the ghost layer
1057  RedistributeElements(new_ranks, target_elements, false);
1058 
1059  // *** STEP 2: derefine now, communication similar to Refine() ***
1060 
1062 
1063  // create derefinement messages to all neighbors (NOTE: some may be empty)
1064  Array<int> neighbors;
1065  NeighborProcessors(neighbors);
1066  for (int i = 0; i < neighbors.Size(); i++)
1067  {
1068  send_deref[neighbors[i]].SetNCMesh(this);
1069  }
1070 
1071  // derefinements that occur next to the processor boundary need to be sent
1072  // to the adjoining neighbors to keep their ghost layers in sync
1073  Array<int> ranks;
1074  ranks.Reserve(64);
1075  for (int i = 0; i < derefs.Size(); i++)
1076  {
1077  const int* fine = derefinements.GetRow(derefs[i]);
1078  Element* parent = old_elements[fine[0]]->parent;
1079 
1080  // send derefinement to neighbors
1081  ElementNeighborProcessors(parent, ranks);
1082  for (int j = 0; j < ranks.Size(); j++)
1083  {
1084  send_deref[ranks[j]].AddDerefinement(parent, new_ranks[fine[0]]);
1085  }
1086  }
1088 
1089  // restore old (pre-redistribution) element indices, for SetDerefMatrixCodes
1090  for (int i = 0; i < leaf_elements.Size(); i++)
1091  {
1092  leaf_elements[i]->index = -1;
1093  }
1094  for (int i = 0; i < old_elements.Size(); i++)
1095  {
1096  old_elements[i]->index = i;
1097  }
1098 
1099  // do local derefinements
1100  Array<Element*> coarse;
1101  old_elements.Copy(coarse);
1102  for (int i = 0; i < derefs.Size(); i++)
1103  {
1104  const int* fine = derefinements.GetRow(derefs[i]);
1105  Element* parent = old_elements[fine[0]]->parent;
1106 
1107  // record the relation of the fine elements to their parent
1108  SetDerefMatrixCodes(parent, coarse);
1109 
1110  NCMesh::DerefineElement(parent);
1111  }
1112 
1113  // receive ghost layer derefinements from all neighbors
1114  for (int j = 0; j < neighbors.Size(); j++)
1115  {
1116  int rank, size;
1118 
1120  msg.SetNCMesh(this);
1121  msg.Recv(rank, size, MyComm);
1122 
1123  // do the ghost derefinements
1124  for (int i = 0; i < msg.Size(); i++)
1125  {
1126  Element* elem = msg.elements[i];
1127  if (elem->ref_type)
1128  {
1129  SetDerefMatrixCodes(elem, coarse);
1131  }
1132  elem->rank = msg.values[i];
1133  }
1134  }
1135 
1136  // update leaf_elements, Element::index etc.
1137  Update();
1138 
1139  UpdateLayers();
1140 
1141  // link old fine elements to the new coarse elements
1142  for (int i = 0; i < coarse.Size(); i++)
1143  {
1144  int index = coarse[i]->index;
1145  if (element_type[index] == 0)
1146  {
1147  // this coarse element will get pruned, encode who owns it now
1148  index = -1 - coarse[i]->rank;
1149  }
1150  transforms.embeddings[i].parent = index;
1151  }
1152 
1153  leaf_elements.Copy(old_elements);
1154 
1155  Prune();
1156 
1157  // renumber coarse element indices after pruning
1158  for (int i = 0; i < coarse.Size(); i++)
1159  {
1160  int &index = transforms.embeddings[i].parent;
1161  if (index >= 0)
1162  {
1163  index = old_elements[index]->index;
1164  }
1165  }
1166 
1167  // make sure we can delete all send buffers
1169 }
1170 
1171 
1172 template<typename Type>
1174  const Table &deref_table)
1175 {
1176  const MPI_Datatype datatype = MPITypeMap<Type>::mpi_type;
1177 
1178  Array<MPI_Request*> requests;
1179  Array<int> neigh;
1180 
1181  requests.Reserve(64);
1182  neigh.Reserve(8);
1183 
1184  // make room for ghost values (indices beyond NumElements)
1185  elem_data.SetSize(leaf_elements.Size(), 0);
1186 
1187  for (int i = 0; i < deref_table.Size(); i++)
1188  {
1189  const int* fine = deref_table.GetRow(i);
1190  int size = deref_table.RowSize(i);
1191  MFEM_ASSERT(size <= 8, "");
1192 
1193  int ranks[8], min_rank = INT_MAX, max_rank = INT_MIN;
1194  for (int j = 0; j < size; j++)
1195  {
1196  ranks[j] = leaf_elements[fine[j]]->rank;
1197  min_rank = std::min(min_rank, ranks[j]);
1198  max_rank = std::max(max_rank, ranks[j]);
1199  }
1200 
1201  // exchange values for derefinements that straddle processor boundaries
1202  if (min_rank != max_rank)
1203  {
1204  neigh.SetSize(0);
1205  for (int j = 0; j < size; j++)
1206  {
1207  if (ranks[j] != MyRank) { neigh.Append(ranks[j]); }
1208  }
1209  neigh.Sort();
1210  neigh.Unique();
1211 
1212  for (int j = 0; j < size; j++/*pass*/)
1213  {
1214  Type *data = &elem_data[fine[j]];
1215 
1216  int rnk = ranks[j], len = 1; /*j;
1217  do { j++; } while (j < size && ranks[j] == rnk);
1218  len = j - len;*/
1219 
1220  if (rnk == MyRank)
1221  {
1222  for (int k = 0; k < neigh.Size(); k++)
1223  {
1224  MPI_Request* req = new MPI_Request;
1225  MPI_Isend(data, len, datatype, neigh[k], 292, MyComm, req);
1226  requests.Append(req);
1227  }
1228  }
1229  else
1230  {
1231  MPI_Request* req = new MPI_Request;
1232  MPI_Irecv(data, len, datatype, rnk, 292, MyComm, req);
1233  requests.Append(req);
1234  }
1235  }
1236  }
1237  }
1238 
1239  for (int i = 0; i < requests.Size(); i++)
1240  {
1241  MPI_Wait(requests[i], MPI_STATUS_IGNORE);
1242  delete requests[i];
1243  }
1244 }
1245 
1246 // instantiate SynchronizeDerefinementData for int and double
1247 template void
1248 ParNCMesh::SynchronizeDerefinementData<int>(Array<int> &, const Table &);
1249 template void
1250 ParNCMesh::SynchronizeDerefinementData<double>(Array<double> &, const Table &);
1251 
1252 
1254  Array<int> &level_ok, int max_nc_level)
1255 {
1256  Array<int> leaf_ok(leaf_elements.Size());
1257  leaf_ok = 1;
1258 
1259  // check elements that we own
1260  for (int i = 0; i < deref_table.Size(); i++)
1261  {
1262  const int* fine = deref_table.GetRow(i),
1263  size = deref_table.RowSize(i);
1264 
1265  Element* parent = leaf_elements[fine[0]]->parent;
1266  for (int j = 0; j < size; j++)
1267  {
1268  Element* child = leaf_elements[fine[j]];
1269  if (child->rank == MyRank)
1270  {
1271  int splits[3];
1272  CountSplits(child, splits);
1273 
1274  for (int k = 0; k < Dim; k++)
1275  {
1276  if ((parent->ref_type & (1 << k)) &&
1277  splits[k] >= max_nc_level)
1278  {
1279  leaf_ok[fine[j]] = 0; break;
1280  }
1281  }
1282  }
1283  }
1284  }
1285 
1286  SynchronizeDerefinementData(leaf_ok, deref_table);
1287 
1288  level_ok.SetSize(deref_table.Size());
1289  level_ok = 1;
1290 
1291  for (int i = 0; i < deref_table.Size(); i++)
1292  {
1293  const int* fine = deref_table.GetRow(i),
1294  size = deref_table.RowSize(i);
1295 
1296  for (int j = 0; j < size; j++)
1297  {
1298  if (!leaf_ok[fine[j]])
1299  {
1300  level_ok[i] = 0; break;
1301  }
1302  }
1303  }
1304 }
1305 
1306 
1308 
1310 {
1311  send_rebalance_dofs.clear();
1312  recv_rebalance_dofs.clear();
1313 
1314  Array<Element*> old_elements;
1315  leaf_elements.GetSubArray(0, NElements, old_elements);
1316 
1317  // figure out new assignments for Element::rank
1318  long local_elems = NElements, total_elems = 0;
1319  MPI_Allreduce(&local_elems, &total_elems, 1, MPI_LONG, MPI_SUM, MyComm);
1320 
1321  long first_elem_global = 0;
1322  MPI_Scan(&local_elems, &first_elem_global, 1, MPI_LONG, MPI_SUM, MyComm);
1323  first_elem_global -= local_elems;
1324 
1325  Array<int> new_ranks(leaf_elements.Size());
1326  new_ranks = -1;
1327 
1328  for (int i = 0, j = 0; i < leaf_elements.Size(); i++)
1329  {
1330  if (leaf_elements[i]->rank == MyRank)
1331  {
1332  new_ranks[i] = Partition(first_elem_global + (j++), total_elems);
1333  }
1334  }
1335 
1336  int target_elements = PartitionFirstIndex(MyRank+1, total_elems)
1337  - PartitionFirstIndex(MyRank, total_elems);
1338 
1339  // assign the new ranks and send elements (plus ghosts) to new owners
1340  RedistributeElements(new_ranks, target_elements, true);
1341 
1342  // set up the old index array
1344  old_index_or_rank = -1;
1345  for (int i = 0; i < old_elements.Size(); i++)
1346  {
1347  Element* e = old_elements[i];
1348  if (e->rank == MyRank) { old_index_or_rank[e->index] = i; }
1349  }
1350 
1351  // get rid of elements beyond the new ghost layer
1352  Prune();
1353 }
1354 
1355 
1356 bool ParNCMesh::compare_ranks(const Element* a, const Element* b)
1357 {
1358  return a->rank < b->rank;
1359 }
1360 
1361 void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
1362  bool record_comm)
1363 {
1364  UpdateLayers();
1365 
1366  // *** STEP 1: communicate new rank assignments for the ghost layer ***
1367 
1368  NeighborElementRankMessage::Map send_ghost_ranks, recv_ghost_ranks;
1369 
1370  ghost_layer.Sort(compare_ranks);
1371  {
1372  Array<Element*> rank_neighbors;
1373 
1374  // loop over neighbor ranks and their elements
1375  int begin = 0, end = 0;
1376  while (end < ghost_layer.Size())
1377  {
1378  // find range of elements belonging to one rank
1379  int rank = ghost_layer[begin]->rank;
1380  while (end < ghost_layer.Size() &&
1381  ghost_layer[end]->rank == rank) { end++; }
1382 
1383  Array<Element*> rank_elems;
1384  rank_elems.MakeRef(&ghost_layer[begin], end - begin);
1385 
1386  // find elements within boundary_layer that are neighbors to 'rank'
1387  rank_neighbors.SetSize(0);
1388  NeighborExpand(rank_elems, rank_neighbors, &boundary_layer);
1389 
1390  // send a message with new rank assignments within 'rank_neighbors'
1391  NeighborElementRankMessage& msg = send_ghost_ranks[rank];
1392  msg.SetNCMesh(this);
1393 
1394  msg.Reserve(rank_neighbors.Size());
1395  for (int i = 0; i < rank_neighbors.Size(); i++)
1396  {
1397  Element* e = rank_neighbors[i];
1398  msg.AddElementRank(e, new_ranks[e->index]);
1399  }
1400 
1401  msg.Isend(rank, MyComm);
1402 
1403  // prepare to receive a message from the neighbor too, these will
1404  // be new the new rank assignments for our ghost layer
1405  recv_ghost_ranks[rank].SetNCMesh(this);
1406 
1407  begin = end;
1408  }
1409  }
1410 
1411  NeighborElementRankMessage::RecvAll(recv_ghost_ranks, MyComm);
1412 
1413  // read new ranks for the ghost layer from messages received
1414  NeighborElementRankMessage::Map::iterator it;
1415  for (it = recv_ghost_ranks.begin(); it != recv_ghost_ranks.end(); ++it)
1416  {
1417  NeighborElementRankMessage &msg = it->second;
1418  for (int i = 0; i < msg.Size(); i++)
1419  {
1420  int ghost_index = msg.elements[i]->index;
1421  MFEM_ASSERT(element_type[ghost_index] == 2, "");
1422  new_ranks[ghost_index] = msg.values[i];
1423  }
1424  }
1425 
1426  recv_ghost_ranks.clear();
1427 
1428  // *** STEP 2: send elements that no longer belong to us to new assignees ***
1429 
1430  /* The result thus far is just the array 'new_ranks' containing new owners
1431  for elements that we currently own plus new owners for the ghost layer.
1432  Next we keep elements that still belong to us and send ElementSets with
1433  the remaining elements to their new owners. Each batch of elements needs
1434  to be sent together with their neighbors so the receiver also gets a
1435  ghost layer that is up to date (this is why we needed Step 1). */
1436 
1437  int received_elements = 0;
1438  for (int i = 0; i < leaf_elements.Size(); i++)
1439  {
1440  Element* e = leaf_elements[i];
1441  if (e->rank == MyRank && new_ranks[i] == MyRank)
1442  {
1443  received_elements++; // initialize to number of elements we're keeping
1444  }
1445  e->rank = new_ranks[i];
1446  }
1447 
1448  RebalanceMessage::Map send_elems;
1449  {
1450  // sort elements we own by the new rank
1451  Array<Element*> owned_elements;
1452  owned_elements.MakeRef(leaf_elements.GetData(), NElements);
1453  owned_elements.Sort(compare_ranks);
1454 
1455  Array<Element*> batch;
1456  batch.Reserve(1024);
1457 
1458  // send elements to new owners
1459  int begin = 0, end = 0;
1460  while (end < NElements)
1461  {
1462  // find range of elements belonging to one rank
1463  int rank = owned_elements[begin]->rank;
1464  while (end < owned_elements.Size() &&
1465  owned_elements[end]->rank == rank) { end++; }
1466 
1467  if (rank != MyRank)
1468  {
1469  Array<Element*> rank_elems;
1470  rank_elems.MakeRef(&owned_elements[begin], end - begin);
1471 
1472  // expand the 'rank_elems' set by its neighbor elements (ghosts)
1473  batch.SetSize(0);
1474  NeighborExpand(rank_elems, batch);
1475 
1476  // send the batch
1477  RebalanceMessage &msg = send_elems[rank];
1478  msg.SetNCMesh(this);
1479 
1480  msg.Reserve(batch.Size());
1481  for (int i = 0; i < batch.Size(); i++)
1482  {
1483  Element* e = batch[i];
1484  if ((element_type[e->index] & 1) || e->rank != rank)
1485  {
1486  msg.AddElementRank(e, e->rank);
1487  }
1488  // NOTE: we skip 'ghosts' that are of the receiver's rank because
1489  // they are not really ghosts and would get sent multiple times,
1490  // disrupting the termination mechanism in Step 4.
1491  }
1492 
1493  msg.Isend(rank, MyComm);
1494 
1495  // also: record what elements we sent (excluding the ghosts)
1496  // so that SendRebalanceDofs can later send data for them
1497  if (record_comm)
1498  {
1499  send_rebalance_dofs[rank].SetElements(rank_elems, this);
1500  }
1501  }
1502 
1503  begin = end;
1504  }
1505  }
1506 
1507  // *** STEP 3: receive elements from others ***
1508 
1509  /* We don't know from whom we're going to receive so we need to probe.
1510  Fortunately, we do know how many elements we're going to own eventually
1511  so the termination condition is easy. */
1512 
1513  RebalanceMessage msg;
1514  msg.SetNCMesh(this);
1515 
1516  while (received_elements < target_elements)
1517  {
1518  int rank, size;
1519  RebalanceMessage::Probe(rank, size, MyComm);
1520 
1521  // receive message; note: elements are created as the message is decoded
1522  msg.Recv(rank, size, MyComm);
1523 
1524  for (int i = 0; i < msg.Size(); i++)
1525  {
1526  int elem_rank = msg.values[i];
1527  msg.elements[i]->rank = elem_rank;
1528 
1529  if (elem_rank == MyRank) { received_elements++; }
1530  }
1531 
1532  // save the ranks we received from, for later use in RecvRebalanceDofs
1533  if (record_comm)
1534  {
1535  recv_rebalance_dofs[rank].SetNCMesh(this);
1536  }
1537  }
1538 
1539  Update();
1540 
1541  // make sure we can delete all send buffers
1542  NeighborElementRankMessage::WaitAllSent(send_ghost_ranks);
1544 }
1545 
1546 
1547 void ParNCMesh::SendRebalanceDofs(int old_ndofs,
1548  const Table &old_element_dofs,
1549  long old_global_offset,
1550  FiniteElementSpace *space)
1551 {
1552  Array<int> dofs;
1553  int vdim = space->GetVDim();
1554 
1555  // fill messages (prepared by Rebalance) with element DOFs
1556  RebalanceDofMessage::Map::iterator it;
1557  for (it = send_rebalance_dofs.begin(); it != send_rebalance_dofs.end(); ++it)
1558  {
1559  RebalanceDofMessage &msg = it->second;
1560  msg.dofs.clear();
1561  int ne = msg.elem_ids.size();
1562  if (ne)
1563  {
1564  msg.dofs.reserve(old_element_dofs.RowSize(msg.elem_ids[0]) * ne * vdim);
1565  }
1566  for (int i = 0; i < ne; i++)
1567  {
1568  old_element_dofs.GetRow(msg.elem_ids[i], dofs);
1569  space->DofsToVDofs(dofs, old_ndofs);
1570  msg.dofs.insert(msg.dofs.end(), dofs.begin(), dofs.end());
1571  }
1572  msg.dof_offset = old_global_offset;
1573  }
1574 
1575  // send the DOFs to element recipients from last Rebalance()
1577 }
1578 
1579 
1581 {
1582  // receive from the same ranks as in last Rebalance()
1584 
1585  // count the size of the result
1586  int ne = 0, nd = 0;
1587  RebalanceDofMessage::Map::iterator it;
1588  for (it = recv_rebalance_dofs.begin(); it != recv_rebalance_dofs.end(); ++it)
1589  {
1590  RebalanceDofMessage &msg = it->second;
1591  ne += msg.elem_ids.size();
1592  nd += msg.dofs.size();
1593  }
1594 
1595  elements.SetSize(ne);
1596  dofs.SetSize(nd);
1597 
1598  // copy element indices and their DOFs
1599  ne = nd = 0;
1600  for (it = recv_rebalance_dofs.begin(); it != recv_rebalance_dofs.end(); ++it)
1601  {
1602  RebalanceDofMessage &msg = it->second;
1603  for (unsigned i = 0; i < msg.elem_ids.size(); i++)
1604  {
1605  elements[ne++] = msg.elem_ids[i];
1606  }
1607  for (unsigned i = 0; i < msg.dofs.size(); i++)
1608  {
1609  dofs[nd++] = msg.dof_offset + msg.dofs[i];
1610  }
1611  }
1612 
1614 }
1615 
1616 
1618 
1620  : ncmesh(other.ncmesh), include_ref_types(other.include_ref_types)
1621 {
1622  other.data.Copy(data);
1623 }
1624 
1626 {
1627  // helper to put an int to the data array
1628  data.Append(value & 0xff);
1629  data.Append((value >> 8) & 0xff);
1630  data.Append((value >> 16) & 0xff);
1631  data.Append((value >> 24) & 0xff);
1632 }
1633 
1635 {
1636  // helper to get an int from the data array
1637  return (int) data[pos] +
1638  ((int) data[pos+1] << 8) +
1639  ((int) data[pos+2] << 16) +
1640  ((int) data[pos+3] << 24);
1641 }
1642 
1644  char flag)
1645 {
1646  for (int i = 0; i < elements.Size(); i++)
1647  {
1648  Element* e = elements[i];
1649  while (e && e->flag != flag)
1650  {
1651  e->flag = flag;
1652  e = e->parent;
1653  }
1654  }
1655 }
1656 
1658 {
1659  if (!elem->ref_type)
1660  {
1661  // we reached a leaf, mark this as zero child mask
1662  data.Append(0);
1663  }
1664  else
1665  {
1666  // check which subtrees contain marked elements
1667  int mask = 0;
1668  for (int i = 0; i < 8; i++)
1669  {
1670  if (elem->child[i] && elem->child[i]->flag)
1671  {
1672  mask |= 1 << i;
1673  }
1674  }
1675 
1676  // write the bit mask and visit the subtrees
1677  data.Append(mask);
1678  if (include_ref_types)
1679  {
1680  data.Append(elem->ref_type);
1681  }
1682 
1683  for (int i = 0; i < 8; i++)
1684  {
1685  if (mask & (1 << i))
1686  {
1687  EncodeTree(elem->child[i]);
1688  }
1689  }
1690  }
1691 }
1692 
1694 {
1695  FlagElements(elements, 1);
1696 
1697  // Each refinement tree that contains at least one element from the set
1698  // is encoded as HEADER + TREE, where HEADER is the root element number and
1699  // TREE is the output of EncodeTree().
1700  Array<Element*> &roots = ncmesh->root_elements;
1701  for (int i = 0; i < roots.Size(); i++)
1702  {
1703  if (roots[i]->flag)
1704  {
1705  WriteInt(i);
1706  EncodeTree(roots[i]);
1707  }
1708  }
1709  WriteInt(-1); // mark end of data
1710 
1711  FlagElements(elements, 0);
1712 }
1713 
1715  Array<Element*> &elements) const
1716 {
1717  int mask = data[pos++];
1718  if (!mask)
1719  {
1720  elements.Append(elem);
1721  }
1722  else
1723  {
1724  if (include_ref_types)
1725  {
1726  int ref_type = data[pos++];
1727  if (!elem->ref_type)
1728  {
1729  ncmesh->RefineElement(elem, ref_type);
1730  }
1731  else { MFEM_ASSERT(ref_type == elem->ref_type, "") }
1732  }
1733  else { MFEM_ASSERT(elem->ref_type != 0, ""); }
1734 
1735  for (int i = 0; i < 8; i++)
1736  {
1737  if (mask & (1 << i))
1738  {
1739  DecodeTree(elem->child[i], pos, elements);
1740  }
1741  }
1742  }
1743 }
1744 
1746 {
1747  int root, pos = 0;
1748  while ((root = GetInt(pos)) >= 0)
1749  {
1750  pos += 4;
1751  DecodeTree(ncmesh->root_elements[root], pos, elements);
1752  }
1753 }
1754 
1755 template<typename T>
1756 static inline void write(std::ostream& os, T value)
1757 {
1758  os.write((char*) &value, sizeof(T));
1759 }
1760 
1761 template<typename T>
1762 static inline T read(std::istream& is)
1763 {
1764  T value;
1765  is.read((char*) &value, sizeof(T));
1766  return value;
1767 }
1768 
1769 void ParNCMesh::ElementSet::Dump(std::ostream &os) const
1770 {
1771  write<int>(os, data.Size());
1772  os.write((const char*) data.GetData(), data.Size());
1773 }
1774 
1775 void ParNCMesh::ElementSet::Load(std::istream &is)
1776 {
1777  data.SetSize(read<int>(is));
1778  is.read((char*) data.GetData(), data.Size());
1779 }
1780 
1781 
1783 
1784 void ParNCMesh::EncodeMeshIds(std::ostream &os, Array<MeshId> ids[])
1785 {
1786  std::map<Element*, int> element_id;
1787 
1788  // get a list of elements involved, dump them to 'os' and create the mapping
1789  // element_id: (Element* -> stream ID)
1790  {
1791  Array<Element*> elements;
1792  for (int type = 0; type < 3; type++)
1793  {
1794  for (int i = 0; i < ids[type].Size(); i++)
1795  {
1796  elements.Append(ids[type][i].element);
1797  }
1798  }
1799 
1800  ElementSet eset(this);
1801  eset.Encode(elements);
1802  eset.Dump(os);
1803 
1804  Array<Element*> decoded;
1805  eset.Decode(decoded);
1806 
1807  for (int i = 0; i < decoded.Size(); i++)
1808  {
1809  element_id[decoded[i]] = i;
1810  }
1811  }
1812 
1813  // write the IDs as element/local pairs
1814  for (int type = 0; type < 3; type++)
1815  {
1816  write<int>(os, ids[type].Size());
1817  for (int i = 0; i < ids[type].Size(); i++)
1818  {
1819  const MeshId& id = ids[type][i];
1820  write<int>(os, element_id[id.element]); // TODO: variable 1-4 bytes
1821  write<char>(os, id.local);
1822  }
1823  }
1824 }
1825 
1826 void ParNCMesh::DecodeMeshIds(std::istream &is, Array<MeshId> ids[])
1827 {
1828  // read the list of elements
1829  ElementSet eset(this);
1830  eset.Load(is);
1831 
1832  Array<Element*> elements;
1833  eset.Decode(elements);
1834 
1835  // read vertex/edge/face IDs
1836  for (int type = 0; type < 3; type++)
1837  {
1838  int ne = read<int>(is);
1839  ids[type].SetSize(ne);
1840 
1841  for (int i = 0; i < ne; i++)
1842  {
1843  int el_num = read<int>(is);
1844  Element* elem = elements[el_num];
1845  MFEM_VERIFY(!elem->ref_type, "not a leaf element: " << el_num);
1846 
1847  MeshId &id = ids[type][i];
1848  id.element = elem;
1849  id.local = read<char>(is);
1850 
1851  // find vertex/edge/face index
1852  GeomInfo &gi = GI[(int) elem->geom];
1853  switch (type)
1854  {
1855  case 0:
1856  {
1857  id.index = elem->node[id.local]->vertex->index;
1858  break;
1859  }
1860  case 1:
1861  {
1862  const int* ev = gi.edges[id.local];
1863  Node* node = nodes.Peek(elem->node[ev[0]], elem->node[ev[1]]);
1864  MFEM_ASSERT(node && node->edge, "edge not found.");
1865  id.index = node->edge->index;
1866  break;
1867  }
1868  default:
1869  {
1870  const int* fv = gi.faces[id.local];
1871  Face* face = faces.Peek(elem->node[fv[0]], elem->node[fv[1]],
1872  elem->node[fv[2]], elem->node[fv[3]]);
1873  MFEM_ASSERT(face, "face not found.");
1874  id.index = face->index;
1875  }
1876  }
1877  }
1878  }
1879 }
1880 
1881 
1883 
1885  const Array<int> &dofs)
1886 {
1887  MFEM_ASSERT(type >= 0 && type < 3, "");
1888  id_dofs[type][id].assign(dofs.GetData(), dofs.GetData() + dofs.Size());
1889 }
1890 
1892  Array<int>& dofs, int &ndofs)
1893 {
1894  MFEM_ASSERT(type >= 0 && type < 3, "");
1895 #ifdef MFEM_DEBUG
1896  if (id_dofs[type].find(id) == id_dofs[type].end())
1897  {
1898  MFEM_ABORT("type/ID " << type << "/" << id.index << " not found in "
1899  "neighbor message. Ghost layers out of sync?");
1900  }
1901 #endif
1902  std::vector<int> &vec = id_dofs[type][id];
1903  dofs.SetSize(vec.size());
1904  dofs.Assign(vec.data());
1905  ndofs = this->ndofs;
1906 }
1907 
1909  std::vector<int> &dofs)
1910 {
1911  // Reorder the DOFs into/from a neutral ordering, independent of local
1912  // edge orientation. The processor neutral edge orientation is given by
1913  // the element local vertex numbering, not the mesh vertex numbering.
1914 
1915  const int *ev = NCMesh::GI[(int) id.element->geom].edges[id.local];
1916  int v0 = id.element->node[ev[0]]->vertex->index;
1917  int v1 = id.element->node[ev[1]]->vertex->index;
1918 
1919  if ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1]))
1920  {
1921  std::vector<int> tmp(dofs);
1922 
1923  int nv = fec->DofForGeometry(Geometry::POINT);
1924  int ne = fec->DofForGeometry(Geometry::SEGMENT);
1925  MFEM_ASSERT((int) dofs.size() == 2*nv + ne, "");
1926 
1927  // swap the two vertex DOFs
1928  for (int i = 0; i < 2; i++)
1929  {
1930  for (int k = 0; k < nv; k++)
1931  {
1932  dofs[nv*i + k] = tmp[nv*(1-i) + k];
1933  }
1934  }
1935 
1936  // reorder the edge DOFs
1937  int* ind = fec->DofOrderForOrientation(Geometry::SEGMENT, 0);
1938  for (int i = 0; i < ne; i++)
1939  {
1940  dofs[2*nv + i] = (ind[i] >= 0) ? tmp[2*nv + ind[i]]
1941  /* */ : -1 - tmp[2*nv + (-1 - ind[i])];
1942  }
1943  }
1944 }
1945 
1946 static void write_dofs(std::ostream &os, const std::vector<int> &dofs)
1947 {
1948  write<int>(os, dofs.size());
1949  // TODO: we should compress the ints, mostly they are contiguous ranges
1950  os.write((const char*) dofs.data(), dofs.size() * sizeof(int));
1951 }
1952 
1953 static void read_dofs(std::istream &is, std::vector<int> &dofs)
1954 {
1955  dofs.resize(read<int>(is));
1956  is.read((char*) dofs.data(), dofs.size() * sizeof(int));
1957 }
1958 
1960 {
1961  IdToDofs::iterator it;
1962 
1963  // collect vertex/edge/face IDs
1964  Array<NCMesh::MeshId> ids[3];
1965  for (int type = 0; type < 3; type++)
1966  {
1967  ids[type].Reserve(id_dofs[type].size());
1968  for (it = id_dofs[type].begin(); it != id_dofs[type].end(); ++it)
1969  {
1970  ids[type].Append(it->first);
1971  }
1972  }
1973 
1974  // encode the IDs
1975  std::ostringstream stream;
1976  pncmesh->EncodeMeshIds(stream, ids);
1977 
1978  // dump the DOFs
1979  for (int type = 0; type < 3; type++)
1980  {
1981  for (it = id_dofs[type].begin(); it != id_dofs[type].end(); ++it)
1982  {
1983  if (type == 1) { ReorderEdgeDofs(it->first, it->second); }
1984  write_dofs(stream, it->second);
1985  }
1986 
1987  // no longer need the original data
1988  id_dofs[type].clear();
1989  }
1990 
1991  write<int>(stream, ndofs);
1992 
1993  stream.str().swap(data);
1994 }
1995 
1997 {
1998  std::istringstream stream(data);
1999 
2000  // decode vertex/edge/face IDs
2001  Array<NCMesh::MeshId> ids[3];
2002  pncmesh->DecodeMeshIds(stream, ids);
2003 
2004  // load DOFs
2005  for (int type = 0; type < 3; type++)
2006  {
2007  id_dofs[type].clear();
2008  for (int i = 0; i < ids[type].Size(); i++)
2009  {
2010  const NCMesh::MeshId &id = ids[type][i];
2011  read_dofs(stream, id_dofs[type][id]);
2012  if (type == 1) { ReorderEdgeDofs(id, id_dofs[type][id]); }
2013  }
2014  }
2015 
2016  ndofs = read<int>(stream);
2017 
2018  // no longer need the raw data
2019  data.clear();
2020 }
2021 
2023 {
2024  std::ostringstream stream;
2025 
2026  // write the int set to the stream
2027  write<int>(stream, rows.size());
2028  for (std::set<int>::iterator it = rows.begin(); it != rows.end(); ++it)
2029  {
2030  write<int>(stream, *it);
2031  }
2032 
2033  rows.clear();
2034  stream.str().swap(data);
2035 }
2036 
2038 {
2039  std::istringstream stream(data);
2040 
2041  // read the int set from the stream
2042  rows.clear();
2043  int size = read<int>(stream);
2044  for (int i = 0; i < size; i++)
2045  {
2046  rows.insert(rows.end(), read<int>(stream));
2047  }
2048 
2049  data.clear();
2050 }
2051 
2052 void NeighborRowReply::AddRow(int row, const Array<int> &cols,
2053  const Vector &srow)
2054 {
2055  MFEM_ASSERT(rows.find(row) == rows.end(), "");
2056  Row& row_data = rows[row];
2057  row_data.cols.assign(cols.GetData(), cols.GetData() + cols.Size());
2058  row_data.srow = srow;
2059 }
2060 
2061 void NeighborRowReply::GetRow(int row, Array<int> &cols, Vector &srow)
2062 {
2063  MFEM_ASSERT(rows.find(row) != rows.end(),
2064  "row " << row << " not found in neighbor message.");
2065  Row& row_data = rows[row];
2066  cols.SetSize(row_data.cols.size());
2067  cols.Assign(row_data.cols.data());
2068  srow = row_data.srow;
2069 }
2070 
2072 {
2073  std::ostringstream stream;
2074 
2075  // dump the rows to the stream
2076  write<int>(stream, rows.size());
2077  for (std::map<int, Row>::iterator it = rows.begin(); it != rows.end(); ++it)
2078  {
2079  write<int>(stream, it->first); // row number
2080  Row& row_data = it->second;
2081  MFEM_ASSERT((int) row_data.cols.size() == row_data.srow.Size(), "");
2082  write_dofs(stream, row_data.cols);
2083  stream.write((const char*) row_data.srow.GetData(),
2084  sizeof(double) * row_data.srow.Size());
2085  }
2086 
2087  rows.clear();
2088  stream.str().swap(data);
2089 }
2090 
2092 {
2093  std::istringstream stream(data); // stream makes a copy of data
2094 
2095  // NOTE: there is no rows.clear() since a row reply can be received
2096  // repeatedly and the received rows accumulate.
2097 
2098  // read the rows
2099  int size = read<int>(stream);
2100  for (int i = 0; i < size; i++)
2101  {
2102  Row& row_data = rows[read<int>(stream)];
2103  read_dofs(stream, row_data.cols);
2104  row_data.srow.SetSize(row_data.cols.size());
2105  stream.read((char*) row_data.srow.GetData(),
2106  sizeof(double) * row_data.srow.Size());
2107  }
2108 
2109  data.clear();
2110 }
2111 
2112 template<class ValueType, bool RefTypes, int Tag>
2114 {
2115  std::ostringstream ostream;
2116 
2117  Array<Element*> tmp_elements;
2118  tmp_elements.MakeRef(elements.data(), elements.size());
2119 
2120  ElementSet eset(pncmesh, RefTypes);
2121  eset.Encode(tmp_elements);
2122  eset.Dump(ostream);
2123 
2124  // decode the element set to obtain a local numbering of elements
2125  Array<Element*> decoded;
2126  eset.Decode(decoded);
2127 
2128  std::map<Element*, int> element_index;
2129  for (int i = 0; i < decoded.Size(); i++)
2130  {
2131  element_index[decoded[i]] = i;
2132  }
2133 
2134  write<int>(ostream, values.size());
2135  MFEM_ASSERT(elements.size() == values.size(), "");
2136 
2137  for (unsigned i = 0; i < values.size(); i++)
2138  {
2139  write<int>(ostream, element_index[elements[i]]); // element number
2140  write<ValueType>(ostream, values[i]);
2141  }
2142 
2143  ostream.str().swap(data);
2144 }
2145 
2146 template<class ValueType, bool RefTypes, int Tag>
2148 {
2149  std::istringstream istream(data);
2150 
2151  ElementSet eset(pncmesh, RefTypes);
2152  eset.Load(istream);
2153 
2154  Array<Element*> tmp_elements;
2155  eset.Decode(tmp_elements);
2156 
2157  Element** el = tmp_elements.GetData();
2158  elements.assign(el, el + tmp_elements.Size());
2159  values.resize(elements.size());
2160 
2161  int count = read<int>(istream);
2162  for (int i = 0; i < count; i++)
2163  {
2164  int index = read<int>(istream);
2165  MFEM_ASSERT(index >= 0 && (size_t)index < values.size(), "");
2166  values[index] = read<ValueType>(istream);
2167  }
2168 
2169  // no longer need the raw data
2170  data.clear();
2171 }
2172 
2174  NCMesh *ncmesh)
2175 {
2176  eset.SetNCMesh(ncmesh);
2177  eset.Encode(elems);
2178 
2179  Array<Element*> decoded;
2180  decoded.Reserve(elems.Size());
2181  eset.Decode(decoded);
2182 
2183  elem_ids.resize(decoded.Size());
2184  for (int i = 0; i < decoded.Size(); i++)
2185  {
2186  elem_ids[i] = decoded[i]->index;
2187  }
2188 }
2189 
2191 {
2192  std::ostringstream stream;
2193 
2194  eset.Dump(stream);
2195  write<long>(stream, dof_offset);
2196  write_dofs(stream, dofs);
2197 
2198  stream.str().swap(data);
2199 }
2200 
2202 {
2203  std::istringstream stream(data);
2204 
2205  eset.Load(stream);
2206  dof_offset = read<long>(stream);
2207  read_dofs(stream, dofs);
2208 
2209  data.clear();
2210 
2211  Array<Element*> elems;
2212  eset.Decode(elems);
2213 
2214  elem_ids.resize(elems.Size());
2215  for (int i = 0; i < elems.Size(); i++)
2216  {
2217  elem_ids[i] = elems[i]->index;
2218  }
2219 }
2220 
2221 
2223 
2224 void ParNCMesh::GetDebugMesh(Mesh &debug_mesh) const
2225 {
2226  // create a serial NCMesh containing all our elements (ghosts and all)
2227  NCMesh* copy = new NCMesh(*this);
2228 
2229  Array<Element*> &cle = copy->leaf_elements;
2230  for (int i = 0; i < cle.Size(); i++)
2231  {
2232  cle[i]->attribute = cle[i]->rank + 1;
2233  }
2234 
2235  debug_mesh.InitFromNCMesh(*copy);
2236  debug_mesh.SetAttributes();
2237  debug_mesh.ncmesh = copy;
2238 }
2239 
2240 
2241 } // namespace mfem
2242 
2243 #endif // MFEM_USE_MPI
NCList face_list
lazy-initialized list of faces, see GetFaceList
Definition: ncmesh.hpp:434
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
Definition: ncmesh.hpp:435
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
Definition: ncmesh.cpp:68
NCList shared_edges
Definition: pncmesh.hpp:245
int Partition(long index, long total_elements) const
Return the processor number for a global element number.
Definition: pncmesh.hpp:274
ElementSet(NCMesh *ncmesh=NULL, bool include_ref_types=false)
Definition: pncmesh.hpp:323
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:1547
void DecodeMeshIds(std::istream &is, Array< MeshId > ids[])
Read from &#39;is&#39; a processor-independent encoding of vertex/edge/face IDs.
Definition: pncmesh.cpp:1826
Array< Element * > ghost_layer
list of elements whose &#39;element_type&#39; == 2.
Definition: pncmesh.hpp:268
void AddElementRank(Element *elem, int rank)
Definition: pncmesh.hpp:442
Array< char > element_type
Definition: pncmesh.hpp:266
int Size() const
Logical size of the array.
Definition: array.hpp:109
std::vector< ValueType > values
Definition: pncmesh.hpp:386
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:282
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
Definition: ncmesh.hpp:668
void FindNeighbors(const Element *elem, Array< Element * > &neighbors, const Array< Element * > *search_set=NULL)
Definition: ncmesh.cpp:2280
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:36
void CalcFaceOrientations()
Definition: pncmesh.cpp:439
void Unique()
Definition: array.hpp:197
int index
edge number in the Mesh
Definition: ncmesh.hpp:322
Helper struct to convert a C++ type to an MPI type.
void AddDofs(int type, const NCMesh::MeshId &id, const Array< int > &dofs)
Add vertex/edge/face DOFs to an outgoing message.
Definition: pncmesh.cpp:1884
char flag
generic flag/marker, can be used by algorithms
Definition: ncmesh.hpp:395
Array< char > face_orient
Definition: pncmesh.hpp:258
void WriteInt(int value)
Definition: pncmesh.cpp:1625
Iterator over items contained in the HashTable.
Definition: hash.hpp:160
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
MPI_Comm MyComm
Definition: pncmesh.hpp:235
virtual void BuildFaceList()
Definition: pncmesh.cpp:238
void Dump(std::ostream &os) const
Definition: pncmesh.cpp:1769
void SetSize(int s)
Resize the vector if the new size is different.
Definition: vector.hpp:263
int index
element number in the Mesh, -1 if refined
Definition: ncmesh.hpp:396
Array< Element * > face_nbr_elements
Definition: pmesh.hpp:107
Array< int > face_owner
Definition: pncmesh.hpp:251
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:167
void SetNCMesh(ParNCMesh *pncmesh)
Set pointer to ParNCMesh (needed to encode the message).
Definition: pncmesh.hpp:395
std::map< int, NeighborElementRankMessage > Map
Definition: pncmesh.hpp:432
void BuildSharedVertices()
Definition: pncmesh.cpp:366
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:160
T * GetData()
Returns the data.
Definition: array.hpp:91
void Load(std::istream &is)
Definition: pncmesh.cpp:1775
Table derefinements
possible derefinements, see GetDerefinementTable
Definition: ncmesh.hpp:467
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
void InitDerefTransforms()
Definition: ncmesh.cpp:1436
int Size() const
Returns the size of the vector.
Definition: vector.hpp:86
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
char geom
Geometry::Type of the element.
Definition: ncmesh.hpp:393
static int find_node(Element *elem, Node *node)
Definition: ncmesh.cpp:1745
void DecodeTree(Element *elem, int &pos, Array< Element * > &elements) const
Definition: pncmesh.cpp:1714
void SetDerefMatrixCodes(Element *parent, Array< Element * > &coarse)
Definition: ncmesh.cpp:1451
void SynchronizeDerefinementData(Array< Type > &elem_data, const Table &deref_table)
Definition: pncmesh.cpp:1173
int index
vertex number in the Mesh
Definition: ncmesh.hpp:309
void NeighborExpand(const Array< Element * > &elements, Array< Element * > &expanded, const Array< Element * > *search_set=NULL)
Definition: ncmesh.cpp:2342
int GetInt(int pos) const
Definition: pncmesh.cpp:1634
virtual void OnMeshUpdated(Mesh *mesh)
Definition: pncmesh.cpp:147
virtual void Derefine(const Array< int > &derefs)
Definition: pncmesh.cpp:1001
std::vector< MeshId > conforming
Definition: ncmesh.hpp:169
T * end() const
Definition: array.hpp:219
Array< int > face_nbr_group
Definition: pmesh.hpp:104
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:3212
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
Definition: ncmesh.cpp:3075
T * begin() const
Definition: array.hpp:218
Element * parent
parent element, NULL if this is a root element
Definition: ncmesh.hpp:404
double * GetData() const
Definition: vector.hpp:90
virtual void AssignLeafIndices()
Definition: pncmesh.cpp:67
const NCList & GetSharedEdges()
Definition: pncmesh.hpp:106
void MakeShared(const Table &groups, const NCList &list, NCList &shared)
Definition: pncmesh.cpp:338
Array< int > vertex_nodeId
Definition: ncmesh.hpp:432
int index
Mesh element number.
Definition: ncmesh.hpp:35
int master
master number (in Mesh numbering)
Definition: ncmesh.hpp:155
Element * element
NCMesh::Element containing this vertex/edge/face.
Definition: ncmesh.hpp:135
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:181
Array< NCFaceInfo > nc_faces_info
Definition: mesh.hpp:99
virtual void OnMeshUpdated(Mesh *mesh)
Definition: ncmesh.cpp:1671
virtual void Decode()
Definition: pncmesh.cpp:2037
std::vector< int > cols
Definition: pncmesh.hpp:568
ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh)
Definition: pncmesh.cpp:27
Array< int > vertex_owner
Definition: pncmesh.hpp:249
RebalanceDofMessage::Map send_rebalance_dofs
Definition: pncmesh.hpp:475
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3039
virtual void Update()
Definition: pncmesh.cpp:54
virtual void Decode()
Definition: pncmesh.cpp:1996
Table face_group
Definition: pncmesh.hpp:256
Array< int > face_nbr_elements_offset
Definition: pmesh.hpp:105
bool PruneTree(Element *elem)
Internal. Recursive part of Prune().
Definition: pncmesh.cpp:841
int index
face number in the Mesh
Definition: ncmesh.hpp:369
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:264
std::vector< Element * > elements
Definition: pncmesh.hpp:385
Array< DenseMatrix * > aux_pm_store
Stores modified point matrices created by GetFaceNeighbors.
Definition: pncmesh.hpp:484
int local
local number within &#39;element&#39;
Definition: ncmesh.hpp:134
void SetAttributes()
Definition: mesh.cpp:768
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:376
void RedistributeElements(Array< int > &new_ranks, int target_elements, bool record_comm)
Definition: pncmesh.cpp:1361
Element * elem[2]
up to 2 elements sharing the face
Definition: ncmesh.hpp:370
std::vector< Master > masters
Definition: ncmesh.hpp:170
void FlagElements(const Array< Element * > &elements, char flag)
Definition: pncmesh.cpp:1643
Array< int > old_index_or_rank
Definition: pncmesh.hpp:481
RebalanceDofMessage::Map recv_rebalance_dofs
Definition: pncmesh.hpp:476
int InitialPartition(int index) const
Helper to get the partition when the serial mesh is being split initially.
Definition: pncmesh.hpp:278
void EncodeTree(Element *elem)
Definition: pncmesh.cpp:1657
void RefineElement(Element *elem, char ref_type)
Definition: ncmesh.cpp:756
Array< Element * > boundary_layer
list of type 3 elements
Definition: pncmesh.hpp:269
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:122
Array< Embedding > embeddings
fine element positions in their parents
Definition: ncmesh.hpp:54
NCList shared_vertices
Definition: pncmesh.hpp:244
virtual void BuildEdgeList()
Definition: ncmesh.cpp:1987
Identifies a vertex/edge/face in both Mesh and NCMesh.
Definition: ncmesh.hpp:131
void Assign(const T *)
Copy data from a pointer. Size() elements are copied.
Definition: array.hpp:518
std::map< int, NeighborRefinementMessage > Map
Definition: pncmesh.hpp:413
static void IsendAll(MapT &rank_msg, MPI_Comm comm)
Helper to send all messages in a rank-to-message map container.
void Sort()
Sorts the array. This requires operator&lt; to be defined for T.
Definition: array.hpp:189
Table edge_group
Definition: pncmesh.hpp:255
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:151
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
Definition: ncmesh.hpp:79
void ClearAuxPM()
Definition: pncmesh.cpp:829
NCList shared_faces
Definition: pncmesh.hpp:246
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:2932
void SetElements(const Array< Element * > &elems, NCMesh *ncmesh)
Definition: pncmesh.cpp:2173
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: pncmesh.cpp:463
void Rebalance()
Definition: pncmesh.cpp:1309
static GeomInfo GI[Geometry::NumGeom]
Definition: ncmesh.hpp:712
int slaves_end
slave faces
Definition: ncmesh.hpp:145
Vertex * vertex
Definition: ncmesh.hpp:342
bool Iso
true if the mesh only contains isotropic refinements
Definition: ncmesh.hpp:277
Table vertex_group
Definition: pncmesh.hpp:254
void DerefineElement(Element *elem)
Definition: ncmesh.cpp:1211
std::map< int, NeighborDerefinementMessage > Map
Definition: pncmesh.hpp:422
virtual void BuildEdgeList()
Definition: pncmesh.cpp:214
static bool compare_ranks(const Element *a, const Element *b)
Definition: pncmesh.cpp:1356
void FindSetNeighbors(const Array< char > &elem_set, Array< Element * > *neighbors, Array< char > *neighbor_set=NULL)
Definition: ncmesh.cpp:2196
Array< Vertex > face_nbr_vertices
Definition: pmesh.hpp:108
std::vector< Slave > slaves
Definition: ncmesh.hpp:171
virtual void BuildFaceList()
Definition: ncmesh.cpp:1876
void GetDebugMesh(Mesh &debug_mesh) const
Definition: pncmesh.cpp:2224
int Find(const T &el) const
Return the first index where &#39;el&#39; is found; return -1 if not found.
Definition: array.hpp:437
void GetDofs(int type, const NCMesh::MeshId &id, Array< int > &dofs, int &ndofs)
Definition: pncmesh.cpp:1891
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
Element * child[8]
2-8 children (if ref_type != 0)
Definition: ncmesh.hpp:402
Face * GetFace(Element *elem, int face_no)
Definition: ncmesh.cpp:391
Array< Element * > leaf_elements
Definition: ncmesh.hpp:430
mfem::Element * NewMeshElement(int geom) const
Definition: ncmesh.cpp:1583
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:24
virtual void Encode()
Definition: pncmesh.cpp:2071
void UpdateLayers()
Definition: pncmesh.cpp:488
void ElementNeighborProcessors(Element *elem, Array< int > &ranks)
Definition: pncmesh.cpp:543
virtual void ElementSharesFace(Element *elem, Face *face)
Definition: pncmesh.cpp:204
virtual void ElementSharesEdge(Element *elem, Edge *edge)
Definition: pncmesh.cpp:192
HashTable< Node > nodes
Definition: ncmesh.hpp:416
void Decode(Array< Element * > &elements) const
Definition: pncmesh.cpp:1745
bool CheckElementType(Element *elem, int type)
Definition: pncmesh.cpp:526
virtual void Decode()
Definition: pncmesh.cpp:2091
virtual void Refine(const Array< Refinement > &refinements)
Definition: pncmesh.cpp:911
virtual ~ParNCMesh()
Definition: pncmesh.cpp:49
static bool compare_ranks_indices(const Element *a, const Element *b)
Definition: pncmesh.cpp:588
virtual void Encode()
Definition: pncmesh.cpp:1959
void RecvRebalanceDofs(Array< int > &elements, Array< long > &dofs)
Receive element DOFs sent by SendRebalanceDofs().
Definition: pncmesh.cpp:1580
void GetRow(int row, Array< int > &cols, Vector &srow)
Definition: pncmesh.cpp:2061
void Isend(int rank, MPI_Comm comm)
Non-blocking send to processor &#39;rank&#39;.
Array< Connection > index_rank
Definition: pncmesh.hpp:309
Array< int > edge_owner
Definition: pncmesh.hpp:250
Node * node[8]
element corners (if ref_type == 0)
Definition: ncmesh.hpp:401
Array< FaceInfo > faces_info
Definition: mesh.hpp:98
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:53
virtual void AssignLeafIndices()
Definition: ncmesh.cpp:1574
static int find_hex_face(int a, int b, int c)
Definition: ncmesh.cpp:1779
NCMesh * ncmesh
Definition: mesh.hpp:144
T & Last()
Return the last element in the array.
Definition: array.hpp:411
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
Definition: pncmesh.cpp:1253
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:502
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:474
Array< Element * > root_elements
Definition: ncmesh.hpp:414
void AddElementRank(Element *elem, int rank)
Definition: pncmesh.hpp:431
virtual void Update()
Definition: ncmesh.cpp:187
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[])
Write to &#39;os&#39; a processor-independent encoding of vertex/edge/face IDs.
Definition: pncmesh.cpp:1784
Vector data type.
Definition: vector.hpp:33
void AddRow(int row, const Array< int > &cols, const Vector &srow)
Definition: pncmesh.cpp:2052
void AddMasterSlaveRanks(int nitems, const NCList &list)
Definition: pncmesh.cpp:271
void Encode(const Array< Element * > &elements)
Definition: pncmesh.cpp:1693
const NCList & GetSharedFaces()
Definition: pncmesh.hpp:114
static int find_element_edge(Element *elem, int v0, int v1)
Definition: ncmesh.cpp:1763
std::map< int, RebalanceMessage > Map
Definition: pncmesh.hpp:443
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition: ncmesh.hpp:188
void ReorderEdgeDofs(const NCMesh::MeshId &id, std::vector< int > &dofs)
Definition: pncmesh.cpp:1908
Array< Element * > tmp_neighbors
Definition: pncmesh.hpp:356
void NeighborProcessors(Array< int > &neighbors)
Definition: pncmesh.cpp:576
int RowSize(int i) const
Definition: table.hpp:102
void CountSplits(Element *elem, int splits[3]) const
Definition: ncmesh.cpp:3028
Table send_face_nbr_elements
Definition: pmesh.hpp:110
static int get_face_orientation(Face *face, Element *e1, Element *e2, int local[2]=NULL)
Definition: pncmesh.cpp:416
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
Definition: ncmesh.hpp:394
virtual void Encode()
Definition: pncmesh.cpp:2022
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
Definition: mesh.cpp:5864
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
Definition: pncmesh.cpp:101
int index
Mesh number.
Definition: ncmesh.hpp:133
Class for parallel meshes.
Definition: pmesh.hpp:28
Abstract data type element.
Definition: element.hpp:27
int rank
processor number (ParNCMesh)
Definition: ncmesh.hpp:397
virtual void LimitNCLevel(int max_nc_level)
Parallel version of NCMesh::LimitNCLevel.
Definition: pncmesh.cpp:983
Array< unsigned char > data
encoded refinement (sub-)trees
Definition: pncmesh.hpp:336
DenseMatrix point_matrix
position within the master edge/face
Definition: ncmesh.hpp:157
static void Probe(int &rank, int &size, MPI_Comm comm)
HashTable< Face > faces
Definition: ncmesh.hpp:417
void GetFaceNeighbors(ParMesh &pmesh)
Definition: pncmesh.cpp:594
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:67