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