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