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