MFEM  v3.1
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 
20 #include <map>
21 #include <limits>
22 
23 namespace mfem
24 {
25 
26 ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh)
27  : NCMesh(ncmesh)
28 {
29  MyComm = comm;
30  MPI_Comm_size(MyComm, &NRanks);
31  MPI_Comm_rank(MyComm, &MyRank);
32 
33  // assign leaf elements to the 'NRanks' processors
34  for (int i = 0; i < leaf_elements.Size(); i++)
35  {
36  leaf_elements[i]->rank = InitialPartition(i);
37  }
38 
41 }
42 
44 {
46 
50 
52  ghost_layer.SetSize(0);
53 }
54 
56 {
57  // This is an override of NCMesh::AssignLeafIndices(). The difference is
58  // that we shift all elements we own to the beginning of the array
59  // 'leaf_elements' and assign all ghost elements indices >= NElements. This
60  // will make the ghosts skipped in NCMesh::GetMeshComponents.
61 
63  for (int i = 0; i < leaf_elements.Size(); i++)
64  {
65  if (leaf_elements[i]->rank == MyRank)
66  {
67  std::swap(leaf_elements[NElements++], leaf_elements[i]);
68  }
69  else
70  {
72  }
73  }
74 
75  for (int i = 0; i < leaf_elements.Size(); i++)
76  {
77  leaf_elements[i]->index = i;
78  }
79 }
80 
82 {
83  // This is an override of NCMesh::UpdateVertices. This version first
84  // assigns Vertex::index to vertices of elements of our rank. Only these
85  // vertices then make it to the Mesh in NCMesh::GetMeshComponents.
86  // The remaining (ghost) vertices are assigned indices greater or equal to
87  // Mesh::GetNV().
88 
89  for (HashTable<Node>::Iterator it(nodes); it; ++it)
90  {
91  if (it->vertex) { it->vertex->index = -1; }
92  }
93 
94  NVertices = 0;
95  for (int i = 0; i < leaf_elements.Size(); i++)
96  {
97  Element* elem = leaf_elements[i];
98  if (elem->rank == MyRank)
99  {
100  for (int j = 0; j < GI[(int) elem->geom].nv; j++)
101  {
102  int &vindex = elem->node[j]->vertex->index;
103  if (vindex < 0) { vindex = NVertices++; }
104  }
105  }
106  }
107 
109  for (HashTable<Node>::Iterator it(nodes); it; ++it)
110  {
111  if (it->vertex && it->vertex->index >= 0)
112  {
113  vertex_nodeId[it->vertex->index] = it->id;
114  }
115  }
116 
117  NGhostVertices = 0;
118  for (HashTable<Node>::Iterator it(nodes); it; ++it)
119  {
120  if (it->vertex && it->vertex->index < 0)
121  {
122  it->vertex->index = NVertices + (NGhostVertices++);
123  }
124  }
125 }
126 
128 {
129  // This is an override (or extension of) NCMesh::OnMeshUpdated().
130  // In addition to getting edge/face indices from 'mesh', we also
131  // assign indices to ghost edges/faces that don't exist in the 'mesh'.
132 
133  // clear Edge:: and Face::index
134  for (HashTable<Node>::Iterator it(nodes); it; ++it)
135  {
136  if (it->edge) { it->edge->index = -1; }
137  }
138  for (HashTable<Face>::Iterator it(faces); it; ++it) { it->index = -1; }
139 
140  // go assign existing edge/face indices
141  NCMesh::OnMeshUpdated(mesh);
142 
143  // assign ghost edge indices
144  NEdges = mesh->GetNEdges();
145  NGhostEdges = 0;
146  for (HashTable<Node>::Iterator it(nodes); it; ++it)
147  {
148  if (it->edge && it->edge->index < 0)
149  {
150  it->edge->index = NEdges + (NGhostEdges++);
151  }
152  }
153 
154  // assign ghost face indices
155  NFaces = mesh->GetNFaces();
156  NGhostFaces = 0;
157  for (HashTable<Face>::Iterator it(faces); it; ++it)
158  {
159  if (it->index < 0) { it->index = NFaces + (NGhostFaces++); }
160  }
161 }
162 
164 {
165  // Called by NCMesh::BuildEdgeList when an edge is visited in a leaf element.
166  // This allows us to determine edge ownership and processors that share it
167  // without duplicating all the HashTable lookups in NCMesh::BuildEdgeList().
168 
169  int &owner = edge_owner[edge->index];
170  owner = std::min(owner, elem->rank);
171 
172  index_rank.Append(Connection(edge->index, elem->rank));
173 }
174 
176 {
177  // Analogous to ElementHasEdge.
178 
179  int &owner = face_owner[face->index];
180  owner = std::min(owner, elem->rank);
181 
182  index_rank.Append(Connection(face->index, elem->rank));
183 }
184 
186 {
187  // This is an extension of NCMesh::BuildEdgeList() which also determines
188  // edge ownership, creates edge processor groups and lists shared edges.
189 
190  int nedges = NEdges + NGhostEdges;
191  edge_owner.SetSize(nedges);
192  edge_owner = std::numeric_limits<int>::max();
193 
194  index_rank.SetSize(12*leaf_elements.Size() * 3/2);
195  index_rank.SetSize(0);
196 
198 
200 
201  index_rank.Sort();
202  index_rank.Unique();
204  index_rank.DeleteAll();
205 
207 }
208 
210 {
211  // This is an extension of NCMesh::BuildFaceList() which also determines
212  // face ownership, creates face processor groups and lists shared faces.
213 
214  int nfaces = NFaces + NGhostFaces;
215  face_owner.SetSize(nfaces);
216  face_owner = std::numeric_limits<int>::max();
217 
218  index_rank.SetSize(6*leaf_elements.Size() * 3/2);
219  index_rank.SetSize(0);
220 
222 
224 
225  index_rank.Sort();
226  index_rank.Unique();
228  index_rank.DeleteAll();
229 
231 
233 }
234 
235 struct MasterSlaveInfo
236 {
237  int master; // master index if this is a slave
238  int slaves_begin, slaves_end; // slave list if this is a master
239  MasterSlaveInfo() : master(-1), slaves_begin(0), slaves_end(0) {}
240 };
241 
242 void ParNCMesh::AddMasterSlaveRanks(int nitems, const NCList& list)
243 {
244  // create an auxiliary structure for each edge/face
245  std::vector<MasterSlaveInfo> info(nitems);
246 
247  for (unsigned i = 0; i < list.masters.size(); i++)
248  {
249  const Master &mf = list.masters[i];
250  info[mf.index].slaves_begin = mf.slaves_begin;
251  info[mf.index].slaves_end = mf.slaves_end;
252  }
253  for (unsigned i = 0; i < list.slaves.size(); i++)
254  {
255  const Slave& sf = list.slaves[i];
256  info[sf.index].master = sf.master;
257  }
258 
259  // We need the processor groups of master edges/faces to contain the ranks of
260  // their slaves (so that master DOFs get sent to those who share the slaves).
261  // Conversely, we need the groups of slave edges/faces to contain the ranks
262  // of their masters. Both can be done by appending more items to the
263  // 'index_rank' array, before it is sorted and converted to the group table.
264  // (Note that a master/slave edge can be shared by more than one processor.)
265 
266  int size = index_rank.Size();
267  for (int i = 0; i < size; i++)
268  {
269  int index = index_rank[i].from;
270  int rank = index_rank[i].to;
271 
272  const MasterSlaveInfo &msi = info[index];
273  if (msi.master >= 0)
274  {
275  // 'index' is a slave, add its rank to the master's group
276  index_rank.Append(Connection(msi.master, rank));
277  }
278  else
279  {
280  for (int j = msi.slaves_begin; j < msi.slaves_end; j++)
281  {
282  // 'index' is a master, add its rank to the groups of the slaves
283  index_rank.Append(Connection(list.slaves[j].index, rank));
284  }
285  }
286  }
287 }
288 
289 static bool is_shared(const Table& groups, int index, int MyRank)
290 {
291  // A vertex/edge/face is shared if its group contains more than one processor
292  // and at the same time one of them is ourselves.
293 
294  int size = groups.RowSize(index);
295  if (size <= 1)
296  {
297  return false;
298  }
299 
300  const int* group = groups.GetRow(index);
301  for (int i = 0; i < size; i++)
302  {
303  if (group[i] == MyRank) { return true; }
304  }
305 
306  return false;
307 }
308 
309 void ParNCMesh::MakeShared(const Table &groups, const NCList &list,
310  NCList &shared)
311 {
312  shared.Clear();
313 
314  for (unsigned i = 0; i < list.conforming.size(); i++)
315  {
316  if (is_shared(groups, list.conforming[i].index, MyRank))
317  {
318  shared.conforming.push_back(list.conforming[i]);
319  }
320  }
321  for (unsigned i = 0; i < list.masters.size(); i++)
322  {
323  if (is_shared(groups, list.masters[i].index, MyRank))
324  {
325  shared.masters.push_back(list.masters[i]);
326  }
327  }
328  for (unsigned i = 0; i < list.slaves.size(); i++)
329  {
330  if (is_shared(groups, list.slaves[i].index, MyRank))
331  {
332  shared.slaves.push_back(list.slaves[i]);
333  }
334  }
335 }
336 
338 {
339  int nvertices = NVertices + NGhostVertices;
340  vertex_owner.SetSize(nvertices);
341  vertex_owner = std::numeric_limits<int>::max();
342 
343  index_rank.SetSize(8*leaf_elements.Size());
344  index_rank.SetSize(0);
345 
346  Array<MeshId> vertex_id(nvertices);
347 
348  // similarly to edges/faces, we loop over the vertices of all leaf elements
349  // to determine which processors share each vertex
350  for (int i = 0; i < leaf_elements.Size(); i++)
351  {
352  Element* elem = leaf_elements[i];
353  for (int j = 0; j < GI[(int) elem->geom].nv; j++)
354  {
355  Node* node = elem->node[j];
356  int index = node->vertex->index;
357 
358  int &owner = vertex_owner[index];
359  owner = std::min(owner, elem->rank);
360 
361  index_rank.Append(Connection(index, elem->rank));
362 
363  MeshId &id = vertex_id[index];
364  id.index = (node->edge ? -1 : index);
365  id.element = elem;
366  id.local = j;
367  }
368  }
369 
370  index_rank.Sort();
371  index_rank.Unique();
373  index_rank.DeleteAll();
374 
375  // create a list of shared vertices, skip obviously slave vertices
376  // (for simplicity, we don't guarantee to skip all slave vertices)
378  for (int i = 0; i < nvertices; i++)
379  {
380  if (is_shared(vertex_group, i, MyRank) && vertex_id[i].index >= 0)
381  {
382  shared_vertices.conforming.push_back(vertex_id[i]);
383  }
384  }
385 }
386 
388 {
389  // Calculate orientation of shared conforming faces.
390  // NOTE: face orientation is calculated relative to its lower rank element.
391  // Thanks to the ghost layer this can be done locally, without communication.
392 
394  face_orient = 0;
395 
396  for (HashTable<Face>::Iterator it(faces); it; ++it)
397  {
398  if (it->ref_count == 2 && it->index < NFaces)
399  {
400  Element* e[2] = { it->elem[0], it->elem[1] };
401  if (e[0]->rank == e[1]->rank) { continue; }
402  if (e[0]->rank > e[1]->rank) { std::swap(e[0], e[1]); }
403 
404  int ids[2][4];
405  for (int i = 0; i < 2; i++)
406  {
407  int f = find_hex_face(find_node(e[i], it->p1),
408  find_node(e[i], it->p2),
409  find_node(e[i], it->p3));
410 
411  // get node IDs for the face as seen from e[i]
412  const int* fv = GI[Geometry::CUBE].faces[f];
413  for (int j = 0; j < 4; j++)
414  {
415  ids[i][j] = e[i]->node[fv[j]]->id;
416  }
417  }
418 
419  face_orient[it->index] = Mesh::GetQuadOrientation(ids[0], ids[1]);
420  }
421  }
422 }
423 
424 void ParNCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
425  Array<int> &bdr_vertices,
426  Array<int> &bdr_edges)
427 {
428  NCMesh::GetBoundaryClosure(bdr_attr_is_ess, bdr_vertices, bdr_edges);
429 
430  int i, j;
431  // filter out ghost vertices
432  for (i = j = 0; i < bdr_vertices.Size(); i++)
433  {
434  if (bdr_vertices[i] < NVertices) { bdr_vertices[j++] = bdr_vertices[i]; }
435  }
436  bdr_vertices.SetSize(j);
437 
438  // filter out ghost edges
439  for (i = j = 0; i < bdr_edges.Size(); i++)
440  {
441  if (bdr_edges[i] < NEdges) { bdr_edges[j++] = bdr_edges[i]; }
442  }
443  bdr_edges.SetSize(j);
444 }
445 
447 
449 {
450  if (element_type.Size()) { return; }
451 
452  int nleaves = leaf_elements.Size();
453 
454  element_type.SetSize(nleaves);
455  for (int i = 0; i < nleaves; i++)
456  {
457  element_type[i] = (leaf_elements[i]->rank == MyRank) ? 1 : 0;
458  }
459 
460  // determine the ghost layer
461  Array<char> ghost_set;
462  FindSetNeighbors(element_type, NULL, &ghost_set);
463 
464  // find the neighbors of the ghost layer
465  Array<char> boundary_set;
466  FindSetNeighbors(ghost_set, NULL, &boundary_set);
467 
468  ghost_layer.SetSize(0);
469  for (int i = 0; i < nleaves; i++)
470  {
471  if (ghost_set[i])
472  {
473  element_type[i] = 2;
474  ghost_layer.Append(leaf_elements[i]);
475  }
476  else if (boundary_set[i] && element_type[i])
477  {
478  element_type[i] = 3;
479  }
480  }
481 }
482 
484  Array<int> &ranks)
485 {
486  MFEM_ASSERT(!elem->ref_type, "not a leaf.");
487 
488  ranks.SetSize(0); // preserve capacity
489 
490  // big shortcut: there are no neighbors if element_type == 1
491  if (element_type[elem->index] == 1) { return; }
492 
493  // ok, we do need to look for neighbors;
494  // at least we can only search in the ghost layer
495  tmp_neighbors.SetSize(0);
497 
498  // return a list of processors
499  for (int i = 0; i < tmp_neighbors.Size(); i++)
500  {
501  ranks.Append(tmp_neighbors[i]->rank);
502  }
503  ranks.Sort();
504  ranks.Unique();
505 }
506 
508 {
509  UpdateLayers();
510 
511  std::set<int> ranks;
512  for (int i = 0; i < ghost_layer.Size(); i++)
513  {
514  ranks.insert(ghost_layer[i]->rank);
515  }
516 
517  neighbors.DeleteAll();
518  neighbors.Reserve(ranks.size());
519 
520  std::set<int>::iterator it;
521  for (it = ranks.begin(); it != ranks.end(); ++it)
522  {
523  neighbors.Append(*it);
524  }
525 }
526 
527 
529 
531 {
532  if (elem->ref_type)
533  {
534  bool remove[8];
535  bool removeAll = true;
536 
537  // determine which subtrees can be removed (and whether it's all of them)
538  for (int i = 0; i < 8; i++)
539  {
540  remove[i] = false;
541  if (elem->child[i])
542  {
543  remove[i] = PruneTree(elem->child[i]);
544  if (!remove[i]) { removeAll = false; }
545  }
546  }
547 
548  // all children can be removed, let the (maybe indirect) parent do it
549  if (removeAll) { return true; }
550 
551  // not all children can be removed, but remove those that can be
552  for (int i = 0; i < 8; i++)
553  {
554  if (remove[i]) { DerefineElement(elem->child[i]); }
555  }
556 
557  return false; // need to keep this element and up
558  }
559  else
560  {
561  // return true if this leaf can be removed
562  return element_type[elem->index] == 0;
563  }
564 }
565 
567 {
568  if (!Iso)
569  {
570  MFEM_WARNING("Can't prune aniso meshes yet.");
571  return;
572  }
573 
574  UpdateLayers();
575 
576  // derefine subtrees whose leaves are all unneeded
577  for (int i = 0; i < root_elements.Size(); i++)
578  {
579  if (PruneTree(root_elements[i]))
580  {
582  }
583  }
584 
585  Update();
586 }
587 
588 void ParNCMesh::Refine(const Array<Refinement> &refinements)
589 {
590  for (int i = 0; i < refinements.Size(); i++)
591  {
592  const Refinement &ref = refinements[i];
593  MFEM_VERIFY(ref.ref_type == 7 || Dim < 3,
594  "anisotropic parallel refinement not supported yet in 3D.");
595  }
596  MFEM_VERIFY(Iso || Dim < 3,
597  "parallel refinement of 3D aniso meshes not supported yet.");
598 
600 
601  // create refinement messages to all neighbors (NOTE: some may be empty)
602  Array<int> neighbors;
603  NeighborProcessors(neighbors);
604  for (int i = 0; i < neighbors.Size(); i++)
605  {
606  send_ref[neighbors[i]].SetNCMesh(this);
607  }
608 
609  // populate messages: all refinements that occur next to the processor
610  // boundary need to be sent to the adjoining neighbors so they can keep
611  // their ghost layer up to date
612  Array<int> ranks;
613  ranks.Reserve(64);
614  for (int i = 0; i < refinements.Size(); i++)
615  {
616  const Refinement &ref = refinements[i];
617  MFEM_ASSERT(ref.index < NElements, "");
618  Element* elem = leaf_elements[ref.index];
619  ElementNeighborProcessors(elem, ranks);
620  for (int j = 0; j < ranks.Size(); j++)
621  {
622  send_ref[ranks[j]].AddRefinement(elem, ref.ref_type);
623  }
624  }
625 
626  // send the messages (overlap with local refinements)
628 
629  // do local refinements
630  for (int i = 0; i < refinements.Size(); i++)
631  {
632  const Refinement &ref = refinements[i];
634  }
635 
636  // receive (ghost layer) refinements from all neighbors
637  for (int j = 0; j < neighbors.Size(); j++)
638  {
639  int rank, size;
641 
643  msg.SetNCMesh(this);
644  msg.Recv(rank, size, MyComm);
645 
646  // do the ghost refinements
647  for (unsigned i = 0; i < msg.refinements.size(); i++)
648  {
649  const ElemRefType &ref = msg.refinements[i];
651  }
652  }
653 
654  Update();
655 
656  // make sure we can delete the send buffers
658 }
659 
660 
661 void ParNCMesh::LimitNCLevel(int max_level)
662 {
663  if (NRanks > 1)
664  {
665  MFEM_ABORT("not implemented in parallel yet.");
666  }
667  NCMesh::LimitNCLevel(max_level);
668 }
669 
670 
672 
673 void ParNCMesh::ElementSet::SetInt(int pos, int value)
674 {
675  // helper to put an int to the data array
676  data[pos] = value & 0xff;
677  data[pos+1] = (value >> 8) & 0xff;
678  data[pos+2] = (value >> 16) & 0xff;
679  data[pos+3] = (value >> 24) & 0xff;
680 }
681 
683 {
684  // helper to get an int from the data array
685  return (int) data[pos] +
686  ((int) data[pos+1] << 8) +
687  ((int) data[pos+2] << 16) +
688  ((int) data[pos+3] << 24);
689 }
690 
692 (Element* elem, const std::set<Element*> &elements)
693 {
694  // is 'elem' in the set?
695  if (elements.find(elem) != elements.end())
696  {
697  // we reached a 'leaf' of our subtree, mark this as zero child mask
698  data.Append(0);
699  return true;
700  }
701  else if (elem->ref_type)
702  {
703  // write a bit mask telling what subtrees contain elements from the set
704  int mpos = data.Size();
705  data.Append(0);
706 
707  // check the subtrees
708  int mask = 0;
709  for (int i = 0; i < 8; i++)
710  {
711  if (elem->child[i])
712  {
713  if (EncodeTree(elem->child[i], elements))
714  {
715  mask |= (unsigned char) 1 << i;
716  }
717  }
718  }
719 
720  if (mask)
721  {
722  data[mpos] = mask;
723  }
724  else
725  {
726  data.DeleteLast();
727  }
728 
729  return mask != 0;
730  }
731  return false;
732 }
733 
734 ParNCMesh::ElementSet::ElementSet(const std::set<Element*> &elements,
735  const Array<Element*> &ncmesh_roots)
736 {
737  int header_pos = 0;
738  data.SetSize(4);
739 
740  // Each refinement tree that contains at least one element from the set
741  // is encoded as HEADER + TREE, where HEADER is the root element number and
742  // TREE is the output of EncodeTree().
743  for (int i = 0; i < ncmesh_roots.Size(); i++)
744  {
745  if (EncodeTree(ncmesh_roots[i], elements))
746  {
747  SetInt(header_pos, i);
748 
749  // make room for the next header
750  header_pos = data.Size();
751  data.SetSize(header_pos + 4);
752  }
753  }
754 
755  // mark end of data
756  SetInt(header_pos, -1);
757 }
758 
760  Array<Element*> &elements) const
761 {
762  int mask = data[pos++];
763  if (!mask)
764  {
765  elements.Append(elem);
766  }
767  else
768  {
769  for (int i = 0; i < 8; i++)
770  {
771  if (mask & (1 << i))
772  {
773  DecodeTree(elem->child[i], pos, elements);
774  }
775  }
776  }
777 }
778 
780  const Array<Element*> &ncmesh_roots) const
781 {
782  int root, pos = 0;
783  while ((root = GetInt(pos)) >= 0)
784  {
785  pos += 4;
786  DecodeTree(ncmesh_roots[root], pos, elements);
787  }
788 }
789 
790 template<typename T>
791 static inline void write(std::ostream& os, T value)
792 {
793  os.write((char*) &value, sizeof(T));
794 }
795 
796 template<typename T>
797 static inline T read(std::istream& is)
798 {
799  T value;
800  is.read((char*) &value, sizeof(T));
801  return value;
802 }
803 
804 void ParNCMesh::ElementSet::Dump(std::ostream &os) const
805 {
806  write<int>(os, data.Size());
807  os.write((const char*) data.GetData(), data.Size());
808 }
809 
810 void ParNCMesh::ElementSet::Load(std::istream &is)
811 {
812  data.SetSize(read<int>(is));
813  is.read((char*) data.GetData(), data.Size());
814 }
815 
816 
818 
819 void ParNCMesh::EncodeMeshIds(std::ostream &os, Array<MeshId> ids[],
820  int dim) const
821 {
822  std::map<Element*, int> element_id;
823 
824  // get a list of elements involved, dump them to 'os' and create the mapping
825  // element_id: (Element* -> stream ID)
826  {
827  std::set<Element*> elements;
828  for (int type = 0; type < dim; type++)
829  {
830  for (int i = 0; i < ids[type].Size(); i++)
831  {
832  elements.insert(ids[type][i].element);
833  }
834  }
835 
836  ElementSet eset(elements, root_elements); // encodes 'elements'
837  eset.Dump(os);
838 
839  Array<Element*> decoded;
840  eset.Decode(decoded, root_elements);
841 
842  for (int i = 0; i < decoded.Size(); i++)
843  {
844  element_id[decoded[i]] = i;
845  }
846  }
847 
848  // write the IDs as element/local pairs
849  for (int type = 0; type < dim; type++)
850  {
851  write<int>(os, ids[type].Size());
852  for (int i = 0; i < ids[type].Size(); i++)
853  {
854  const MeshId& id = ids[type][i];
855  write<int>(os, element_id[id.element]); // TODO: variable 1-4 bytes
856  write<char>(os, id.local);
857  }
858  }
859 }
860 
861 void ParNCMesh::DecodeMeshIds(std::istream &is, Array<MeshId> ids[], int dim,
862  bool decode_indices) const
863 {
864  // read the list of elements
865  ElementSet eset(is);
866 
867  Array<Element*> elements;
868  eset.Decode(elements, root_elements);
869 
870  // read vertex/edge/face IDs
871  for (int type = 0; type < dim; type++)
872  {
873  int ne = read<int>(is);
874  ids[type].SetSize(ne);
875 
876  for (int i = 0; i < ne; i++)
877  {
878  int el_num = read<int>(is);
879  Element* elem = elements[el_num];
880  MFEM_VERIFY(!elem->ref_type, "not a leaf element: " << el_num);
881 
882  MeshId &id = ids[type][i];
883  id.element = elem;
884  id.local = read<char>(is);
885 
886  if (!decode_indices) { continue; }
887 
888  // find vertex/edge/face index
889  GeomInfo &gi = GI[(int) elem->geom];
890  switch (type)
891  {
892  case 0:
893  {
894  id.index = elem->node[id.local]->vertex->index;
895  break;
896  }
897  case 1:
898  {
899  const int* ev = gi.edges[id.local];
900  Node* node = nodes.Peek(elem->node[ev[0]], elem->node[ev[1]]);
901  MFEM_ASSERT(node && node->edge, "edge not found.");
902  id.index = node->edge->index;
903  break;
904  }
905  default:
906  {
907  const int* fv = gi.faces[id.local];
908  Face* face = faces.Peek(elem->node[fv[0]], elem->node[fv[1]],
909  elem->node[fv[2]], elem->node[fv[3]]);
910  MFEM_ASSERT(face, "face not found.");
911  id.index = face->index;
912  }
913  }
914  }
915  }
916 }
917 
919 
921  const Array<int> &dofs)
922 {
923  MFEM_ASSERT(type >= 0 && type < 3, "");
924  id_dofs[type][id].assign(dofs.GetData(), dofs.GetData() + dofs.Size());
925 }
926 
928  Array<int>& dofs, int &ndofs)
929 {
930  MFEM_ASSERT(type >= 0 && type < 3, "");
931 #ifdef MFEM_DEBUG
932  if (id_dofs[type].find(id) == id_dofs[type].end())
933  {
934  MFEM_ABORT("type/ID " << type << "/" << id.index << " not found in "
935  "neighbor message. Ghost layers out of sync?");
936  }
937 #endif
938  std::vector<int> &vec = id_dofs[type][id];
939  dofs.SetSize(vec.size());
940  dofs.Assign(vec.data());
941  ndofs = this->ndofs;
942 }
943 
945  std::vector<int> &dofs)
946 {
947  // Reorder the DOFs into/from a neutral ordering, independent of local
948  // edge orientation. The processor neutral edge orientation is given by
949  // the element local vertex numbering, not the mesh vertex numbering.
950 
951  const int *ev = NCMesh::GI[(int) id.element->geom].edges[id.local];
952  int v0 = id.element->node[ev[0]]->vertex->index;
953  int v1 = id.element->node[ev[1]]->vertex->index;
954 
955  if ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1]))
956  {
957  std::vector<int> tmp(dofs);
958 
959  int nv = fec->DofForGeometry(Geometry::POINT);
960  int ne = fec->DofForGeometry(Geometry::SEGMENT);
961  MFEM_ASSERT((int) dofs.size() == 2*nv + ne, "");
962 
963  // swap the two vertex DOFs
964  for (int i = 0; i < 2; i++)
965  {
966  for (int k = 0; k < nv; k++)
967  {
968  dofs[nv*i + k] = tmp[nv*(1-i) + k];
969  }
970  }
971 
972  // reorder the edge DOFs
973  int* ind = fec->DofOrderForOrientation(Geometry::SEGMENT, 0);
974  for (int i = 0; i < ne; i++)
975  {
976  dofs[2*nv + i] = (ind[i] >= 0) ? tmp[2*nv + ind[i]]
977  /* */ : -1 - tmp[2*nv + (-1 - ind[i])];
978  }
979  }
980 }
981 
982 static void write_dofs(std::ostream &os, const std::vector<int> &dofs)
983 {
984  write<int>(os, dofs.size());
985  // TODO: we should compress the ints, mostly they are contiguous ranges
986  os.write((const char*) dofs.data(), dofs.size() * sizeof(int));
987 }
988 
989 static void read_dofs(std::istream &is, std::vector<int> &dofs)
990 {
991  dofs.resize(read<int>(is));
992  is.read((char*) dofs.data(), dofs.size() * sizeof(int));
993 }
994 
996 {
997  IdToDofs::iterator it;
998 
999  // collect vertex/edge/face IDs
1000  Array<NCMesh::MeshId> ids[3];
1001  for (int type = 0; type < 3; type++)
1002  {
1003  ids[type].Reserve(id_dofs[type].size());
1004  for (it = id_dofs[type].begin(); it != id_dofs[type].end(); ++it)
1005  {
1006  ids[type].Append(it->first);
1007  }
1008  }
1009 
1010  // encode the IDs
1011  std::ostringstream stream;
1012  pncmesh->EncodeMeshIds(stream, ids, 3);
1013 
1014  // dump the DOFs
1015  for (int type = 0; type < 3; type++)
1016  {
1017  for (it = id_dofs[type].begin(); it != id_dofs[type].end(); ++it)
1018  {
1019  if (type == 1) { ReorderEdgeDofs(it->first, it->second); }
1020  write_dofs(stream, it->second);
1021  }
1022 
1023  // no longer need the original data
1024  id_dofs[type].clear();
1025  }
1026 
1027  write<int>(stream, ndofs);
1028 
1029  stream.str().swap(data);
1030 }
1031 
1033 {
1034  std::istringstream stream(data);
1035 
1036  // decode vertex/edge/face IDs
1037  Array<NCMesh::MeshId> ids[3];
1038  pncmesh->DecodeMeshIds(stream, ids, 3, true);
1039 
1040  // load DOFs
1041  for (int type = 0; type < 3; type++)
1042  {
1043  id_dofs[type].clear();
1044  for (int i = 0; i < ids[type].Size(); i++)
1045  {
1046  const NCMesh::MeshId &id = ids[type][i];
1047  read_dofs(stream, id_dofs[type][id]);
1048  if (type == 1) { ReorderEdgeDofs(id, id_dofs[type][id]); }
1049  }
1050  }
1051 
1052  ndofs = read<int>(stream);
1053 
1054  // no longer need the raw data
1055  data.clear();
1056 }
1057 
1059 {
1060  std::ostringstream stream;
1061 
1062  // write the int set to the stream
1063  write<int>(stream, rows.size());
1064  for (std::set<int>::iterator it = rows.begin(); it != rows.end(); ++it)
1065  {
1066  write<int>(stream, *it);
1067  }
1068 
1069  rows.clear();
1070  stream.str().swap(data);
1071 }
1072 
1074 {
1075  std::istringstream stream(data);
1076 
1077  // read the int set from the stream
1078  rows.clear();
1079  int size = read<int>(stream);
1080  for (int i = 0; i < size; i++)
1081  {
1082  rows.insert(rows.end(), read<int>(stream));
1083  }
1084 
1085  data.clear();
1086 }
1087 
1088 void NeighborRowReply::AddRow(int row, const Array<int> &cols,
1089  const Vector &srow)
1090 {
1091  MFEM_ASSERT(rows.find(row) == rows.end(), "");
1092  Row& row_data = rows[row];
1093  row_data.cols.assign(cols.GetData(), cols.GetData() + cols.Size());
1094  row_data.srow = srow;
1095 }
1096 
1097 void NeighborRowReply::GetRow(int row, Array<int> &cols, Vector &srow)
1098 {
1099  MFEM_ASSERT(rows.find(row) != rows.end(),
1100  "row " << row << " not found in neighbor message.");
1101  Row& row_data = rows[row];
1102  cols.SetSize(row_data.cols.size());
1103  cols.Assign(row_data.cols.data());
1104  srow = row_data.srow;
1105 }
1106 
1108 {
1109  std::ostringstream stream;
1110 
1111  // dump the rows to the stream
1112  write<int>(stream, rows.size());
1113  for (std::map<int, Row>::iterator it = rows.begin(); it != rows.end(); ++it)
1114  {
1115  write<int>(stream, it->first); // row number
1116  Row& row_data = it->second;
1117  MFEM_ASSERT((int) row_data.cols.size() == row_data.srow.Size(), "");
1118  write_dofs(stream, row_data.cols);
1119  stream.write((const char*) row_data.srow.GetData(),
1120  sizeof(double) * row_data.srow.Size());
1121  }
1122 
1123  rows.clear();
1124  stream.str().swap(data);
1125 }
1126 
1128 {
1129  std::istringstream stream(data); // stream makes a copy of data
1130 
1131  // NOTE: there is no rows.clear() since a row reply can be received
1132  // repeatedly and the received rows accumulate.
1133 
1134  // read the rows
1135  int size = read<int>(stream);
1136  for (int i = 0; i < size; i++)
1137  {
1138  Row& row_data = rows[read<int>(stream)];
1139  read_dofs(stream, row_data.cols);
1140  row_data.srow.SetSize(row_data.cols.size());
1141  stream.read((char*) row_data.srow.GetData(),
1142  sizeof(double) * row_data.srow.Size());
1143  }
1144 
1145  data.clear();
1146 }
1147 
1149 {
1150  Array<MeshId> ids;
1151 
1152  // abuse EncodeMeshIds() to encode the list of refinements
1153  ids.Reserve(refinements.size());
1154  for (unsigned i = 0; i < refinements.size(); i++)
1155  {
1156  const ElemRefType &ref = refinements[i];
1157  ids.Append(MeshId(-1, ref.elem, ref.ref_type));
1158  }
1159 
1160  std::ostringstream stream;
1161  pncmesh->EncodeMeshIds(stream, &ids, 1);
1162 
1163  stream.str().swap(data);
1164 }
1165 
1167 {
1169 
1170  // inverse abuse to Encode()
1171  std::istringstream stream(data);
1172  pncmesh->DecodeMeshIds(stream, &ids, 1, false);
1173 
1174  refinements.clear();
1175  refinements.reserve(ids.Size());
1176  for (int i = 0; i < ids.Size(); i++)
1177  {
1178  AddRefinement(ids[i].element, ids[i].local);
1179  }
1180 
1181  data.clear();
1182 }
1183 
1184 
1186 
1187 void ParNCMesh::GetDebugMesh(Mesh &debug_mesh) const
1188 {
1189  // create a serial NCMesh containing all our elements (ghosts and all)
1190  NCMesh* copy = new NCMesh(*this);
1191 
1192  Array<Element*> &cle = copy->leaf_elements;
1193  for (int i = 0; i < cle.Size(); i++)
1194  {
1195  cle[i]->attribute = (cle[i]->rank == MyRank) ? 1 : 2;
1196  }
1197 
1198  debug_mesh.InitFromNCMesh(*copy);
1199  debug_mesh.SetAttributes();
1200  debug_mesh.ncmesh = copy;
1201 }
1202 
1203 
1204 } // namespace mfem
1205 
1206 #endif // MFEM_USE_MPI
NCList face_list
lazy-initialized list of faces, see GetFaceList
Definition: ncmesh.hpp:383
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
Definition: ncmesh.hpp:384
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
Definition: ncmesh.cpp:55
NCList shared_edges
Definition: pncmesh.hpp:194
virtual void LimitNCLevel(int max_level)
Definition: ncmesh.cpp:2530
Array< Element * > ghost_layer
list of elements whose &#39;element_type&#39; == 2.
Definition: pncmesh.hpp:217
Array< char > element_type
Definition: pncmesh.hpp:215
void EncodeMeshIds(std::ostream &os, Array< MeshId > ids[], int dim) const
Write to &#39;os&#39; a processor-independent encoding of vertex/edge/face IDs.
Definition: pncmesh.cpp:819
int Size() const
Logical size of the array.
Definition: array.hpp:109
void FindNeighbors(const Element *elem, Array< Element * > &neighbors, const Array< Element * > *search_set=NULL)
Definition: ncmesh.cpp:1965
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:35
void CalcFaceOrientations()
Definition: pncmesh.cpp:387
void Unique()
Definition: array.hpp:193
bool EncodeTree(Element *elem, const std::set< Element * > &elements)
Definition: pncmesh.cpp:692
int index
edge number in the Mesh
Definition: ncmesh.hpp:278
void AddDofs(int type, const NCMesh::MeshId &id, const Array< int > &dofs)
Add vertex/edge/face DOFs to an outgoing message.
Definition: pncmesh.cpp:920
Array< char > face_orient
Definition: pncmesh.hpp:207
Iterator over items contained in the HashTable.
Definition: hash.hpp:160
MPI_Comm MyComm
Definition: pncmesh.hpp:184
virtual void BuildFaceList()
Definition: pncmesh.cpp:209
void Dump(std::ostream &os) const
Definition: pncmesh.cpp:804
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
Array< int > face_owner
Definition: pncmesh.hpp:200
void BuildSharedVertices()
Definition: pncmesh.cpp:337
T * GetData()
Returns the data.
Definition: array.hpp:91
void Load(std::istream &is)
Definition: pncmesh.cpp:810
int Size() const
Returns the size of the vector.
Definition: vector.hpp:84
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
static int find_node(Element *elem, Node *node)
Definition: ncmesh.cpp:1479
void DecodeTree(Element *elem, int &pos, Array< Element * > &elements) const
Definition: pncmesh.cpp:759
int index
vertex number in the Mesh
Definition: ncmesh.hpp:265
int GetInt(int pos) const
Definition: pncmesh.cpp:682
virtual void OnMeshUpdated(Mesh *mesh)
Definition: pncmesh.cpp:127
std::vector< MeshId > conforming
Definition: ncmesh.hpp:126
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:3924
double * GetData() const
Definition: vector.hpp:88
virtual void AssignLeafIndices()
Definition: pncmesh.cpp:55
void MakeShared(const Table &groups, const NCList &list, NCList &shared)
Definition: pncmesh.cpp:309
Array< int > vertex_nodeId
Definition: ncmesh.hpp:381
void DeleteAll()
Delete whole array.
Definition: array.hpp:455
int index
Mesh element number.
Definition: ncmesh.hpp:34
int master
master number (in Mesh numbering)
Definition: ncmesh.hpp:117
void Recv(int rank, int size, MPI_Comm comm)
Post-probe receive from processor &#39;rank&#39; of message size &#39;size&#39;.
virtual void OnMeshUpdated(Mesh *mesh)
Definition: ncmesh.cpp:1426
virtual void Decode()
Definition: pncmesh.cpp:1073
std::vector< int > cols
Definition: pncmesh.hpp:412
ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh)
Definition: pncmesh.cpp:26
Array< int > vertex_owner
Definition: pncmesh.hpp:198
virtual void Update()
Definition: pncmesh.cpp:43
virtual void Decode()
Definition: pncmesh.cpp:1032
Table face_group
Definition: pncmesh.hpp:205
void DecodeMeshIds(std::istream &is, Array< MeshId > ids[], int dim, bool decode_indices) const
Read from &#39;is&#39; a processor-independent encoding of vertex/edge/face IDs.
Definition: pncmesh.cpp:861
int dim
Definition: ex3.cpp:48
bool PruneTree(Element *elem)
Internal. Recursive part of Prune().
Definition: pncmesh.cpp:530
int index
face number in the Mesh
Definition: ncmesh.hpp:323
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:260
void SetAttributes()
Definition: mesh.cpp:796
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:368
int InitialPartition(int index) const
Assigns elements to processors at the initial stage (ParMesh creation).
Definition: pncmesh.hpp:222
void RefineElement(Element *elem, char ref_type)
Definition: ncmesh.cpp:725
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:122
NCList shared_vertices
Definition: pncmesh.hpp:193
virtual void BuildEdgeList()
Definition: ncmesh.cpp:1693
Identifies a vertex/edge/face in both Mesh and NCMesh.
Definition: ncmesh.hpp:93
void Assign(const T *)
Copy data from a pointer. Size() elements are copied.
Definition: array.hpp:510
std::map< int, NeighborRefinementMessage > Map
Definition: pncmesh.hpp:322
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
void SetInt(int pos, int value)
Definition: pncmesh.cpp:673
Table edge_group
Definition: pncmesh.hpp:204
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
Definition: ncmesh.hpp:63
NCList shared_faces
Definition: pncmesh.hpp:195
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:2387
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: pncmesh.cpp:424
static GeomInfo GI[Geometry::NumGeom]
Definition: ncmesh.hpp:632
int slaves_end
slave faces
Definition: ncmesh.hpp:107
Vertex * vertex
Definition: ncmesh.hpp:297
bool Iso
true if the mesh only contains isotropic refinements
Definition: ncmesh.hpp:233
Table vertex_group
Definition: pncmesh.hpp:203
void DerefineElement(Element *elem)
Definition: ncmesh.cpp:1179
virtual void BuildEdgeList()
Definition: pncmesh.cpp:185
void FindSetNeighbors(const Array< char > &elem_set, Array< Element * > *neighbors, Array< char > *neighbor_set=NULL)
Definition: ncmesh.cpp:1881
virtual void BuildFaceList()
Definition: ncmesh.cpp:1591
void GetDebugMesh(Mesh &debug_mesh) const
Definition: pncmesh.cpp:1187
void GetDofs(int type, const NCMesh::MeshId &id, Array< int > &dofs, int &ndofs)
Definition: pncmesh.cpp:927
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:323
Element * child[8]
2-8 children (if ref_type != 0)
Definition: ncmesh.hpp:355
Array< Element * > leaf_elements
Definition: ncmesh.hpp:379
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:24
void Decode(Array< Element * > &elements, const Array< Element * > &ncmesh_roots) const
Definition: pncmesh.cpp:779
virtual void Encode()
Definition: pncmesh.cpp:1107
void UpdateLayers()
Definition: pncmesh.cpp:448
void ElementNeighborProcessors(Element *elem, Array< int > &ranks)
Definition: pncmesh.cpp:483
virtual void ElementSharesFace(Element *elem, Face *face)
Definition: pncmesh.cpp:175
virtual void ElementSharesEdge(Element *elem, Edge *edge)
Definition: pncmesh.cpp:163
HashTable< Node > nodes
Definition: ncmesh.hpp:365
virtual void Decode()
Definition: pncmesh.cpp:1127
virtual void Refine(const Array< Refinement > &refinements)
Definition: pncmesh.cpp:588
virtual void Encode()
Definition: pncmesh.cpp:995
void GetRow(int row, Array< int > &cols, Vector &srow)
Definition: pncmesh.cpp:1097
std::vector< ElemRefType > refinements
Definition: pncmesh.hpp:314
Array< Connection > index_rank
Definition: pncmesh.hpp:247
Array< int > edge_owner
Definition: pncmesh.hpp:199
static int find_hex_face(int a, int b, int c)
Definition: ncmesh.cpp:1497
void SetNCMesh(ParNCMesh *pncmesh)
Set pointer to ParNCMesh (needed to encode the message).
Definition: pncmesh.hpp:320
NCMesh * ncmesh
Definition: mesh.hpp:143
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:460
Array< Element * > root_elements
Definition: ncmesh.hpp:363
virtual void Update()
Definition: ncmesh.cpp:174
static void WaitAllSent(MapT &rank_msg)
Helper to wait for all messages in a map container to be sent.
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:463
Vector data type.
Definition: vector.hpp:33
void AddRow(int row, const Array< int > &cols, const Vector &srow)
Definition: pncmesh.cpp:1088
void AddMasterSlaveRanks(int nitems, const NCList &list)
Definition: pncmesh.cpp:242
void ReorderEdgeDofs(const NCMesh::MeshId &id, std::vector< int > &dofs)
Definition: pncmesh.cpp:944
Array< Element * > tmp_neighbors
Definition: pncmesh.hpp:287
void NeighborProcessors(Array< int > &neighbors)
Definition: pncmesh.cpp:507
int RowSize(int i) const
Definition: table.hpp:102
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
Definition: ncmesh.hpp:348
virtual void Encode()
Definition: pncmesh.cpp:1058
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
Definition: mesh.cpp:6741
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
Definition: pncmesh.cpp:81
int index
Mesh number.
Definition: ncmesh.hpp:95
Abstract data type element.
Definition: element.hpp:27
virtual void LimitNCLevel(int max_level)
To be implemented.
Definition: pncmesh.cpp:661
Array< unsigned char > data
encoded refinement (sub-)trees
Definition: pncmesh.hpp:271
static void Probe(int &rank, int &size, MPI_Comm comm)
HashTable< Face > faces
Definition: ncmesh.hpp:366