MFEM  v4.6.0
Finite element discretization library
pmesh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "mesh_headers.hpp"
17 #include "../fem/fem.hpp"
18 #include "../general/sets.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../general/text.hpp"
21 #include "../general/globals.hpp"
22 
23 #include <iostream>
24 #include <fstream>
25 
26 using namespace std;
27 
28 namespace mfem
29 {
30 
31 ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
32  : Mesh(pmesh, false),
33  group_svert(pmesh.group_svert),
34  group_sedge(pmesh.group_sedge),
35  group_stria(pmesh.group_stria),
36  group_squad(pmesh.group_squad),
37  face_nbr_el_to_face(NULL),
38  glob_elem_offset(-1),
39  glob_offset_sequence(-1),
40  gtopo(pmesh.gtopo)
41 {
42  MyComm = pmesh.MyComm;
43  NRanks = pmesh.NRanks;
44  MyRank = pmesh.MyRank;
45 
46  // Duplicate the shared_edges
47  shared_edges.SetSize(pmesh.shared_edges.Size());
48  for (int i = 0; i < shared_edges.Size(); i++)
49  {
50  shared_edges[i] = pmesh.shared_edges[i]->Duplicate(this);
51  }
52 
53  shared_trias = pmesh.shared_trias;
54  shared_quads = pmesh.shared_quads;
55 
56  // Copy the shared-to-local index Arrays
59  sface_lface = pmesh.sface_lface;
60 
61  // Do not copy face-neighbor data (can be generated if needed)
62  have_face_nbr_data = false;
63 
64  // If pmesh has a ParNURBSExtension, it was copied by the Mesh copy ctor, so
65  // there is no need to do anything here.
66 
67  // Copy ParNCMesh, if present
68  if (pmesh.pncmesh)
69  {
70  pncmesh = new ParNCMesh(*pmesh.pncmesh);
71  pncmesh->OnMeshUpdated(this);
72  }
73  else
74  {
75  pncmesh = NULL;
76  }
77  ncmesh = pncmesh;
78 
79  // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
80  // and the FiniteElementSpace (as a ParFiniteElementSpace)
81  if (pmesh.Nodes && copy_nodes)
82  {
83  FiniteElementSpace *fes = pmesh.Nodes->FESpace();
84  const FiniteElementCollection *fec = fes->FEColl();
85  FiniteElementCollection *fec_copy =
87  ParFiniteElementSpace *pfes_copy =
88  new ParFiniteElementSpace(*fes, *this, fec_copy);
89  Nodes = new ParGridFunction(pfes_copy);
90  Nodes->MakeOwner(fec_copy);
91  *Nodes = *pmesh.Nodes;
92  own_nodes = 1;
93  }
94 }
95 
97 {
98  Swap(mesh);
99 }
100 
102 {
103  Swap(mesh);
104  return *this;
105 }
106 
107 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
108  int part_method)
109  : face_nbr_el_to_face(NULL)
110  , glob_elem_offset(-1)
111  , glob_offset_sequence(-1)
112  , gtopo(comm)
113 {
114  int *partitioning = NULL;
115  Array<bool> activeBdrElem;
116 
117  MyComm = comm;
118  MPI_Comm_size(MyComm, &NRanks);
119  MPI_Comm_rank(MyComm, &MyRank);
120 
121  if (mesh.Nonconforming())
122  {
123  if (partitioning_)
124  {
125  partitioning = partitioning_;
126  }
127  ncmesh = pncmesh = new ParNCMesh(comm, *mesh.ncmesh, partitioning);
128  if (!partitioning)
129  {
130  partitioning = new int[mesh.GetNE()];
131  for (int i = 0; i < mesh.GetNE(); i++)
132  {
133  partitioning[i] = pncmesh->InitialPartition(i);
134  }
135  }
136 
137  pncmesh->Prune();
138 
140  pncmesh->OnMeshUpdated(this);
141 
143 
144  // SetMeshGen(); // called by Mesh::InitFromNCMesh(...) above
145  meshgen = mesh.meshgen; // copy the global 'meshgen'
146 
147  mesh.attributes.Copy(attributes);
149 
151  }
152  else // mesh.Conforming()
153  {
154  Dim = mesh.Dim;
155  spaceDim = mesh.spaceDim;
156 
157  ncmesh = pncmesh = NULL;
158 
159  if (partitioning_)
160  {
161  partitioning = partitioning_;
162  }
163  else
164  {
165  partitioning = mesh.GeneratePartitioning(NRanks, part_method);
166  }
167 
168  // re-enumerate the partitions to better map to actual processor
169  // interconnect topology !?
170 
171  Array<int> vert_global_local;
172  NumOfVertices = BuildLocalVertices(mesh, partitioning, vert_global_local);
173  NumOfElements = BuildLocalElements(mesh, partitioning, vert_global_local);
174 
175  Table *edge_element = NULL;
176  NumOfBdrElements = BuildLocalBoundary(mesh, partitioning,
177  vert_global_local,
178  activeBdrElem, edge_element);
179 
180  SetMeshGen();
181  meshgen = mesh.meshgen; // copy the global 'meshgen'
182 
183  mesh.attributes.Copy(attributes);
185 
186  NumOfEdges = NumOfFaces = 0;
187 
188  if (Dim > 1)
189  {
190  el_to_edge = new Table;
192  }
193 
194  STable3D *faces_tbl = NULL;
195  if (Dim == 3)
196  {
197  faces_tbl = GetElementToFaceTable(1);
198  }
199 
200  GenerateFaces();
201 
202  ListOfIntegerSets groups;
203  {
204  // the first group is the local one
205  IntegerSet group;
206  group.Recreate(1, &MyRank);
207  groups.Insert(group);
208  }
209 
210  MFEM_ASSERT(mesh.GetNFaces() == 0 || Dim >= 3, "");
211 
212  Array<int> face_group(mesh.GetNFaces());
213  Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
214 
215  FindSharedFaces(mesh, partitioning, face_group, groups);
216  int nsedges = FindSharedEdges(mesh, partitioning, edge_element, groups);
217  int nsvert = FindSharedVertices(partitioning, vert_element, groups);
218 
219  // build the group communication topology
220  gtopo.Create(groups, 822);
221 
222  // fill out group_sface, group_sedge, group_svert
223  int ngroups = groups.Size()-1, nstris, nsquads;
224  BuildFaceGroup(ngroups, mesh, face_group, nstris, nsquads);
225  BuildEdgeGroup(ngroups, *edge_element);
226  BuildVertexGroup(ngroups, *vert_element);
227 
228  // build shared_faces and sface_lface mapping
229  BuildSharedFaceElems(nstris, nsquads, mesh, partitioning, faces_tbl,
230  face_group, vert_global_local);
231  delete faces_tbl;
232 
233  // build shared_edges and sedge_ledge mapping
234  BuildSharedEdgeElems(nsedges, mesh, vert_global_local, edge_element);
235  delete edge_element;
236 
237  // build svert_lvert mapping
238  BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
239  delete vert_element;
240 
241  SetMeshGen();
242  meshgen = mesh.meshgen; // copy the global 'meshgen'
243  }
244 
245  if (mesh.NURBSext)
246  {
247  MFEM_ASSERT(mesh.GetNodes() &&
248  mesh.GetNodes()->FESpace()->GetNURBSext() == mesh.NURBSext,
249  "invalid NURBS mesh");
250  NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
251  activeBdrElem);
252  }
253 
254  if (mesh.GetNodes()) // curved mesh
255  {
256  if (!NURBSext)
257  {
258  Nodes = new ParGridFunction(this, mesh.GetNodes());
259  }
260  else
261  {
262  const FiniteElementSpace *glob_fes = mesh.GetNodes()->FESpace();
264  FiniteElementCollection::New(glob_fes->FEColl()->Name());
265  ParFiniteElementSpace *pfes =
266  new ParFiniteElementSpace(this, nfec, glob_fes->GetVDim(),
267  glob_fes->GetOrdering());
268  Nodes = new ParGridFunction(pfes);
269  Nodes->MakeOwner(nfec); // Nodes will own nfec and pfes
270  }
271  own_nodes = 1;
272 
273  Array<int> gvdofs, lvdofs;
274  Vector lnodes;
275  int element_counter = 0;
276  for (int i = 0; i < mesh.GetNE(); i++)
277  {
278  if (partitioning[i] == MyRank)
279  {
280  Nodes->FESpace()->GetElementVDofs(element_counter, lvdofs);
281  mesh.GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
282  mesh.GetNodes()->GetSubVector(gvdofs, lnodes);
283  Nodes->SetSubVector(lvdofs, lnodes);
284  element_counter++;
285  }
286  }
287 
288  // set meaningful values to 'vertices' even though we have Nodes,
289  // for compatibility (e.g., Mesh::GetVertex())
291  }
292 
293  if (partitioning != partitioning_)
294  {
295  delete [] partitioning;
296  }
297 
298  have_face_nbr_data = false;
299 }
300 
301 
303  const int* partitioning,
304  Array<int> &vert_global_local)
305 {
306  vert_global_local.SetSize(mesh.GetNV());
307  vert_global_local = -1;
308 
309  int vert_counter = 0;
310  for (int i = 0; i < mesh.GetNE(); i++)
311  {
312  if (partitioning[i] == MyRank)
313  {
314  Array<int> vert;
315  mesh.GetElementVertices(i, vert);
316  for (int j = 0; j < vert.Size(); j++)
317  {
318  if (vert_global_local[vert[j]] < 0)
319  {
320  vert_global_local[vert[j]] = vert_counter++;
321  }
322  }
323  }
324  }
325 
326  // re-enumerate the local vertices to preserve the global ordering
327  vert_counter = 0;
328  for (int i = 0; i < vert_global_local.Size(); i++)
329  {
330  if (vert_global_local[i] >= 0)
331  {
332  vert_global_local[i] = vert_counter++;
333  }
334  }
335 
336  vertices.SetSize(vert_counter);
337 
338  for (int i = 0; i < vert_global_local.Size(); i++)
339  {
340  if (vert_global_local[i] >= 0)
341  {
342  vertices[vert_global_local[i]].SetCoords(mesh.SpaceDimension(),
343  mesh.GetVertex(i));
344  }
345  }
346 
347  return vert_counter;
348 }
349 
350 int ParMesh::BuildLocalElements(const Mesh& mesh, const int* partitioning,
351  const Array<int>& vert_global_local)
352 {
353  int nelems = 0;
354  for (int i = 0; i < mesh.GetNE(); i++)
355  {
356  if (partitioning[i] == MyRank) { nelems++; }
357  }
358 
359  elements.SetSize(nelems);
360 
361  // Determine elements, enumerating the local elements to preserve the global
362  // order. This is used, e.g. by the ParGridFunction ctor that takes a global
363  // GridFunction as input parameter.
364  int element_counter = 0;
365  for (int i = 0; i < mesh.GetNE(); i++)
366  {
367  if (partitioning[i] == MyRank)
368  {
369  elements[element_counter] = mesh.GetElement(i)->Duplicate(this);
370  int *v = elements[element_counter]->GetVertices();
371  int nv = elements[element_counter]->GetNVertices();
372  for (int j = 0; j < nv; j++)
373  {
374  v[j] = vert_global_local[v[j]];
375  }
376  element_counter++;
377  }
378  }
379 
380  return element_counter;
381 }
382 
383 int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
384  const Array<int>& vert_global_local,
385  Array<bool>& activeBdrElem,
386  Table*& edge_element)
387 {
388  int nbdry = 0;
389 
390  if (mesh.NURBSext)
391  {
392  activeBdrElem.SetSize(mesh.GetNBE());
393  activeBdrElem = false;
394  }
395  // build boundary elements
396  if (Dim == 3)
397  {
398  for (int i = 0; i < mesh.GetNBE(); i++)
399  {
400  int face, o, el1, el2;
401  mesh.GetBdrElementFace(i, &face, &o);
402  mesh.GetFaceElements(face, &el1, &el2);
403  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
404  {
405  nbdry++;
406  if (mesh.NURBSext)
407  {
408  activeBdrElem[i] = true;
409  }
410  }
411  }
412 
413  int bdrelem_counter = 0;
414  boundary.SetSize(nbdry);
415  for (int i = 0; i < mesh.GetNBE(); i++)
416  {
417  int face, o, el1, el2;
418  mesh.GetBdrElementFace(i, &face, &o);
419  mesh.GetFaceElements(face, &el1, &el2);
420  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
421  {
422  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
423  int *v = boundary[bdrelem_counter]->GetVertices();
424  int nv = boundary[bdrelem_counter]->GetNVertices();
425  for (int j = 0; j < nv; j++)
426  {
427  v[j] = vert_global_local[v[j]];
428  }
429  bdrelem_counter++;
430  }
431  }
432  }
433  else if (Dim == 2)
434  {
435  edge_element = new Table;
436  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
437 
438  for (int i = 0; i < mesh.GetNBE(); i++)
439  {
440  int edge = mesh.GetBdrElementEdgeIndex(i);
441  int el1 = edge_element->GetRow(edge)[0];
442  if (partitioning[el1] == MyRank)
443  {
444  nbdry++;
445  if (mesh.NURBSext)
446  {
447  activeBdrElem[i] = true;
448  }
449  }
450  }
451 
452  int bdrelem_counter = 0;
453  boundary.SetSize(nbdry);
454  for (int i = 0; i < mesh.GetNBE(); i++)
455  {
456  int edge = mesh.GetBdrElementEdgeIndex(i);
457  int el1 = edge_element->GetRow(edge)[0];
458  if (partitioning[el1] == MyRank)
459  {
460  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
461  int *v = boundary[bdrelem_counter]->GetVertices();
462  int nv = boundary[bdrelem_counter]->GetNVertices();
463  for (int j = 0; j < nv; j++)
464  {
465  v[j] = vert_global_local[v[j]];
466  }
467  bdrelem_counter++;
468  }
469  }
470  }
471  else if (Dim == 1)
472  {
473  for (int i = 0; i < mesh.GetNBE(); i++)
474  {
475  int vert = mesh.boundary[i]->GetVertices()[0];
476  int el1, el2;
477  mesh.GetFaceElements(vert, &el1, &el2);
478  if (partitioning[el1] == MyRank)
479  {
480  nbdry++;
481  }
482  }
483 
484  int bdrelem_counter = 0;
485  boundary.SetSize(nbdry);
486  for (int i = 0; i < mesh.GetNBE(); i++)
487  {
488  int vert = mesh.boundary[i]->GetVertices()[0];
489  int el1, el2;
490  mesh.GetFaceElements(vert, &el1, &el2);
491  if (partitioning[el1] == MyRank)
492  {
493  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
494  int *v = boundary[bdrelem_counter]->GetVertices();
495  v[0] = vert_global_local[v[0]];
496  bdrelem_counter++;
497  }
498  }
499  }
500 
501  return nbdry;
502 }
503 
504 void ParMesh::FindSharedFaces(const Mesh &mesh, const int *partitioning,
505  Array<int> &face_group,
506  ListOfIntegerSets &groups)
507 {
508  IntegerSet group;
509 
510  // determine shared faces
511  face_group.SetSize(mesh.GetNFaces());
512  for (int i = 0; i < face_group.Size(); i++)
513  {
514  int el[2];
515  face_group[i] = -1;
516  mesh.GetFaceElements(i, &el[0], &el[1]);
517  if (el[1] >= 0)
518  {
519  el[0] = partitioning[el[0]];
520  el[1] = partitioning[el[1]];
521  if ((el[0] == MyRank && el[1] != MyRank) ||
522  (el[0] != MyRank && el[1] == MyRank))
523  {
524  group.Recreate(2, el);
525  face_group[i] = groups.Insert(group) - 1;
526  }
527  }
528  }
529 }
530 
531 int ParMesh::FindSharedEdges(const Mesh &mesh, const int *partitioning,
532  Table*& edge_element,
533  ListOfIntegerSets& groups)
534 {
535  IntegerSet group;
536 
537  // determine shared edges
538  int sedge_counter = 0;
539  if (!edge_element)
540  {
541  edge_element = new Table;
542  if (Dim == 1)
543  {
544  edge_element->SetDims(0,0);
545  }
546  else
547  {
548  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
549  }
550  }
551 
552  for (int i = 0; i < edge_element->Size(); i++)
553  {
554  int me = 0, others = 0;
555  for (int j = edge_element->GetI()[i]; j < edge_element->GetI()[i+1]; j++)
556  {
557  int k = edge_element->GetJ()[j];
558  int rank = partitioning[k];
559  edge_element->GetJ()[j] = rank;
560  if (rank == MyRank)
561  {
562  me = 1;
563  }
564  else
565  {
566  others = 1;
567  }
568  }
569 
570  if (me && others)
571  {
572  sedge_counter++;
573  group.Recreate(edge_element->RowSize(i), edge_element->GetRow(i));
574  edge_element->GetRow(i)[0] = groups.Insert(group) - 1;
575  }
576  else
577  {
578  edge_element->GetRow(i)[0] = -1;
579  }
580  }
581 
582  return sedge_counter;
583 }
584 
585 int ParMesh::FindSharedVertices(const int *partitioning, Table *vert_element,
586  ListOfIntegerSets &groups)
587 {
588  IntegerSet group;
589 
590  // determine shared vertices
591  int svert_counter = 0;
592  for (int i = 0; i < vert_element->Size(); i++)
593  {
594  int me = 0, others = 0;
595  for (int j = vert_element->GetI()[i]; j < vert_element->GetI()[i+1]; j++)
596  {
597  vert_element->GetJ()[j] = partitioning[vert_element->GetJ()[j]];
598  if (vert_element->GetJ()[j] == MyRank)
599  {
600  me = 1;
601  }
602  else
603  {
604  others = 1;
605  }
606  }
607 
608  if (me && others)
609  {
610  svert_counter++;
611  group.Recreate(vert_element->RowSize(i), vert_element->GetRow(i));
612  vert_element->GetI()[i] = groups.Insert(group) - 1;
613  }
614  else
615  {
616  vert_element->GetI()[i] = -1;
617  }
618  }
619  return svert_counter;
620 }
621 
622 void ParMesh::BuildFaceGroup(int ngroups, const Mesh &mesh,
623  const Array<int> &face_group,
624  int &nstria, int &nsquad)
625 {
626  // build group_stria and group_squad
627  group_stria.MakeI(ngroups);
628  group_squad.MakeI(ngroups);
629 
630  for (int i = 0; i < face_group.Size(); i++)
631  {
632  if (face_group[i] >= 0)
633  {
634  if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
635  {
636  group_stria.AddAColumnInRow(face_group[i]);
637  }
638  else
639  {
640  group_squad.AddAColumnInRow(face_group[i]);
641  }
642  }
643  }
644 
645  group_stria.MakeJ();
646  group_squad.MakeJ();
647 
648  nstria = nsquad = 0;
649  for (int i = 0; i < face_group.Size(); i++)
650  {
651  if (face_group[i] >= 0)
652  {
653  if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
654  {
655  group_stria.AddConnection(face_group[i], nstria++);
656  }
657  else
658  {
659  group_squad.AddConnection(face_group[i], nsquad++);
660  }
661  }
662  }
663 
666 }
667 
668 void ParMesh::BuildEdgeGroup(int ngroups, const Table &edge_element)
669 {
670  group_sedge.MakeI(ngroups);
671 
672  for (int i = 0; i < edge_element.Size(); i++)
673  {
674  if (edge_element.GetRow(i)[0] >= 0)
675  {
676  group_sedge.AddAColumnInRow(edge_element.GetRow(i)[0]);
677  }
678  }
679 
680  group_sedge.MakeJ();
681 
682  int sedge_counter = 0;
683  for (int i = 0; i < edge_element.Size(); i++)
684  {
685  if (edge_element.GetRow(i)[0] >= 0)
686  {
687  group_sedge.AddConnection(edge_element.GetRow(i)[0], sedge_counter++);
688  }
689  }
690 
692 }
693 
694 void ParMesh::BuildVertexGroup(int ngroups, const Table &vert_element)
695 {
696  group_svert.MakeI(ngroups);
697 
698  for (int i = 0; i < vert_element.Size(); i++)
699  {
700  if (vert_element.GetI()[i] >= 0)
701  {
702  group_svert.AddAColumnInRow(vert_element.GetI()[i]);
703  }
704  }
705 
706  group_svert.MakeJ();
707 
708  int svert_counter = 0;
709  for (int i = 0; i < vert_element.Size(); i++)
710  {
711  if (vert_element.GetI()[i] >= 0)
712  {
713  group_svert.AddConnection(vert_element.GetI()[i], svert_counter++);
714  }
715  }
716 
718 }
719 
720 void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
721  const Mesh& mesh, int *partitioning,
722  const STable3D *faces_tbl,
723  const Array<int> &face_group,
724  const Array<int> &vert_global_local)
725 {
726  shared_trias.SetSize(ntri_faces);
727  shared_quads.SetSize(nquad_faces);
728  sface_lface. SetSize(ntri_faces + nquad_faces);
729 
730  if (Dim < 3) { return; }
731 
732  int stria_counter = 0;
733  int squad_counter = 0;
734  for (int i = 0; i < face_group.Size(); i++)
735  {
736  if (face_group[i] < 0) { continue; }
737 
738  const Element *face = mesh.GetFace(i);
739  const int *fv = face->GetVertices();
740  switch (face->GetType())
741  {
742  case Element::TRIANGLE:
743  {
744  shared_trias[stria_counter].Set(fv);
745  int *v = shared_trias[stria_counter].v;
746  for (int j = 0; j < 3; j++)
747  {
748  v[j] = vert_global_local[v[j]];
749  }
750  const int lface = (*faces_tbl)(v[0], v[1], v[2]);
751  sface_lface[stria_counter] = lface;
752  if (meshgen == 1) // Tet-only mesh
753  {
754  Tetrahedron *tet = dynamic_cast<Tetrahedron *>
755  (elements[faces_info[lface].Elem1No]);
756  // mark the shared face for refinement by reorienting
757  // it according to the refinement flag in the tetrahedron
758  // to which this shared face belongs to.
759  if (tet->GetRefinementFlag())
760  {
761  tet->GetMarkedFace(faces_info[lface].Elem1Inf/64, v);
762  // flip the shared face in the processor that owns the
763  // second element (in 'mesh')
764  int gl_el1, gl_el2;
765  mesh.GetFaceElements(i, &gl_el1, &gl_el2);
766  if (MyRank == partitioning[gl_el2])
767  {
768  std::swap(v[0], v[1]);
769  }
770  }
771  }
772  stria_counter++;
773  break;
774  }
775 
777  {
778  shared_quads[squad_counter].Set(fv);
779  int *v = shared_quads[squad_counter].v;
780  for (int j = 0; j < 4; j++)
781  {
782  v[j] = vert_global_local[v[j]];
783  }
784  sface_lface[shared_trias.Size() + squad_counter] =
785  (*faces_tbl)(v[0], v[1], v[2], v[3]);
786  squad_counter++;
787  break;
788  }
789 
790  default:
791  MFEM_ABORT("unknown face type: " << face->GetType());
792  break;
793  }
794  }
795 }
796 
797 void ParMesh::BuildSharedEdgeElems(int nedges, Mesh& mesh,
798  const Array<int>& vert_global_local,
799  const Table* edge_element)
800 {
801  // The passed in mesh is still the global mesh. "this" mesh is the
802  // local partitioned mesh.
803 
804  shared_edges.SetSize(nedges);
805  sedge_ledge. SetSize(nedges);
806 
807  {
808  DSTable v_to_v(NumOfVertices);
809  GetVertexToVertexTable(v_to_v);
810 
811  int sedge_counter = 0;
812  for (int i = 0; i < edge_element->Size(); i++)
813  {
814  if (edge_element->GetRow(i)[0] >= 0)
815  {
816  Array<int> vert;
817  mesh.GetEdgeVertices(i, vert);
818 
819  shared_edges[sedge_counter] =
820  new Segment(vert_global_local[vert[0]],
821  vert_global_local[vert[1]], 1);
822 
823  sedge_ledge[sedge_counter] = v_to_v(vert_global_local[vert[0]],
824  vert_global_local[vert[1]]);
825 
826  MFEM_VERIFY(sedge_ledge[sedge_counter] >= 0, "Error in v_to_v.");
827 
828  sedge_counter++;
829  }
830  }
831  }
832 }
833 
835  const mfem::Table *vert_element,
836  const Array<int> &vert_global_local)
837 {
838  // build svert_lvert
839  svert_lvert.SetSize(nvert);
840 
841  int svert_counter = 0;
842  for (int i = 0; i < vert_element->Size(); i++)
843  {
844  if (vert_element->GetI()[i] >= 0)
845  {
846  svert_lvert[svert_counter++] = vert_global_local[i];
847  }
848  }
849 }
850 
851 
852 // protected method, used by Nonconforming(De)Refinement and Rebalance
854  : MyComm(pncmesh.MyComm)
855  , NRanks(pncmesh.NRanks)
856  , MyRank(pncmesh.MyRank)
857  , face_nbr_el_to_face(NULL)
858  , glob_elem_offset(-1)
859  , glob_offset_sequence(-1)
860  , gtopo(MyComm)
861  , pncmesh(NULL)
862 {
864  ReduceMeshGen();
865  have_face_nbr_data = false;
866 }
867 
869 {
870  if (glob_offset_sequence != sequence) // mesh has changed
871  {
872  long long local_elems = NumOfElements;
873  MPI_Scan(&local_elems, &glob_elem_offset, 1, MPI_LONG_LONG, MPI_SUM,
874  MyComm);
875  glob_elem_offset -= local_elems;
876 
877  glob_offset_sequence = sequence; // don't recalculate until refinement etc.
878  }
879 }
880 
882 {
883  int loc_meshgen = meshgen;
884  MPI_Allreduce(&loc_meshgen, &meshgen, 1, MPI_INT, MPI_BOR, MyComm);
885 }
886 
888 {
889  // Determine sedge_ledge
891  if (shared_edges.Size())
892  {
893  DSTable v_to_v(NumOfVertices);
894  GetVertexToVertexTable(v_to_v);
895  for (int se = 0; se < shared_edges.Size(); se++)
896  {
897  const int *v = shared_edges[se]->GetVertices();
898  const int l_edge = v_to_v(v[0], v[1]);
899  MFEM_ASSERT(l_edge >= 0, "invalid shared edge");
900  sedge_ledge[se] = l_edge;
901  }
902  }
903 
904  // Determine sface_lface
905  const int nst = shared_trias.Size();
906  sface_lface.SetSize(nst + shared_quads.Size());
907  if (sface_lface.Size())
908  {
909  STable3D *faces_tbl = GetFacesTable();
910  for (int st = 0; st < nst; st++)
911  {
912  const int *v = shared_trias[st].v;
913  sface_lface[st] = (*faces_tbl)(v[0], v[1], v[2]);
914  }
915  for (int sq = 0; sq < shared_quads.Size(); sq++)
916  {
917  const int *v = shared_quads[sq].v;
918  sface_lface[nst+sq] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
919  }
920  delete faces_tbl;
921  }
922 }
923 
924 ParMesh::ParMesh(MPI_Comm comm, istream &input, bool refine)
925  : face_nbr_el_to_face(NULL)
926  , glob_elem_offset(-1)
927  , glob_offset_sequence(-1)
928  , gtopo(comm)
929 {
930  MyComm = comm;
931  MPI_Comm_size(MyComm, &NRanks);
932  MPI_Comm_rank(MyComm, &MyRank);
933 
934  have_face_nbr_data = false;
935  pncmesh = NULL;
936 
937  const int gen_edges = 1;
938 
939  Load(input, gen_edges, refine, true);
940 }
941 
942 void ParMesh::Load(istream &input, int generate_edges, int refine,
943  bool fix_orientation)
944 {
946 
947  // Tell Loader() to read up to 'mfem_serial_mesh_end' instead of
948  // 'mfem_mesh_end', as we have additional parallel mesh data to load in from
949  // the stream.
950  Loader(input, generate_edges, "mfem_serial_mesh_end");
951 
952  ReduceMeshGen(); // determine the global 'meshgen'
953 
954  if (Conforming())
955  {
956  LoadSharedEntities(input);
957  }
958  else
959  {
960  // the ParNCMesh instance was already constructed in 'Loader'
961  pncmesh = dynamic_cast<ParNCMesh*>(ncmesh);
962  MFEM_ASSERT(pncmesh, "internal error");
963 
964  // in the NC case we don't need to load extra data from the file,
965  // as the shared entities can be constructed from the ghost layer
967  }
968 
969  Finalize(refine, fix_orientation);
970 
971  EnsureParNodes();
972 
973  // note: attributes and bdr_attributes are local lists
974 
975  // TODO: NURBS meshes?
976 }
977 
978 void ParMesh::LoadSharedEntities(istream &input)
979 {
980  string ident;
981  skip_comment_lines(input, '#');
982 
983  // read the group topology
984  input >> ident;
985  MFEM_VERIFY(ident == "communication_groups",
986  "input stream is not a parallel MFEM mesh");
987  gtopo.Load(input);
988 
989  skip_comment_lines(input, '#');
990 
991  // read and set the sizes of svert_lvert, group_svert
992  {
993  int num_sverts;
994  input >> ident >> num_sverts;
995  MFEM_VERIFY(ident == "total_shared_vertices", "invalid mesh file");
996  svert_lvert.SetSize(num_sverts);
997  group_svert.SetDims(GetNGroups()-1, num_sverts);
998  }
999  // read and set the sizes of sedge_ledge, group_sedge
1000  if (Dim >= 2)
1001  {
1002  skip_comment_lines(input, '#');
1003  int num_sedges;
1004  input >> ident >> num_sedges;
1005  MFEM_VERIFY(ident == "total_shared_edges", "invalid mesh file");
1006  sedge_ledge.SetSize(num_sedges);
1007  shared_edges.SetSize(num_sedges);
1008  group_sedge.SetDims(GetNGroups()-1, num_sedges);
1009  }
1010  else
1011  {
1012  group_sedge.SetSize(GetNGroups()-1, 0); // create empty group_sedge
1013  }
1014  // read and set the sizes of sface_lface, group_{stria,squad}
1015  if (Dim >= 3)
1016  {
1017  skip_comment_lines(input, '#');
1018  int num_sface;
1019  input >> ident >> num_sface;
1020  MFEM_VERIFY(ident == "total_shared_faces", "invalid mesh file");
1021  sface_lface.SetSize(num_sface);
1024  }
1025  else
1026  {
1027  group_stria.SetSize(GetNGroups()-1, 0); // create empty group_stria
1028  group_squad.SetSize(GetNGroups()-1, 0); // create empty group_squad
1029  }
1030 
1031  // read, group by group, the contents of group_svert, svert_lvert,
1032  // group_sedge, shared_edges, group_{stria,squad}, shared_{trias,quads}
1033  int svert_counter = 0, sedge_counter = 0;
1034  for (int gr = 1; gr < GetNGroups(); gr++)
1035  {
1036  skip_comment_lines(input, '#');
1037 #if 0
1038  // implementation prior to prism-dev merge
1039  int g;
1040  input >> ident >> g; // group
1041  if (g != gr)
1042  {
1043  mfem::err << "ParMesh::ParMesh : expecting group " << gr
1044  << ", read group " << g << endl;
1045  mfem_error();
1046  }
1047 #endif
1048  {
1049  int nv;
1050  input >> ident >> nv; // shared_vertices (in this group)
1051  MFEM_VERIFY(ident == "shared_vertices", "invalid mesh file");
1052  nv += svert_counter;
1053  MFEM_VERIFY(nv <= group_svert.Size_of_connections(),
1054  "incorrect number of total_shared_vertices");
1055  group_svert.GetI()[gr] = nv;
1056  for ( ; svert_counter < nv; svert_counter++)
1057  {
1058  group_svert.GetJ()[svert_counter] = svert_counter;
1059  input >> svert_lvert[svert_counter];
1060  }
1061  }
1062  if (Dim >= 2)
1063  {
1064  int ne, v[2];
1065  input >> ident >> ne; // shared_edges (in this group)
1066  MFEM_VERIFY(ident == "shared_edges", "invalid mesh file");
1067  ne += sedge_counter;
1068  MFEM_VERIFY(ne <= group_sedge.Size_of_connections(),
1069  "incorrect number of total_shared_edges");
1070  group_sedge.GetI()[gr] = ne;
1071  for ( ; sedge_counter < ne; sedge_counter++)
1072  {
1073  group_sedge.GetJ()[sedge_counter] = sedge_counter;
1074  input >> v[0] >> v[1];
1075  shared_edges[sedge_counter] = new Segment(v[0], v[1], 1);
1076  }
1077  }
1078  if (Dim >= 3)
1079  {
1080  int nf, tstart = shared_trias.Size(), qstart = shared_quads.Size();
1081  input >> ident >> nf; // shared_faces (in this group)
1082  for (int i = 0; i < nf; i++)
1083  {
1084  int geom, *v;
1085  input >> geom;
1086  switch (geom)
1087  {
1088  case Geometry::TRIANGLE:
1089  shared_trias.SetSize(shared_trias.Size()+1);
1090  v = shared_trias.Last().v;
1091  for (int ii = 0; ii < 3; ii++) { input >> v[ii]; }
1092  break;
1093  case Geometry::SQUARE:
1094  shared_quads.SetSize(shared_quads.Size()+1);
1095  v = shared_quads.Last().v;
1096  for (int ii = 0; ii < 4; ii++) { input >> v[ii]; }
1097  break;
1098  default:
1099  MFEM_ABORT("invalid shared face geometry: " << geom);
1100  }
1101  }
1102  group_stria.AddColumnsInRow(gr-1, shared_trias.Size()-tstart);
1103  group_squad.AddColumnsInRow(gr-1, shared_quads.Size()-qstart);
1104  }
1105  }
1106  if (Dim >= 3)
1107  {
1108  MFEM_VERIFY(shared_trias.Size() + shared_quads.Size()
1109  == sface_lface.Size(),
1110  "incorrect number of total_shared_faces");
1111  // Define the J arrays of group_stria and group_squad -- they just contain
1112  // consecutive numbers starting from 0 up to shared_trias.Size()-1 and
1113  // shared_quads.Size()-1, respectively.
1114  group_stria.MakeJ();
1115  for (int i = 0; i < shared_trias.Size(); i++)
1116  {
1117  group_stria.GetJ()[i] = i;
1118  }
1119  group_squad.MakeJ();
1120  for (int i = 0; i < shared_quads.Size(); i++)
1121  {
1122  group_squad.GetJ()[i] = i;
1123  }
1124  }
1125 }
1126 
1127 ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
1128 {
1129  MakeRefined_(*orig_mesh, ref_factor, ref_type);
1130 }
1131 
1132 void ParMesh::MakeRefined_(ParMesh &orig_mesh, int ref_factor, int ref_type)
1133 {
1134  MyComm = orig_mesh.GetComm();
1135  NRanks = orig_mesh.GetNRanks();
1136  MyRank = orig_mesh.GetMyRank();
1137  face_nbr_el_to_face = NULL;
1138  glob_elem_offset = -1;
1139  glob_offset_sequence = -1;
1140  gtopo = orig_mesh.gtopo;
1141  have_face_nbr_data = false;
1142  pncmesh = NULL;
1143 
1144  Array<int> ref_factors(orig_mesh.GetNE());
1145  ref_factors = ref_factor;
1146  Mesh::MakeRefined_(orig_mesh, ref_factors, ref_type);
1147 
1148  // Need to initialize:
1149  // - shared_edges, shared_{trias,quads}
1150  // - group_svert, group_sedge, group_{stria,squad}
1151  // - svert_lvert, sedge_ledge, sface_lface
1152 
1153  meshgen = orig_mesh.meshgen; // copy the global 'meshgen'
1154 
1155  H1_FECollection rfec(ref_factor, Dim, ref_type);
1156  ParFiniteElementSpace rfes(&orig_mesh, &rfec);
1157 
1158  // count the number of entries in each row of group_s{vert,edge,face}
1159  group_svert.MakeI(GetNGroups()-1); // exclude the local group 0
1163  for (int gr = 1; gr < GetNGroups(); gr++)
1164  {
1165  // orig vertex -> vertex
1166  group_svert.AddColumnsInRow(gr-1, orig_mesh.GroupNVertices(gr));
1167  // orig edge -> (ref_factor-1) vertices and (ref_factor) edges
1168  const int orig_ne = orig_mesh.GroupNEdges(gr);
1169  group_svert.AddColumnsInRow(gr-1, (ref_factor-1)*orig_ne);
1170  group_sedge.AddColumnsInRow(gr-1, ref_factor*orig_ne);
1171  // orig face -> (?) vertices, (?) edges, and (?) faces
1172  const int orig_nt = orig_mesh.GroupNTriangles(gr);
1173  if (orig_nt > 0)
1174  {
1175  const Geometry::Type geom = Geometry::TRIANGLE;
1176  const int nvert = Geometry::NumVerts[geom];
1177  RefinedGeometry &RG =
1178  *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1179 
1180  // count internal vertices
1181  group_svert.AddColumnsInRow(gr-1, orig_nt*rfec.DofForGeometry(geom));
1182  // count internal edges
1183  group_sedge.AddColumnsInRow(gr-1, orig_nt*(RG.RefEdges.Size()/2-
1184  RG.NumBdrEdges));
1185  // count refined faces
1186  group_stria.AddColumnsInRow(gr-1, orig_nt*(RG.RefGeoms.Size()/nvert));
1187  }
1188  const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1189  if (orig_nq > 0)
1190  {
1191  const Geometry::Type geom = Geometry::SQUARE;
1192  const int nvert = Geometry::NumVerts[geom];
1193  RefinedGeometry &RG =
1194  *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1195 
1196  // count internal vertices
1197  group_svert.AddColumnsInRow(gr-1, orig_nq*rfec.DofForGeometry(geom));
1198  // count internal edges
1199  group_sedge.AddColumnsInRow(gr-1, orig_nq*(RG.RefEdges.Size()/2-
1200  RG.NumBdrEdges));
1201  // count refined faces
1202  group_squad.AddColumnsInRow(gr-1, orig_nq*(RG.RefGeoms.Size()/nvert));
1203  }
1204  }
1205 
1206  group_svert.MakeJ();
1208 
1209  group_sedge.MakeJ();
1212 
1213  group_stria.MakeJ();
1214  group_squad.MakeJ();
1217  sface_lface.SetSize(shared_trias.Size() + shared_quads.Size());
1218 
1219  Array<int> rdofs;
1220  for (int gr = 1; gr < GetNGroups(); gr++)
1221  {
1222  // add shared vertices from original shared vertices
1223  const int orig_n_verts = orig_mesh.GroupNVertices(gr);
1224  for (int j = 0; j < orig_n_verts; j++)
1225  {
1226  rfes.GetVertexDofs(orig_mesh.GroupVertex(gr, j), rdofs);
1227  group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[0])-1);
1228  }
1229 
1230  // add refined shared edges; add shared vertices from refined shared edges
1231  const int orig_n_edges = orig_mesh.GroupNEdges(gr);
1232  if (orig_n_edges > 0)
1233  {
1234  const Geometry::Type geom = Geometry::SEGMENT;
1235  const int nvert = Geometry::NumVerts[geom];
1236  RefinedGeometry &RG = *GlobGeometryRefiner.Refine(geom, ref_factor);
1237  const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1238 
1239  for (int e = 0; e < orig_n_edges; e++)
1240  {
1241  rfes.GetSharedEdgeDofs(gr, e, rdofs);
1242  MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1243  // add the internal edge 'rdofs' as shared vertices
1244  for (int j = 2; j < rdofs.Size(); j++)
1245  {
1246  group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1247  }
1248  for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1249  {
1250  Element *elem = NewElement(geom);
1251  int *v = elem->GetVertices();
1252  for (int k = 0; k < nvert; k++)
1253  {
1254  int cid = RG.RefGeoms[j+k]; // local Cartesian index
1255  v[k] = rdofs[c2h_map[cid]];
1256  }
1257  group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1258  }
1259  }
1260  }
1261  // add refined shared faces; add shared edges and shared vertices from
1262  // refined shared faces
1263  const int orig_nt = orig_mesh.GroupNTriangles(gr);
1264  if (orig_nt > 0)
1265  {
1266  const Geometry::Type geom = Geometry::TRIANGLE;
1267  const int nvert = Geometry::NumVerts[geom];
1268  RefinedGeometry &RG =
1269  *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1270  const int num_int_verts = rfec.DofForGeometry(geom);
1271  const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1272 
1273  for (int f = 0; f < orig_nt; f++)
1274  {
1275  rfes.GetSharedTriangleDofs(gr, f, rdofs);
1276  MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1277  // add the internal face 'rdofs' as shared vertices
1278  for (int j = rdofs.Size()-num_int_verts; j < rdofs.Size(); j++)
1279  {
1280  group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1281  }
1282  // add the internal (for the shared face) edges as shared edges
1283  for (int j = 2*RG.NumBdrEdges; j < RG.RefEdges.Size(); j += 2)
1284  {
1286  int *v = elem->GetVertices();
1287  for (int k = 0; k < 2; k++)
1288  {
1289  v[k] = rdofs[c2h_map[RG.RefEdges[j+k]]];
1290  }
1291  group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1292  }
1293  // add refined shared faces
1294  for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1295  {
1296  shared_trias.SetSize(shared_trias.Size()+1);
1297  int *v = shared_trias.Last().v;
1298  for (int k = 0; k < nvert; k++)
1299  {
1300  int cid = RG.RefGeoms[j+k]; // local Cartesian index
1301  v[k] = rdofs[c2h_map[cid]];
1302  }
1303  group_stria.AddConnection(gr-1, shared_trias.Size()-1);
1304  }
1305  }
1306  }
1307  const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1308  if (orig_nq > 0)
1309  {
1310  const Geometry::Type geom = Geometry::SQUARE;
1311  const int nvert = Geometry::NumVerts[geom];
1312  RefinedGeometry &RG =
1313  *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1314  const int num_int_verts = rfec.DofForGeometry(geom);
1315  const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1316 
1317  for (int f = 0; f < orig_nq; f++)
1318  {
1319  rfes.GetSharedQuadrilateralDofs(gr, f, rdofs);
1320  MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1321  // add the internal face 'rdofs' as shared vertices
1322  for (int j = rdofs.Size()-num_int_verts; j < rdofs.Size(); j++)
1323  {
1324  group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1325  }
1326  // add the internal (for the shared face) edges as shared edges
1327  for (int j = 2*RG.NumBdrEdges; j < RG.RefEdges.Size(); j += 2)
1328  {
1330  int *v = elem->GetVertices();
1331  for (int k = 0; k < 2; k++)
1332  {
1333  v[k] = rdofs[c2h_map[RG.RefEdges[j+k]]];
1334  }
1335  group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1336  }
1337  // add refined shared faces
1338  for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1339  {
1340  shared_quads.SetSize(shared_quads.Size()+1);
1341  int *v = shared_quads.Last().v;
1342  for (int k = 0; k < nvert; k++)
1343  {
1344  int cid = RG.RefGeoms[j+k]; // local Cartesian index
1345  v[k] = rdofs[c2h_map[cid]];
1346  }
1347  group_squad.AddConnection(gr-1, shared_quads.Size()-1);
1348  }
1349  }
1350  }
1351  }
1356 
1357  FinalizeParTopo();
1358 
1359  if (Nodes != NULL)
1360  {
1361  // This call will turn the Nodes into a ParGridFunction
1362  SetCurvature(1, GetNodalFESpace()->IsDGSpace(), spaceDim,
1363  GetNodalFESpace()->GetOrdering());
1364  }
1365 }
1366 
1367 ParMesh ParMesh::MakeRefined(ParMesh &orig_mesh, int ref_factor, int ref_type)
1368 {
1369  ParMesh mesh;
1370  mesh.MakeRefined_(orig_mesh, ref_factor, ref_type);
1371  return mesh;
1372 }
1373 
1375 {
1376  ParMesh mesh;
1377 
1378  mesh.MyComm = orig_mesh.GetComm();
1379  mesh.NRanks = orig_mesh.GetNRanks();
1380  mesh.MyRank = orig_mesh.GetMyRank();
1381  mesh.glob_elem_offset = -1;
1382  mesh.glob_offset_sequence = -1;
1383  mesh.gtopo = orig_mesh.gtopo;
1384  mesh.have_face_nbr_data = false;
1385  mesh.pncmesh = NULL;
1386  mesh.meshgen = orig_mesh.meshgen;
1387 
1388  H1_FECollection fec(1, orig_mesh.Dimension());
1389  ParFiniteElementSpace fes(&orig_mesh, &fec);
1390 
1391  Array<int> vglobal(orig_mesh.GetNV());
1392  for (int iv=0; iv<orig_mesh.GetNV(); ++iv)
1393  {
1394  vglobal[iv] = fes.GetGlobalTDofNumber(iv);
1395  }
1396  mesh.MakeSimplicial_(orig_mesh, vglobal);
1397 
1398  // count the number of entries in each row of group_s{vert,edge,face}
1399  mesh.group_svert.MakeI(mesh.GetNGroups()-1); // exclude the local group 0
1400  mesh.group_sedge.MakeI(mesh.GetNGroups()-1);
1401  mesh.group_stria.MakeI(mesh.GetNGroups()-1);
1402  mesh.group_squad.MakeI(mesh.GetNGroups()-1);
1403  for (int gr = 1; gr < mesh.GetNGroups(); gr++)
1404  {
1405  mesh.group_svert.AddColumnsInRow(gr-1, orig_mesh.GroupNVertices(gr));
1406  mesh.group_sedge.AddColumnsInRow(gr-1, orig_mesh.GroupNEdges(gr));
1407  // Every quad gives an extra edge
1408  const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1409  mesh.group_sedge.AddColumnsInRow(gr-1, orig_nq);
1410  // Every quad is subdivided into two triangles
1411  mesh.group_stria.AddColumnsInRow(gr-1, 2*orig_nq);
1412  // Existing triangles remain unchanged
1413  const int orig_nt = orig_mesh.GroupNTriangles(gr);
1414  mesh.group_stria.AddColumnsInRow(gr-1, orig_nt);
1415  }
1416  mesh.group_svert.MakeJ();
1418 
1419  mesh.group_sedge.MakeJ();
1420  mesh.shared_edges.Reserve(mesh.group_sedge.Size_of_connections());
1422 
1423  mesh.group_stria.MakeJ();
1424  mesh.shared_trias.Reserve(mesh.group_stria.Size_of_connections());
1425  mesh.sface_lface.SetSize(mesh.shared_trias.Size());
1426 
1427  mesh.group_squad.MakeJ();
1428 
1429  constexpr int ntris = 2, nv_tri = 3, nv_quad = 4;
1430 
1431  Array<int> dofs;
1432  for (int gr = 1; gr < mesh.GetNGroups(); gr++)
1433  {
1434  // add shared vertices from original shared vertices
1435  const int orig_n_verts = orig_mesh.GroupNVertices(gr);
1436  for (int j = 0; j < orig_n_verts; j++)
1437  {
1438  fes.GetVertexDofs(orig_mesh.GroupVertex(gr, j), dofs);
1439  mesh.group_svert.AddConnection(gr-1, mesh.svert_lvert.Append(dofs[0])-1);
1440  }
1441 
1442  // add original shared edges
1443  const int orig_n_edges = orig_mesh.GroupNEdges(gr);
1444  for (int e = 0; e < orig_n_edges; e++)
1445  {
1446  int iedge, o;
1447  orig_mesh.GroupEdge(gr, e, iedge, o);
1448  Element *elem = mesh.NewElement(Geometry::SEGMENT);
1449  Array<int> edge_verts;
1450  orig_mesh.GetEdgeVertices(iedge, edge_verts);
1451  elem->SetVertices(edge_verts);
1452  mesh.group_sedge.AddConnection(gr-1, mesh.shared_edges.Append(elem)-1);
1453  }
1454  // add original shared triangles
1455  const int orig_nt = orig_mesh.GroupNTriangles(gr);
1456  for (int e = 0; e < orig_nt; e++)
1457  {
1458  int itri, o;
1459  orig_mesh.GroupTriangle(gr, e, itri, o);
1460  const int *v = orig_mesh.GetFace(itri)->GetVertices();
1461  mesh.shared_trias.SetSize(mesh.shared_trias.Size()+1);
1462  int *v2 = mesh.shared_trias.Last().v;
1463  for (int iv=0; iv<nv_tri; ++iv) { v2[iv] = v[iv]; }
1464  mesh.group_stria.AddConnection(gr-1, mesh.shared_trias.Size()-1);
1465  }
1466  // add triangles from split quads and add resulting diagonal edge
1467  const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1468  if (orig_nq > 0)
1469  {
1470  static const int trimap[12] =
1471  {
1472  0, 0, 0, 1,
1473  1, 2, 1, 2,
1474  2, 3, 3, 3
1475  };
1476  static const int diagmap[4] = { 0, 2, 1, 3 };
1477  for (int f = 0; f < orig_nq; ++f)
1478  {
1479  int iquad, o;
1480  orig_mesh.GroupQuadrilateral(gr, f, iquad, o);
1481  const int *v = orig_mesh.GetFace(iquad)->GetVertices();
1482  // Split quad according the smallest (global) vertex
1483  int vg[nv_quad];
1484  for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
1485  int iv_min = std::min_element(vg, vg+nv_quad) - vg;
1486  int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
1487  // Add diagonal
1488  Element *diag = mesh.NewElement(Geometry::SEGMENT);
1489  int *v_diag = diag->GetVertices();
1490  v_diag[0] = v[diagmap[0 + isplit*2]];
1491  v_diag[1] = v[diagmap[1 + isplit*2]];
1492  mesh.group_sedge.AddConnection(gr-1, mesh.shared_edges.Append(diag)-1);
1493  // Add two new triangles
1494  for (int itri=0; itri<ntris; ++itri)
1495  {
1496  mesh.shared_trias.SetSize(mesh.shared_trias.Size()+1);
1497  int *v2 = mesh.shared_trias.Last().v;
1498  for (int iv=0; iv<nv_tri; ++iv)
1499  {
1500  v2[iv] = v[trimap[itri + isplit*2 + iv*ntris*2]];
1501  }
1502  mesh.group_stria.AddConnection(gr-1, mesh.shared_trias.Size()-1);
1503  }
1504  }
1505  }
1506  }
1507  mesh.group_svert.ShiftUpI();
1508  mesh.group_sedge.ShiftUpI();
1509  mesh.group_stria.ShiftUpI();
1510 
1511  mesh.FinalizeParTopo();
1512 
1513  return mesh;
1514 }
1515 
1516 void ParMesh::Finalize(bool refine, bool fix_orientation)
1517 {
1518  const int meshgen_save = meshgen; // Mesh::Finalize() may call SetMeshGen()
1519 
1520  Mesh::Finalize(refine, fix_orientation);
1521 
1522  meshgen = meshgen_save;
1523  // Note: if Mesh::Finalize() calls MarkTetMeshForRefinement() then the
1524  // shared_trias have been rotated as necessary.
1525 
1526  // Setup secondary parallel mesh data: sedge_ledge, sface_lface
1527  FinalizeParTopo();
1528 }
1529 
1530 int ParMesh::GetLocalElementNum(long long global_element_num) const
1531 {
1533  long long local = global_element_num - glob_elem_offset;
1534  if (local < 0 || local >= NumOfElements) { return -1; }
1535  return local;
1536 }
1537 
1538 long long ParMesh::GetGlobalElementNum(int local_element_num) const
1539 {
1541  return glob_elem_offset + local_element_num;
1542 }
1543 
1545 {
1546  // Determine the largest attribute number across all processors
1547  int max_attr = attr.Size() ? attr.Max() : 1 /*allow empty ranks*/;
1548  int glb_max_attr = -1;
1549  MPI_Allreduce(&max_attr, &glb_max_attr, 1, MPI_INT, MPI_MAX, MyComm);
1550 
1551  // Create marker arrays to indicate which attributes are present
1552  // assuming attribute numbers are in the range [1,glb_max_attr].
1553  bool *attr_marker = new bool[glb_max_attr];
1554  bool *glb_attr_marker = new bool[glb_max_attr];
1555  for (int i = 0; i < glb_max_attr; i++)
1556  {
1557  attr_marker[i] = false;
1558  }
1559  for (int i = 0; i < attr.Size(); i++)
1560  {
1561  attr_marker[attr[i] - 1] = true;
1562  }
1563  MPI_Allreduce(attr_marker, glb_attr_marker, glb_max_attr,
1564  MPI_C_BOOL, MPI_LOR, MyComm);
1565  delete [] attr_marker;
1566 
1567  // Translate from the marker array to a unique, sorted list of attributes
1568  attr.SetSize(0);
1569  attr.Reserve(glb_max_attr);
1570  for (int i = 0; i < glb_max_attr; i++)
1571  {
1572  if (glb_attr_marker[i])
1573  {
1574  attr.Append(i + 1);
1575  }
1576  }
1577  delete [] glb_attr_marker;
1578 }
1579 
1581 {
1582  // Determine the attributes occurring in local interior and boundary elements
1584 
1586  if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1587  {
1588  MFEM_WARNING("Non-positive boundary element attributes found!");
1589  }
1590 
1592  if (attributes.Size() > 0 && attributes[0] <= 0)
1593  {
1594  MFEM_WARNING("Non-positive element attributes found!");
1595  }
1596 }
1597 
1599 {
1600  // maximum number of boundary elements over all ranks
1601  int maxNumOfBdrElements;
1602  MPI_Allreduce(&NumOfBdrElements, &maxNumOfBdrElements, 1,
1603  MPI_INT, MPI_MAX, MyComm);
1604  return (maxNumOfBdrElements > 0);
1605 }
1606 
1607 void ParMesh::GroupEdge(int group, int i, int &edge, int &o) const
1608 {
1609  int sedge = group_sedge.GetRow(group-1)[i];
1610  edge = sedge_ledge[sedge];
1611  int *v = shared_edges[sedge]->GetVertices();
1612  o = (v[0] < v[1]) ? (+1) : (-1);
1613 }
1614 
1615 void ParMesh::GroupTriangle(int group, int i, int &face, int &o) const
1616 {
1617  int stria = group_stria.GetRow(group-1)[i];
1618  face = sface_lface[stria];
1619  // face gives the base orientation
1620  MFEM_ASSERT(faces[face]->GetType() == Element::TRIANGLE,
1621  "Expecting a triangular face.");
1622 
1623  o = GetTriOrientation(faces[face]->GetVertices(), shared_trias[stria].v);
1624 }
1625 
1626 void ParMesh::GroupQuadrilateral(int group, int i, int &face, int &o) const
1627 {
1628  int squad = group_squad.GetRow(group-1)[i];
1629  face = sface_lface[shared_trias.Size()+squad];
1630  // face gives the base orientation
1631  MFEM_ASSERT(faces[face]->GetType() == Element::QUADRILATERAL,
1632  "Expecting a quadrilateral face.");
1633 
1634  o = GetQuadOrientation(faces[face]->GetVertices(), shared_quads[squad].v);
1635 }
1636 
1638  GroupCommunicator& sedge_comm) const
1639 {
1640  Table &gr_sedge = sedge_comm.GroupLDofTable();
1641  gr_sedge.SetDims(GetNGroups(), shared_edges.Size());
1642  gr_sedge.GetI()[0] = 0;
1643  for (int gr = 1; gr <= GetNGroups(); gr++)
1644  {
1645  gr_sedge.GetI()[gr] = group_sedge.GetI()[gr-1];
1646  }
1647  for (int k = 0; k < shared_edges.Size(); k++)
1648  {
1649  if (ordering == 1)
1650  {
1651  gr_sedge.GetJ()[k] =k;
1652  }
1653  else
1654  {
1655  gr_sedge.GetJ()[k] = group_sedge.GetJ()[k];
1656  }
1657  }
1658  sedge_comm.Finalize();
1659 }
1660 
1662  GroupCommunicator& svert_comm) const
1663 {
1664  Table &gr_svert = svert_comm.GroupLDofTable();
1665  gr_svert.SetDims(GetNGroups(), svert_lvert.Size());
1666  gr_svert.GetI()[0] = 0;
1667  for (int gr = 1; gr <= GetNGroups(); gr++)
1668  {
1669  gr_svert.GetI()[gr] = group_svert.GetI()[gr-1];
1670  }
1671  for (int k = 0; k < svert_lvert.Size(); k++)
1672  {
1673  if (ordering == 1)
1674  {
1675  gr_svert.GetJ()[k] = k;
1676  }
1677  else
1678  {
1679  gr_svert.GetJ()[k] = group_svert.GetJ()[k];
1680  }
1681  }
1682  svert_comm.Finalize();
1683 }
1684 
1686  GroupCommunicator& squad_comm) const
1687 {
1688  Table &gr_squad = squad_comm.GroupLDofTable();
1689  gr_squad.SetDims(GetNGroups(), shared_quads.Size());
1690  gr_squad.GetI()[0] = 0;
1691  for (int gr = 1; gr <= GetNGroups(); gr++)
1692  {
1693  gr_squad.GetI()[gr] = group_squad.GetI()[gr-1];
1694  }
1695  for (int k = 0; k < shared_quads.Size(); k++)
1696  {
1697  if (ordering == 1)
1698  {
1699  gr_squad.GetJ()[k] = k;
1700  }
1701  else
1702  {
1703  gr_squad.GetJ()[k] = group_squad.GetJ()[k];
1704  }
1705  }
1706  squad_comm.Finalize();
1707 }
1708 
1710  GroupCommunicator& stria_comm) const
1711 {
1712  Table &gr_stria = stria_comm.GroupLDofTable();
1713  gr_stria.SetDims(GetNGroups(), shared_trias.Size());
1714  gr_stria.GetI()[0] = 0;
1715  for (int gr = 1; gr <= GetNGroups(); gr++)
1716  {
1717  gr_stria.GetI()[gr] = group_stria.GetI()[gr-1];
1718  }
1719  for (int k = 0; k < shared_trias.Size(); k++)
1720  {
1721  if (ordering == 1)
1722  {
1723  gr_stria.GetJ()[k] = k;
1724  }
1725  else
1726  {
1727  gr_stria.GetJ()[k] = group_stria.GetJ()[k];
1728  }
1729  }
1730  stria_comm.Finalize();
1731 }
1732 
1734 {
1735  Array<int> order;
1736  GetEdgeOrdering(v_to_v, order); // local edge ordering
1737 
1738  // create a GroupCommunicator on the shared edges
1739  GroupCommunicator sedge_comm(gtopo);
1740  GetSharedEdgeCommunicator(0, sedge_comm);
1741 
1742  Array<int> sedge_ord(shared_edges.Size());
1743  Array<Pair<int,int> > sedge_ord_map(shared_edges.Size());
1744  for (int k = 0; k < shared_edges.Size(); k++)
1745  {
1746  // sedge_ledge may be undefined -- use shared_edges and v_to_v instead
1747  const int sedge = group_sedge.GetJ()[k];
1748  const int *v = shared_edges[sedge]->GetVertices();
1749  sedge_ord[k] = order[v_to_v(v[0], v[1])];
1750  }
1751 
1752  sedge_comm.Bcast<int>(sedge_ord, 1);
1753 
1754  for (int k = 0, gr = 1; gr < GetNGroups(); gr++)
1755  {
1756  const int n = group_sedge.RowSize(gr-1);
1757  if (n == 0) { continue; }
1758  sedge_ord_map.SetSize(n);
1759  for (int j = 0; j < n; j++)
1760  {
1761  sedge_ord_map[j].one = sedge_ord[k+j];
1762  sedge_ord_map[j].two = j;
1763  }
1764  SortPairs<int, int>(sedge_ord_map, n);
1765  for (int j = 0; j < n; j++)
1766  {
1767  const int sedge_from = group_sedge.GetJ()[k+j];
1768  const int *v = shared_edges[sedge_from]->GetVertices();
1769  sedge_ord[k+j] = order[v_to_v(v[0], v[1])];
1770  }
1771  std::sort(&sedge_ord[k], &sedge_ord[k] + n);
1772  for (int j = 0; j < n; j++)
1773  {
1774  const int sedge_to = group_sedge.GetJ()[k+sedge_ord_map[j].two];
1775  const int *v = shared_edges[sedge_to]->GetVertices();
1776  order[v_to_v(v[0], v[1])] = sedge_ord[k+j];
1777  }
1778  k += n;
1779  }
1780 
1781 #ifdef MFEM_DEBUG
1782  {
1783  Array<Pair<int, double> > ilen_len(order.Size());
1784 
1785  for (int i = 0; i < NumOfVertices; i++)
1786  {
1787  for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
1788  {
1789  int j = it.Index();
1790  ilen_len[j].one = order[j];
1791  ilen_len[j].two = GetLength(i, it.Column());
1792  }
1793  }
1794 
1795  SortPairs<int, double>(ilen_len, order.Size());
1796 
1797  double d_max = 0.;
1798  for (int i = 1; i < order.Size(); i++)
1799  {
1800  d_max = std::max(d_max, ilen_len[i-1].two-ilen_len[i].two);
1801  }
1802 
1803 #if 0
1804  // Debug message from every MPI rank.
1805  mfem::out << "proc. " << MyRank << '/' << NRanks << ": d_max = " << d_max
1806  << endl;
1807 #else
1808  // Debug message just from rank 0.
1809  double glob_d_max;
1810  MPI_Reduce(&d_max, &glob_d_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
1811  if (MyRank == 0)
1812  {
1813  mfem::out << "glob_d_max = " << glob_d_max << endl;
1814  }
1815 #endif
1816  }
1817 #endif
1818 
1819  // use 'order' to mark the tets, the boundary triangles, and the shared
1820  // triangle faces
1821  for (int i = 0; i < NumOfElements; i++)
1822  {
1823  if (elements[i]->GetType() == Element::TETRAHEDRON)
1824  {
1825  elements[i]->MarkEdge(v_to_v, order);
1826  }
1827  }
1828 
1829  for (int i = 0; i < NumOfBdrElements; i++)
1830  {
1831  if (boundary[i]->GetType() == Element::TRIANGLE)
1832  {
1833  boundary[i]->MarkEdge(v_to_v, order);
1834  }
1835  }
1836 
1837  for (int i = 0; i < shared_trias.Size(); i++)
1838  {
1839  Triangle::MarkEdge(shared_trias[i].v, v_to_v, order);
1840  }
1841 }
1842 
1843 // For a line segment with vertices v[0] and v[1], return a number with
1844 // the following meaning:
1845 // 0 - the edge was not refined
1846 // 1 - the edge e was refined once by splitting v[0],v[1]
1847 int ParMesh::GetEdgeSplittings(Element *edge, const DSTable &v_to_v,
1848  int *middle)
1849 {
1850  int m, *v = edge->GetVertices();
1851 
1852  if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
1853  {
1854  return 1;
1855  }
1856  else
1857  {
1858  return 0;
1859  }
1860 }
1861 
1862 void ParMesh::GetFaceSplittings(const int *fv, const HashTable<Hashed2> &v_to_v,
1863  Array<unsigned> &codes)
1864 {
1865  typedef Triple<int,int,int> face_t;
1866  Array<face_t> face_stack;
1867 
1868  unsigned code = 0;
1869  face_stack.Append(face_t(fv[0], fv[1], fv[2]));
1870  for (unsigned bit = 0; face_stack.Size() > 0; bit++)
1871  {
1872  if (bit == 8*sizeof(unsigned))
1873  {
1874  codes.Append(code);
1875  code = bit = 0;
1876  }
1877 
1878  const face_t &f = face_stack.Last();
1879  int mid = v_to_v.FindId(f.one, f.two);
1880  if (mid == -1)
1881  {
1882  // leave a 0 at bit 'bit'
1883  face_stack.DeleteLast();
1884  }
1885  else
1886  {
1887  code += (1 << bit); // set bit 'bit' to 1
1888  mid += NumOfVertices;
1889  face_stack.Append(face_t(f.three, f.one, mid));
1890  face_t &r = face_stack[face_stack.Size()-2];
1891  r = face_t(r.two, r.three, mid);
1892  }
1893  }
1894  codes.Append(code);
1895 }
1896 
1898  const Array<unsigned> &codes, int &pos)
1899 {
1900  typedef Triple<int,int,int> face_t;
1901  Array<face_t> face_stack;
1902 
1903  bool need_refinement = 0;
1904  face_stack.Append(face_t(v[0], v[1], v[2]));
1905  for (unsigned bit = 0, code = codes[pos++]; face_stack.Size() > 0; bit++)
1906  {
1907  if (bit == 8*sizeof(unsigned))
1908  {
1909  code = codes[pos++];
1910  bit = 0;
1911  }
1912 
1913  if ((code & (1 << bit)) == 0) { face_stack.DeleteLast(); continue; }
1914 
1915  const face_t &f = face_stack.Last();
1916  int mid = v_to_v.FindId(f.one, f.two);
1917  if (mid == -1)
1918  {
1919  mid = v_to_v.GetId(f.one, f.two);
1920  int ind[2] = { f.one, f.two };
1921  vertices.Append(Vertex());
1922  AverageVertices(ind, 2, vertices.Size()-1);
1923  need_refinement = 1;
1924  }
1925  mid += NumOfVertices;
1926  face_stack.Append(face_t(f.three, f.one, mid));
1927  face_t &r = face_stack[face_stack.Size()-2];
1928  r = face_t(r.two, r.three, mid);
1929  }
1930  return need_refinement;
1931 }
1932 
1933 void ParMesh::GenerateOffsets(int N, HYPRE_BigInt loc_sizes[],
1934  Array<HYPRE_BigInt> *offsets[]) const
1935 {
1936  if (HYPRE_AssumedPartitionCheck())
1937  {
1938  Array<HYPRE_BigInt> temp(N);
1939  MPI_Scan(loc_sizes, temp.GetData(), N, HYPRE_MPI_BIG_INT, MPI_SUM, MyComm);
1940  for (int i = 0; i < N; i++)
1941  {
1942  offsets[i]->SetSize(3);
1943  (*offsets[i])[0] = temp[i] - loc_sizes[i];
1944  (*offsets[i])[1] = temp[i];
1945  }
1946  MPI_Bcast(temp.GetData(), N, HYPRE_MPI_BIG_INT, NRanks-1, MyComm);
1947  for (int i = 0; i < N; i++)
1948  {
1949  (*offsets[i])[2] = temp[i];
1950  // check for overflow
1951  MFEM_VERIFY((*offsets[i])[0] >= 0 && (*offsets[i])[1] >= 0,
1952  "overflow in offsets");
1953  }
1954  }
1955  else
1956  {
1957  Array<HYPRE_BigInt> temp(N*NRanks);
1958  MPI_Allgather(loc_sizes, N, HYPRE_MPI_BIG_INT, temp.GetData(), N,
1959  HYPRE_MPI_BIG_INT, MyComm);
1960  for (int i = 0; i < N; i++)
1961  {
1962  Array<HYPRE_BigInt> &offs = *offsets[i];
1963  offs.SetSize(NRanks+1);
1964  offs[0] = 0;
1965  for (int j = 0; j < NRanks; j++)
1966  {
1967  offs[j+1] = offs[j] + temp[i+N*j];
1968  }
1969  // Check for overflow
1970  MFEM_VERIFY(offs[MyRank] >= 0 && offs[MyRank+1] >= 0,
1971  "overflow in offsets");
1972  }
1973  }
1974 }
1975 
1977  int i, IsoparametricTransformation *ElTr)
1978 {
1979  DenseMatrix &pointmat = ElTr->GetPointMat();
1980  Element *elem = face_nbr_elements[i];
1981 
1982  ElTr->Attribute = elem->GetAttribute();
1983  ElTr->ElementNo = NumOfElements + i;
1985  ElTr->mesh = this;
1986  ElTr->Reset();
1987 
1988  if (Nodes == NULL)
1989  {
1990  const int nv = elem->GetNVertices();
1991  const int *v = elem->GetVertices();
1992 
1993  pointmat.SetSize(spaceDim, nv);
1994  for (int k = 0; k < spaceDim; k++)
1995  {
1996  for (int j = 0; j < nv; j++)
1997  {
1998  pointmat(k, j) = face_nbr_vertices[v[j]](k);
1999  }
2000  }
2001 
2003  }
2004  else
2005  {
2006  Array<int> vdofs;
2007  ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
2008  if (pNodes)
2009  {
2010  pNodes->ParFESpace()->GetFaceNbrElementVDofs(i, vdofs);
2011  int n = vdofs.Size()/spaceDim;
2012  pointmat.SetSize(spaceDim, n);
2013  for (int k = 0; k < spaceDim; k++)
2014  {
2015  for (int j = 0; j < n; j++)
2016  {
2017  pointmat(k,j) = (pNodes->FaceNbrData())(vdofs[n*k+j]);
2018  }
2019  }
2020 
2021  ElTr->SetFE(pNodes->ParFESpace()->GetFaceNbrFE(i));
2022  }
2023  else
2024  {
2025  MFEM_ABORT("Nodes are not ParGridFunction!");
2026  }
2027  }
2028 }
2029 
2030 double ParMesh::GetFaceNbrElementSize(int i, int type)
2031 {
2033 }
2034 
2036 {
2037  if (!have_face_nbr_data)
2038  {
2039  return;
2040  }
2041 
2042  have_face_nbr_data = false;
2046  for (int i = 0; i < face_nbr_elements.Size(); i++)
2047  {
2049  }
2050  face_nbr_elements.DeleteAll();
2051  face_nbr_vertices.DeleteAll();
2054 }
2055 
2056 void ParMesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
2057 {
2058  space_dim = (space_dim == -1) ? spaceDim : space_dim;
2060  if (discont)
2061  {
2062  nfec = new L2_FECollection(order, Dim, BasisType::GaussLobatto);
2063  }
2064  else
2065  {
2066  nfec = new H1_FECollection(order, Dim);
2067  }
2068  ParFiniteElementSpace* nfes = new ParFiniteElementSpace(this, nfec, space_dim,
2069  ordering);
2070  auto pnodes = new ParGridFunction(nfes);
2071  GetNodes(*pnodes);
2072  NewNodes(*pnodes, true);
2073  Nodes->MakeOwner(nfec);
2074 }
2075 
2077 {
2078  ParFiniteElementSpace *npfes = dynamic_cast<ParFiniteElementSpace*>(nfes);
2079  if (npfes)
2080  {
2081  SetNodalFESpace(npfes);
2082  }
2083  else
2084  {
2085  Mesh::SetNodalFESpace(nfes);
2086  }
2087 }
2088 
2090 {
2091  ParGridFunction *nodes = new ParGridFunction(npfes);
2092  SetNodalGridFunction(nodes, true);
2093 }
2094 
2096 {
2097  if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
2098  {
2099  ParFiniteElementSpace *pfes =
2100  new ParFiniteElementSpace(*Nodes->FESpace(), *this);
2101  ParGridFunction *new_nodes = new ParGridFunction(pfes);
2102 
2103  *new_nodes = *Nodes;
2104 
2105  if (Nodes->OwnFEC())
2106  {
2107  new_nodes->MakeOwner(Nodes->OwnFEC());
2108  Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
2109  delete Nodes->FESpace();
2110  }
2111 
2112  delete Nodes;
2113  Nodes = new_nodes;
2114  }
2115 }
2116 
2118 {
2119  if (have_face_nbr_data)
2120  {
2121  return;
2122  }
2123 
2124  if (Nonconforming())
2125  {
2126  // with ParNCMesh we can set up face neighbors without communication
2127  pncmesh->GetFaceNeighbors(*this);
2128  have_face_nbr_data = true;
2129 
2131  return;
2132  }
2133 
2134  Table *gr_sface;
2135  int *s2l_face;
2136  bool del_tables = false;
2137  if (Dim == 1)
2138  {
2139  gr_sface = &group_svert;
2140  s2l_face = svert_lvert;
2141  }
2142  else if (Dim == 2)
2143  {
2144  gr_sface = &group_sedge;
2145  s2l_face = sedge_ledge;
2146  }
2147  else
2148  {
2149  s2l_face = sface_lface;
2150  if (shared_trias.Size() == sface_lface.Size())
2151  {
2152  // All shared faces are Triangular
2153  gr_sface = &group_stria;
2154  }
2155  else if (shared_quads.Size() == sface_lface.Size())
2156  {
2157  // All shared faced are Quadrilateral
2158  gr_sface = &group_squad;
2159  }
2160  else
2161  {
2162  // Shared faces contain a mixture of triangles and quads
2163  gr_sface = new Table;
2164  del_tables = true;
2165 
2166  // Merge the Tables group_stria and group_squad
2167  gr_sface->MakeI(group_stria.Size());
2168  for (int gr=0; gr<group_stria.Size(); gr++)
2169  {
2170  gr_sface->AddColumnsInRow(gr,
2171  group_stria.RowSize(gr) +
2172  group_squad.RowSize(gr));
2173  }
2174  gr_sface->MakeJ();
2175  const int nst = shared_trias.Size();
2176  for (int gr=0; gr<group_stria.Size(); gr++)
2177  {
2178  gr_sface->AddConnections(gr,
2179  group_stria.GetRow(gr),
2180  group_stria.RowSize(gr));
2181  for (int c=0; c<group_squad.RowSize(gr); c++)
2182  {
2183  gr_sface->AddConnection(gr,
2184  nst + group_squad.GetRow(gr)[c]);
2185  }
2186  }
2187  gr_sface->ShiftUpI();
2188  }
2189  }
2190 
2191  ExchangeFaceNbrData(gr_sface, s2l_face);
2192 
2193  if (Dim == 3)
2194  {
2196  }
2197 
2198  if (del_tables) { delete gr_sface; }
2199 
2200  if ( have_face_nbr_data ) { return; }
2201 
2202  have_face_nbr_data = true;
2203 
2205 }
2206 
2207 void ParMesh::ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
2208 {
2209  int num_face_nbrs = 0;
2210  for (int g = 1; g < GetNGroups(); g++)
2211  {
2212  if (gr_sface->RowSize(g-1) > 0)
2213  {
2214  num_face_nbrs++;
2215  }
2216  }
2217 
2218  face_nbr_group.SetSize(num_face_nbrs);
2219 
2220  if (num_face_nbrs == 0)
2221  {
2222  have_face_nbr_data = true;
2223  return;
2224  }
2225 
2226  {
2227  // sort face-neighbors by processor rank
2228  Array<Pair<int, int> > rank_group(num_face_nbrs);
2229 
2230  for (int g = 1, counter = 0; g < GetNGroups(); g++)
2231  {
2232  if (gr_sface->RowSize(g-1) > 0)
2233  {
2234  MFEM_ASSERT(gtopo.GetGroupSize(g) == 2, "group size is not 2!");
2235 
2236  const int *nbs = gtopo.GetGroup(g);
2237  int lproc = (nbs[0]) ? nbs[0] : nbs[1];
2238  rank_group[counter].one = gtopo.GetNeighborRank(lproc);
2239  rank_group[counter].two = g;
2240  counter++;
2241  }
2242  }
2243 
2244  SortPairs<int, int>(rank_group, rank_group.Size());
2245 
2246  for (int fn = 0; fn < num_face_nbrs; fn++)
2247  {
2248  face_nbr_group[fn] = rank_group[fn].two;
2249  }
2250  }
2251 
2252  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
2253  MPI_Request *send_requests = requests;
2254  MPI_Request *recv_requests = requests + num_face_nbrs;
2255  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
2256 
2257  int *nbr_data = new int[6*num_face_nbrs];
2258  int *nbr_send_data = nbr_data;
2259  int *nbr_recv_data = nbr_data + 3*num_face_nbrs;
2260 
2261  Array<int> el_marker(GetNE());
2262  Array<int> vertex_marker(GetNV());
2263  el_marker = -1;
2264  vertex_marker = -1;
2265 
2266  Array<int> fcs, cor;
2267 
2268  Table send_face_nbr_elemdata, send_face_nbr_facedata;
2269 
2270  send_face_nbr_elements.MakeI(num_face_nbrs);
2271  send_face_nbr_vertices.MakeI(num_face_nbrs);
2272  send_face_nbr_elemdata.MakeI(num_face_nbrs);
2273  send_face_nbr_facedata.MakeI(num_face_nbrs);
2274  for (int fn = 0; fn < num_face_nbrs; fn++)
2275  {
2276  int nbr_group = face_nbr_group[fn];
2277  int num_sfaces = gr_sface->RowSize(nbr_group-1);
2278  int *sface = gr_sface->GetRow(nbr_group-1);
2279  for (int i = 0; i < num_sfaces; i++)
2280  {
2281  int lface = s2l_face[sface[i]];
2282  int el = faces_info[lface].Elem1No;
2283  if (el_marker[el] != fn)
2284  {
2285  el_marker[el] = fn;
2287 
2288  const int nv = elements[el]->GetNVertices();
2289  const int *v = elements[el]->GetVertices();
2290  for (int j = 0; j < nv; j++)
2291  if (vertex_marker[v[j]] != fn)
2292  {
2293  vertex_marker[v[j]] = fn;
2295  }
2296 
2297  const int nf = elements[el]->GetNFaces();
2298 
2299  send_face_nbr_elemdata.AddColumnsInRow(fn, nv + nf + 2);
2300  }
2301  }
2302  send_face_nbr_facedata.AddColumnsInRow(fn, 2*num_sfaces);
2303 
2304  nbr_send_data[3*fn ] = send_face_nbr_elements.GetI()[fn];
2305  nbr_send_data[3*fn+1] = send_face_nbr_vertices.GetI()[fn];
2306  nbr_send_data[3*fn+2] = send_face_nbr_elemdata.GetI()[fn];
2307 
2308  int nbr_rank = GetFaceNbrRank(fn);
2309  int tag = 0;
2310 
2311  MPI_Isend(&nbr_send_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
2312  &send_requests[fn]);
2313  MPI_Irecv(&nbr_recv_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
2314  &recv_requests[fn]);
2315  }
2318  send_face_nbr_elemdata.MakeJ();
2319  send_face_nbr_facedata.MakeJ();
2320  el_marker = -1;
2321  vertex_marker = -1;
2322  const int nst = shared_trias.Size();
2323  for (int fn = 0; fn < num_face_nbrs; fn++)
2324  {
2325  int nbr_group = face_nbr_group[fn];
2326  int num_sfaces = gr_sface->RowSize(nbr_group-1);
2327  int *sface = gr_sface->GetRow(nbr_group-1);
2328  for (int i = 0; i < num_sfaces; i++)
2329  {
2330  const int sf = sface[i];
2331  int lface = s2l_face[sf];
2332  int el = faces_info[lface].Elem1No;
2333  if (el_marker[el] != fn)
2334  {
2335  el_marker[el] = fn;
2337 
2338  const int nv = elements[el]->GetNVertices();
2339  const int *v = elements[el]->GetVertices();
2340  for (int j = 0; j < nv; j++)
2341  if (vertex_marker[v[j]] != fn)
2342  {
2343  vertex_marker[v[j]] = fn;
2345  }
2346 
2347  send_face_nbr_elemdata.AddConnection(fn, GetAttribute(el));
2348  send_face_nbr_elemdata.AddConnection(
2349  fn, GetElementBaseGeometry(el));
2350  send_face_nbr_elemdata.AddConnections(fn, v, nv);
2351 
2352  if (Dim == 3)
2353  {
2354  const int nf = elements[el]->GetNFaces();
2355  GetElementFaces(el, fcs, cor);
2356  send_face_nbr_elemdata.AddConnections(fn, cor, nf);
2357  }
2358  }
2359  send_face_nbr_facedata.AddConnection(fn, el);
2360  int info = faces_info[lface].Elem1Inf;
2361  // change the orientation in info to be relative to the shared face
2362  // in 1D and 2D keep the orientation equal to 0
2363  if (Dim == 3)
2364  {
2365  const int *lf_v = faces[lface]->GetVertices();
2366  if (sf < nst) // triangle shared face
2367  {
2368  info += GetTriOrientation(shared_trias[sf].v, lf_v);
2369  }
2370  else // quad shared face
2371  {
2372  info += GetQuadOrientation(shared_quads[sf-nst].v, lf_v);
2373  }
2374  }
2375  send_face_nbr_facedata.AddConnection(fn, info);
2376  }
2377  }
2380  send_face_nbr_elemdata.ShiftUpI();
2381  send_face_nbr_facedata.ShiftUpI();
2382 
2383  // convert the vertex indices in send_face_nbr_elemdata
2384  // convert the element indices in send_face_nbr_facedata
2385  for (int fn = 0; fn < num_face_nbrs; fn++)
2386  {
2387  int num_elems = send_face_nbr_elements.RowSize(fn);
2388  int *elems = send_face_nbr_elements.GetRow(fn);
2389  int num_verts = send_face_nbr_vertices.RowSize(fn);
2390  int *verts = send_face_nbr_vertices.GetRow(fn);
2391  int *elemdata = send_face_nbr_elemdata.GetRow(fn);
2392  int num_sfaces = send_face_nbr_facedata.RowSize(fn)/2;
2393  int *facedata = send_face_nbr_facedata.GetRow(fn);
2394 
2395  for (int i = 0; i < num_verts; i++)
2396  {
2397  vertex_marker[verts[i]] = i;
2398  }
2399 
2400  for (int el = 0; el < num_elems; el++)
2401  {
2402  const int nv = elements[elems[el]]->GetNVertices();
2403  const int nf = (Dim == 3) ? elements[elems[el]]->GetNFaces() : 0;
2404  elemdata += 2; // skip the attribute and the geometry type
2405  for (int j = 0; j < nv; j++)
2406  {
2407  elemdata[j] = vertex_marker[elemdata[j]];
2408  }
2409  elemdata += nv + nf;
2410 
2411  el_marker[elems[el]] = el;
2412  }
2413 
2414  for (int i = 0; i < num_sfaces; i++)
2415  {
2416  facedata[2*i] = el_marker[facedata[2*i]];
2417  }
2418  }
2419 
2420  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2421 
2422  Array<int> recv_face_nbr_facedata;
2423  Table recv_face_nbr_elemdata;
2424 
2425  // fill-in face_nbr_elements_offset, face_nbr_vertices_offset
2426  face_nbr_elements_offset.SetSize(num_face_nbrs + 1);
2427  face_nbr_vertices_offset.SetSize(num_face_nbrs + 1);
2428  recv_face_nbr_elemdata.MakeI(num_face_nbrs);
2429  face_nbr_elements_offset[0] = 0;
2430  face_nbr_vertices_offset[0] = 0;
2431  for (int fn = 0; fn < num_face_nbrs; fn++)
2432  {
2433  face_nbr_elements_offset[fn+1] =
2434  face_nbr_elements_offset[fn] + nbr_recv_data[3*fn];
2435  face_nbr_vertices_offset[fn+1] =
2436  face_nbr_vertices_offset[fn] + nbr_recv_data[3*fn+1];
2437  recv_face_nbr_elemdata.AddColumnsInRow(fn, nbr_recv_data[3*fn+2]);
2438  }
2439  recv_face_nbr_elemdata.MakeJ();
2440 
2441  MPI_Waitall(num_face_nbrs, send_requests, statuses);
2442 
2443  // send and receive the element data
2444  for (int fn = 0; fn < num_face_nbrs; fn++)
2445  {
2446  int nbr_rank = GetFaceNbrRank(fn);
2447  int tag = 0;
2448 
2449  MPI_Isend(send_face_nbr_elemdata.GetRow(fn),
2450  send_face_nbr_elemdata.RowSize(fn),
2451  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
2452 
2453  MPI_Irecv(recv_face_nbr_elemdata.GetRow(fn),
2454  recv_face_nbr_elemdata.RowSize(fn),
2455  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
2456  }
2457 
2458  // convert the element data into face_nbr_elements
2459  face_nbr_elements.SetSize(face_nbr_elements_offset[num_face_nbrs]);
2462  while (true)
2463  {
2464  int fn;
2465  MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2466 
2467  if (fn == MPI_UNDEFINED)
2468  {
2469  break;
2470  }
2471 
2472  int vert_off = face_nbr_vertices_offset[fn];
2473  int elem_off = face_nbr_elements_offset[fn];
2474  int num_elems = face_nbr_elements_offset[fn+1] - elem_off;
2475  int *recv_elemdata = recv_face_nbr_elemdata.GetRow(fn);
2476 
2477  for (int i = 0; i < num_elems; i++)
2478  {
2479  Element *el = NewElement(recv_elemdata[1]);
2480  el->SetAttribute(recv_elemdata[0]);
2481  recv_elemdata += 2;
2482  int nv = el->GetNVertices();
2483  for (int j = 0; j < nv; j++)
2484  {
2485  recv_elemdata[j] += vert_off;
2486  }
2487  el->SetVertices(recv_elemdata);
2488  recv_elemdata += nv;
2489  if (Dim == 3)
2490  {
2491  int nf = el->GetNFaces();
2492  int * fn_ori = face_nbr_el_ori.GetRow(elem_off);
2493  for (int j = 0; j < nf; j++)
2494  {
2495  fn_ori[j] = recv_elemdata[j];
2496  }
2497  recv_elemdata += nf;
2498  }
2499  face_nbr_elements[elem_off++] = el;
2500  }
2501  }
2503 
2504  MPI_Waitall(num_face_nbrs, send_requests, statuses);
2505 
2506  // send and receive the face data
2507  recv_face_nbr_facedata.SetSize(
2508  send_face_nbr_facedata.Size_of_connections());
2509  for (int fn = 0; fn < num_face_nbrs; fn++)
2510  {
2511  int nbr_rank = GetFaceNbrRank(fn);
2512  int tag = 0;
2513 
2514  MPI_Isend(send_face_nbr_facedata.GetRow(fn),
2515  send_face_nbr_facedata.RowSize(fn),
2516  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
2517 
2518  // the size of the send and receive face data is the same
2519  MPI_Irecv(&recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]],
2520  send_face_nbr_facedata.RowSize(fn),
2521  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
2522  }
2523 
2524  // transfer the received face data into faces_info
2525  while (true)
2526  {
2527  int fn;
2528  MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2529 
2530  if (fn == MPI_UNDEFINED)
2531  {
2532  break;
2533  }
2534 
2535  int elem_off = face_nbr_elements_offset[fn];
2536  int nbr_group = face_nbr_group[fn];
2537  int num_sfaces = gr_sface->RowSize(nbr_group-1);
2538  int *sface = gr_sface->GetRow(nbr_group-1);
2539  int *facedata =
2540  &recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]];
2541 
2542  for (int i = 0; i < num_sfaces; i++)
2543  {
2544  const int sf = sface[i];
2545  int lface = s2l_face[sf];
2546  FaceInfo &face_info = faces_info[lface];
2547  face_info.Elem2No = -1 - (facedata[2*i] + elem_off);
2548  int info = facedata[2*i+1];
2549  // change the orientation in info to be relative to the local face
2550  if (Dim < 3)
2551  {
2552  info++; // orientation 0 --> orientation 1
2553  }
2554  else
2555  {
2556  int nbr_ori = info%64, nbr_v[4];
2557  const int *lf_v = faces[lface]->GetVertices();
2558 
2559  if (sf < nst) // triangle shared face
2560  {
2561  // apply the nbr_ori to sf_v to get nbr_v
2562  const int *perm = tri_t::Orient[nbr_ori];
2563  const int *sf_v = shared_trias[sf].v;
2564  for (int j = 0; j < 3; j++)
2565  {
2566  nbr_v[perm[j]] = sf_v[j];
2567  }
2568  // get the orientation of nbr_v w.r.t. the local face
2569  nbr_ori = GetTriOrientation(lf_v, nbr_v);
2570  }
2571  else // quad shared face
2572  {
2573  // apply the nbr_ori to sf_v to get nbr_v
2574  const int *perm = quad_t::Orient[nbr_ori];
2575  const int *sf_v = shared_quads[sf-nst].v;
2576  for (int j = 0; j < 4; j++)
2577  {
2578  nbr_v[perm[j]] = sf_v[j];
2579  }
2580  // get the orientation of nbr_v w.r.t. the local face
2581  nbr_ori = GetQuadOrientation(lf_v, nbr_v);
2582  }
2583 
2584  info = 64*(info/64) + nbr_ori;
2585  }
2586  face_info.Elem2Inf = info;
2587  }
2588  }
2589 
2590  MPI_Waitall(num_face_nbrs, send_requests, statuses);
2591 
2592  // allocate the face_nbr_vertices
2593  face_nbr_vertices.SetSize(face_nbr_vertices_offset[num_face_nbrs]);
2594 
2595  delete [] nbr_data;
2596 
2597  delete [] statuses;
2598  delete [] requests;
2599 }
2600 
2602 {
2603  if (!have_face_nbr_data)
2604  {
2605  ExchangeFaceNbrData(); // calls this method at the end
2606  }
2607  else if (Nodes == NULL)
2608  {
2609  if (Nonconforming())
2610  {
2611  // with ParNCMesh we already have the vertices
2612  return;
2613  }
2614 
2615  int num_face_nbrs = GetNFaceNeighbors();
2616 
2617  if (!num_face_nbrs) { return; }
2618 
2619  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
2620  MPI_Request *send_requests = requests;
2621  MPI_Request *recv_requests = requests + num_face_nbrs;
2622  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
2623 
2624  // allocate buffer and copy the vertices to be sent
2626  for (int i = 0; i < send_vertices.Size(); i++)
2627  {
2628  send_vertices[i] = vertices[send_face_nbr_vertices.GetJ()[i]];
2629  }
2630 
2631  // send and receive the vertices
2632  for (int fn = 0; fn < num_face_nbrs; fn++)
2633  {
2634  int nbr_rank = GetFaceNbrRank(fn);
2635  int tag = 0;
2636 
2637  MPI_Isend(send_vertices[send_face_nbr_vertices.GetI()[fn]](),
2639  MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
2640 
2642  3*(face_nbr_vertices_offset[fn+1] -
2644  MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
2645  }
2646 
2647  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2648  MPI_Waitall(num_face_nbrs, send_requests, statuses);
2649 
2650  delete [] statuses;
2651  delete [] requests;
2652  }
2653  else
2654  {
2655  ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
2656  MFEM_VERIFY(pNodes != NULL, "Nodes are not ParGridFunction!");
2657  pNodes->ExchangeFaceNbrData();
2658  }
2659 }
2660 
2662 {
2663  STable3D *sfaces_tbl = new STable3D(face_nbr_vertices.Size());
2664  for (int i = 0; i < face_nbr_elements.Size(); i++)
2665  {
2666  const int *v = face_nbr_elements[i]->GetVertices();
2667  switch (face_nbr_elements[i]->GetType())
2668  {
2669  case Element::TETRAHEDRON:
2670  {
2671  for (int j = 0; j < 4; j++)
2672  {
2673  const int *fv = tet_t::FaceVert[j];
2674  sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2675  }
2676  break;
2677  }
2678  case Element::WEDGE:
2679  {
2680  for (int j = 0; j < 2; j++)
2681  {
2682  const int *fv = pri_t::FaceVert[j];
2683  sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2684  }
2685  for (int j = 2; j < 5; j++)
2686  {
2687  const int *fv = pri_t::FaceVert[j];
2688  sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2689  }
2690  break;
2691  }
2692  case Element::PYRAMID:
2693  {
2694  for (int j = 0; j < 1; j++)
2695  {
2696  const int *fv = pyr_t::FaceVert[j];
2697  sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2698  }
2699  for (int j = 1; j < 5; j++)
2700  {
2701  const int *fv = pyr_t::FaceVert[j];
2702  sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2703  }
2704  break;
2705  }
2706  case Element::HEXAHEDRON:
2707  {
2708  // find the face by the vertices with the smallest 3 numbers
2709  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
2710  for (int j = 0; j < 6; j++)
2711  {
2712  const int *fv = hex_t::FaceVert[j];
2713  sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2714  }
2715  break;
2716  }
2717  default:
2718  MFEM_ABORT("Unexpected type of Element.");
2719  }
2720  }
2721  return sfaces_tbl;
2722 }
2723 
2725 {
2726  int i, *v;
2727  STable3D * faces_tbl = GetFacesTable();
2728  STable3D * sfaces_tbl = GetSharedFacesTable();
2729 
2730  if (face_nbr_el_to_face != NULL)
2731  {
2732  delete face_nbr_el_to_face;
2733  }
2734  face_nbr_el_to_face = new Table(face_nbr_elements.Size(), 6);
2735  for (i = 0; i < face_nbr_elements.Size(); i++)
2736  {
2737  v = face_nbr_elements[i]->GetVertices();
2738  switch (face_nbr_elements[i]->GetType())
2739  {
2740  case Element::TETRAHEDRON:
2741  {
2742  for (int j = 0; j < 4; j++)
2743  {
2744  const int *fv = tet_t::FaceVert[j];
2745  int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2746  if (lf < 0)
2747  {
2748  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2749  if (lf >= 0)
2750  {
2751  lf += NumOfFaces;
2752  }
2753  }
2754  face_nbr_el_to_face->Push(i, lf);
2755  }
2756  break;
2757  }
2758  case Element::WEDGE:
2759  {
2760  for (int j = 0; j < 2; j++)
2761  {
2762  const int *fv = pri_t::FaceVert[j];
2763  int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2764  if (lf < 0)
2765  {
2766  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2767  if (lf >= 0)
2768  {
2769  lf += NumOfFaces;
2770  }
2771  }
2772  face_nbr_el_to_face->Push(i, lf);
2773  }
2774  for (int j = 2; j < 5; j++)
2775  {
2776  const int *fv = pri_t::FaceVert[j];
2777  int k = 0;
2778  int max = v[fv[0]];
2779 
2780  if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2781  if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2782  if (max < v[fv[3]]) { k = 3; }
2783 
2784  int v0 = -1, v1 = -1, v2 = -1;
2785  switch (k)
2786  {
2787  case 0:
2788  v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2789  break;
2790  case 1:
2791  v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2792  break;
2793  case 2:
2794  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2795  break;
2796  case 3:
2797  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2798  break;
2799  }
2800  int lf = faces_tbl->Index(v0, v1, v2);
2801  if (lf < 0)
2802  {
2803  lf = sfaces_tbl->Index(v0, v1, v2);
2804  if (lf >= 0)
2805  {
2806  lf += NumOfFaces;
2807  }
2808  }
2809  face_nbr_el_to_face->Push(i, lf);
2810  }
2811  break;
2812  }
2813  case Element::PYRAMID:
2814  {
2815  for (int j = 0; j < 1; j++)
2816  {
2817  const int *fv = pyr_t::FaceVert[j];
2818  int k = 0;
2819  int max = v[fv[0]];
2820 
2821  if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2822  if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2823  if (max < v[fv[3]]) { k = 3; }
2824 
2825  int v0 = -1, v1 = -1, v2 = -1;
2826  switch (k)
2827  {
2828  case 0:
2829  v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2830  break;
2831  case 1:
2832  v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2833  break;
2834  case 2:
2835  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2836  break;
2837  case 3:
2838  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2839  break;
2840  }
2841  int lf = faces_tbl->Index(v0, v1, v2);
2842  if (lf < 0)
2843  {
2844  lf = sfaces_tbl->Index(v0, v1, v2);
2845  if (lf >= 0)
2846  {
2847  lf += NumOfFaces;
2848  }
2849  }
2850  face_nbr_el_to_face->Push(i, lf);
2851  }
2852  for (int j = 1; j < 5; j++)
2853  {
2854  const int *fv = pyr_t::FaceVert[j];
2855  int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2856  if (lf < 0)
2857  {
2858  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2859  if (lf >= 0)
2860  {
2861  lf += NumOfFaces;
2862  }
2863  }
2864  face_nbr_el_to_face->Push(i, lf);
2865  }
2866  break;
2867  }
2868  case Element::HEXAHEDRON:
2869  {
2870  // find the face by the vertices with the smallest 3 numbers
2871  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
2872  for (int j = 0; j < 6; j++)
2873  {
2874  const int *fv = hex_t::FaceVert[j];
2875  int k = 0;
2876  int max = v[fv[0]];
2877 
2878  if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2879  if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2880  if (max < v[fv[3]]) { k = 3; }
2881 
2882  int v0 = -1, v1 = -1, v2 = -1;
2883  switch (k)
2884  {
2885  case 0:
2886  v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2887  break;
2888  case 1:
2889  v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2890  break;
2891  case 2:
2892  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2893  break;
2894  case 3:
2895  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2896  break;
2897  }
2898  int lf = faces_tbl->Index(v0, v1, v2);
2899  if (lf < 0)
2900  {
2901  lf = sfaces_tbl->Index(v0, v1, v2);
2902  if (lf >= 0)
2903  {
2904  lf += NumOfFaces;
2905  }
2906  }
2907  face_nbr_el_to_face->Push(i, lf);
2908  }
2909  break;
2910  }
2911  default:
2912  MFEM_ABORT("Unexpected type of Element.");
2913  }
2914  }
2916 
2917  delete sfaces_tbl;
2918  if (ret_ftbl)
2919  {
2920  return faces_tbl;
2921  }
2922  delete faces_tbl;
2923  return NULL;
2924 }
2925 
2926 int ParMesh::GetFaceNbrRank(int fn) const
2927 {
2928  if (Conforming())
2929  {
2930  int nbr_group = face_nbr_group[fn];
2931  const int *nbs = gtopo.GetGroup(nbr_group);
2932  int nbr_lproc = (nbs[0]) ? nbs[0] : nbs[1];
2933  int nbr_rank = gtopo.GetNeighborRank(nbr_lproc);
2934  return nbr_rank;
2935  }
2936  else
2937  {
2938  // NC: simplified handling of face neighbor ranks
2939  return face_nbr_group[fn];
2940  }
2941 }
2942 
2943 void
2945 {
2946  int n, j;
2947  int el_nbr = i - GetNE();
2948  if (face_nbr_el_to_face)
2949  {
2950  face_nbr_el_to_face->GetRow(el_nbr, fcs);
2951  }
2952  else
2953  {
2954  MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2955  "face_nbr_el_to_face not generated.");
2956  }
2957  if (el_nbr < face_nbr_el_ori.Size())
2958  {
2959  const int * row = face_nbr_el_ori.GetRow(el_nbr);
2960  n = fcs.Size();
2961  cor.SetSize(n);
2962  for (j=0; j<n; j++)
2963  {
2964  cor[j] = row[j];
2965  }
2966  }
2967  else
2968  {
2969  MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2970  "face_nbr_el_to_face not generated.");
2971  }
2972 }
2973 
2975 {
2976  const Array<int> *s2l_face;
2977  if (Dim == 1)
2978  {
2979  s2l_face = &svert_lvert;
2980  }
2981  else if (Dim == 2)
2982  {
2983  s2l_face = &sedge_ledge;
2984  }
2985  else
2986  {
2987  s2l_face = &sface_lface;
2988  }
2989 
2990  Table *face_elem = new Table;
2991 
2992  face_elem->MakeI(faces_info.Size());
2993 
2994  for (int i = 0; i < faces_info.Size(); i++)
2995  {
2996  if (faces_info[i].Elem2No >= 0)
2997  {
2998  face_elem->AddColumnsInRow(i, 2);
2999  }
3000  else
3001  {
3002  face_elem->AddAColumnInRow(i);
3003  }
3004  }
3005  for (int i = 0; i < s2l_face->Size(); i++)
3006  {
3007  face_elem->AddAColumnInRow((*s2l_face)[i]);
3008  }
3009 
3010  face_elem->MakeJ();
3011 
3012  for (int i = 0; i < faces_info.Size(); i++)
3013  {
3014  face_elem->AddConnection(i, faces_info[i].Elem1No);
3015  if (faces_info[i].Elem2No >= 0)
3016  {
3017  face_elem->AddConnection(i, faces_info[i].Elem2No);
3018  }
3019  }
3020  for (int i = 0; i < s2l_face->Size(); i++)
3021  {
3022  int lface = (*s2l_face)[i];
3023  int nbr_elem_idx = -1 - faces_info[lface].Elem2No;
3024  face_elem->AddConnection(lface, NumOfElements + nbr_elem_idx);
3025  }
3026 
3027  face_elem->ShiftUpI();
3028 
3029  return face_elem;
3030 }
3031 
3033  FaceElementTransformations* FETr, Element::Type face_type,
3034  Geometry::Type face_geom)
3035 {
3036  // calculate composition of FETr->Loc1 and FETr->Elem1
3037  DenseMatrix &face_pm = FETr->GetPointMat();
3038  FETr->Reset();
3039  if (Nodes == NULL)
3040  {
3041  FETr->Elem1->Transform(FETr->Loc1.Transf.GetPointMat(), face_pm);
3042  FETr->SetFE(GetTransformationFEforElementType(face_type));
3043  }
3044  else
3045  {
3046  const FiniteElement* face_el =
3047  Nodes->FESpace()->GetTraceElement(FETr->Elem1No, face_geom);
3048  MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
3049  "Mesh requires nodal Finite Element.");
3050 
3051 #if 0 // TODO: handle the case of non-interpolatory Nodes
3052  DenseMatrix I;
3053  face_el->Project(Transformation.GetFE(), FETr->Loc1.Transf, I);
3054  MultABt(Transformation.GetPointMat(), I, pm_face);
3055 #else
3056  IntegrationRule eir(face_el->GetDof());
3057  FETr->Loc1.Transform(face_el->GetNodes(), eir);
3058  Nodes->GetVectorValues(*FETr->Elem1, eir, face_pm);
3059 #endif
3060  FETr->SetFE(face_el);
3061  }
3062 }
3063 
3065  int FaceNo,
3066  int mask)
3067 {
3068  if (FaceNo < GetNumFaces())
3069  {
3070  return Mesh::GetFaceElementTransformations(FaceNo, mask);
3071  }
3072  else
3073  {
3074  const bool fill2 = mask & 10; // Elem2 and/or Loc2
3075  return GetSharedFaceTransformationsByLocalIndex(FaceNo, fill2);
3076  }
3077 }
3078 
3080 GetSharedFaceTransformations(int sf, bool fill2)
3081 {
3082  int FaceNo = GetSharedFace(sf);
3083 
3084  return GetSharedFaceTransformationsByLocalIndex(FaceNo, fill2);
3085 }
3086 
3089 {
3090  FaceInfo &face_info = faces_info[FaceNo];
3091  MFEM_VERIFY(face_info.Elem2Inf >= 0, "The face must be shared.");
3092 
3093  bool is_slave = Nonconforming() && IsSlaveFace(face_info);
3094  bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
3095 
3096  int mask = 0;
3098  FaceElemTr.Elem1 = NULL;
3099  FaceElemTr.Elem2 = NULL;
3100 
3101  NCFaceInfo* nc_info = NULL;
3102  if (is_slave) { nc_info = &nc_faces_info[face_info.NCFace]; }
3103 
3104  int local_face = is_ghost ? nc_info->MasterFace : FaceNo;
3105  Element::Type face_type = GetFaceElementType(local_face);
3106  Geometry::Type face_geom = GetFaceGeometry(local_face);
3107 
3108  // setup the transformation for the first element
3109  FaceElemTr.Elem1No = face_info.Elem1No;
3113 
3114  // setup the transformation for the second (neighbor) element
3115  int Elem2NbrNo;
3116  if (fill2)
3117  {
3118  Elem2NbrNo = -1 - face_info.Elem2No;
3119  // Store the "shifted index" for element 2 in FaceElemTr.Elem2No.
3120  // `Elem2NbrNo` is the index of the face neighbor (starting from 0),
3121  // and `FaceElemTr.Elem2No` will be offset by the number of (local)
3122  // elements in the mesh.
3123  FaceElemTr.Elem2No = NumOfElements + Elem2NbrNo;
3127  }
3128  else
3129  {
3130  FaceElemTr.Elem2No = -1;
3131  }
3132 
3133  // setup the face transformation if the face is not a ghost
3134  if (!is_ghost)
3135  {
3137  // NOTE: The above call overwrites FaceElemTr.Loc1
3139  }
3140  else
3141  {
3142  FaceElemTr.SetGeometryType(face_geom);
3143  }
3144 
3145  // setup Loc1 & Loc2
3146  int elem_type = GetElementType(face_info.Elem1No);
3147  GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc1.Transf,
3148  face_info.Elem1Inf);
3150 
3151  if (fill2)
3152  {
3153  elem_type = face_nbr_elements[Elem2NbrNo]->GetType();
3154  GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc2.Transf,
3155  face_info.Elem2Inf);
3157  }
3158 
3159  // adjust Loc1 or Loc2 of the master face if this is a slave face
3160  if (is_slave)
3161  {
3162  if (is_ghost || fill2)
3163  {
3164  // is_ghost -> modify side 1, otherwise -> modify side 2:
3165  ApplyLocalSlaveTransformation(FaceElemTr, face_info, is_ghost);
3166  }
3167  }
3168 
3169  // for ghost faces we need a special version of GetFaceTransformation
3170  if (is_ghost)
3171  {
3172  GetGhostFaceTransformation(&FaceElemTr, face_type, face_geom);
3174  }
3175 
3177 
3178  // This check can be useful for internal debugging, however it will fail on
3179  // periodic boundary faces, so we keep it disabled in general.
3180 #if 0
3181 #ifdef MFEM_DEBUG
3182  double dist = FaceElemTr.CheckConsistency();
3183  if (dist >= 1e-12)
3184  {
3185  mfem::out << "\nInternal error: face id = " << FaceNo
3186  << ", dist = " << dist << ", rank = " << MyRank << '\n';
3187  FaceElemTr.CheckConsistency(1); // print coordinates
3188  MFEM_ABORT("internal error");
3189  }
3190 #endif
3191 #endif
3192 
3193  return &FaceElemTr;
3194 }
3195 
3197 {
3198  if (Conforming())
3199  {
3200  switch (Dim)
3201  {
3202  case 1: return svert_lvert.Size();
3203  case 2: return sedge_ledge.Size();
3204  default: return sface_lface.Size();
3205  }
3206  }
3207  else
3208  {
3209  MFEM_ASSERT(Dim > 1, "");
3210  const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3211  return shared.conforming.Size() + shared.slaves.Size();
3212  }
3213 }
3214 
3215 int ParMesh::GetSharedFace(int sface) const
3216 {
3217  if (Conforming())
3218  {
3219  switch (Dim)
3220  {
3221  case 1: return svert_lvert[sface];
3222  case 2: return sedge_ledge[sface];
3223  default: return sface_lface[sface];
3224  }
3225  }
3226  else
3227  {
3228  MFEM_ASSERT(Dim > 1, "");
3229  const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3230  int csize = (int) shared.conforming.Size();
3231  return sface < csize
3232  ? shared.conforming[sface].index
3233  : shared.slaves[sface - csize].index;
3234  }
3235 }
3236 
3238 {
3239  const_cast<ParMesh*>(this)->ExchangeFaceNbrData();
3240  return Mesh::GetNFbyType(type);
3241 }
3242 
3243 // shift cyclically 3 integers a, b, c, so that the smallest of
3244 // order[a], order[b], order[c] is first
3245 static inline
3246 void Rotate3Indirect(int &a, int &b, int &c,
3247  const Array<std::int64_t> &order)
3248 {
3249  if (order[a] < order[b])
3250  {
3251  if (order[a] > order[c])
3252  {
3253  ShiftRight(a, b, c);
3254  }
3255  }
3256  else
3257  {
3258  if (order[b] < order[c])
3259  {
3260  ShiftRight(c, b, a);
3261  }
3262  else
3263  {
3264  ShiftRight(a, b, c);
3265  }
3266  }
3267 }
3268 
3270 {
3271  if (Dim != 3 || !(meshgen & 1))
3272  {
3273  return;
3274  }
3275 
3276  ResetLazyData();
3277 
3278  DSTable *old_v_to_v = NULL;
3279  Table *old_elem_vert = NULL;
3280 
3281  if (Nodes)
3282  {
3283  PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3284  }
3285 
3286  // create a GroupCommunicator over shared vertices
3287  GroupCommunicator svert_comm(gtopo);
3288  GetSharedVertexCommunicator(0, svert_comm);
3289 
3290  // communicate the local index of each shared vertex from the group master to
3291  // other ranks in the group
3292  Array<int> svert_master_rank(svert_lvert.Size());
3293  Array<int> svert_master_index(svert_lvert);
3294  {
3295  for (int i = 0; i < group_svert.Size(); i++)
3296  {
3297  int rank = gtopo.GetGroupMasterRank(i+1);
3298  for (int j = 0; j < group_svert.RowSize(i); j++)
3299  {
3300  svert_master_rank[group_svert.GetRow(i)[j]] = rank;
3301  }
3302  }
3303  svert_comm.Bcast(svert_master_index);
3304  }
3305 
3306  // the pairs (master rank, master local index) define a globally consistent
3307  // vertex ordering
3308  Array<std::int64_t> glob_vert_order(vertices.Size());
3309  {
3310  Array<int> lvert_svert(vertices.Size());
3311  lvert_svert = -1;
3312  for (int i = 0; i < svert_lvert.Size(); i++)
3313  {
3314  lvert_svert[svert_lvert[i]] = i;
3315  }
3316 
3317  for (int i = 0; i < vertices.Size(); i++)
3318  {
3319  int s = lvert_svert[i];
3320  if (s >= 0)
3321  {
3322  glob_vert_order[i] =
3323  (std::int64_t(svert_master_rank[s]) << 32) + svert_master_index[s];
3324  }
3325  else
3326  {
3327  glob_vert_order[i] = (std::int64_t(MyRank) << 32) + i;
3328  }
3329  }
3330  }
3331 
3332  // rotate tetrahedra so that vertex zero is the lowest (global) index vertex,
3333  // vertex 1 is the second lowest (global) index and vertices 2 and 3 preserve
3334  // positive orientation of the element
3335  for (int i = 0; i < NumOfElements; i++)
3336  {
3338  {
3339  int *v = elements[i]->GetVertices();
3340 
3341  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3342 
3343  if (glob_vert_order[v[0]] < glob_vert_order[v[3]])
3344  {
3345  Rotate3Indirect(v[1], v[2], v[3], glob_vert_order);
3346  }
3347  else
3348  {
3349  ShiftRight(v[0], v[1], v[3]);
3350  }
3351  }
3352  }
3353 
3354  // rotate also boundary triangles
3355  for (int i = 0; i < NumOfBdrElements; i++)
3356  {
3358  {
3359  int *v = boundary[i]->GetVertices();
3360 
3361  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3362  }
3363  }
3364 
3365  const bool check_consistency = true;
3366  if (check_consistency)
3367  {
3368  // create a GroupCommunicator on the shared triangles
3369  GroupCommunicator stria_comm(gtopo);
3370  GetSharedTriCommunicator(0, stria_comm);
3371 
3372  Array<int> stria_flag(shared_trias.Size());
3373  for (int i = 0; i < stria_flag.Size(); i++)
3374  {
3375  const int *v = shared_trias[i].v;
3376  if (glob_vert_order[v[0]] < glob_vert_order[v[1]])
3377  {
3378  stria_flag[i] = (glob_vert_order[v[0]] < glob_vert_order[v[2]]) ? 0 : 2;
3379  }
3380  else // v[1] < v[0]
3381  {
3382  stria_flag[i] = (glob_vert_order[v[1]] < glob_vert_order[v[2]]) ? 1 : 2;
3383  }
3384  }
3385 
3386  Array<int> stria_master_flag(stria_flag);
3387  stria_comm.Bcast(stria_master_flag);
3388  for (int i = 0; i < stria_flag.Size(); i++)
3389  {
3390  const int *v = shared_trias[i].v;
3391  MFEM_VERIFY(stria_flag[i] == stria_master_flag[i],
3392  "inconsistent vertex ordering found, shared triangle "
3393  << i << ": ("
3394  << v[0] << ", " << v[1] << ", " << v[2] << "), "
3395  << "local flag: " << stria_flag[i]
3396  << ", master flag: " << stria_master_flag[i]);
3397  }
3398  }
3399 
3400  // rotate shared triangle faces
3401  for (int i = 0; i < shared_trias.Size(); i++)
3402  {
3403  int *v = shared_trias[i].v;
3404 
3405  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3406  }
3407 
3408  // finalize
3409  if (!Nodes)
3410  {
3412  GenerateFaces();
3413  if (el_to_edge)
3414  {
3416  }
3417  }
3418  else
3419  {
3420  DoNodeReorder(old_v_to_v, old_elem_vert);
3421  delete old_elem_vert;
3422  delete old_v_to_v;
3423  }
3424 
3425  // the local edge and face numbering is changed therefore we need to
3426  // update sedge_ledge and sface_lface.
3427  FinalizeParTopo();
3428 }
3429 
3430 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
3431 {
3432  if (pncmesh)
3433  {
3434  MFEM_ABORT("Local and nonconforming refinements cannot be mixed.");
3435  }
3436 
3438 
3440 
3441  if (Dim == 3)
3442  {
3443  int uniform_refinement = 0;
3444  if (type < 0)
3445  {
3446  type = -type;
3447  uniform_refinement = 1;
3448  }
3449 
3450  // 1. Hash table of vertex to vertex connections corresponding to refined
3451  // edges.
3452  HashTable<Hashed2> v_to_v;
3453 
3454  // 2. Do the red refinement.
3455  switch (type)
3456  {
3457  case 1:
3458  for (int i = 0; i < marked_el.Size(); i++)
3459  {
3460  Bisection(marked_el[i], v_to_v);
3461  }
3462  break;
3463  case 2:
3464  for (int i = 0; i < marked_el.Size(); i++)
3465  {
3466  Bisection(marked_el[i], v_to_v);
3467 
3468  Bisection(NumOfElements - 1, v_to_v);
3469  Bisection(marked_el[i], v_to_v);
3470  }
3471  break;
3472  case 3:
3473  for (int i = 0; i < marked_el.Size(); i++)
3474  {
3475  Bisection(marked_el[i], v_to_v);
3476 
3477  int j = NumOfElements - 1;
3478  Bisection(j, v_to_v);
3479  Bisection(NumOfElements - 1, v_to_v);
3480  Bisection(j, v_to_v);
3481 
3482  Bisection(marked_el[i], v_to_v);
3483  Bisection(NumOfElements-1, v_to_v);
3484  Bisection(marked_el[i], v_to_v);
3485  }
3486  break;
3487  }
3488 
3489  // 3. Do the green refinement (to get conforming mesh).
3490  int need_refinement;
3491  int max_faces_in_group = 0;
3492  // face_splittings identify how the shared faces have been split
3493  Array<unsigned> *face_splittings = new Array<unsigned>[GetNGroups()-1];
3494  for (int i = 0; i < GetNGroups()-1; i++)
3495  {
3496  const int faces_in_group = GroupNTriangles(i+1);
3497  face_splittings[i].Reserve(faces_in_group);
3498  if (faces_in_group > max_faces_in_group)
3499  {
3500  max_faces_in_group = faces_in_group;
3501  }
3502  }
3503  int neighbor;
3504  Array<unsigned> iBuf(max_faces_in_group);
3505 
3506  MPI_Request *requests = new MPI_Request[GetNGroups()-1];
3507  MPI_Status status;
3508 
3509 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3510  int ref_loops_all = 0, ref_loops_par = 0;
3511 #endif
3512  do
3513  {
3514  need_refinement = 0;
3515  for (int i = 0; i < NumOfElements; i++)
3516  {
3517  if (elements[i]->NeedRefinement(v_to_v))
3518  {
3519  need_refinement = 1;
3520  Bisection(i, v_to_v);
3521  }
3522  }
3523 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3524  ref_loops_all++;
3525 #endif
3526 
3527  if (uniform_refinement)
3528  {
3529  continue;
3530  }
3531 
3532  // if the mesh is locally conforming start making it globally
3533  // conforming
3534  if (need_refinement == 0)
3535  {
3536 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3537  ref_loops_par++;
3538 #endif
3539  // MPI_Barrier(MyComm);
3540  const int tag = 293;
3541 
3542  // (a) send the type of interface splitting
3543  int req_count = 0;
3544  for (int i = 0; i < GetNGroups()-1; i++)
3545  {
3546  const int *group_faces = group_stria.GetRow(i);
3547  const int faces_in_group = group_stria.RowSize(i);
3548  // it is enough to communicate through the faces
3549  if (faces_in_group == 0) { continue; }
3550 
3551  face_splittings[i].SetSize(0);
3552  for (int j = 0; j < faces_in_group; j++)
3553  {
3554  GetFaceSplittings(shared_trias[group_faces[j]].v, v_to_v,
3555  face_splittings[i]);
3556  }
3557  const int *nbs = gtopo.GetGroup(i+1);
3558  neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3559  MPI_Isend(face_splittings[i], face_splittings[i].Size(),
3560  MPI_UNSIGNED, neighbor, tag, MyComm,
3561  &requests[req_count++]);
3562  }
3563 
3564  // (b) receive the type of interface splitting
3565  for (int i = 0; i < GetNGroups()-1; i++)
3566  {
3567  const int *group_faces = group_stria.GetRow(i);
3568  const int faces_in_group = group_stria.RowSize(i);
3569  if (faces_in_group == 0) { continue; }
3570 
3571  const int *nbs = gtopo.GetGroup(i+1);
3572  neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3573  MPI_Probe(neighbor, tag, MyComm, &status);
3574  int count;
3575  MPI_Get_count(&status, MPI_UNSIGNED, &count);
3576  iBuf.SetSize(count);
3577  MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag, MyComm,
3578  MPI_STATUS_IGNORE);
3579 
3580  for (int j = 0, pos = 0; j < faces_in_group; j++)
3581  {
3582  const int *v = shared_trias[group_faces[j]].v;
3583  need_refinement |= DecodeFaceSplittings(v_to_v, v, iBuf, pos);
3584  }
3585  }
3586 
3587  int nr = need_refinement;
3588  MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3589 
3590  MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
3591  }
3592  }
3593  while (need_refinement == 1);
3594 
3595 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3596  {
3597  int i = ref_loops_all;
3598  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3599  if (MyRank == 0)
3600  {
3601  mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3602  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3603  << '\n' << endl;
3604  }
3605  }
3606 #endif
3607 
3608  delete [] requests;
3609  iBuf.DeleteAll();
3610  delete [] face_splittings;
3611 
3612  // 4. Update the boundary elements.
3613  do
3614  {
3615  need_refinement = 0;
3616  for (int i = 0; i < NumOfBdrElements; i++)
3617  {
3618  if (boundary[i]->NeedRefinement(v_to_v))
3619  {
3620  need_refinement = 1;
3621  BdrBisection(i, v_to_v);
3622  }
3623  }
3624  }
3625  while (need_refinement == 1);
3626 
3627  if (NumOfBdrElements != boundary.Size())
3628  {
3629  mfem_error("ParMesh::LocalRefinement :"
3630  " (NumOfBdrElements != boundary.Size())");
3631  }
3632 
3633  ResetLazyData();
3634 
3635  const int old_nv = NumOfVertices;
3636  NumOfVertices = vertices.Size();
3637 
3638  RefineGroups(old_nv, v_to_v);
3639 
3640  // 5. Update the groups after refinement.
3641  if (el_to_face != NULL)
3642  {
3644  GenerateFaces();
3645  }
3646 
3647  // 6. Update element-to-edge relations.
3648  if (el_to_edge != NULL)
3649  {
3651  }
3652  } // 'if (Dim == 3)'
3653 
3654 
3655  if (Dim == 2)
3656  {
3657  int uniform_refinement = 0;
3658  if (type < 0)
3659  {
3660  // type = -type; // not used
3661  uniform_refinement = 1;
3662  }
3663 
3664  // 1. Get table of vertex to vertex connections.
3665  DSTable v_to_v(NumOfVertices);
3666  GetVertexToVertexTable(v_to_v);
3667 
3668  // 2. Get edge to element connections in arrays edge1 and edge2
3669  int nedges = v_to_v.NumberOfEntries();
3670  int *edge1 = new int[nedges];
3671  int *edge2 = new int[nedges];
3672  int *middle = new int[nedges];
3673 
3674  for (int i = 0; i < nedges; i++)
3675  {
3676  edge1[i] = edge2[i] = middle[i] = -1;
3677  }
3678 
3679  for (int i = 0; i < NumOfElements; i++)
3680  {
3681  int *v = elements[i]->GetVertices();
3682  for (int j = 0; j < 3; j++)
3683  {
3684  int ind = v_to_v(v[j], v[(j+1)%3]);
3685  (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
3686  }
3687  }
3688 
3689  // 3. Do the red refinement.
3690  for (int i = 0; i < marked_el.Size(); i++)
3691  {
3692  RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
3693  }
3694 
3695  // 4. Do the green refinement (to get conforming mesh).
3696  int need_refinement;
3697  int edges_in_group, max_edges_in_group = 0;
3698  // edge_splittings identify how the shared edges have been split
3699  int **edge_splittings = new int*[GetNGroups()-1];
3700  for (int i = 0; i < GetNGroups()-1; i++)
3701  {
3702  edges_in_group = GroupNEdges(i+1);
3703  edge_splittings[i] = new int[edges_in_group];
3704  if (edges_in_group > max_edges_in_group)
3705  {
3706  max_edges_in_group = edges_in_group;
3707  }
3708  }
3709  int neighbor, *iBuf = new int[max_edges_in_group];
3710 
3711  Array<int> group_edges;
3712 
3713  MPI_Request request;
3714  MPI_Status status;
3715  Vertex V;
3716  V(2) = 0.0;
3717 
3718 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3719  int ref_loops_all = 0, ref_loops_par = 0;
3720 #endif
3721  do
3722  {
3723  need_refinement = 0;
3724  for (int i = 0; i < nedges; i++)
3725  {
3726  if (middle[i] != -1 && edge1[i] != -1)
3727  {
3728  need_refinement = 1;
3729  GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
3730  }
3731  }
3732 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3733  ref_loops_all++;
3734 #endif
3735 
3736  if (uniform_refinement)
3737  {
3738  continue;
3739  }
3740 
3741  // if the mesh is locally conforming start making it globally
3742  // conforming
3743  if (need_refinement == 0)
3744  {
3745 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3746  ref_loops_par++;
3747 #endif
3748  // MPI_Barrier(MyComm);
3749 
3750  // (a) send the type of interface splitting
3751  for (int i = 0; i < GetNGroups()-1; i++)
3752  {
3753  group_sedge.GetRow(i, group_edges);
3754  edges_in_group = group_edges.Size();
3755  // it is enough to communicate through the edges
3756  if (edges_in_group != 0)
3757  {
3758  for (int j = 0; j < edges_in_group; j++)
3759  {
3760  edge_splittings[i][j] =
3761  GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
3762  middle);
3763  }
3764  const int *nbs = gtopo.GetGroup(i+1);
3765  if (nbs[0] == 0)
3766  {
3767  neighbor = gtopo.GetNeighborRank(nbs[1]);
3768  }
3769  else
3770  {
3771  neighbor = gtopo.GetNeighborRank(nbs[0]);
3772  }
3773  MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
3774  neighbor, 0, MyComm, &request);
3775  }
3776  }
3777 
3778  // (b) receive the type of interface splitting
3779  for (int i = 0; i < GetNGroups()-1; i++)
3780  {
3781  group_sedge.GetRow(i, group_edges);
3782  edges_in_group = group_edges.Size();
3783  if (edges_in_group != 0)
3784  {
3785  const int *nbs = gtopo.GetGroup(i+1);
3786  if (nbs[0] == 0)
3787  {
3788  neighbor = gtopo.GetNeighborRank(nbs[1]);
3789  }
3790  else
3791  {
3792  neighbor = gtopo.GetNeighborRank(nbs[0]);
3793  }
3794  MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
3795  MPI_ANY_TAG, MyComm, &status);
3796 
3797  for (int j = 0; j < edges_in_group; j++)
3798  {
3799  if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
3800  {
3801  int *v = shared_edges[group_edges[j]]->GetVertices();
3802  int ii = v_to_v(v[0], v[1]);
3803 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3804  if (middle[ii] != -1)
3805  {
3806  mfem_error("ParMesh::LocalRefinement (triangles) : "
3807  "Oops!");
3808  }
3809 #endif
3810  need_refinement = 1;
3811  middle[ii] = NumOfVertices++;
3812  for (int c = 0; c < 2; c++)
3813  {
3814  V(c) = 0.5 * (vertices[v[0]](c) + vertices[v[1]](c));
3815  }
3816  vertices.Append(V);
3817  }
3818  }
3819  }
3820  }
3821 
3822  int nr = need_refinement;
3823  MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3824  }
3825  }
3826  while (need_refinement == 1);
3827 
3828 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3829  {
3830  int i = ref_loops_all;
3831  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3832  if (MyRank == 0)
3833  {
3834  mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3835  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3836  << '\n' << endl;
3837  }
3838  }
3839 #endif
3840 
3841  for (int i = 0; i < GetNGroups()-1; i++)
3842  {
3843  delete [] edge_splittings[i];
3844  }
3845  delete [] edge_splittings;
3846 
3847  delete [] iBuf;
3848 
3849  // 5. Update the boundary elements.
3850  int v1[2], v2[2], bisect, temp;
3851  temp = NumOfBdrElements;
3852  for (int i = 0; i < temp; i++)
3853  {
3854  int *v = boundary[i]->GetVertices();
3855  bisect = v_to_v(v[0], v[1]);
3856  if (middle[bisect] != -1)
3857  {
3858  // the element was refined (needs updating)
3859  if (boundary[i]->GetType() == Element::SEGMENT)
3860  {
3861  v1[0] = v[0]; v1[1] = middle[bisect];
3862  v2[0] = middle[bisect]; v2[1] = v[1];
3863 
3864  boundary[i]->SetVertices(v1);
3865  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
3866  }
3867  else
3868  {
3869  mfem_error("Only bisection of segment is implemented for bdr"
3870  " elem.");
3871  }
3872  }
3873  }
3874  NumOfBdrElements = boundary.Size();
3875 
3876  ResetLazyData();
3877 
3878  // 5a. Update the groups after refinement.
3879  RefineGroups(v_to_v, middle);
3880 
3881  // 6. Free the allocated memory.
3882  delete [] edge1;
3883  delete [] edge2;
3884  delete [] middle;
3885 
3886  if (el_to_edge != NULL)
3887  {
3889  GenerateFaces();
3890  }
3891  } // 'if (Dim == 2)'
3892 
3893  if (Dim == 1) // --------------------------------------------------------
3894  {
3895  int cne = NumOfElements, cnv = NumOfVertices;
3896  NumOfVertices += marked_el.Size();
3897  NumOfElements += marked_el.Size();
3898  vertices.SetSize(NumOfVertices);
3899  elements.SetSize(NumOfElements);
3901 
3902  for (int j = 0; j < marked_el.Size(); j++)
3903  {
3904  int i = marked_el[j];
3905  Segment *c_seg = (Segment *)elements[i];
3906  int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
3907  int new_v = cnv + j, new_e = cne + j;
3908  AverageVertices(vert, 2, new_v);
3909  elements[new_e] = new Segment(new_v, vert[1], attr);
3910  vert[1] = new_v;
3911 
3914  }
3915 
3916  static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
3918  UseExternalData(seg_children, 1, 2, 3);
3919 
3920  GenerateFaces();
3921  } // end of 'if (Dim == 1)'
3922 
3924  sequence++;
3925 
3926  UpdateNodes();
3927 
3928 #ifdef MFEM_DEBUG
3929  CheckElementOrientation(false);
3931 #endif
3932 }
3933 
3935  int nc_limit)
3936 {
3937  if (NURBSext)
3938  {
3939  MFEM_ABORT("NURBS meshes are not supported. Please project the "
3940  "NURBS to Nodes first with SetCurvature().");
3941  }
3942 
3943  if (!pncmesh)
3944  {
3945  MFEM_ABORT("Can't convert conforming ParMesh to nonconforming ParMesh "
3946  "(you need to initialize the ParMesh from a nonconforming "
3947  "serial Mesh)");
3948  }
3949 
3951 
3952  // NOTE: no check of !refinements.Size(), in parallel we would have to reduce
3953 
3954  // do the refinements
3956  pncmesh->Refine(refinements);
3957 
3958  if (nc_limit > 0)
3959  {
3960  pncmesh->LimitNCLevel(nc_limit);
3961  }
3962 
3963  // create a second mesh containing the finest elements from 'pncmesh'
3964  ParMesh* pmesh2 = new ParMesh(*pncmesh);
3965  pncmesh->OnMeshUpdated(pmesh2);
3966 
3967  attributes.Copy(pmesh2->attributes);
3969 
3970  // now swap the meshes, the second mesh will become the old coarse mesh
3971  // and this mesh will be the new fine mesh
3972  Mesh::Swap(*pmesh2, false);
3973 
3974  delete pmesh2; // NOTE: old face neighbors destroyed here
3975 
3977 
3979 
3981  sequence++;
3982 
3983  UpdateNodes();
3984 }
3985 
3987  double threshold, int nc_limit, int op)
3988 {
3989  MFEM_VERIFY(pncmesh, "Only supported for non-conforming meshes.");
3990  MFEM_VERIFY(!NURBSext, "Derefinement of NURBS meshes is not supported. "
3991  "Project the NURBS to Nodes first.");
3992 
3993  const Table &dt = pncmesh->GetDerefinementTable();
3994 
3995  pncmesh->SynchronizeDerefinementData(elem_error, dt);
3996 
3997  Array<int> level_ok;
3998  if (nc_limit > 0)
3999  {
4000  pncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
4001  }
4002 
4003  Array<int> derefs;
4004  for (int i = 0; i < dt.Size(); i++)
4005  {
4006  if (nc_limit > 0 && !level_ok[i]) { continue; }
4007 
4008  double error =
4009  AggregateError(elem_error, dt.GetRow(i), dt.RowSize(i), op);
4010 
4011  if (error < threshold) { derefs.Append(i); }
4012  }
4013 
4014  long long glob_size = ReduceInt(derefs.Size());
4015  if (!glob_size) { return false; }
4016 
4017  // Destroy face-neighbor data only when actually de-refining.
4019 
4020  pncmesh->Derefine(derefs);
4021 
4022  ParMesh* mesh2 = new ParMesh(*pncmesh);
4023  pncmesh->OnMeshUpdated(mesh2);
4024 
4025  attributes.Copy(mesh2->attributes);
4027 
4028  Mesh::Swap(*mesh2, false);
4029  delete mesh2;
4030 
4032 
4034 
4036  sequence++;
4037 
4038  UpdateNodes();
4039 
4040  return true;
4041 }
4042 
4043 
4045 {
4046  RebalanceImpl(NULL); // default SFC-based partition
4047 }
4048 
4049 void ParMesh::Rebalance(const Array<int> &partition)
4050 {
4051  RebalanceImpl(&partition);
4052 }
4053 
4054 void ParMesh::RebalanceImpl(const Array<int> *partition)
4055 {
4056  if (Conforming())
4057  {
4058  MFEM_ABORT("Load balancing is currently not supported for conforming"
4059  " meshes.");
4060  }
4061 
4062  if (Nodes)
4063  {
4064  // check that Nodes use a parallel FE space, so we can call UpdateNodes()
4065  MFEM_VERIFY(dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace())
4066  != NULL, "internal error");
4067  }
4068 
4070 
4071  pncmesh->Rebalance(partition);
4072 
4073  ParMesh* pmesh2 = new ParMesh(*pncmesh);
4074  pncmesh->OnMeshUpdated(pmesh2);
4075 
4076  attributes.Copy(pmesh2->attributes);
4078 
4079  Mesh::Swap(*pmesh2, false);
4080  delete pmesh2;
4081 
4083 
4085 
4087  sequence++;
4088 
4089  UpdateNodes();
4090 }
4091 
4092 void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
4093 {
4094  // Refine groups after LocalRefinement in 2D (triangle meshes)
4095 
4096  MFEM_ASSERT(Dim == 2 && meshgen == 1, "internal error");
4097 
4098  Array<int> group_verts, group_edges;
4099 
4100  // To update the groups after a refinement, we observe that:
4101  // - every (new and old) vertex, edge and face belongs to exactly one group
4102  // - the refinement does not create new groups
4103  // - a new vertex appears only as the middle of a refined edge
4104  // - a face can be refined 2, 3 or 4 times producing new edges and faces
4105 
4106  int *I_group_svert, *J_group_svert;
4107  int *I_group_sedge, *J_group_sedge;
4108 
4109  I_group_svert = Memory<int>(GetNGroups()+1);
4110  I_group_sedge = Memory<int>(GetNGroups()+1);
4111 
4112  I_group_svert[0] = I_group_svert[1] = 0;
4113  I_group_sedge[0] = I_group_sedge[1] = 0;
4114 
4115  // overestimate the size of the J arrays
4116  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4118  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4119 
4120  for (int group = 0; group < GetNGroups()-1; group++)
4121  {
4122  // Get the group shared objects
4123  group_svert.GetRow(group, group_verts);
4124  group_sedge.GetRow(group, group_edges);
4125 
4126  // Check which edges have been refined
4127  for (int i = 0; i < group_sedge.RowSize(group); i++)
4128  {
4129  int *v = shared_edges[group_edges[i]]->GetVertices();
4130  const int ind = middle[v_to_v(v[0], v[1])];
4131  if (ind != -1)
4132  {
4133  // add a vertex
4134  group_verts.Append(svert_lvert.Append(ind)-1);
4135  // update the edges
4136  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4137  shared_edges.Append(new Segment(v[1], ind, attr));
4138  group_edges.Append(sedge_ledge.Append(-1)-1);
4139  v[1] = ind;
4140  }
4141  }
4142 
4143  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4144  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4145 
4146  int *J;
4147  J = J_group_svert+I_group_svert[group];
4148  for (int i = 0; i < group_verts.Size(); i++)
4149  {
4150  J[i] = group_verts[i];
4151  }
4152  J = J_group_sedge+I_group_sedge[group];
4153  for (int i = 0; i < group_edges.Size(); i++)
4154  {
4155  J[i] = group_edges[i];
4156  }
4157  }
4158 
4159  FinalizeParTopo();
4160 
4161  group_svert.SetIJ(I_group_svert, J_group_svert);
4162  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4163 }
4164 
4165 void ParMesh::RefineGroups(int old_nv, const HashTable<Hashed2> &v_to_v)
4166 {
4167  // Refine groups after LocalRefinement in 3D (tetrahedral meshes)
4168 
4169  MFEM_ASSERT(Dim == 3 && meshgen == 1, "internal error");
4170 
4171  Array<int> group_verts, group_edges, group_trias;
4172 
4173  // To update the groups after a refinement, we observe that:
4174  // - every (new and old) vertex, edge and face belongs to exactly one group
4175  // - the refinement does not create new groups
4176  // - a new vertex appears only as the middle of a refined edge
4177  // - a face can be refined multiple times producing new edges and faces
4178 
4179  Array<Segment *> sedge_stack;
4180  Array<Vert3> sface_stack;
4181 
4182  Array<int> I_group_svert, J_group_svert;
4183  Array<int> I_group_sedge, J_group_sedge;
4184  Array<int> I_group_stria, J_group_stria;
4185 
4186  I_group_svert.SetSize(GetNGroups());
4187  I_group_sedge.SetSize(GetNGroups());
4188  I_group_stria.SetSize(GetNGroups());
4189 
4190  I_group_svert[0] = 0;
4191  I_group_sedge[0] = 0;
4192  I_group_stria[0] = 0;
4193 
4194  for (int group = 0; group < GetNGroups()-1; group++)
4195  {
4196  // Get the group shared objects
4197  group_svert.GetRow(group, group_verts);
4198  group_sedge.GetRow(group, group_edges);
4199  group_stria.GetRow(group, group_trias);
4200 
4201  // Check which edges have been refined
4202  for (int i = 0; i < group_sedge.RowSize(group); i++)
4203  {
4204  int *v = shared_edges[group_edges[i]]->GetVertices();
4205  int ind = v_to_v.FindId(v[0], v[1]);
4206  if (ind == -1) { continue; }
4207 
4208  // This shared edge is refined: walk the whole refinement tree
4209  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4210  do
4211  {
4212  ind += old_nv;
4213  // Add new shared vertex
4214  group_verts.Append(svert_lvert.Append(ind)-1);
4215  // Put the right sub-edge on top of the stack
4216  sedge_stack.Append(new Segment(ind, v[1], attr));
4217  // The left sub-edge replaces the original edge
4218  v[1] = ind;
4219  ind = v_to_v.FindId(v[0], ind);
4220  }
4221  while (ind != -1);
4222  // Process all edges in the edge stack
4223  do
4224  {
4225  Segment *se = sedge_stack.Last();
4226  v = se->GetVertices();
4227  ind = v_to_v.FindId(v[0], v[1]);
4228  if (ind == -1)
4229  {
4230  // The edge 'se' is not refined
4231  sedge_stack.DeleteLast();
4232  // Add new shared edge
4233  shared_edges.Append(se);
4234  group_edges.Append(sedge_ledge.Append(-1)-1);
4235  }
4236  else
4237  {
4238  // The edge 'se' is refined
4239  ind += old_nv;
4240  // Add new shared vertex
4241  group_verts.Append(svert_lvert.Append(ind)-1);
4242  // Put the left sub-edge on top of the stack
4243  sedge_stack.Append(new Segment(v[0], ind, attr));
4244  // The right sub-edge replaces the original edge
4245  v[0] = ind;
4246  }
4247  }
4248  while (sedge_stack.Size() > 0);
4249  }
4250 
4251  // Check which triangles have been refined
4252  for (int i = 0; i < group_stria.RowSize(group); i++)
4253  {
4254  int *v = shared_trias[group_trias[i]].v;
4255  int ind = v_to_v.FindId(v[0], v[1]);
4256  if (ind == -1) { continue; }
4257 
4258  // This shared face is refined: walk the whole refinement tree
4259  const int edge_attr = 1;
4260  do
4261  {
4262  ind += old_nv;
4263  // Add the refinement edge to the edge stack
4264  sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4265  // Put the right sub-triangle on top of the face stack
4266  sface_stack.Append(Vert3(v[1], v[2], ind));
4267  // The left sub-triangle replaces the original one
4268  v[1] = v[0]; v[0] = v[2]; v[2] = ind;
4269  ind = v_to_v.FindId(v[0], v[1]);
4270  }
4271  while (ind != -1);
4272  // Process all faces (triangles) in the face stack
4273  do
4274  {
4275  Vert3 &st = sface_stack.Last();
4276  v = st.v;
4277  ind = v_to_v.FindId(v[0], v[1]);
4278  if (ind == -1)
4279  {
4280  // The triangle 'st' is not refined
4281  // Add new shared face
4282  shared_trias.Append(st);
4283  group_trias.Append(sface_lface.Append(-1)-1);
4284  sface_stack.DeleteLast();
4285  }
4286  else
4287  {
4288  // The triangle 'st' is refined
4289  ind += old_nv;
4290  // Add the refinement edge to the edge stack
4291  sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4292  // Put the left sub-triangle on top of the face stack
4293  sface_stack.Append(Vert3(v[2], v[0], ind));
4294  // Note that the above Append() may invalidate 'v'
4295  v = sface_stack[sface_stack.Size()-2].v;
4296  // The right sub-triangle replaces the original one
4297  v[0] = v[1]; v[1] = v[2]; v[2] = ind;
4298  }
4299  }
4300  while (sface_stack.Size() > 0);
4301  // Process all edges in the edge stack (same code as above)
4302  do
4303  {
4304  Segment *se = sedge_stack.Last();
4305  v = se->GetVertices();
4306  ind = v_to_v.FindId(v[0], v[1]);
4307  if (ind == -1)
4308  {
4309  // The edge 'se' is not refined
4310  sedge_stack.DeleteLast();
4311  // Add new shared edge
4312  shared_edges.Append(se);
4313  group_edges.Append(sedge_ledge.Append(-1)-1);
4314  }
4315  else
4316  {
4317  // The edge 'se' is refined
4318  ind += old_nv;
4319  // Add new shared vertex
4320  group_verts.Append(svert_lvert.Append(ind)-1);
4321  // Put the left sub-edge on top of the stack
4322  sedge_stack.Append(new Segment(v[0], ind, edge_attr));
4323  // The right sub-edge replaces the original edge
4324  v[0] = ind;
4325  }
4326  }
4327  while (sedge_stack.Size() > 0);
4328  }
4329 
4330  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4331  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4332  I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4333 
4334  J_group_svert.Append(group_verts);
4335  J_group_sedge.Append(group_edges);
4336  J_group_stria.Append(group_trias);
4337  }
4338 
4339  FinalizeParTopo();
4340 
4341  group_svert.SetIJ(I_group_svert, J_group_svert);
4342  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4343  group_stria.SetIJ(I_group_stria, J_group_stria);
4344  I_group_svert.LoseData(); J_group_svert.LoseData();
4345  I_group_sedge.LoseData(); J_group_sedge.LoseData();
4346  I_group_stria.LoseData(); J_group_stria.LoseData();
4347 }
4348 
4350 {
4351  Array<int> sverts, sedges;
4352 
4353  int *I_group_svert, *J_group_svert;
4354  int *I_group_sedge, *J_group_sedge;
4355 
4356  I_group_svert = Memory<int>(GetNGroups());
4357  I_group_sedge = Memory<int>(GetNGroups());
4358 
4359  I_group_svert[0] = 0;
4360  I_group_sedge[0] = 0;
4361 
4362  // compute the size of the J arrays
4363  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4365  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4366 
4367  for (int group = 0; group < GetNGroups()-1; group++)
4368  {
4369  // Get the group shared objects
4370  group_svert.GetRow(group, sverts);
4371  group_sedge.GetRow(group, sedges);
4372 
4373  // Process all the edges
4374  for (int i = 0; i < group_sedge.RowSize(group); i++)
4375  {
4376  int *v = shared_edges[sedges[i]]->GetVertices();
4377  const int ind = old_nv + sedge_ledge[sedges[i]];
4378  // add a vertex
4379  sverts.Append(svert_lvert.Append(ind)-1);
4380  // update the edges
4381  const int attr = shared_edges[sedges[i]]->GetAttribute();
4382  shared_edges.Append(new Segment(v[1], ind, attr));
4383  sedges.Append(sedge_ledge.Append(-1)-1);
4384  v[1] = ind;
4385  }
4386 
4387  I_group_svert[group+1] = I_group_svert[group] + sverts.Size();
4388  I_group_sedge[group+1] = I_group_sedge[group] + sedges.Size();
4389 
4390  sverts.CopyTo(J_group_svert + I_group_svert[group]);
4391  sedges.CopyTo(J_group_sedge + I_group_sedge[group]);
4392  }
4393 
4394  FinalizeParTopo();
4395 
4396  group_svert.SetIJ(I_group_svert, J_group_svert);
4397  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4398 }
4399 
4400 void ParMesh::UniformRefineGroups3D(int old_nv, int old_nedges,
4401  const DSTable &old_v_to_v,
4402  const STable3D &old_faces,
4403  Array<int> *f2qf)
4404 {
4405  // f2qf can be NULL if all faces are quads or there are no quad faces
4406 
4407  Array<int> group_verts, group_edges, group_trias, group_quads;
4408 
4409  int *I_group_svert, *J_group_svert;
4410  int *I_group_sedge, *J_group_sedge;
4411  int *I_group_stria, *J_group_stria;
4412  int *I_group_squad, *J_group_squad;
4413 
4414  I_group_svert = Memory<int>(GetNGroups());
4415  I_group_sedge = Memory<int>(GetNGroups());
4416  I_group_stria = Memory<int>(GetNGroups());
4417  I_group_squad = Memory<int>(GetNGroups());
4418 
4419  I_group_svert[0] = 0;
4420  I_group_sedge[0] = 0;
4421  I_group_stria[0] = 0;
4422  I_group_squad[0] = 0;
4423 
4424  // compute the size of the J arrays
4425  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4428  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections() +
4431  J_group_stria = Memory<int>(4*group_stria.Size_of_connections());
4432  J_group_squad = Memory<int>(4*group_squad.Size_of_connections());
4433 
4434  const int oface = old_nv + old_nedges;
4435 
4436  for (int group = 0; group < GetNGroups()-1; group++)
4437  {
4438  // Get the group shared objects
4439  group_svert.GetRow(group, group_verts);
4440  group_sedge.GetRow(group, group_edges);
4441  group_stria.GetRow(group, group_trias);
4442  group_squad.GetRow(group, group_quads);
4443 
4444  // Process the edges that have been refined
4445  for (int i = 0; i < group_sedge.RowSize(group); i++)
4446  {
4447  int *v = shared_edges[group_edges[i]]->GetVertices();
4448  const int ind = old_nv + old_v_to_v(v[0], v[1]);
4449  // add a vertex
4450  group_verts.Append(svert_lvert.Append(ind)-1);
4451  // update the edges
4452  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4453  shared_edges.Append(new Segment(v[1], ind, attr));
4454  group_edges.Append(sedge_ledge.Append(-1)-1);
4455  v[1] = ind; // v[0] remains the same
4456  }
4457 
4458  // Process the triangles that have been refined
4459  for (int i = 0; i < group_stria.RowSize(group); i++)
4460  {
4461  int m[3];
4462  const int stria = group_trias[i];
4463  int *v = shared_trias[stria].v;
4464  // add the refinement edges
4465  m[0] = old_nv + old_v_to_v(v[0], v[1]);
4466  m[1] = old_nv + old_v_to_v(v[1], v[2]);
4467  m[2] = old_nv + old_v_to_v(v[2], v[0]);
4468  const int edge_attr = 1;
4469  shared_edges.Append(new Segment(m[0], m[1], edge_attr));
4470  group_edges.Append(sedge_ledge.Append(-1)-1);
4471  shared_edges.Append(new Segment(m[1], m[2], edge_attr));
4472  group_edges.Append(sedge_ledge.Append(-1)-1);
4473  shared_edges.Append(new Segment(m[0], m[2], edge_attr));
4474  group_edges.Append(sedge_ledge.Append(-1)-1);
4475  // update faces
4476  const int nst = shared_trias.Size();
4477  shared_trias.SetSize(nst+3);
4478  // The above SetSize() may invalidate 'v'
4479  v = shared_trias[stria].v;
4480  shared_trias[nst+0].Set(m[1],m[2],m[0]);
4481  shared_trias[nst+1].Set(m[0],v[1],m[1]);
4482  shared_trias[nst+2].Set(m[2],m[1],v[2]);
4483  v[1] = m[0]; v[2] = m[2]; // v[0] remains the same
4484  group_trias.Append(nst+0);
4485  group_trias.Append(nst+1);
4486  group_trias.Append(nst+2);
4487  // sface_lface is set later
4488  }
4489 
4490  // Process the quads that have been refined
4491  for (int i = 0; i < group_squad.RowSize(group); i++)
4492  {
4493  int m[5];
4494  const int squad = group_quads[i];
4495  int *v = shared_quads[squad].v;
4496  const int olf = old_faces(v[0], v[1], v[2], v[3]);
4497  // f2qf can be NULL if all faces are quads
4498  m[0] = oface + (f2qf ? (*f2qf)[olf] : olf);
4499  // add a vertex
4500  group_verts.Append(svert_lvert.Append(m[0])-1);
4501  // add the refinement edges
4502  m[1] = old_nv + old_v_to_v(v[0], v[1]);
4503  m[2] = old_nv + old_v_to_v(v[1], v[2]);
4504  m[3] = old_nv + old_v_to_v(v[2], v[3]);
4505  m[4] = old_nv + old_v_to_v(v[3], v[0]);
4506  const int edge_attr = 1;
4507  shared_edges.Append(new Segment(m[1], m[0], edge_attr));
4508  group_edges.Append(sedge_ledge.Append(-1)-1);
4509  shared_edges.Append(new Segment(m[2], m[0], edge_attr));
4510  group_edges.Append(sedge_ledge.Append(-1)-1);
4511  shared_edges.Append(new Segment(m[3], m[0], edge_attr));
4512  group_edges.Append(sedge_ledge.Append(-1)-1);
4513  shared_edges.Append(new Segment(m[4], m[0], edge_attr));
4514  group_edges.Append(sedge_ledge.Append(-1)-1);
4515  // update faces
4516  const int nsq = shared_quads.Size();
4517  shared_quads.SetSize(nsq+3);
4518  // The above SetSize() may invalidate 'v'
4519  v = shared_quads[squad].v;
4520  shared_quads[nsq+0].Set(m[1],v[1],m[2],m[0]);
4521  shared_quads[nsq+1].Set(m[0],m[2],v[2],m[3]);
4522  shared_quads[nsq+2].Set(m[4],m[0],m[3],v[3]);
4523  v[1] = m[1]; v[2] = m[0]; v[3] = m[4]; // v[0] remains the same
4524  group_quads.Append(nsq+0);
4525  group_quads.Append(nsq+1);
4526  group_quads.Append(nsq+2);
4527  // sface_lface is set later
4528  }
4529 
4530  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4531  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4532  I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4533  I_group_squad[group+1] = I_group_squad[group] + group_quads.Size();
4534 
4535  group_verts.CopyTo(J_group_svert + I_group_svert[group]);
4536  group_edges.CopyTo(J_group_sedge + I_group_sedge[group]);
4537  group_trias.CopyTo(J_group_stria + I_group_stria[group]);
4538  group_quads.CopyTo(J_group_squad + I_group_squad[group]);
4539  }
4540 
4541  FinalizeParTopo();
4542 
4543  group_svert.SetIJ(I_group_svert, J_group_svert);
4544  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4545  group_stria.SetIJ(I_group_stria, J_group_stria);
4546  group_squad.SetIJ(I_group_squad, J_group_squad);
4547 }
4548 
4550 {
4552 
4553  const int old_nv = NumOfVertices;
4554 
4555  // call Mesh::UniformRefinement2D so that it won't update the nodes
4556  {
4557  const bool update_nodes = false;
4558  Mesh::UniformRefinement2D_base(update_nodes);
4559  }
4560 
4561  // update the groups
4562  UniformRefineGroups2D(old_nv);
4563 
4564  UpdateNodes();
4565 
4566 #ifdef MFEM_DEBUG
4567  // If there are no Nodes, the orientation is checked in the call to
4568  // UniformRefinement2D_base() above.
4569  if (Nodes) { CheckElementOrientation(false); }
4570 #endif
4571 }
4572 
4574 {
4576 
4577  const int old_nv = NumOfVertices;
4578  const int old_nedges = NumOfEdges;
4579 
4580  DSTable v_to_v(NumOfVertices);
4581  GetVertexToVertexTable(v_to_v);
4582  STable3D *faces_tbl = GetFacesTable();
4583 
4584  // call Mesh::UniformRefinement3D_base so that it won't update the nodes
4585  Array<int> f2qf;
4586  {
4587  const bool update_nodes = false;
4588  UniformRefinement3D_base(&f2qf, &v_to_v, update_nodes);
4589  // Note: for meshes that have triangular faces, v_to_v is modified by the
4590  // above call to return different edge indices - this is used when
4591  // updating the groups. This is needed by ReorientTetMesh().
4592  }
4593 
4594  // update the groups
4595  UniformRefineGroups3D(old_nv, old_nedges, v_to_v, *faces_tbl,
4596  f2qf.Size() ? &f2qf : NULL);
4597  delete faces_tbl;
4598 
4599  UpdateNodes();
4600 }
4601 
4603 {
4604  if (MyRank == 0)
4605  {
4606  mfem::out << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
4607  }
4608 }
4609 
4610 void ParMesh::PrintXG(std::ostream &os) const
4611 {
4612  MFEM_ASSERT(Dim == spaceDim, "2D manifolds not supported");
4613  if (Dim == 3 && meshgen == 1)
4614  {
4615  int i, j, nv;
4616  const int *ind;
4617 
4618  os << "NETGEN_Neutral_Format\n";
4619  // print the vertices
4620  os << NumOfVertices << '\n';
4621  for (i = 0; i < NumOfVertices; i++)
4622  {
4623  for (j = 0; j < Dim; j++)
4624  {
4625  os << " " << vertices[i](j);
4626  }
4627  os << '\n';
4628  }
4629 
4630  // print the elements
4631  os << NumOfElements << '\n';
4632  for (i = 0; i < NumOfElements; i++)
4633  {
4634  nv = elements[i]->GetNVertices();
4635  ind = elements[i]->GetVertices();
4636  os << elements[i]->GetAttribute();
4637  for (j = 0; j < nv; j++)
4638  {
4639  os << " " << ind[j]+1;
4640  }
4641  os << '\n';
4642  }
4643 
4644  // print the boundary + shared faces information
4645  os << NumOfBdrElements + sface_lface.Size() << '\n';
4646  // boundary
4647  for (i = 0; i < NumOfBdrElements; i++)
4648  {
4649  nv = boundary[i]->GetNVertices();
4650  ind = boundary[i]->GetVertices();
4651  os << boundary[i]->GetAttribute();
4652  for (j = 0; j < nv; j++)
4653  {
4654  os << " " << ind[j]+1;
4655  }
4656  os << '\n';
4657  }
4658  // shared faces
4659  const int sf_attr =
4660  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4661  for (i = 0; i < shared_trias.Size(); i++)
4662  {
4663  ind = shared_trias[i].v;
4664  os << sf_attr;
4665  for (j = 0; j < 3; j++)
4666  {
4667  os << ' ' << ind[j]+1;
4668  }
4669  os << '\n';
4670  }
4671  // There are no quad shared faces
4672  }
4673 
4674  if (Dim == 3 && meshgen == 2)
4675  {
4676  int i, j, nv;
4677  const int *ind;
4678 
4679  os << "TrueGrid\n"
4680  << "1 " << NumOfVertices << " " << NumOfElements
4681  << " 0 0 0 0 0 0 0\n"
4682  << "0 0 0 1 0 0 0 0 0 0 0\n"
4683  << "0 0 " << NumOfBdrElements+sface_lface.Size()
4684  << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
4685  << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
4686  << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
4687 
4688  // print the vertices
4689  for (i = 0; i < NumOfVertices; i++)
4690  {
4691  os << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
4692  << " " << vertices[i](2) << " 0.0\n";
4693  }
4694 
4695  // print the elements
4696  for (i = 0; i < NumOfElements; i++)
4697  {
4698  nv = elements[i]->GetNVertices();
4699  ind = elements[i]->GetVertices();
4700  os << i+1 << " " << elements[i]->GetAttribute();
4701  for (j = 0; j < nv; j++)
4702  {
4703  os << " " << ind[j]+1;
4704  }
4705  os << '\n';
4706  }
4707 
4708  // print the boundary information
4709  for (i = 0; i < NumOfBdrElements; i++)
4710  {
4711  nv = boundary[i]->GetNVertices();
4712  ind = boundary[i]->GetVertices();
4713  os << boundary[i]->GetAttribute();
4714  for (j = 0; j < nv; j++)
4715  {
4716  os << " " << ind[j]+1;
4717  }
4718  os << " 1.0 1.0 1.0 1.0\n";
4719  }
4720 
4721  // print the shared faces information
4722  const int sf_attr =
4723  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4724  // There are no shared triangle faces
4725  for (i = 0; i < shared_quads.Size(); i++)
4726  {
4727  ind = shared_quads[i].v;
4728  os << sf_attr;
4729  for (j = 0; j < 4; j++)
4730  {
4731  os << ' ' << ind[j]+1;
4732  }
4733  os << " 1.0 1.0 1.0 1.0\n";
4734  }
4735  }
4736 
4737  if (Dim == 2)
4738  {
4739  int i, j, attr;
4740  Array<int> v;
4741 
4742  os << "areamesh2\n\n";
4743 
4744  // print the boundary + shared edges information
4745  os << NumOfBdrElements + shared_edges.Size() << '\n';
4746  // boundary
4747  for (i = 0; i < NumOfBdrElements; i++)
4748  {
4749  attr = boundary[i]->GetAttribute();
4750  boundary[i]->GetVertices(v);
4751  os << attr << " ";
4752  for (j = 0; j < v.Size(); j++)
4753  {
4754  os << v[j] + 1 << " ";
4755  }
4756  os << '\n';
4757  }
4758  // shared edges
4759  for (i = 0; i < shared_edges.Size(); i++)
4760  {
4761  attr = shared_edges[i]->GetAttribute();
4762  shared_edges[i]->GetVertices(v);
4763  os << attr << " ";
4764  for (j = 0; j < v.Size(); j++)
4765  {
4766  os << v[j] + 1 << " ";
4767  }
4768  os << '\n';
4769  }
4770 
4771  // print the elements
4772  os << NumOfElements << '\n';
4773  for (i = 0; i < NumOfElements; i++)
4774  {
4775  attr = elements[i]->GetAttribute();
4776  elements[i]->GetVertices(v);
4777 
4778  os << attr << " ";
4779  if ((j = GetElementType(i)) == Element::TRIANGLE)
4780  {
4781  os << 3 << " ";
4782  }
4783  else if (j == Element::QUADRILATERAL)
4784  {
4785  os << 4 << " ";
4786  }
4787  else if (j == Element::SEGMENT)
4788  {
4789  os << 2 << " ";
4790  }
4791  for (j = 0; j < v.Size(); j++)
4792  {
4793  os << v[j] + 1 << " ";
4794  }
4795  os << '\n';
4796  }
4797 
4798  // print the vertices
4799  os << NumOfVertices << '\n';
4800  for (i = 0; i < NumOfVertices; i++)
4801  {
4802  for (j = 0; j < Dim; j++)
4803  {
4804  os << vertices[i](j) << " ";
4805  }
4806  os << '\n';
4807  }
4808  }
4809  os.flush();
4810 }
4811 
4813 {
4814  // In 2D, this is a workaround for a CPU boundary rendering artifact. We need
4815  // to skip a shared master edge if one of its slaves has the same rank.
4816 
4817  const NCMesh::NCList &list = pncmesh->GetEdgeList();
4818  for (int i = master.slaves_begin; i < master.slaves_end; i++)
4819  {
4820  if (!pncmesh->IsGhost(1, list.slaves[i].index)) { return true; }
4821  }
4822  return false;
4823 }
4824 
4825 void ParMesh::Print(std::ostream &os) const
4826 {
4827  int shared_bdr_attr;
4828  Array<int> nc_shared_faces;
4829 
4830  if (NURBSext)
4831  {
4832  Printer(os); // does not print shared boundary
4833  return;
4834  }
4835 
4836  const Array<int>* s2l_face;
4837  if (!pncmesh)
4838  {
4839  s2l_face = ((Dim == 1) ? &svert_lvert :
4840  ((Dim == 2) ? &sedge_ledge : &sface_lface));
4841  }
4842  else
4843  {
4844  s2l_face = &nc_shared_faces;
4845  if (Dim >= 2)
4846  {
4847  // get a list of all shared non-ghost faces
4848  const NCMesh::NCList& sfaces =
4849  (Dim == 3) ? pncmesh->GetSharedFaces() : pncmesh->GetSharedEdges();
4850  const int nfaces = GetNumFaces();
4851  for (int i = 0; i < sfaces.conforming.Size(); i++)
4852  {
4853  int index = sfaces.conforming[i].index;
4854  if (index < nfaces) { nc_shared_faces.Append(index); }
4855  }
4856  for (int i = 0; i < sfaces.masters.Size(); i++)
4857  {
4858  if (Dim == 2 && WantSkipSharedMaster(sfaces.masters[i])) { continue; }
4859  int index = sfaces.masters[i].index;
4860  if (index < nfaces) { nc_shared_faces.Append(index); }
4861  }
4862  for (int i = 0; i < sfaces.slaves.Size(); i++)
4863  {
4864  int index = sfaces.slaves[i].index;
4865  if (index < nfaces) { nc_shared_faces.Append(index); }
4866  }
4867  }
4868  }
4869 
4870  os << "MFEM mesh v1.0\n";
4871 
4872  // optional
4873  os <<
4874  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4875  "# POINT = 0\n"
4876  "# SEGMENT = 1\n"
4877  "# TRIANGLE = 2\n"
4878  "# SQUARE = 3\n"
4879  "# TETRAHEDRON = 4\n"
4880  "# CUBE = 5\n"
4881  "# PRISM = 6\n"
4882  "#\n";
4883 
4884  os << "\ndimension\n" << Dim
4885  << "\n\nelements\n" << NumOfElements << '\n';
4886  for (int i = 0; i < NumOfElements; i++)
4887  {
4888  PrintElement(elements[i], os);
4889  }
4890 
4891  int num_bdr_elems = NumOfBdrElements;
4892  if (print_shared && Dim > 1)
4893  {
4894  num_bdr_elems += s2l_face->Size();
4895  }
4896  os << "\nboundary\n" << num_bdr_elems << '\n';
4897  for (int i = 0; i < NumOfBdrElements; i++)
4898  {
4899  PrintElement(boundary[i], os);
4900  }
4901 
4902  if (print_shared && Dim > 1)
4903  {
4904  if (bdr_attributes.Size())
4905  {
4906  shared_bdr_attr = bdr_attributes.Max() + MyRank + 1;
4907  }
4908  else
4909  {
4910  shared_bdr_attr = MyRank + 1;
4911  }
4912  for (int i = 0; i < s2l_face->Size(); i++)
4913  {
4914  // Modify the attributes of the faces (not used otherwise?)
4915  faces[(*s2l_face)[i]]->SetAttribute(shared_bdr_attr);
4916  PrintElement(faces[(*s2l_face)[i]], os);
4917  }
4918  }
4919  os << "\nvertices\n" << NumOfVertices << '\n';
4920  if (Nodes == NULL)
4921  {
4922  os << spaceDim << '\n';
4923  for (int i = 0; i < NumOfVertices; i++)
4924  {
4925  os << vertices[i](0);
4926  for (int j = 1; j < spaceDim; j++)
4927  {
4928  os << ' ' << vertices[i](j);
4929  }
4930  os << '\n';
4931  }
4932  os.flush();
4933  }
4934  else
4935  {
4936  os << "\nnodes\n";
4937  Nodes->Save(os);
4938  }
4939 }
4940 
4941 void ParMesh::Save(const std::string &fname, int precision) const
4942 {
4943  ostringstream fname_with_suffix;
4944  fname_with_suffix << fname << "." << setfill('0') << setw(6) << MyRank;
4945  ofstream ofs(fname_with_suffix.str().c_str());
4946  ofs.precision(precision);
4947  Print(ofs);
4948 }
4949 
4950 #ifdef MFEM_USE_ADIOS2
4952 {
4953  Mesh::Print(os);
4954 }
4955 #endif
4956 
4957 static void dump_element(const Element* elem, Array<int> &data)
4958 {
4959  data.Append(elem->GetGeometryType());
4960 
4961  int nv = elem->GetNVertices();
4962  const int *v = elem->GetVertices();
4963  for (int i = 0; i < nv; i++)
4964  {
4965  data.Append(v[i]);
4966  }
4967 }
4968 
4969 void ParMesh::PrintAsOne(std::ostream &os) const
4970 {
4971  int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
4972  const int *v;
4973  MPI_Status status;
4974  Array<double> vert;
4975  Array<int> ints;
4976 
4977  if (MyRank == 0)
4978  {
4979  os << "MFEM mesh v1.0\n";
4980 
4981  // optional
4982  os <<
4983  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4984  "# POINT = 0\n"
4985  "# SEGMENT = 1\n"
4986  "# TRIANGLE = 2\n"
4987  "# SQUARE = 3\n"
4988  "# TETRAHEDRON = 4\n"
4989  "# CUBE = 5\n"
4990  "# PRISM = 6\n"
4991  "#\n";
4992 
4993  os << "\ndimension\n" << Dim;
4994  }
4995 
4996  nv = NumOfElements;
4997  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
4998  if (MyRank == 0)
4999  {
5000  os << "\n\nelements\n" << ne << '\n';
5001  for (i = 0; i < NumOfElements; i++)
5002  {
5003  // processor number + 1 as attribute and geometry type
5004  os << 1 << ' ' << elements[i]->GetGeometryType();
5005  // vertices
5006  nv = elements[i]->GetNVertices();
5007  v = elements[i]->GetVertices();
5008  for (j = 0; j < nv; j++)
5009  {
5010  os << ' ' << v[j];
5011  }
5012  os << '\n';
5013  }
5014  vc = NumOfVertices;
5015  for (p = 1; p < NRanks; p++)
5016  {
5017  MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
5018  ints.SetSize(ne);
5019  if (ne)
5020  {
5021  MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
5022  }
5023  for (i = 0; i < ne; )
5024  {
5025  // processor number + 1 as attribute and geometry type
5026  os << p+1 << ' ' << ints[i];
5027  // vertices
5028  k = Geometries.GetVertices(ints[i++])->GetNPoints();
5029  for (j = 0; j < k; j++)
5030  {
5031  os << ' ' << vc + ints[i++];
5032  }
5033  os << '\n';
5034  }
5035  vc += nv;
5036  }
5037  }
5038  else
5039  {
5040  // for each element send its geometry type and its vertices
5041  ne = 0;
5042  for (i = 0; i < NumOfElements; i++)
5043  {
5044  ne += 1 + elements[i]->GetNVertices();
5045  }
5046  nv = NumOfVertices;
5047  MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
5048 
5049  ints.Reserve(ne);
5050  ints.SetSize(0);
5051  for (i = 0; i < NumOfElements; i++)
5052  {
5053  dump_element(elements[i], ints);
5054  }
5055  MFEM_ASSERT(ints.Size() == ne, "");
5056  if (ne)
5057  {
5058  MPI_Send(&ints[0], ne, MPI_INT, 0, 445, MyComm);
5059  }
5060  }
5061 
5062  // boundary + shared boundary
5063  ne = NumOfBdrElements;
5064  if (!pncmesh)
5065  {
5066  ne += GetNSharedFaces();
5067  }
5068  else if (Dim > 1)
5069  {
5070  const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5071  ne += list.conforming.Size() + list.masters.Size() + list.slaves.Size();
5072  // In addition to the number returned by GetNSharedFaces(), include the
5073  // the master shared faces as well.
5074  }
5075  ints.Reserve(ne * (1 + 2*(Dim-1))); // just an upper bound
5076  ints.SetSize(0);
5077 
5078  // for each boundary and shared boundary element send its geometry type
5079  // and its vertices
5080  ne = 0;
5081  for (i = j = 0; i < NumOfBdrElements; i++)
5082  {
5083  dump_element(boundary[i], ints); ne++;
5084  }
5085  if (!pncmesh)
5086  {
5087  switch (Dim)
5088  {
5089  case 1:
5090  for (i = 0; i < svert_lvert.Size(); i++)
5091  {
5092  ints.Append(Geometry::POINT);
5093  ints.Append(svert_lvert[i]);
5094  ne++;
5095  }
5096  break;
5097 
5098  case 2:
5099  for (i = 0; i < shared_edges.Size(); i++)
5100  {
5101  dump_element(shared_edges[i], ints); ne++;
5102  }
5103  break;
5104 
5105  case 3:
5106  for (i = 0; i < shared_trias.Size(); i++)
5107  {
5108  ints.Append(Geometry::TRIANGLE);
5109  ints.Append(shared_trias[i].v, 3);
5110  ne++;
5111  }
5112  for (i = 0; i < shared_quads.Size(); i++)
5113  {
5114  ints.Append(Geometry::SQUARE);
5115  ints.Append(shared_quads[i].v, 4);
5116  ne++;
5117  }
5118  break;
5119 
5120  default:
5121  MFEM_ABORT("invalid dimension: " << Dim);
5122  }
5123  }
5124  else if (Dim > 1)
5125  {
5126  const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5127  const int nfaces = GetNumFaces();
5128  for (i = 0; i < list.conforming.Size(); i++)
5129  {
5130  int index = list.conforming[i].index;
5131  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5132  }
5133  for (i = 0; i < list.masters.Size(); i++)
5134  {
5135  int index = list.masters[i].index;
5136  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5137  }
5138  for (i = 0; i < list.slaves.Size(); i++)
5139  {
5140  int index = list.slaves[i].index;
5141  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5142  }
5143  }
5144 
5145  MPI_Reduce(&ne, &k, 1, MPI_INT, MPI_SUM, 0, MyComm);
5146  if (MyRank == 0)
5147  {
5148  os << "\nboundary\n" << k << '\n';
5149  vc = 0;
5150  for (p = 0; p < NRanks; p++)
5151  {
5152  if (p)
5153  {
5154  MPI_Recv(nv_ne, 2, MPI_INT, p, 446, MyComm, &status);
5155  ints.SetSize(ne);
5156  if (ne)
5157  {
5158  MPI_Recv(ints.GetData(), ne, MPI_INT, p, 447, MyComm, &status);
5159  }
5160  }
5161  else
5162  {
5163  ne = ints.Size();
5164  nv = NumOfVertices;
5165  }
5166  for (i = 0; i < ne; )
5167  {
5168  // processor number + 1 as bdr. attr. and bdr. geometry type
5169  os << p+1 << ' ' << ints[i];
5170  k = Geometries.NumVerts[ints[i++]];
5171  // vertices
5172  for (j = 0; j < k; j++)
5173  {
5174  os << ' ' << vc + ints[i++];
5175  }
5176  os << '\n';
5177  }
5178  vc += nv;
5179  }
5180  }
5181  else
5182  {
5183  nv = NumOfVertices;
5184  ne = ints.Size();
5185  MPI_Send(nv_ne, 2, MPI_INT, 0, 446, MyComm);
5186  if (ne)
5187  {
5188  MPI_Send(ints.GetData(), ne, MPI_INT, 0, 447, MyComm);
5189  }
5190  }
5191 
5192  // vertices / nodes
5193  MPI_Reduce(&NumOfVertices, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5194  if (MyRank == 0)
5195  {
5196  os << "\nvertices\n" << nv << '\n';
5197  }
5198  if (Nodes == NULL)
5199  {
5200  if (MyRank == 0)
5201  {
5202  os << spaceDim << '\n';
5203  for (i = 0; i < NumOfVertices; i++)
5204  {
5205  os << vertices[i](0);
5206  for (j = 1; j < spaceDim; j++)
5207  {
5208  os << ' ' << vertices[i](j);
5209  }
5210  os << '\n';
5211  }
5212  for (p = 1; p < NRanks; p++)
5213  {
5214  MPI_Recv(&nv, 1, MPI_INT, p, 448, MyComm, &status);
5215  vert.SetSize(nv*spaceDim);
5216  if (nv)
5217  {
5218  MPI_Recv(&vert[0], nv*spaceDim, MPI_DOUBLE, p, 449, MyComm, &status);
5219  }
5220  for (i = 0; i < nv; i++)
5221  {
5222  os << vert[i*spaceDim];
5223  for (j = 1; j < spaceDim; j++)
5224  {
5225  os << ' ' << vert[i*spaceDim+j];
5226  }
5227  os << '\n';
5228  }
5229  }
5230  os.flush();
5231  }
5232  else
5233  {
5234  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 448, MyComm);
5236  for (i = 0; i < NumOfVertices; i++)
5237  {
5238  for (j = 0; j < spaceDim; j++)
5239  {
5240  vert[i*spaceDim+j] = vertices[i](j);
5241  }
5242  }
5243  if (NumOfVertices)
5244  {
5245  MPI_Send(&vert[0], NumOfVertices*spaceDim, MPI_DOUBLE, 0, 449, MyComm);
5246  }
5247  }
5248  }
5249  else
5250  {
5251  if (MyRank == 0)
5252  {
5253  os << "\nnodes\n";
5254  }
5255  ParGridFunction *pnodes = dynamic_cast<ParGridFunction *>(Nodes);
5256  if (pnodes)
5257  {
5258  pnodes->SaveAsOne(os);
5259  }
5260  else
5261  {
5262  ParFiniteElementSpace *pfes =
5263  dynamic_cast<ParFiniteElementSpace *>(Nodes->FESpace());
5264  if (pfes)
5265  {
5266  // create a wrapper ParGridFunction
5267  ParGridFunction ParNodes(pfes, Nodes);
5268  ParNodes.SaveAsOne(os);
5269  }
5270  else
5271  {
5272  mfem_error("ParMesh::PrintAsOne : Nodes have no parallel info!");
5273  }
5274  }
5275  }
5276 }
5277 
5278 void ParMesh::PrintAsSerial(std::ostream &os) const
5279 {
5280  int save_rank = 0;
5281  Mesh serialmesh = GetSerialMesh(save_rank);
5282  if (MyRank == save_rank)
5283  {
5284  serialmesh.Printer(os);
5285  }
5286  MPI_Barrier(MyComm);
5287 }
5288 
5289 Mesh ParMesh::GetSerialMesh(int save_rank) const
5290 {
5291  if (pncmesh || NURBSext)
5292  {
5293  MFEM_ABORT("Nonconforming meshes and NURBS meshes are not yet supported.");
5294  }
5295 
5296  // Define linear H1 space for vertex numbering
5297  H1_FECollection fec_linear(1, Dim);
5298  ParMesh *pm = const_cast<ParMesh *>(this);
5299  ParFiniteElementSpace pfespace_linear(pm, &fec_linear);
5300 
5301  long long ne_glob_l = GetGlobalNE(); // needs to be called by all ranks
5302  MFEM_VERIFY(int(ne_glob_l) == ne_glob_l,
5303  "overflow in the number of elements!");
5304  int ne_glob = (save_rank == MyRank) ? int(ne_glob_l) : 0;
5305 
5306  long long nvertices = pfespace_linear.GetTrueVSize();
5307  long long nvertices_glob_l = 0;
5308  MPI_Reduce(&nvertices, &nvertices_glob_l, 1, MPI_LONG_LONG, MPI_SUM,
5309  save_rank, MyComm);
5310  int nvertices_glob = int(nvertices_glob_l);
5311  MFEM_VERIFY(nvertices_glob == nvertices_glob_l,
5312  "overflow in the number of vertices!");
5313 
5314  long long nbe = NumOfBdrElements;
5315  long long nbe_glob_l = 0;
5316  MPI_Reduce(&nbe, &nbe_glob_l, 1, MPI_LONG_LONG, MPI_SUM, save_rank, MyComm);
5317  int nbe_glob = int(nbe_glob_l);
5318  MFEM_VERIFY(nbe_glob == nbe_glob_l,
5319  "overflow in the number of boundary elements!");
5320 
5321  // On ranks other than save_rank, the *_glob variables are 0, so the serial
5322  // mesh is empty.
5323  Mesh serialmesh(Dim, nvertices_glob, ne_glob, nbe_glob, spaceDim);
5324 
5325  int n_send_recv;
5326  MPI_Status status;
5327  Array<double> vert;
5328  Array<int> ints, dofs;
5329 
5330  // First set the connectivity of serial mesh using the True Dofs from
5331  // the linear H1 space.
5332  if (MyRank == save_rank)
5333  {
5334  for (int e = 0; e < NumOfElements; e++)
5335  {
5336  const int attr = elements[e]->GetAttribute();
5337  const int geom_type = elements[e]->GetGeometryType();
5338  pfespace_linear.GetElementDofs(e, dofs);
5339  for (int j = 0; j < dofs.Size(); j++)
5340  {
5341  dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
5342  }
5343  Element *elem = serialmesh.NewElement(geom_type);
5344  elem->SetAttribute(attr);
5345  elem->SetVertices(dofs);
5346  serialmesh.AddElement(elem);
5347  }
5348 
5349  for (int p = 0; p < NRanks; p++)
5350  {
5351  if (p == save_rank) { continue; }
5352  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 444, MyComm, &status);
5353  ints.SetSize(n_send_recv);
5354  if (n_send_recv)
5355  {
5356  MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 445, MyComm, &status);
5357  }
5358  for (int i = 0; i < n_send_recv; )
5359  {
5360  int attr = ints[i++];
5361  int geom_type = ints[i++];
5362  Element *elem = serialmesh.NewElement(geom_type);
5363  elem->SetAttribute(attr);
5364  elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
5365  serialmesh.AddElement(elem);
5366  }
5367  }
5368  }
5369  else
5370  {
5371  n_send_recv = 0;
5372  for (int e = 0; e < NumOfElements; e++)
5373  {
5374  n_send_recv += 2 + elements[e]->GetNVertices();
5375  }
5376  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 444, MyComm);
5377  ints.Reserve(n_send_recv);
5378  ints.SetSize(0);
5379  for (int e = 0; e < NumOfElements; e++)
5380  {
5381  const int attr = elements[e]->GetAttribute();
5382  const int geom_type = elements[e]->GetGeometryType();;
5383  ints.Append(attr);
5384  ints.Append(geom_type);
5385  pfespace_linear.GetElementDofs(e, dofs);
5386  for (int j = 0; j < dofs.Size(); j++)
5387  {
5388  ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
5389  }
5390  }
5391  if (n_send_recv)
5392  {
5393  MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 445, MyComm);
5394  }
5395  }
5396 
5397  // Write out boundary elements
5398  if (MyRank == save_rank)
5399  {
5400  for (int e = 0; e < NumOfBdrElements; e++)
5401  {
5402  const int attr = boundary[e]->GetAttribute();
5403  const int geom_type = boundary[e]->GetGeometryType();
5404  pfespace_linear.GetBdrElementDofs(e, dofs);
5405  for (int j = 0; j < dofs.Size(); j++)
5406  {
5407  dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
5408  }
5409  Element *elem = serialmesh.NewElement(geom_type);
5410  elem->SetAttribute(attr);
5411  elem->SetVertices(dofs);
5412  serialmesh.AddBdrElement(elem);
5413  }
5414 
5415  for (int p = 0; p < NRanks; p++)
5416  {
5417  if (p == save_rank) { continue; }
5418  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 446, MyComm, &status);
5419  ints.SetSize(n_send_recv);
5420  if (n_send_recv)
5421  {
5422  MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 447, MyComm, &status);
5423  }
5424  for (int i = 0; i < n_send_recv; )
5425  {
5426  int attr = ints[i++];
5427  int geom_type = ints[i++];
5428  Element *elem = serialmesh.NewElement(geom_type);
5429  elem->SetAttribute(attr);
5430  elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
5431  serialmesh.AddBdrElement(elem);
5432  }
5433  }
5434  } // MyRank == save_rank
5435  else
5436  {
5437  n_send_recv = 0;
5438  for (int e = 0; e < NumOfBdrElements; e++)
5439  {
5440  n_send_recv += 2 + GetBdrElement(e)->GetNVertices();
5441  }
5442  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 446, MyComm);
5443  ints.Reserve(n_send_recv);
5444  ints.SetSize(0);
5445  for (int e = 0; e < NumOfBdrElements; e++)
5446  {
5447  const int attr = boundary[e]->GetAttribute();
5448  const int geom_type = boundary[e]->GetGeometryType();
5449  ints.Append(attr);
5450  ints.Append(geom_type);
5451  pfespace_linear.GetBdrElementDofs(e, dofs);
5452  for (int j = 0; j < dofs.Size(); j++)
5453  {
5454  ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
5455  }
5456  }
5457  if (n_send_recv)
5458  {
5459  MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 447, MyComm);
5460  }
5461  } // MyRank != save_rank
5462 
5463  if (MyRank == save_rank)
5464  {
5465  for (int v = 0; v < nvertices_glob; v++)
5466  {
5467  serialmesh.AddVertex(0.0); // all other coordinates are 0 by default
5468  }
5469  serialmesh.FinalizeTopology();
5470  }
5471 
5472  // From each processor, we send element-wise vertex/dof locations and
5473  // overwrite the vertex/dof locations of the serial mesh.
5474  if (MyRank == save_rank && Nodes)
5475  {
5476  FiniteElementSpace *fespace_serial = NULL;
5477  // Duplicate the FE collection to make sure the serial mesh is completely
5478  // independent of the parallel mesh:
5479  auto fec_serial = FiniteElementCollection::New(
5480  GetNodalFESpace()->FEColl()->Name());
5481  fespace_serial = new FiniteElementSpace(&serialmesh,
5482  fec_serial,
5483  spaceDim,
5484  GetNodalFESpace()->GetOrdering());
5485  serialmesh.SetNodalFESpace(fespace_serial);
5486  serialmesh.GetNodes()->MakeOwner(fec_serial);
5487  // The serial mesh owns its Nodes and they, in turn, own fec_serial and
5488  // fespace_serial.
5489  }
5490 
5491  int elem_count = 0; // To keep track of element count in serial mesh
5492  if (MyRank == save_rank)
5493  {
5494  Vector nodeloc;
5495  Array<int> ints_serial;
5496  for (int e = 0; e < NumOfElements; e++)
5497  {
5498  if (Nodes)
5499  {
5500  Nodes->GetElementDofValues(e, nodeloc);
5501  serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
5502  serialmesh.GetNodes()->SetSubVector(dofs, nodeloc);
5503  }
5504  else
5505  {
5506  GetElementVertices(e, ints);
5507  serialmesh.GetElementVertices(elem_count++, ints_serial);
5508  for (int i = 0; i < ints.Size(); i++)
5509  {
5510  const double *vdata = GetVertex(ints[i]);
5511  double *vdata_serial = serialmesh.GetVertex(ints_serial[i]);
5512  for (int d = 0; d < spaceDim; d++)
5513  {
5514  vdata_serial[d] = vdata[d];
5515  }
5516  }
5517  }
5518  }
5519 
5520  for (int p = 0; p < NRanks; p++)
5521  {
5522  if (p == save_rank) { continue; }
5523  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, MyComm, &status);
5524  vert.SetSize(n_send_recv);
5525  if (n_send_recv)
5526  {
5527  MPI_Recv(&vert[0], n_send_recv, MPI_DOUBLE, p, 449, MyComm, &status);
5528  }
5529  for (int i = 0; i < n_send_recv; )
5530  {
5531  if (Nodes)
5532  {
5533  serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
5534  serialmesh.GetNodes()->SetSubVector(dofs, &vert[i]);
5535  i += dofs.Size();
5536  }
5537  else
5538  {
5539  serialmesh.GetElementVertices(elem_count++, ints_serial);
5540  for (int j = 0; j < ints_serial.Size(); j++)
5541  {
5542  double *vdata_serial = serialmesh.GetVertex(ints_serial[j]);
5543  for (int d = 0; d < spaceDim; d++)
5544  {
5545  vdata_serial[d] = vert[i++];
5546  }
5547  }
5548  }
5549  }
5550  }
5551  } // MyRank == save_rank
5552  else
5553  {
5554  n_send_recv = 0;
5555  Vector nodeloc;
5556  for (int e = 0; e < NumOfElements; e++)
5557  {
5558  if (Nodes)
5559  {
5560  const FiniteElement *fe = Nodes->FESpace()->GetFE(e);
5561  n_send_recv += spaceDim*fe->GetDof();
5562  }
5563  else
5564  {
5565  n_send_recv += elements[e]->GetNVertices()*spaceDim;
5566  }
5567  }
5568  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, MyComm);
5569  vert.Reserve(n_send_recv);
5570  vert.SetSize(0);
5571  for (int e = 0; e < NumOfElements; e++)
5572  {
5573  if (Nodes)
5574  {
5575  Nodes->GetElementDofValues(e, nodeloc);
5576  for (int j = 0; j < nodeloc.Size(); j++)
5577  {
5578  vert.Append(nodeloc(j));
5579  }
5580  }
5581  else
5582  {
5583  GetElementVertices(e, ints);
5584  for (int i = 0; i < ints.Size(); i++)
5585  {
5586  const double *vdata = GetVertex(ints[i]);
5587  for (int d = 0; d < spaceDim; d++)
5588  {
5589  vert.Append(vdata[d]);
5590  }
5591  }
5592  }
5593  }
5594  if (n_send_recv)
5595  {
5596  MPI_Send(&vert[0], n_send_recv, MPI_DOUBLE, save_rank, 449, MyComm);
5597  }
5598  }
5599 
5600  MPI_Barrier(MyComm);
5601  return serialmesh;
5602 }
5603 
5604 void ParMesh::SaveAsOne(const std::string &fname, int precision) const
5605 {
5606  ofstream ofs;
5607  if (MyRank == 0)
5608  {
5609  ofs.open(fname);
5610  ofs.precision(precision);
5611  }
5612  PrintAsOne(ofs);
5613 }
5614 
5615 void ParMesh::PrintAsOneXG(std::ostream &os)
5616 {
5617  MFEM_ASSERT(Dim == spaceDim, "2D Manifolds not supported.");
5618  if (Dim == 3 && meshgen == 1)
5619  {
5620  int i, j, k, nv, ne, p;
5621  const int *ind, *v;
5622  MPI_Status status;
5623  Array<double> vert;
5624  Array<int> ints;
5625 
5626  if (MyRank == 0)
5627  {
5628  os << "NETGEN_Neutral_Format\n";
5629  // print the vertices
5630  ne = NumOfVertices;
5631  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5632  os << nv << '\n';
5633  for (i = 0; i < NumOfVertices; i++)
5634  {
5635  for (j = 0; j < Dim; j++)
5636  {
5637  os << " " << vertices[i](j);
5638  }
5639  os << '\n';
5640  }
5641  for (p = 1; p < NRanks; p++)
5642  {
5643  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5644  vert.SetSize(Dim*nv);
5645  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
5646  for (i = 0; i < nv; i++)
5647  {
5648  for (j = 0; j < Dim; j++)
5649  {
5650  os << " " << vert[Dim*i+j];
5651  }
5652  os << '\n';
5653  }
5654  }
5655 
5656  // print the elements
5657  nv = NumOfElements;
5658  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5659  os << ne << '\n';
5660  for (i = 0; i < NumOfElements; i++)
5661  {
5662  nv = elements[i]->GetNVertices();
5663  ind = elements[i]->GetVertices();
5664  os << 1;
5665  for (j = 0; j < nv; j++)
5666  {
5667  os << " " << ind[j]+1;
5668  }
5669  os << '\n';
5670  }
5671  k = NumOfVertices;
5672  for (p = 1; p < NRanks; p++)
5673  {
5674  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5675  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5676  ints.SetSize(4*ne);
5677  MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
5678  for (i = 0; i < ne; i++)
5679  {
5680  os << p+1;
5681  for (j = 0; j < 4; j++)
5682  {
5683  os << " " << k+ints[i*4+j]+1;
5684  }
5685  os << '\n';
5686  }
5687  k += nv;
5688  }
5689  // print the boundary + shared faces information
5691  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5692  os << ne << '\n';
5693  // boundary
5694  for (i = 0; i < NumOfBdrElements; i++)
5695  {
5696  nv = boundary[i]->GetNVertices();
5697  ind = boundary[i]->GetVertices();
5698  os << 1;
5699  for (j = 0; j < nv; j++)
5700  {
5701  os << " " << ind[j]+1;
5702  }
5703  os << '\n';
5704  }
5705  // shared faces
5706  const int sf_attr =
5707  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
5708  for (i = 0; i < shared_trias.Size(); i++)
5709  {
5710  ind = shared_trias[i].v;
5711  os << sf_attr;
5712  for (j = 0; j < 3; j++)
5713  {
5714  os << ' ' << ind[j]+1;
5715  }
5716  os << '\n';
5717  }
5718  // There are no quad shared faces
5719  k = NumOfVertices;
5720  for (p = 1; p < NRanks; p++)
5721  {
5722  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5723  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5724  ints.SetSize(3*ne);
5725  MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
5726  for (i = 0; i < ne; i++)
5727  {
5728  os << p+1;
5729  for (j = 0; j < 3; j++)
5730  {
5731  os << ' ' << k+ints[i*3+j]+1;
5732  }
5733  os << '\n';
5734  }
5735  k += nv;
5736  }
5737  }
5738  else
5739  {
5740  ne = NumOfVertices;
5741  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5742  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5743  vert.SetSize(Dim*NumOfVertices);
5744  for (i = 0; i < NumOfVertices; i++)
5745  for (j = 0; j < Dim; j++)
5746  {
5747  vert[Dim*i+j] = vertices[i](j);
5748  }
5749  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
5750  0, 445, MyComm);
5751  // elements
5752  ne = NumOfElements;
5753  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5754  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5755  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
5756  ints.SetSize(NumOfElements*4);
5757  for (i = 0; i < NumOfElements; i++)
5758  {
5759  v = elements[i]->GetVertices();
5760  for (j = 0; j < 4; j++)
5761  {
5762  ints[4*i+j] = v[j];
5763  }
5764  }
5765  MPI_Send(&ints[0], 4*NumOfElements, MPI_INT, 0, 447, MyComm);
5766  // boundary + shared faces
5768  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5769  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5771  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
5772  ints.SetSize(3*ne);
5773  for (i = 0; i < NumOfBdrElements; i++)
5774  {
5775  v = boundary[i]->GetVertices();
5776  for (j = 0; j < 3; j++)
5777  {
5778  ints[3*i+j] = v[j];
5779  }
5780  }
5781  for ( ; i < ne; i++)
5782  {
5783  v = shared_trias[i-NumOfBdrElements].v; // tet mesh
5784  for (j = 0; j < 3; j++)
5785  {
5786  ints[3*i+j] = v[j];
5787  }
5788  }
5789  MPI_Send(&ints[0], 3*ne, MPI_INT, 0, 447, MyComm);
5790  }
5791  }
5792 
5793  if (Dim == 3 && meshgen == 2)
5794  {
5795  int i, j, k, nv, ne, p;
5796  const int *ind, *v;
5797  MPI_Status status;
5798  Array<double> vert;
5799  Array<int> ints;
5800 
5801  int TG_nv, TG_ne, TG_nbe;
5802 
5803  if (MyRank == 0)
5804  {
5805  MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5806  MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5808  MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
5809 
5810  os << "TrueGrid\n"
5811  << "1 " << TG_nv << " " << TG_ne << " 0 0 0 0 0 0 0\n"