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