MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pmesh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "mesh_headers.hpp"
17 #include "../fem/fem.hpp"
18 #include "../general/sets.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../general/text.hpp"
21 #include "../general/globals.hpp"
22 
23 #include <iostream>
24 #include <fstream>
25 
26 using namespace std;
27 
28 namespace mfem
29 {
30 
31 ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
32  : Mesh(pmesh, false),
33  group_svert(pmesh.group_svert),
34  group_sedge(pmesh.group_sedge),
35  group_stria(pmesh.group_stria),
36  group_squad(pmesh.group_squad),
37  face_nbr_el_to_face(NULL),
38  glob_elem_offset(-1),
39  glob_offset_sequence(-1),
40  gtopo(pmesh.gtopo)
41 {
42  MyComm = pmesh.MyComm;
43  NRanks = pmesh.NRanks;
44  MyRank = pmesh.MyRank;
45 
46  // Duplicate the shared_edges
47  shared_edges.SetSize(pmesh.shared_edges.Size());
48  for (int i = 0; i < shared_edges.Size(); i++)
49  {
50  shared_edges[i] = pmesh.shared_edges[i]->Duplicate(this);
51  }
52 
53  shared_trias = pmesh.shared_trias;
54  shared_quads = pmesh.shared_quads;
55 
56  // Copy the shared-to-local index Arrays
59  sface_lface = pmesh.sface_lface;
60 
61  // Do not copy face-neighbor data (can be generated if needed)
62  have_face_nbr_data = false;
63 
64  // If pmesh has a ParNURBSExtension, it was copied by the Mesh copy ctor, so
65  // there is no need to do anything here.
66 
67  // Copy ParNCMesh, if present
68  if (pmesh.pncmesh)
69  {
70  pncmesh = new ParNCMesh(*pmesh.pncmesh);
71  pncmesh->OnMeshUpdated(this);
72  }
73  else
74  {
75  pncmesh = NULL;
76  }
77  ncmesh = pncmesh;
78 
79  // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
80  // and the FiniteElementSpace (as a ParFiniteElementSpace)
81  if (pmesh.Nodes && copy_nodes)
82  {
83  FiniteElementSpace *fes = pmesh.Nodes->FESpace();
84  const FiniteElementCollection *fec = fes->FEColl();
85  FiniteElementCollection *fec_copy =
87  ParFiniteElementSpace *pfes_copy =
88  new ParFiniteElementSpace(*fes, *this, fec_copy);
89  Nodes = new ParGridFunction(pfes_copy);
90  Nodes->MakeOwner(fec_copy);
91  *Nodes = *pmesh.Nodes;
92  own_nodes = 1;
93  }
94 }
95 
97 {
98  Swap(mesh);
99 }
100 
102 {
103  Swap(mesh);
104  return *this;
105 }
106 
107 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
108  int part_method)
109  : face_nbr_el_to_face(NULL)
110  , glob_elem_offset(-1)
111  , glob_offset_sequence(-1)
112  , gtopo(comm)
113 {
114  int *partitioning = NULL;
115  Array<bool> activeBdrElem;
116 
117  MyComm = comm;
118  MPI_Comm_size(MyComm, &NRanks);
119  MPI_Comm_rank(MyComm, &MyRank);
120 
121  if (mesh.Nonconforming())
122  {
123  if (partitioning_)
124  {
125  partitioning = partitioning_;
126  }
127  ncmesh = pncmesh = new ParNCMesh(comm, *mesh.ncmesh, partitioning);
128  if (!partitioning)
129  {
130  partitioning = new int[mesh.GetNE()];
131  for (int i = 0; i < mesh.GetNE(); i++)
132  {
133  partitioning[i] = pncmesh->InitialPartition(i);
134  }
135  }
136 
137  pncmesh->Prune();
138 
140  pncmesh->OnMeshUpdated(this);
141 
143 
144  // SetMeshGen(); // called by Mesh::InitFromNCMesh(...) above
145  meshgen = mesh.meshgen; // copy the global 'meshgen'
146 
147  mesh.attributes.Copy(attributes);
149 
151  }
152  else // mesh.Conforming()
153  {
154  Dim = mesh.Dim;
155  spaceDim = mesh.spaceDim;
156 
157  ncmesh = pncmesh = NULL;
158 
159  if (partitioning_)
160  {
161  partitioning = partitioning_;
162  }
163  else
164  {
165  partitioning = mesh.GeneratePartitioning(NRanks, part_method);
166  }
167 
168  // re-enumerate the partitions to better map to actual processor
169  // interconnect topology !?
170 
171  Array<int> vert_global_local;
172  NumOfVertices = BuildLocalVertices(mesh, partitioning, vert_global_local);
173  NumOfElements = BuildLocalElements(mesh, partitioning, vert_global_local);
174 
175  Table *edge_element = NULL;
176  NumOfBdrElements = BuildLocalBoundary(mesh, partitioning,
177  vert_global_local,
178  activeBdrElem, edge_element);
179 
180  SetMeshGen();
181  meshgen = mesh.meshgen; // copy the global 'meshgen'
182 
183  mesh.attributes.Copy(attributes);
185 
186  NumOfEdges = NumOfFaces = 0;
187 
188  if (Dim > 1)
189  {
190  el_to_edge = new Table;
192  }
193 
194  STable3D *faces_tbl = NULL;
195  if (Dim == 3)
196  {
197  faces_tbl = GetElementToFaceTable(1);
198  }
199 
200  GenerateFaces();
201 
202  ListOfIntegerSets groups;
203  {
204  // the first group is the local one
205  IntegerSet group;
206  group.Recreate(1, &MyRank);
207  groups.Insert(group);
208  }
209 
210  MFEM_ASSERT(mesh.GetNFaces() == 0 || Dim >= 3, "");
211 
212  Array<int> face_group(mesh.GetNFaces());
213  Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
214 
215  FindSharedFaces(mesh, partitioning, face_group, groups);
216  int nsedges = FindSharedEdges(mesh, partitioning, edge_element, groups);
217  int nsvert = FindSharedVertices(partitioning, vert_element, groups);
218 
219  // build the group communication topology
220  gtopo.Create(groups, 822);
221 
222  // fill out group_sface, group_sedge, group_svert
223  int ngroups = groups.Size()-1, nstris, nsquads;
224  BuildFaceGroup(ngroups, mesh, face_group, nstris, nsquads);
225  BuildEdgeGroup(ngroups, *edge_element);
226  BuildVertexGroup(ngroups, *vert_element);
227 
228  // build shared_faces and sface_lface mapping
229  BuildSharedFaceElems(nstris, nsquads, mesh, partitioning, faces_tbl,
230  face_group, vert_global_local);
231  delete faces_tbl;
232 
233  // build shared_edges and sedge_ledge mapping
234  BuildSharedEdgeElems(nsedges, mesh, vert_global_local, edge_element);
235  delete edge_element;
236 
237  // build svert_lvert mapping
238  BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
239  delete vert_element;
240 
241  SetMeshGen();
242  meshgen = mesh.meshgen; // copy the global 'meshgen'
243  }
244 
245  if (mesh.NURBSext)
246  {
247  MFEM_ASSERT(mesh.GetNodes() &&
248  mesh.GetNodes()->FESpace()->GetNURBSext() == mesh.NURBSext,
249  "invalid NURBS mesh");
250  NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
251  activeBdrElem);
252  }
253 
254  if (mesh.GetNodes()) // curved mesh
255  {
256  if (!NURBSext)
257  {
258  Nodes = new ParGridFunction(this, mesh.GetNodes());
259  }
260  else
261  {
262  const FiniteElementSpace *glob_fes = mesh.GetNodes()->FESpace();
264  FiniteElementCollection::New(glob_fes->FEColl()->Name());
265  ParFiniteElementSpace *pfes =
266  new ParFiniteElementSpace(this, nfec, glob_fes->GetVDim(),
267  glob_fes->GetOrdering());
268  Nodes = new ParGridFunction(pfes);
269  Nodes->MakeOwner(nfec); // Nodes will own nfec and pfes
270  }
271  own_nodes = 1;
272 
273  Array<int> gvdofs, lvdofs;
274  Vector lnodes;
275  int element_counter = 0;
276  for (int i = 0; i < mesh.GetNE(); i++)
277  {
278  if (partitioning[i] == MyRank)
279  {
280  Nodes->FESpace()->GetElementVDofs(element_counter, lvdofs);
281  mesh.GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
282  mesh.GetNodes()->GetSubVector(gvdofs, lnodes);
283  Nodes->SetSubVector(lvdofs, lnodes);
284  element_counter++;
285  }
286  }
287 
288  // set meaningful values to 'vertices' even though we have Nodes,
289  // for compatibility (e.g., Mesh::GetVertex())
291  }
292 
293  if (partitioning != partitioning_)
294  {
295  delete [] partitioning;
296  }
297 
298  have_face_nbr_data = false;
299 }
300 
301 
303  const int* partitioning,
304  Array<int> &vert_global_local)
305 {
306  vert_global_local.SetSize(mesh.GetNV());
307  vert_global_local = -1;
308 
309  int vert_counter = 0;
310  for (int i = 0; i < mesh.GetNE(); i++)
311  {
312  if (partitioning[i] == MyRank)
313  {
314  Array<int> vert;
315  mesh.GetElementVertices(i, vert);
316  for (int j = 0; j < vert.Size(); j++)
317  {
318  if (vert_global_local[vert[j]] < 0)
319  {
320  vert_global_local[vert[j]] = vert_counter++;
321  }
322  }
323  }
324  }
325 
326  // re-enumerate the local vertices to preserve the global ordering
327  vert_counter = 0;
328  for (int i = 0; i < vert_global_local.Size(); i++)
329  {
330  if (vert_global_local[i] >= 0)
331  {
332  vert_global_local[i] = vert_counter++;
333  }
334  }
335 
336  vertices.SetSize(vert_counter);
337 
338  for (int i = 0; i < vert_global_local.Size(); i++)
339  {
340  if (vert_global_local[i] >= 0)
341  {
342  vertices[vert_global_local[i]].SetCoords(mesh.SpaceDimension(),
343  mesh.GetVertex(i));
344  }
345  }
346 
347  return vert_counter;
348 }
349 
350 int ParMesh::BuildLocalElements(const Mesh& mesh, const int* partitioning,
351  const Array<int>& vert_global_local)
352 {
353  int nelems = 0;
354  for (int i = 0; i < mesh.GetNE(); i++)
355  {
356  if (partitioning[i] == MyRank) { nelems++; }
357  }
358 
359  elements.SetSize(nelems);
360 
361  // Determine elements, enumerating the local elements to preserve the global
362  // order. This is used, e.g. by the ParGridFunction ctor that takes a global
363  // GridFunction as input parameter.
364  int element_counter = 0;
365  for (int i = 0; i < mesh.GetNE(); i++)
366  {
367  if (partitioning[i] == MyRank)
368  {
369  elements[element_counter] = mesh.GetElement(i)->Duplicate(this);
370  int *v = elements[element_counter]->GetVertices();
371  int nv = elements[element_counter]->GetNVertices();
372  for (int j = 0; j < nv; j++)
373  {
374  v[j] = vert_global_local[v[j]];
375  }
376  element_counter++;
377  }
378  }
379 
380  return element_counter;
381 }
382 
383 int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
384  const Array<int>& vert_global_local,
385  Array<bool>& activeBdrElem,
386  Table*& edge_element)
387 {
388  int nbdry = 0;
389 
390  if (mesh.NURBSext)
391  {
392  activeBdrElem.SetSize(mesh.GetNBE());
393  activeBdrElem = false;
394  }
395  // build boundary elements
396  if (Dim == 3)
397  {
398  for (int i = 0; i < mesh.GetNBE(); i++)
399  {
400  int face, o, el1, el2;
401  mesh.GetBdrElementFace(i, &face, &o);
402  mesh.GetFaceElements(face, &el1, &el2);
403  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
404  {
405  nbdry++;
406  if (mesh.NURBSext)
407  {
408  activeBdrElem[i] = true;
409  }
410  }
411  }
412 
413  int bdrelem_counter = 0;
414  boundary.SetSize(nbdry);
415  for (int i = 0; i < mesh.GetNBE(); i++)
416  {
417  int face, o, el1, el2;
418  mesh.GetBdrElementFace(i, &face, &o);
419  mesh.GetFaceElements(face, &el1, &el2);
420  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
421  {
422  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
423  int *v = boundary[bdrelem_counter]->GetVertices();
424  int nv = boundary[bdrelem_counter]->GetNVertices();
425  for (int j = 0; j < nv; j++)
426  {
427  v[j] = vert_global_local[v[j]];
428  }
429  bdrelem_counter++;
430  }
431  }
432  }
433  else if (Dim == 2)
434  {
435  edge_element = new Table;
436  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
437 
438  for (int i = 0; i < mesh.GetNBE(); i++)
439  {
440  int edge = mesh.GetBdrElementEdgeIndex(i);
441  int el1 = edge_element->GetRow(edge)[0];
442  if (partitioning[el1] == MyRank)
443  {
444  nbdry++;
445  if (mesh.NURBSext)
446  {
447  activeBdrElem[i] = true;
448  }
449  }
450  }
451 
452  int bdrelem_counter = 0;
453  boundary.SetSize(nbdry);
454  for (int i = 0; i < mesh.GetNBE(); i++)
455  {
456  int edge = mesh.GetBdrElementEdgeIndex(i);
457  int el1 = edge_element->GetRow(edge)[0];
458  if (partitioning[el1] == MyRank)
459  {
460  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
461  int *v = boundary[bdrelem_counter]->GetVertices();
462  int nv = boundary[bdrelem_counter]->GetNVertices();
463  for (int j = 0; j < nv; j++)
464  {
465  v[j] = vert_global_local[v[j]];
466  }
467  bdrelem_counter++;
468  }
469  }
470  }
471  else if (Dim == 1)
472  {
473  for (int i = 0; i < mesh.GetNBE(); i++)
474  {
475  int vert = mesh.boundary[i]->GetVertices()[0];
476  int el1, el2;
477  mesh.GetFaceElements(vert, &el1, &el2);
478  if (partitioning[el1] == MyRank)
479  {
480  nbdry++;
481  }
482  }
483 
484  int bdrelem_counter = 0;
485  boundary.SetSize(nbdry);
486  for (int i = 0; i < mesh.GetNBE(); i++)
487  {
488  int vert = mesh.boundary[i]->GetVertices()[0];
489  int el1, el2;
490  mesh.GetFaceElements(vert, &el1, &el2);
491  if (partitioning[el1] == MyRank)
492  {
493  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
494  int *v = boundary[bdrelem_counter]->GetVertices();
495  v[0] = vert_global_local[v[0]];
496  bdrelem_counter++;
497  }
498  }
499  }
500 
501  return nbdry;
502 }
503 
504 void ParMesh::FindSharedFaces(const Mesh &mesh, const int *partitioning,
505  Array<int> &face_group,
506  ListOfIntegerSets &groups)
507 {
508  IntegerSet group;
509 
510  // determine shared faces
511  face_group.SetSize(mesh.GetNFaces());
512  for (int i = 0; i < face_group.Size(); i++)
513  {
514  int el[2];
515  face_group[i] = -1;
516  mesh.GetFaceElements(i, &el[0], &el[1]);
517  if (el[1] >= 0)
518  {
519  el[0] = partitioning[el[0]];
520  el[1] = partitioning[el[1]];
521  if ((el[0] == MyRank && el[1] != MyRank) ||
522  (el[0] != MyRank && el[1] == MyRank))
523  {
524  group.Recreate(2, el);
525  face_group[i] = groups.Insert(group) - 1;
526  }
527  }
528  }
529 }
530 
531 int ParMesh::FindSharedEdges(const Mesh &mesh, const int *partitioning,
532  Table*& edge_element,
533  ListOfIntegerSets& groups)
534 {
535  IntegerSet group;
536 
537  // determine shared edges
538  int sedge_counter = 0;
539  if (!edge_element)
540  {
541  edge_element = new Table;
542  if (Dim == 1)
543  {
544  edge_element->SetDims(0,0);
545  }
546  else
547  {
548  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
549  }
550  }
551 
552  for (int i = 0; i < edge_element->Size(); i++)
553  {
554  int me = 0, others = 0;
555  for (int j = edge_element->GetI()[i]; j < edge_element->GetI()[i+1]; j++)
556  {
557  int k = edge_element->GetJ()[j];
558  int rank = partitioning[k];
559  edge_element->GetJ()[j] = rank;
560  if (rank == MyRank)
561  {
562  me = 1;
563  }
564  else
565  {
566  others = 1;
567  }
568  }
569 
570  if (me && others)
571  {
572  sedge_counter++;
573  group.Recreate(edge_element->RowSize(i), edge_element->GetRow(i));
574  edge_element->GetRow(i)[0] = groups.Insert(group) - 1;
575  }
576  else
577  {
578  edge_element->GetRow(i)[0] = -1;
579  }
580  }
581 
582  return sedge_counter;
583 }
584 
585 int ParMesh::FindSharedVertices(const int *partitioning, Table *vert_element,
586  ListOfIntegerSets &groups)
587 {
588  IntegerSet group;
589 
590  // determine shared vertices
591  int svert_counter = 0;
592  for (int i = 0; i < vert_element->Size(); i++)
593  {
594  int me = 0, others = 0;
595  for (int j = vert_element->GetI()[i]; j < vert_element->GetI()[i+1]; j++)
596  {
597  vert_element->GetJ()[j] = partitioning[vert_element->GetJ()[j]];
598  if (vert_element->GetJ()[j] == MyRank)
599  {
600  me = 1;
601  }
602  else
603  {
604  others = 1;
605  }
606  }
607 
608  if (me && others)
609  {
610  svert_counter++;
611  group.Recreate(vert_element->RowSize(i), vert_element->GetRow(i));
612  vert_element->GetI()[i] = groups.Insert(group) - 1;
613  }
614  else
615  {
616  vert_element->GetI()[i] = -1;
617  }
618  }
619  return svert_counter;
620 }
621 
622 void ParMesh::BuildFaceGroup(int ngroups, const Mesh &mesh,
623  const Array<int> &face_group,
624  int &nstria, int &nsquad)
625 {
626  // build group_stria and group_squad
627  group_stria.MakeI(ngroups);
628  group_squad.MakeI(ngroups);
629 
630  for (int i = 0; i < face_group.Size(); i++)
631  {
632  if (face_group[i] >= 0)
633  {
634  if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
635  {
636  group_stria.AddAColumnInRow(face_group[i]);
637  }
638  else
639  {
640  group_squad.AddAColumnInRow(face_group[i]);
641  }
642  }
643  }
644 
645  group_stria.MakeJ();
646  group_squad.MakeJ();
647 
648  nstria = nsquad = 0;
649  for (int i = 0; i < face_group.Size(); i++)
650  {
651  if (face_group[i] >= 0)
652  {
653  if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
654  {
655  group_stria.AddConnection(face_group[i], nstria++);
656  }
657  else
658  {
659  group_squad.AddConnection(face_group[i], nsquad++);
660  }
661  }
662  }
663 
666 }
667 
668 void ParMesh::BuildEdgeGroup(int ngroups, const Table &edge_element)
669 {
670  group_sedge.MakeI(ngroups);
671 
672  for (int i = 0; i < edge_element.Size(); i++)
673  {
674  if (edge_element.GetRow(i)[0] >= 0)
675  {
676  group_sedge.AddAColumnInRow(edge_element.GetRow(i)[0]);
677  }
678  }
679 
680  group_sedge.MakeJ();
681 
682  int sedge_counter = 0;
683  for (int i = 0; i < edge_element.Size(); i++)
684  {
685  if (edge_element.GetRow(i)[0] >= 0)
686  {
687  group_sedge.AddConnection(edge_element.GetRow(i)[0], sedge_counter++);
688  }
689  }
690 
692 }
693 
694 void ParMesh::BuildVertexGroup(int ngroups, const Table &vert_element)
695 {
696  group_svert.MakeI(ngroups);
697 
698  for (int i = 0; i < vert_element.Size(); i++)
699  {
700  if (vert_element.GetI()[i] >= 0)
701  {
702  group_svert.AddAColumnInRow(vert_element.GetI()[i]);
703  }
704  }
705 
706  group_svert.MakeJ();
707 
708  int svert_counter = 0;
709  for (int i = 0; i < vert_element.Size(); i++)
710  {
711  if (vert_element.GetI()[i] >= 0)
712  {
713  group_svert.AddConnection(vert_element.GetI()[i], svert_counter++);
714  }
715  }
716 
718 }
719 
720 void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
721  const Mesh& mesh, int *partitioning,
722  const STable3D *faces_tbl,
723  const Array<int> &face_group,
724  const Array<int> &vert_global_local)
725 {
726  shared_trias.SetSize(ntri_faces);
727  shared_quads.SetSize(nquad_faces);
728  sface_lface. SetSize(ntri_faces + nquad_faces);
729 
730  if (Dim < 3) { return; }
731 
732  int stria_counter = 0;
733  int squad_counter = 0;
734  for (int i = 0; i < face_group.Size(); i++)
735  {
736  if (face_group[i] < 0) { continue; }
737 
738  const Element *face = mesh.GetFace(i);
739  const int *fv = face->GetVertices();
740  switch (face->GetType())
741  {
742  case Element::TRIANGLE:
743  {
744  shared_trias[stria_counter].Set(fv);
745  int *v = shared_trias[stria_counter].v;
746  for (int j = 0; j < 3; j++)
747  {
748  v[j] = vert_global_local[v[j]];
749  }
750  const int lface = (*faces_tbl)(v[0], v[1], v[2]);
751  sface_lface[stria_counter] = lface;
752  if (meshgen == 1) // Tet-only mesh
753  {
754  Tetrahedron *tet = dynamic_cast<Tetrahedron *>
755  (elements[faces_info[lface].Elem1No]);
756  // mark the shared face for refinement by reorienting
757  // it according to the refinement flag in the tetrahedron
758  // to which this shared face belongs to.
759  if (tet->GetRefinementFlag())
760  {
761  tet->GetMarkedFace(faces_info[lface].Elem1Inf/64, v);
762  // flip the shared face in the processor that owns the
763  // second element (in 'mesh')
764  int gl_el1, gl_el2;
765  mesh.GetFaceElements(i, &gl_el1, &gl_el2);
766  if (MyRank == partitioning[gl_el2])
767  {
768  std::swap(v[0], v[1]);
769  }
770  }
771  }
772  stria_counter++;
773  break;
774  }
775 
777  {
778  shared_quads[squad_counter].Set(fv);
779  int *v = shared_quads[squad_counter].v;
780  for (int j = 0; j < 4; j++)
781  {
782  v[j] = vert_global_local[v[j]];
783  }
784  sface_lface[shared_trias.Size() + squad_counter] =
785  (*faces_tbl)(v[0], v[1], v[2], v[3]);
786  squad_counter++;
787  break;
788  }
789 
790  default:
791  MFEM_ABORT("unknown face type: " << face->GetType());
792  break;
793  }
794  }
795 }
796 
797 void ParMesh::BuildSharedEdgeElems(int nedges, Mesh& mesh,
798  const Array<int>& vert_global_local,
799  const Table* edge_element)
800 {
801  // The passed in mesh is still the global mesh. "this" mesh is the
802  // local partitioned mesh.
803 
804  shared_edges.SetSize(nedges);
805  sedge_ledge. SetSize(nedges);
806 
807  {
808  DSTable v_to_v(NumOfVertices);
809  GetVertexToVertexTable(v_to_v);
810 
811  int sedge_counter = 0;
812  for (int i = 0; i < edge_element->Size(); i++)
813  {
814  if (edge_element->GetRow(i)[0] >= 0)
815  {
816  Array<int> vert;
817  mesh.GetEdgeVertices(i, vert);
818 
819  shared_edges[sedge_counter] =
820  new Segment(vert_global_local[vert[0]],
821  vert_global_local[vert[1]], 1);
822 
823  sedge_ledge[sedge_counter] = v_to_v(vert_global_local[vert[0]],
824  vert_global_local[vert[1]]);
825 
826  MFEM_VERIFY(sedge_ledge[sedge_counter] >= 0, "Error in v_to_v.");
827 
828  sedge_counter++;
829  }
830  }
831  }
832 }
833 
835  const mfem::Table *vert_element,
836  const Array<int> &vert_global_local)
837 {
838  // build svert_lvert
839  svert_lvert.SetSize(nvert);
840 
841  int svert_counter = 0;
842  for (int i = 0; i < vert_element->Size(); i++)
843  {
844  if (vert_element->GetI()[i] >= 0)
845  {
846  svert_lvert[svert_counter++] = vert_global_local[i];
847  }
848  }
849 }
850 
851 
852 // protected method, used by Nonconforming(De)Refinement and Rebalance
854  : MyComm(pncmesh.MyComm)
855  , NRanks(pncmesh.NRanks)
856  , MyRank(pncmesh.MyRank)
857  , face_nbr_el_to_face(NULL)
858  , glob_elem_offset(-1)
859  , glob_offset_sequence(-1)
860  , gtopo(MyComm)
861  , pncmesh(NULL)
862 {
863  Mesh::InitFromNCMesh(pncmesh);
864  ReduceMeshGen();
865  have_face_nbr_data = false;
866 }
867 
869 {
870  if (glob_offset_sequence != sequence) // mesh has changed
871  {
872  long 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(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::HEXAHEDRON:
2693  {
2694  // find the face by the vertices with the smallest 3 numbers
2695  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
2696  for (int j = 0; j < 6; j++)
2697  {
2698  const int *fv = hex_t::FaceVert[j];
2699  sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2700  }
2701  break;
2702  }
2703  default:
2704  MFEM_ABORT("Unexpected type of Element.");
2705  }
2706  }
2707  return sfaces_tbl;
2708 }
2709 
2711 {
2712  int i, *v;
2713  STable3D * faces_tbl = GetFacesTable();
2714  STable3D * sfaces_tbl = GetSharedFacesTable();
2715 
2716  if (face_nbr_el_to_face != NULL)
2717  {
2718  delete face_nbr_el_to_face;
2719  }
2720  face_nbr_el_to_face = new Table(face_nbr_elements.Size(), 6);
2721  for (i = 0; i < face_nbr_elements.Size(); i++)
2722  {
2723  v = face_nbr_elements[i]->GetVertices();
2724  switch (face_nbr_elements[i]->GetType())
2725  {
2726  case Element::TETRAHEDRON:
2727  {
2728  for (int j = 0; j < 4; j++)
2729  {
2730  const int *fv = tet_t::FaceVert[j];
2731  int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2732  if (lf < 0)
2733  {
2734  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2735  if (lf >= 0)
2736  {
2737  lf += NumOfFaces;
2738  }
2739  }
2740  face_nbr_el_to_face->Push(i, lf);
2741  }
2742  break;
2743  }
2744  case Element::WEDGE:
2745  {
2746  for (int j = 0; j < 2; j++)
2747  {
2748  const int *fv = pri_t::FaceVert[j];
2749  int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2750  if (lf < 0)
2751  {
2752  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2753  if (lf >= 0)
2754  {
2755  lf += NumOfFaces;
2756  }
2757  }
2758  face_nbr_el_to_face->Push(i, lf);
2759  }
2760  for (int j = 2; j < 5; j++)
2761  {
2762  const int *fv = pri_t::FaceVert[j];
2763  int k = 0;
2764  int max = v[fv[0]];
2765 
2766  if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2767  if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2768  if (max < v[fv[3]]) { k = 3; }
2769 
2770  int v0 = -1, v1 = -1, v2 = -1;
2771  switch (k)
2772  {
2773  case 0:
2774  v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2775  break;
2776  case 1:
2777  v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2778  break;
2779  case 2:
2780  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2781  break;
2782  case 3:
2783  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2784  break;
2785  }
2786  int lf = faces_tbl->Index(v0, v1, v2);
2787  if (lf < 0)
2788  {
2789  lf = sfaces_tbl->Index(v0, v1, v2);
2790  if (lf >= 0)
2791  {
2792  lf += NumOfFaces;
2793  }
2794  }
2795  face_nbr_el_to_face->Push(i, lf);
2796  }
2797  break;
2798  }
2799  case Element::HEXAHEDRON:
2800  {
2801  // find the face by the vertices with the smallest 3 numbers
2802  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
2803  for (int j = 0; j < 6; j++)
2804  {
2805  const int *fv = hex_t::FaceVert[j];
2806  int k = 0;
2807  int max = v[fv[0]];
2808 
2809  if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2810  if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2811  if (max < v[fv[3]]) { k = 3; }
2812 
2813  int v0 = -1, v1 = -1, v2 = -1;
2814  switch (k)
2815  {
2816  case 0:
2817  v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2818  break;
2819  case 1:
2820  v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2821  break;
2822  case 2:
2823  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2824  break;
2825  case 3:
2826  v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2827  break;
2828  }
2829  int lf = faces_tbl->Index(v0, v1, v2);
2830  if (lf < 0)
2831  {
2832  lf = sfaces_tbl->Index(v0, v1, v2);
2833  if (lf >= 0)
2834  {
2835  lf += NumOfFaces;
2836  }
2837  }
2838  face_nbr_el_to_face->Push(i, lf);
2839  }
2840  break;
2841  }
2842  default:
2843  MFEM_ABORT("Unexpected type of Element.");
2844  }
2845  }
2847 
2848  delete sfaces_tbl;
2849  if (ret_ftbl)
2850  {
2851  return faces_tbl;
2852  }
2853  delete faces_tbl;
2854  return NULL;
2855 }
2856 
2857 int ParMesh::GetFaceNbrRank(int fn) const
2858 {
2859  if (Conforming())
2860  {
2861  int nbr_group = face_nbr_group[fn];
2862  const int *nbs = gtopo.GetGroup(nbr_group);
2863  int nbr_lproc = (nbs[0]) ? nbs[0] : nbs[1];
2864  int nbr_rank = gtopo.GetNeighborRank(nbr_lproc);
2865  return nbr_rank;
2866  }
2867  else
2868  {
2869  // NC: simplified handling of face neighbor ranks
2870  return face_nbr_group[fn];
2871  }
2872 }
2873 
2874 void
2876 {
2877  int n, j;
2878  int el_nbr = i - GetNE();
2879  if (face_nbr_el_to_face)
2880  {
2881  face_nbr_el_to_face->GetRow(el_nbr, fcs);
2882  }
2883  else
2884  {
2885  MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2886  "face_nbr_el_to_face not generated.");
2887  }
2888  if (el_nbr < face_nbr_el_ori.Size())
2889  {
2890  const int * row = face_nbr_el_ori.GetRow(el_nbr);
2891  n = fcs.Size();
2892  cor.SetSize(n);
2893  for (j=0; j<n; j++)
2894  {
2895  cor[j] = row[j];
2896  }
2897  }
2898  else
2899  {
2900  MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2901  "face_nbr_el_to_face not generated.");
2902  }
2903 }
2904 
2906 {
2907  const Array<int> *s2l_face;
2908  if (Dim == 1)
2909  {
2910  s2l_face = &svert_lvert;
2911  }
2912  else if (Dim == 2)
2913  {
2914  s2l_face = &sedge_ledge;
2915  }
2916  else
2917  {
2918  s2l_face = &sface_lface;
2919  }
2920 
2921  Table *face_elem = new Table;
2922 
2923  face_elem->MakeI(faces_info.Size());
2924 
2925  for (int i = 0; i < faces_info.Size(); i++)
2926  {
2927  if (faces_info[i].Elem2No >= 0)
2928  {
2929  face_elem->AddColumnsInRow(i, 2);
2930  }
2931  else
2932  {
2933  face_elem->AddAColumnInRow(i);
2934  }
2935  }
2936  for (int i = 0; i < s2l_face->Size(); i++)
2937  {
2938  face_elem->AddAColumnInRow((*s2l_face)[i]);
2939  }
2940 
2941  face_elem->MakeJ();
2942 
2943  for (int i = 0; i < faces_info.Size(); i++)
2944  {
2945  face_elem->AddConnection(i, faces_info[i].Elem1No);
2946  if (faces_info[i].Elem2No >= 0)
2947  {
2948  face_elem->AddConnection(i, faces_info[i].Elem2No);
2949  }
2950  }
2951  for (int i = 0; i < s2l_face->Size(); i++)
2952  {
2953  int lface = (*s2l_face)[i];
2954  int nbr_elem_idx = -1 - faces_info[lface].Elem2No;
2955  face_elem->AddConnection(lface, NumOfElements + nbr_elem_idx);
2956  }
2957 
2958  face_elem->ShiftUpI();
2959 
2960  return face_elem;
2961 }
2962 
2964  FaceElementTransformations* FETr, Element::Type face_type,
2965  Geometry::Type face_geom)
2966 {
2967  // calculate composition of FETr->Loc1 and FETr->Elem1
2968  DenseMatrix &face_pm = FETr->GetPointMat();
2969  FETr->Reset();
2970  if (Nodes == NULL)
2971  {
2972  FETr->Elem1->Transform(FETr->Loc1.Transf.GetPointMat(), face_pm);
2973  FETr->SetFE(GetTransformationFEforElementType(face_type));
2974  }
2975  else
2976  {
2977  const FiniteElement* face_el =
2978  Nodes->FESpace()->GetTraceElement(FETr->Elem1No, face_geom);
2979  MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
2980  "Mesh requires nodal Finite Element.");
2981 
2982 #if 0 // TODO: handle the case of non-interpolatory Nodes
2983  DenseMatrix I;
2984  face_el->Project(Transformation.GetFE(), FETr->Loc1.Transf, I);
2985  MultABt(Transformation.GetPointMat(), I, pm_face);
2986 #else
2987  IntegrationRule eir(face_el->GetDof());
2988  FETr->Loc1.Transform(face_el->GetNodes(), eir);
2989  Nodes->GetVectorValues(*FETr->Elem1, eir, face_pm);
2990 #endif
2991  FETr->SetFE(face_el);
2992  }
2993 }
2994 
2996  int FaceNo,
2997  int mask)
2998 {
2999  if (FaceNo < GetNumFaces())
3000  {
3001  return Mesh::GetFaceElementTransformations(FaceNo, mask);
3002  }
3003  else
3004  {
3005  const bool fill2 = mask & 10; // Elem2 and/or Loc2
3006  return GetSharedFaceTransformationsByLocalIndex(FaceNo, fill2);
3007  }
3008 }
3009 
3011 GetSharedFaceTransformations(int sf, bool fill2)
3012 {
3013  int FaceNo = GetSharedFace(sf);
3014 
3015  return GetSharedFaceTransformationsByLocalIndex(FaceNo, fill2);
3016 }
3017 
3020 {
3021  FaceInfo &face_info = faces_info[FaceNo];
3022  MFEM_VERIFY(face_info.Elem2Inf >= 0, "The face must be shared.");
3023 
3024  bool is_slave = Nonconforming() && IsSlaveFace(face_info);
3025  bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
3026 
3027  int mask = 0;
3029  FaceElemTr.Elem1 = NULL;
3030  FaceElemTr.Elem2 = NULL;
3031 
3032  NCFaceInfo* nc_info = NULL;
3033  if (is_slave) { nc_info = &nc_faces_info[face_info.NCFace]; }
3034 
3035  int local_face = is_ghost ? nc_info->MasterFace : FaceNo;
3036  Element::Type face_type = GetFaceElementType(local_face);
3037  Geometry::Type face_geom = GetFaceGeometryType(local_face);
3038 
3039  // setup the transformation for the first element
3040  FaceElemTr.Elem1No = face_info.Elem1No;
3044 
3045  // setup the transformation for the second (neighbor) element
3046  int Elem2NbrNo;
3047  if (fill2)
3048  {
3049  Elem2NbrNo = -1 - face_info.Elem2No;
3050  // Store the "shifted index" for element 2 in FaceElemTr.Elem2No.
3051  // `Elem2NbrNo` is the index of the face neighbor (starting from 0),
3052  // and `FaceElemTr.Elem2No` will be offset by the number of (local)
3053  // elements in the mesh.
3054  FaceElemTr.Elem2No = NumOfElements + Elem2NbrNo;
3058  }
3059  else
3060  {
3061  FaceElemTr.Elem2No = -1;
3062  }
3063 
3064  // setup the face transformation if the face is not a ghost
3065  if (!is_ghost)
3066  {
3068  // NOTE: The above call overwrites FaceElemTr.Loc1
3070  }
3071  else
3072  {
3073  FaceElemTr.SetGeometryType(face_geom);
3074  }
3075 
3076  // setup Loc1 & Loc2
3077  int elem_type = GetElementType(face_info.Elem1No);
3078  GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc1.Transf,
3079  face_info.Elem1Inf);
3081 
3082  if (fill2)
3083  {
3084  elem_type = face_nbr_elements[Elem2NbrNo]->GetType();
3085  GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc2.Transf,
3086  face_info.Elem2Inf);
3088  }
3089 
3090  // adjust Loc1 or Loc2 of the master face if this is a slave face
3091  if (is_slave)
3092  {
3093  if (is_ghost || fill2)
3094  {
3095  // is_ghost -> modify side 1, otherwise -> modify side 2:
3096  ApplyLocalSlaveTransformation(FaceElemTr, face_info, is_ghost);
3097  }
3098  }
3099 
3100  // for ghost faces we need a special version of GetFaceTransformation
3101  if (is_ghost)
3102  {
3103  GetGhostFaceTransformation(&FaceElemTr, face_type, face_geom);
3105  }
3106 
3108 
3109  // This check can be useful for internal debugging, however it will fail on
3110  // periodic boundary faces, so we keep it disabled in general.
3111 #if 0
3112 #ifdef MFEM_DEBUG
3113  double dist = FaceElemTr.CheckConsistency();
3114  if (dist >= 1e-12)
3115  {
3116  mfem::out << "\nInternal error: face id = " << FaceNo
3117  << ", dist = " << dist << ", rank = " << MyRank << '\n';
3118  FaceElemTr.CheckConsistency(1); // print coordinates
3119  MFEM_ABORT("internal error");
3120  }
3121 #endif
3122 #endif
3123 
3124  return &FaceElemTr;
3125 }
3126 
3128 {
3129  if (Conforming())
3130  {
3131  switch (Dim)
3132  {
3133  case 1: return svert_lvert.Size();
3134  case 2: return sedge_ledge.Size();
3135  default: return sface_lface.Size();
3136  }
3137  }
3138  else
3139  {
3140  MFEM_ASSERT(Dim > 1, "");
3141  const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3142  return shared.conforming.Size() + shared.slaves.Size();
3143  }
3144 }
3145 
3146 int ParMesh::GetSharedFace(int sface) const
3147 {
3148  if (Conforming())
3149  {
3150  switch (Dim)
3151  {
3152  case 1: return svert_lvert[sface];
3153  case 2: return sedge_ledge[sface];
3154  default: return sface_lface[sface];
3155  }
3156  }
3157  else
3158  {
3159  MFEM_ASSERT(Dim > 1, "");
3160  const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3161  int csize = (int) shared.conforming.Size();
3162  return sface < csize
3163  ? shared.conforming[sface].index
3164  : shared.slaves[sface - csize].index;
3165  }
3166 }
3167 
3169 {
3170  const_cast<ParMesh*>(this)->ExchangeFaceNbrData();
3171  return Mesh::GetNFbyType(type);
3172 }
3173 
3174 // shift cyclically 3 integers a, b, c, so that the smallest of
3175 // order[a], order[b], order[c] is first
3176 static inline
3177 void Rotate3Indirect(int &a, int &b, int &c,
3178  const Array<std::int64_t> &order)
3179 {
3180  if (order[a] < order[b])
3181  {
3182  if (order[a] > order[c])
3183  {
3184  ShiftRight(a, b, c);
3185  }
3186  }
3187  else
3188  {
3189  if (order[b] < order[c])
3190  {
3191  ShiftRight(c, b, a);
3192  }
3193  else
3194  {
3195  ShiftRight(a, b, c);
3196  }
3197  }
3198 }
3199 
3201 {
3202  if (Dim != 3 || !(meshgen & 1))
3203  {
3204  return;
3205  }
3206 
3207  ResetLazyData();
3208 
3209  DSTable *old_v_to_v = NULL;
3210  Table *old_elem_vert = NULL;
3211 
3212  if (Nodes)
3213  {
3214  PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3215  }
3216 
3217  // create a GroupCommunicator over shared vertices
3218  GroupCommunicator svert_comm(gtopo);
3219  GetSharedVertexCommunicator(svert_comm);
3220 
3221  // communicate the local index of each shared vertex from the group master to
3222  // other ranks in the group
3223  Array<int> svert_master_rank(svert_lvert.Size());
3224  Array<int> svert_master_index(svert_lvert);
3225  {
3226  for (int i = 0; i < group_svert.Size(); i++)
3227  {
3228  int rank = gtopo.GetGroupMasterRank(i+1);
3229  for (int j = 0; j < group_svert.RowSize(i); j++)
3230  {
3231  svert_master_rank[group_svert.GetRow(i)[j]] = rank;
3232  }
3233  }
3234  svert_comm.Bcast(svert_master_index);
3235  }
3236 
3237  // the pairs (master rank, master local index) define a globally consistent
3238  // vertex ordering
3239  Array<std::int64_t> glob_vert_order(vertices.Size());
3240  {
3241  Array<int> lvert_svert(vertices.Size());
3242  lvert_svert = -1;
3243  for (int i = 0; i < svert_lvert.Size(); i++)
3244  {
3245  lvert_svert[svert_lvert[i]] = i;
3246  }
3247 
3248  for (int i = 0; i < vertices.Size(); i++)
3249  {
3250  int s = lvert_svert[i];
3251  if (s >= 0)
3252  {
3253  glob_vert_order[i] =
3254  (std::int64_t(svert_master_rank[s]) << 32) + svert_master_index[s];
3255  }
3256  else
3257  {
3258  glob_vert_order[i] = (std::int64_t(MyRank) << 32) + i;
3259  }
3260  }
3261  }
3262 
3263  // rotate tetrahedra so that vertex zero is the lowest (global) index vertex,
3264  // vertex 1 is the second lowest (global) index and vertices 2 and 3 preserve
3265  // positive orientation of the element
3266  for (int i = 0; i < NumOfElements; i++)
3267  {
3269  {
3270  int *v = elements[i]->GetVertices();
3271 
3272  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3273 
3274  if (glob_vert_order[v[0]] < glob_vert_order[v[3]])
3275  {
3276  Rotate3Indirect(v[1], v[2], v[3], glob_vert_order);
3277  }
3278  else
3279  {
3280  ShiftRight(v[0], v[1], v[3]);
3281  }
3282  }
3283  }
3284 
3285  // rotate also boundary triangles
3286  for (int i = 0; i < NumOfBdrElements; i++)
3287  {
3289  {
3290  int *v = boundary[i]->GetVertices();
3291 
3292  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3293  }
3294  }
3295 
3296  const bool check_consistency = true;
3297  if (check_consistency)
3298  {
3299  // create a GroupCommunicator on the shared triangles
3300  GroupCommunicator stria_comm(gtopo);
3301  {
3302  // initialize stria_comm
3303  Table &gr_stria = stria_comm.GroupLDofTable();
3304  // gr_stria differs from group_stria - the latter does not store gr. 0
3305  gr_stria.SetDims(GetNGroups(), shared_trias.Size());
3306  gr_stria.GetI()[0] = 0;
3307  for (int gr = 1; gr <= GetNGroups(); gr++)
3308  {
3309  gr_stria.GetI()[gr] = group_stria.GetI()[gr-1];
3310  }
3311  for (int k = 0; k < shared_trias.Size(); k++)
3312  {
3313  gr_stria.GetJ()[k] = group_stria.GetJ()[k];
3314  }
3315  stria_comm.Finalize();
3316  }
3317  Array<int> stria_flag(shared_trias.Size());
3318  for (int i = 0; i < stria_flag.Size(); i++)
3319  {
3320  const int *v = shared_trias[i].v;
3321  if (glob_vert_order[v[0]] < glob_vert_order[v[1]])
3322  {
3323  stria_flag[i] = (glob_vert_order[v[0]] < glob_vert_order[v[2]]) ? 0 : 2;
3324  }
3325  else // v[1] < v[0]
3326  {
3327  stria_flag[i] = (glob_vert_order[v[1]] < glob_vert_order[v[2]]) ? 1 : 2;
3328  }
3329  }
3330 
3331  Array<int> stria_master_flag(stria_flag);
3332  stria_comm.Bcast(stria_master_flag);
3333  for (int i = 0; i < stria_flag.Size(); i++)
3334  {
3335  const int *v = shared_trias[i].v;
3336  MFEM_VERIFY(stria_flag[i] == stria_master_flag[i],
3337  "inconsistent vertex ordering found, shared triangle "
3338  << i << ": ("
3339  << v[0] << ", " << v[1] << ", " << v[2] << "), "
3340  << "local flag: " << stria_flag[i]
3341  << ", master flag: " << stria_master_flag[i]);
3342  }
3343  }
3344 
3345  // rotate shared triangle faces
3346  for (int i = 0; i < shared_trias.Size(); i++)
3347  {
3348  int *v = shared_trias[i].v;
3349 
3350  Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3351  }
3352 
3353  // finalize
3354  if (!Nodes)
3355  {
3357  GenerateFaces();
3358  if (el_to_edge)
3359  {
3361  }
3362  }
3363  else
3364  {
3365  DoNodeReorder(old_v_to_v, old_elem_vert);
3366  delete old_elem_vert;
3367  delete old_v_to_v;
3368  }
3369 
3370  // the local edge and face numbering is changed therefore we need to
3371  // update sedge_ledge and sface_lface.
3372  FinalizeParTopo();
3373 }
3374 
3375 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
3376 {
3377  if (pncmesh)
3378  {
3379  MFEM_ABORT("Local and nonconforming refinements cannot be mixed.");
3380  }
3381 
3383 
3385 
3386  if (Dim == 3)
3387  {
3388  int uniform_refinement = 0;
3389  if (type < 0)
3390  {
3391  type = -type;
3392  uniform_refinement = 1;
3393  }
3394 
3395  // 1. Hash table of vertex to vertex connections corresponding to refined
3396  // edges.
3397  HashTable<Hashed2> v_to_v;
3398 
3399  // 2. Do the red refinement.
3400  switch (type)
3401  {
3402  case 1:
3403  for (int i = 0; i < marked_el.Size(); i++)
3404  {
3405  Bisection(marked_el[i], v_to_v);
3406  }
3407  break;
3408  case 2:
3409  for (int i = 0; i < marked_el.Size(); i++)
3410  {
3411  Bisection(marked_el[i], v_to_v);
3412 
3413  Bisection(NumOfElements - 1, v_to_v);
3414  Bisection(marked_el[i], v_to_v);
3415  }
3416  break;
3417  case 3:
3418  for (int i = 0; i < marked_el.Size(); i++)
3419  {
3420  Bisection(marked_el[i], v_to_v);
3421 
3422  int j = NumOfElements - 1;
3423  Bisection(j, v_to_v);
3424  Bisection(NumOfElements - 1, v_to_v);
3425  Bisection(j, v_to_v);
3426 
3427  Bisection(marked_el[i], v_to_v);
3428  Bisection(NumOfElements-1, v_to_v);
3429  Bisection(marked_el[i], v_to_v);
3430  }
3431  break;
3432  }
3433 
3434  // 3. Do the green refinement (to get conforming mesh).
3435  int need_refinement;
3436  int max_faces_in_group = 0;
3437  // face_splittings identify how the shared faces have been split
3438  Array<unsigned> *face_splittings = new Array<unsigned>[GetNGroups()-1];
3439  for (int i = 0; i < GetNGroups()-1; i++)
3440  {
3441  const int faces_in_group = GroupNTriangles(i+1);
3442  face_splittings[i].Reserve(faces_in_group);
3443  if (faces_in_group > max_faces_in_group)
3444  {
3445  max_faces_in_group = faces_in_group;
3446  }
3447  }
3448  int neighbor;
3449  Array<unsigned> iBuf(max_faces_in_group);
3450 
3451  MPI_Request *requests = new MPI_Request[GetNGroups()-1];
3452  MPI_Status status;
3453 
3454 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3455  int ref_loops_all = 0, ref_loops_par = 0;
3456 #endif
3457  do
3458  {
3459  need_refinement = 0;
3460  for (int i = 0; i < NumOfElements; i++)
3461  {
3462  if (elements[i]->NeedRefinement(v_to_v))
3463  {
3464  need_refinement = 1;
3465  Bisection(i, v_to_v);
3466  }
3467  }
3468 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3469  ref_loops_all++;
3470 #endif
3471 
3472  if (uniform_refinement)
3473  {
3474  continue;
3475  }
3476 
3477  // if the mesh is locally conforming start making it globally
3478  // conforming
3479  if (need_refinement == 0)
3480  {
3481 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3482  ref_loops_par++;
3483 #endif
3484  // MPI_Barrier(MyComm);
3485  const int tag = 293;
3486 
3487  // (a) send the type of interface splitting
3488  int req_count = 0;
3489  for (int i = 0; i < GetNGroups()-1; i++)
3490  {
3491  const int *group_faces = group_stria.GetRow(i);
3492  const int faces_in_group = group_stria.RowSize(i);
3493  // it is enough to communicate through the faces
3494  if (faces_in_group == 0) { continue; }
3495 
3496  face_splittings[i].SetSize(0);
3497  for (int j = 0; j < faces_in_group; j++)
3498  {
3499  GetFaceSplittings(shared_trias[group_faces[j]].v, v_to_v,
3500  face_splittings[i]);
3501  }
3502  const int *nbs = gtopo.GetGroup(i+1);
3503  neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3504  MPI_Isend(face_splittings[i], face_splittings[i].Size(),
3505  MPI_UNSIGNED, neighbor, tag, MyComm,
3506  &requests[req_count++]);
3507  }
3508 
3509  // (b) receive the type of interface splitting
3510  for (int i = 0; i < GetNGroups()-1; i++)
3511  {
3512  const int *group_faces = group_stria.GetRow(i);
3513  const int faces_in_group = group_stria.RowSize(i);
3514  if (faces_in_group == 0) { continue; }
3515 
3516  const int *nbs = gtopo.GetGroup(i+1);
3517  neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3518  MPI_Probe(neighbor, tag, MyComm, &status);
3519  int count;
3520  MPI_Get_count(&status, MPI_UNSIGNED, &count);
3521  iBuf.SetSize(count);
3522  MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag, MyComm,
3523  MPI_STATUS_IGNORE);
3524 
3525  for (int j = 0, pos = 0; j < faces_in_group; j++)
3526  {
3527  const int *v = shared_trias[group_faces[j]].v;
3528  need_refinement |= DecodeFaceSplittings(v_to_v, v, iBuf, pos);
3529  }
3530  }
3531 
3532  int nr = need_refinement;
3533  MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3534 
3535  MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
3536  }
3537  }
3538  while (need_refinement == 1);
3539 
3540 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3541  {
3542  int i = ref_loops_all;
3543  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3544  if (MyRank == 0)
3545  {
3546  mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3547  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3548  << '\n' << endl;
3549  }
3550  }
3551 #endif
3552 
3553  delete [] requests;
3554  iBuf.DeleteAll();
3555  delete [] face_splittings;
3556 
3557  // 4. Update the boundary elements.
3558  do
3559  {
3560  need_refinement = 0;
3561  for (int i = 0; i < NumOfBdrElements; i++)
3562  {
3563  if (boundary[i]->NeedRefinement(v_to_v))
3564  {
3565  need_refinement = 1;
3566  BdrBisection(i, v_to_v);
3567  }
3568  }
3569  }
3570  while (need_refinement == 1);
3571 
3572  if (NumOfBdrElements != boundary.Size())
3573  {
3574  mfem_error("ParMesh::LocalRefinement :"
3575  " (NumOfBdrElements != boundary.Size())");
3576  }
3577 
3578  ResetLazyData();
3579 
3580  const int old_nv = NumOfVertices;
3581  NumOfVertices = vertices.Size();
3582 
3583  RefineGroups(old_nv, v_to_v);
3584 
3585  // 5. Update the groups after refinement.
3586  if (el_to_face != NULL)
3587  {
3589  GenerateFaces();
3590  }
3591 
3592  // 6. Update element-to-edge relations.
3593  if (el_to_edge != NULL)
3594  {
3596  }
3597  } // 'if (Dim == 3)'
3598 
3599 
3600  if (Dim == 2)
3601  {
3602  int uniform_refinement = 0;
3603  if (type < 0)
3604  {
3605  // type = -type; // not used
3606  uniform_refinement = 1;
3607  }
3608 
3609  // 1. Get table of vertex to vertex connections.
3610  DSTable v_to_v(NumOfVertices);
3611  GetVertexToVertexTable(v_to_v);
3612 
3613  // 2. Get edge to element connections in arrays edge1 and edge2
3614  int nedges = v_to_v.NumberOfEntries();
3615  int *edge1 = new int[nedges];
3616  int *edge2 = new int[nedges];
3617  int *middle = new int[nedges];
3618 
3619  for (int i = 0; i < nedges; i++)
3620  {
3621  edge1[i] = edge2[i] = middle[i] = -1;
3622  }
3623 
3624  for (int i = 0; i < NumOfElements; i++)
3625  {
3626  int *v = elements[i]->GetVertices();
3627  for (int j = 0; j < 3; j++)
3628  {
3629  int ind = v_to_v(v[j], v[(j+1)%3]);
3630  (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
3631  }
3632  }
3633 
3634  // 3. Do the red refinement.
3635  for (int i = 0; i < marked_el.Size(); i++)
3636  {
3637  RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
3638  }
3639 
3640  // 4. Do the green refinement (to get conforming mesh).
3641  int need_refinement;
3642  int edges_in_group, max_edges_in_group = 0;
3643  // edge_splittings identify how the shared edges have been split
3644  int **edge_splittings = new int*[GetNGroups()-1];
3645  for (int i = 0; i < GetNGroups()-1; i++)
3646  {
3647  edges_in_group = GroupNEdges(i+1);
3648  edge_splittings[i] = new int[edges_in_group];
3649  if (edges_in_group > max_edges_in_group)
3650  {
3651  max_edges_in_group = edges_in_group;
3652  }
3653  }
3654  int neighbor, *iBuf = new int[max_edges_in_group];
3655 
3656  Array<int> group_edges;
3657 
3658  MPI_Request request;
3659  MPI_Status status;
3660  Vertex V;
3661  V(2) = 0.0;
3662 
3663 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3664  int ref_loops_all = 0, ref_loops_par = 0;
3665 #endif
3666  do
3667  {
3668  need_refinement = 0;
3669  for (int i = 0; i < nedges; i++)
3670  {
3671  if (middle[i] != -1 && edge1[i] != -1)
3672  {
3673  need_refinement = 1;
3674  GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
3675  }
3676  }
3677 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3678  ref_loops_all++;
3679 #endif
3680 
3681  if (uniform_refinement)
3682  {
3683  continue;
3684  }
3685 
3686  // if the mesh is locally conforming start making it globally
3687  // conforming
3688  if (need_refinement == 0)
3689  {
3690 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3691  ref_loops_par++;
3692 #endif
3693  // MPI_Barrier(MyComm);
3694 
3695  // (a) send the type of interface splitting
3696  for (int i = 0; i < GetNGroups()-1; i++)
3697  {
3698  group_sedge.GetRow(i, group_edges);
3699  edges_in_group = group_edges.Size();
3700  // it is enough to communicate through the edges
3701  if (edges_in_group != 0)
3702  {
3703  for (int j = 0; j < edges_in_group; j++)
3704  {
3705  edge_splittings[i][j] =
3706  GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
3707  middle);
3708  }
3709  const int *nbs = gtopo.GetGroup(i+1);
3710  if (nbs[0] == 0)
3711  {
3712  neighbor = gtopo.GetNeighborRank(nbs[1]);
3713  }
3714  else
3715  {
3716  neighbor = gtopo.GetNeighborRank(nbs[0]);
3717  }
3718  MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
3719  neighbor, 0, MyComm, &request);
3720  }
3721  }
3722 
3723  // (b) receive the type of interface splitting
3724  for (int i = 0; i < GetNGroups()-1; i++)
3725  {
3726  group_sedge.GetRow(i, group_edges);
3727  edges_in_group = group_edges.Size();
3728  if (edges_in_group != 0)
3729  {
3730  const int *nbs = gtopo.GetGroup(i+1);
3731  if (nbs[0] == 0)
3732  {
3733  neighbor = gtopo.GetNeighborRank(nbs[1]);
3734  }
3735  else
3736  {
3737  neighbor = gtopo.GetNeighborRank(nbs[0]);
3738  }
3739  MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
3740  MPI_ANY_TAG, MyComm, &status);
3741 
3742  for (int j = 0; j < edges_in_group; j++)
3743  {
3744  if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
3745  {
3746  int *v = shared_edges[group_edges[j]]->GetVertices();
3747  int ii = v_to_v(v[0], v[1]);
3748 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3749  if (middle[ii] != -1)
3750  {
3751  mfem_error("ParMesh::LocalRefinement (triangles) : "
3752  "Oops!");
3753  }
3754 #endif
3755  need_refinement = 1;
3756  middle[ii] = NumOfVertices++;
3757  for (int c = 0; c < 2; c++)
3758  {
3759  V(c) = 0.5 * (vertices[v[0]](c) + vertices[v[1]](c));
3760  }
3761  vertices.Append(V);
3762  }
3763  }
3764  }
3765  }
3766 
3767  int nr = need_refinement;
3768  MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3769  }
3770  }
3771  while (need_refinement == 1);
3772 
3773 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3774  {
3775  int i = ref_loops_all;
3776  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3777  if (MyRank == 0)
3778  {
3779  mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3780  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3781  << '\n' << endl;
3782  }
3783  }
3784 #endif
3785 
3786  for (int i = 0; i < GetNGroups()-1; i++)
3787  {
3788  delete [] edge_splittings[i];
3789  }
3790  delete [] edge_splittings;
3791 
3792  delete [] iBuf;
3793 
3794  // 5. Update the boundary elements.
3795  int v1[2], v2[2], bisect, temp;
3796  temp = NumOfBdrElements;
3797  for (int i = 0; i < temp; i++)
3798  {
3799  int *v = boundary[i]->GetVertices();
3800  bisect = v_to_v(v[0], v[1]);
3801  if (middle[bisect] != -1)
3802  {
3803  // the element was refined (needs updating)
3804  if (boundary[i]->GetType() == Element::SEGMENT)
3805  {
3806  v1[0] = v[0]; v1[1] = middle[bisect];
3807  v2[0] = middle[bisect]; v2[1] = v[1];
3808 
3809  boundary[i]->SetVertices(v1);
3810  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
3811  }
3812  else
3813  {
3814  mfem_error("Only bisection of segment is implemented for bdr"
3815  " elem.");
3816  }
3817  }
3818  }
3819  NumOfBdrElements = boundary.Size();
3820 
3821  ResetLazyData();
3822 
3823  // 5a. Update the groups after refinement.
3824  RefineGroups(v_to_v, middle);
3825 
3826  // 6. Free the allocated memory.
3827  delete [] edge1;
3828  delete [] edge2;
3829  delete [] middle;
3830 
3831  if (el_to_edge != NULL)
3832  {
3834  GenerateFaces();
3835  }
3836  } // 'if (Dim == 2)'
3837 
3838  if (Dim == 1) // --------------------------------------------------------
3839  {
3840  int cne = NumOfElements, cnv = NumOfVertices;
3841  NumOfVertices += marked_el.Size();
3842  NumOfElements += marked_el.Size();
3843  vertices.SetSize(NumOfVertices);
3844  elements.SetSize(NumOfElements);
3846 
3847  for (int j = 0; j < marked_el.Size(); j++)
3848  {
3849  int i = marked_el[j];
3850  Segment *c_seg = (Segment *)elements[i];
3851  int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
3852  int new_v = cnv + j, new_e = cne + j;
3853  AverageVertices(vert, 2, new_v);
3854  elements[new_e] = new Segment(new_v, vert[1], attr);
3855  vert[1] = new_v;
3856 
3859  }
3860 
3861  static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
3863  UseExternalData(seg_children, 1, 2, 3);
3864 
3865  GenerateFaces();
3866  } // end of 'if (Dim == 1)'
3867 
3869  sequence++;
3870 
3871  UpdateNodes();
3872 
3873 #ifdef MFEM_DEBUG
3874  CheckElementOrientation(false);
3876 #endif
3877 }
3878 
3880  int nc_limit)
3881 {
3882  if (NURBSext)
3883  {
3884  MFEM_ABORT("NURBS meshes are not supported. Please project the "
3885  "NURBS to Nodes first with SetCurvature().");
3886  }
3887 
3888  if (!pncmesh)
3889  {
3890  MFEM_ABORT("Can't convert conforming ParMesh to nonconforming ParMesh "
3891  "(you need to initialize the ParMesh from a nonconforming "
3892  "serial Mesh)");
3893  }
3894 
3896 
3897  // NOTE: no check of !refinements.Size(), in parallel we would have to reduce
3898 
3899  // do the refinements
3901  pncmesh->Refine(refinements);
3902 
3903  if (nc_limit > 0)
3904  {
3905  pncmesh->LimitNCLevel(nc_limit);
3906  }
3907 
3908  // create a second mesh containing the finest elements from 'pncmesh'
3909  ParMesh* pmesh2 = new ParMesh(*pncmesh);
3910  pncmesh->OnMeshUpdated(pmesh2);
3911 
3912  attributes.Copy(pmesh2->attributes);
3914 
3915  // now swap the meshes, the second mesh will become the old coarse mesh
3916  // and this mesh will be the new fine mesh
3917  Mesh::Swap(*pmesh2, false);
3918 
3919  delete pmesh2; // NOTE: old face neighbors destroyed here
3920 
3922 
3924 
3926  sequence++;
3927 
3928  UpdateNodes();
3929 }
3930 
3932  double threshold, int nc_limit, int op)
3933 {
3934  MFEM_VERIFY(pncmesh, "Only supported for non-conforming meshes.");
3935  MFEM_VERIFY(!NURBSext, "Derefinement of NURBS meshes is not supported. "
3936  "Project the NURBS to Nodes first.");
3937 
3938  const Table &dt = pncmesh->GetDerefinementTable();
3939 
3940  pncmesh->SynchronizeDerefinementData(elem_error, dt);
3941 
3942  Array<int> level_ok;
3943  if (nc_limit > 0)
3944  {
3945  pncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
3946  }
3947 
3948  Array<int> derefs;
3949  for (int i = 0; i < dt.Size(); i++)
3950  {
3951  if (nc_limit > 0 && !level_ok[i]) { continue; }
3952 
3953  double error =
3954  AggregateError(elem_error, dt.GetRow(i), dt.RowSize(i), op);
3955 
3956  if (error < threshold) { derefs.Append(i); }
3957  }
3958 
3959  long long glob_size = ReduceInt(derefs.Size());
3960  if (!glob_size) { return false; }
3961 
3962  // Destroy face-neighbor data only when actually de-refining.
3964 
3965  pncmesh->Derefine(derefs);
3966 
3967  ParMesh* mesh2 = new ParMesh(*pncmesh);
3968  pncmesh->OnMeshUpdated(mesh2);
3969 
3970  attributes.Copy(mesh2->attributes);
3972 
3973  Mesh::Swap(*mesh2, false);
3974  delete mesh2;
3975 
3977 
3979 
3981  sequence++;
3982 
3983  UpdateNodes();
3984 
3985  return true;
3986 }
3987 
3988 
3990 {
3991  RebalanceImpl(NULL); // default SFC-based partition
3992 }
3993 
3994 void ParMesh::Rebalance(const Array<int> &partition)
3995 {
3996  RebalanceImpl(&partition);
3997 }
3998 
3999 void ParMesh::RebalanceImpl(const Array<int> *partition)
4000 {
4001  if (Conforming())
4002  {
4003  MFEM_ABORT("Load balancing is currently not supported for conforming"
4004  " meshes.");
4005  }
4006 
4007  if (Nodes)
4008  {
4009  // check that Nodes use a parallel FE space, so we can call UpdateNodes()
4010  MFEM_VERIFY(dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace())
4011  != NULL, "internal error");
4012  }
4013 
4015 
4016  pncmesh->Rebalance(partition);
4017 
4018  ParMesh* pmesh2 = new ParMesh(*pncmesh);
4019  pncmesh->OnMeshUpdated(pmesh2);
4020 
4021  attributes.Copy(pmesh2->attributes);
4023 
4024  Mesh::Swap(*pmesh2, false);
4025  delete pmesh2;
4026 
4028 
4030 
4032  sequence++;
4033 
4034  UpdateNodes();
4035 }
4036 
4037 void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
4038 {
4039  // Refine groups after LocalRefinement in 2D (triangle meshes)
4040 
4041  MFEM_ASSERT(Dim == 2 && meshgen == 1, "internal error");
4042 
4043  Array<int> group_verts, group_edges;
4044 
4045  // To update the groups after a refinement, we observe that:
4046  // - every (new and old) vertex, edge and face belongs to exactly one group
4047  // - the refinement does not create new groups
4048  // - a new vertex appears only as the middle of a refined edge
4049  // - a face can be refined 2, 3 or 4 times producing new edges and faces
4050 
4051  int *I_group_svert, *J_group_svert;
4052  int *I_group_sedge, *J_group_sedge;
4053 
4054  I_group_svert = Memory<int>(GetNGroups()+1);
4055  I_group_sedge = Memory<int>(GetNGroups()+1);
4056 
4057  I_group_svert[0] = I_group_svert[1] = 0;
4058  I_group_sedge[0] = I_group_sedge[1] = 0;
4059 
4060  // overestimate the size of the J arrays
4061  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4063  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4064 
4065  for (int group = 0; group < GetNGroups()-1; group++)
4066  {
4067  // Get the group shared objects
4068  group_svert.GetRow(group, group_verts);
4069  group_sedge.GetRow(group, group_edges);
4070 
4071  // Check which edges have been refined
4072  for (int i = 0; i < group_sedge.RowSize(group); i++)
4073  {
4074  int *v = shared_edges[group_edges[i]]->GetVertices();
4075  const int ind = middle[v_to_v(v[0], v[1])];
4076  if (ind != -1)
4077  {
4078  // add a vertex
4079  group_verts.Append(svert_lvert.Append(ind)-1);
4080  // update the edges
4081  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4082  shared_edges.Append(new Segment(v[1], ind, attr));
4083  group_edges.Append(sedge_ledge.Append(-1)-1);
4084  v[1] = ind;
4085  }
4086  }
4087 
4088  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4089  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4090 
4091  int *J;
4092  J = J_group_svert+I_group_svert[group];
4093  for (int i = 0; i < group_verts.Size(); i++)
4094  {
4095  J[i] = group_verts[i];
4096  }
4097  J = J_group_sedge+I_group_sedge[group];
4098  for (int i = 0; i < group_edges.Size(); i++)
4099  {
4100  J[i] = group_edges[i];
4101  }
4102  }
4103 
4104  FinalizeParTopo();
4105 
4106  group_svert.SetIJ(I_group_svert, J_group_svert);
4107  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4108 }
4109 
4110 void ParMesh::RefineGroups(int old_nv, const HashTable<Hashed2> &v_to_v)
4111 {
4112  // Refine groups after LocalRefinement in 3D (tetrahedral meshes)
4113 
4114  MFEM_ASSERT(Dim == 3 && meshgen == 1, "internal error");
4115 
4116  Array<int> group_verts, group_edges, group_trias;
4117 
4118  // To update the groups after a refinement, we observe that:
4119  // - every (new and old) vertex, edge and face belongs to exactly one group
4120  // - the refinement does not create new groups
4121  // - a new vertex appears only as the middle of a refined edge
4122  // - a face can be refined multiple times producing new edges and faces
4123 
4124  Array<Segment *> sedge_stack;
4125  Array<Vert3> sface_stack;
4126 
4127  Array<int> I_group_svert, J_group_svert;
4128  Array<int> I_group_sedge, J_group_sedge;
4129  Array<int> I_group_stria, J_group_stria;
4130 
4131  I_group_svert.SetSize(GetNGroups());
4132  I_group_sedge.SetSize(GetNGroups());
4133  I_group_stria.SetSize(GetNGroups());
4134 
4135  I_group_svert[0] = 0;
4136  I_group_sedge[0] = 0;
4137  I_group_stria[0] = 0;
4138 
4139  for (int group = 0; group < GetNGroups()-1; group++)
4140  {
4141  // Get the group shared objects
4142  group_svert.GetRow(group, group_verts);
4143  group_sedge.GetRow(group, group_edges);
4144  group_stria.GetRow(group, group_trias);
4145 
4146  // Check which edges have been refined
4147  for (int i = 0; i < group_sedge.RowSize(group); i++)
4148  {
4149  int *v = shared_edges[group_edges[i]]->GetVertices();
4150  int ind = v_to_v.FindId(v[0], v[1]);
4151  if (ind == -1) { continue; }
4152 
4153  // This shared edge is refined: walk the whole refinement tree
4154  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4155  do
4156  {
4157  ind += old_nv;
4158  // Add new shared vertex
4159  group_verts.Append(svert_lvert.Append(ind)-1);
4160  // Put the right sub-edge on top of the stack
4161  sedge_stack.Append(new Segment(ind, v[1], attr));
4162  // The left sub-edge replaces the original edge
4163  v[1] = ind;
4164  ind = v_to_v.FindId(v[0], ind);
4165  }
4166  while (ind != -1);
4167  // Process all edges in the edge stack
4168  do
4169  {
4170  Segment *se = sedge_stack.Last();
4171  v = se->GetVertices();
4172  ind = v_to_v.FindId(v[0], v[1]);
4173  if (ind == -1)
4174  {
4175  // The edge 'se' is not refined
4176  sedge_stack.DeleteLast();
4177  // Add new shared edge
4178  shared_edges.Append(se);
4179  group_edges.Append(sedge_ledge.Append(-1)-1);
4180  }
4181  else
4182  {
4183  // The edge 'se' is refined
4184  ind += old_nv;
4185  // Add new shared vertex
4186  group_verts.Append(svert_lvert.Append(ind)-1);
4187  // Put the left sub-edge on top of the stack
4188  sedge_stack.Append(new Segment(v[0], ind, attr));
4189  // The right sub-edge replaces the original edge
4190  v[0] = ind;
4191  }
4192  }
4193  while (sedge_stack.Size() > 0);
4194  }
4195 
4196  // Check which triangles have been refined
4197  for (int i = 0; i < group_stria.RowSize(group); i++)
4198  {
4199  int *v = shared_trias[group_trias[i]].v;
4200  int ind = v_to_v.FindId(v[0], v[1]);
4201  if (ind == -1) { continue; }
4202 
4203  // This shared face is refined: walk the whole refinement tree
4204  const int edge_attr = 1;
4205  do
4206  {
4207  ind += old_nv;
4208  // Add the refinement edge to the edge stack
4209  sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4210  // Put the right sub-triangle on top of the face stack
4211  sface_stack.Append(Vert3(v[1], v[2], ind));
4212  // The left sub-triangle replaces the original one
4213  v[1] = v[0]; v[0] = v[2]; v[2] = ind;
4214  ind = v_to_v.FindId(v[0], v[1]);
4215  }
4216  while (ind != -1);
4217  // Process all faces (triangles) in the face stack
4218  do
4219  {
4220  Vert3 &st = sface_stack.Last();
4221  v = st.v;
4222  ind = v_to_v.FindId(v[0], v[1]);
4223  if (ind == -1)
4224  {
4225  // The triangle 'st' is not refined
4226  // Add new shared face
4227  shared_trias.Append(st);
4228  group_trias.Append(sface_lface.Append(-1)-1);
4229  sface_stack.DeleteLast();
4230  }
4231  else
4232  {
4233  // The triangle 'st' is refined
4234  ind += old_nv;
4235  // Add the refinement edge to the edge stack
4236  sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4237  // Put the left sub-triangle on top of the face stack
4238  sface_stack.Append(Vert3(v[2], v[0], ind));
4239  // Note that the above Append() may invalidate 'v'
4240  v = sface_stack[sface_stack.Size()-2].v;
4241  // The right sub-triangle replaces the original one
4242  v[0] = v[1]; v[1] = v[2]; v[2] = ind;
4243  }
4244  }
4245  while (sface_stack.Size() > 0);
4246  // Process all edges in the edge stack (same code as above)
4247  do
4248  {
4249  Segment *se = sedge_stack.Last();
4250  v = se->GetVertices();
4251  ind = v_to_v.FindId(v[0], v[1]);
4252  if (ind == -1)
4253  {
4254  // The edge 'se' is not refined
4255  sedge_stack.DeleteLast();
4256  // Add new shared edge
4257  shared_edges.Append(se);
4258  group_edges.Append(sedge_ledge.Append(-1)-1);
4259  }
4260  else
4261  {
4262  // The edge 'se' is refined
4263  ind += old_nv;
4264  // Add new shared vertex
4265  group_verts.Append(svert_lvert.Append(ind)-1);
4266  // Put the left sub-edge on top of the stack
4267  sedge_stack.Append(new Segment(v[0], ind, edge_attr));
4268  // The right sub-edge replaces the original edge
4269  v[0] = ind;
4270  }
4271  }
4272  while (sedge_stack.Size() > 0);
4273  }
4274 
4275  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4276  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4277  I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4278 
4279  J_group_svert.Append(group_verts);
4280  J_group_sedge.Append(group_edges);
4281  J_group_stria.Append(group_trias);
4282  }
4283 
4284  FinalizeParTopo();
4285 
4286  group_svert.SetIJ(I_group_svert, J_group_svert);
4287  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4288  group_stria.SetIJ(I_group_stria, J_group_stria);
4289  I_group_svert.LoseData(); J_group_svert.LoseData();
4290  I_group_sedge.LoseData(); J_group_sedge.LoseData();
4291  I_group_stria.LoseData(); J_group_stria.LoseData();
4292 }
4293 
4295 {
4296  Array<int> sverts, sedges;
4297 
4298  int *I_group_svert, *J_group_svert;
4299  int *I_group_sedge, *J_group_sedge;
4300 
4301  I_group_svert = Memory<int>(GetNGroups());
4302  I_group_sedge = Memory<int>(GetNGroups());
4303 
4304  I_group_svert[0] = 0;
4305  I_group_sedge[0] = 0;
4306 
4307  // compute the size of the J arrays
4308  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4310  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4311 
4312  for (int group = 0; group < GetNGroups()-1; group++)
4313  {
4314  // Get the group shared objects
4315  group_svert.GetRow(group, sverts);
4316  group_sedge.GetRow(group, sedges);
4317 
4318  // Process all the edges
4319  for (int i = 0; i < group_sedge.RowSize(group); i++)
4320  {
4321  int *v = shared_edges[sedges[i]]->GetVertices();
4322  const int ind = old_nv + sedge_ledge[sedges[i]];
4323  // add a vertex
4324  sverts.Append(svert_lvert.Append(ind)-1);
4325  // update the edges
4326  const int attr = shared_edges[sedges[i]]->GetAttribute();
4327  shared_edges.Append(new Segment(v[1], ind, attr));
4328  sedges.Append(sedge_ledge.Append(-1)-1);
4329  v[1] = ind;
4330  }
4331 
4332  I_group_svert[group+1] = I_group_svert[group] + sverts.Size();
4333  I_group_sedge[group+1] = I_group_sedge[group] + sedges.Size();
4334 
4335  sverts.CopyTo(J_group_svert + I_group_svert[group]);
4336  sedges.CopyTo(J_group_sedge + I_group_sedge[group]);
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 }
4344 
4345 void ParMesh::UniformRefineGroups3D(int old_nv, int old_nedges,
4346  const DSTable &old_v_to_v,
4347  const STable3D &old_faces,
4348  Array<int> *f2qf)
4349 {
4350  // f2qf can be NULL if all faces are quads or there are no quad faces
4351 
4352  Array<int> group_verts, group_edges, group_trias, group_quads;
4353 
4354  int *I_group_svert, *J_group_svert;
4355  int *I_group_sedge, *J_group_sedge;
4356  int *I_group_stria, *J_group_stria;
4357  int *I_group_squad, *J_group_squad;
4358 
4359  I_group_svert = Memory<int>(GetNGroups());
4360  I_group_sedge = Memory<int>(GetNGroups());
4361  I_group_stria = Memory<int>(GetNGroups());
4362  I_group_squad = Memory<int>(GetNGroups());
4363 
4364  I_group_svert[0] = 0;
4365  I_group_sedge[0] = 0;
4366  I_group_stria[0] = 0;
4367  I_group_squad[0] = 0;
4368 
4369  // compute the size of the J arrays
4370  J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4373  J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections() +
4376  J_group_stria = Memory<int>(4*group_stria.Size_of_connections());
4377  J_group_squad = Memory<int>(4*group_squad.Size_of_connections());
4378 
4379  const int oface = old_nv + old_nedges;
4380 
4381  for (int group = 0; group < GetNGroups()-1; group++)
4382  {
4383  // Get the group shared objects
4384  group_svert.GetRow(group, group_verts);
4385  group_sedge.GetRow(group, group_edges);
4386  group_stria.GetRow(group, group_trias);
4387  group_squad.GetRow(group, group_quads);
4388 
4389  // Process the edges that have been refined
4390  for (int i = 0; i < group_sedge.RowSize(group); i++)
4391  {
4392  int *v = shared_edges[group_edges[i]]->GetVertices();
4393  const int ind = old_nv + old_v_to_v(v[0], v[1]);
4394  // add a vertex
4395  group_verts.Append(svert_lvert.Append(ind)-1);
4396  // update the edges
4397  const int attr = shared_edges[group_edges[i]]->GetAttribute();
4398  shared_edges.Append(new Segment(v[1], ind, attr));
4399  group_edges.Append(sedge_ledge.Append(-1)-1);
4400  v[1] = ind; // v[0] remains the same
4401  }
4402 
4403  // Process the triangles that have been refined
4404  for (int i = 0; i < group_stria.RowSize(group); i++)
4405  {
4406  int m[3];
4407  const int stria = group_trias[i];
4408  int *v = shared_trias[stria].v;
4409  // add the refinement edges
4410  m[0] = old_nv + old_v_to_v(v[0], v[1]);
4411  m[1] = old_nv + old_v_to_v(v[1], v[2]);
4412  m[2] = old_nv + old_v_to_v(v[2], v[0]);
4413  const int edge_attr = 1;
4414  shared_edges.Append(new Segment(m[0], m[1], edge_attr));
4415  group_edges.Append(sedge_ledge.Append(-1)-1);
4416  shared_edges.Append(new Segment(m[1], m[2], edge_attr));
4417  group_edges.Append(sedge_ledge.Append(-1)-1);
4418  shared_edges.Append(new Segment(m[0], m[2], edge_attr));
4419  group_edges.Append(sedge_ledge.Append(-1)-1);
4420  // update faces
4421  const int nst = shared_trias.Size();
4422  shared_trias.SetSize(nst+3);
4423  // The above SetSize() may invalidate 'v'
4424  v = shared_trias[stria].v;
4425  shared_trias[nst+0].Set(m[1],m[2],m[0]);
4426  shared_trias[nst+1].Set(m[0],v[1],m[1]);
4427  shared_trias[nst+2].Set(m[2],m[1],v[2]);
4428  v[1] = m[0]; v[2] = m[2]; // v[0] remains the same
4429  group_trias.Append(nst+0);
4430  group_trias.Append(nst+1);
4431  group_trias.Append(nst+2);
4432  // sface_lface is set later
4433  }
4434 
4435  // Process the quads that have been refined
4436  for (int i = 0; i < group_squad.RowSize(group); i++)
4437  {
4438  int m[5];
4439  const int squad = group_quads[i];
4440  int *v = shared_quads[squad].v;
4441  const int olf = old_faces(v[0], v[1], v[2], v[3]);
4442  // f2qf can be NULL if all faces are quads
4443  m[0] = oface + (f2qf ? (*f2qf)[olf] : olf);
4444  // add a vertex
4445  group_verts.Append(svert_lvert.Append(m[0])-1);
4446  // add the refinement edges
4447  m[1] = old_nv + old_v_to_v(v[0], v[1]);
4448  m[2] = old_nv + old_v_to_v(v[1], v[2]);
4449  m[3] = old_nv + old_v_to_v(v[2], v[3]);
4450  m[4] = old_nv + old_v_to_v(v[3], v[0]);
4451  const int edge_attr = 1;
4452  shared_edges.Append(new Segment(m[1], m[0], edge_attr));
4453  group_edges.Append(sedge_ledge.Append(-1)-1);
4454  shared_edges.Append(new Segment(m[2], m[0], edge_attr));
4455  group_edges.Append(sedge_ledge.Append(-1)-1);
4456  shared_edges.Append(new Segment(m[3], m[0], edge_attr));
4457  group_edges.Append(sedge_ledge.Append(-1)-1);
4458  shared_edges.Append(new Segment(m[4], m[0], edge_attr));
4459  group_edges.Append(sedge_ledge.Append(-1)-1);
4460  // update faces
4461  const int nsq = shared_quads.Size();
4462  shared_quads.SetSize(nsq+3);
4463  // The above SetSize() may invalidate 'v'
4464  v = shared_quads[squad].v;
4465  shared_quads[nsq+0].Set(m[1],v[1],m[2],m[0]);
4466  shared_quads[nsq+1].Set(m[0],m[2],v[2],m[3]);
4467  shared_quads[nsq+2].Set(m[4],m[0],m[3],v[3]);
4468  v[1] = m[1]; v[2] = m[0]; v[3] = m[4]; // v[0] remains the same
4469  group_quads.Append(nsq+0);
4470  group_quads.Append(nsq+1);
4471  group_quads.Append(nsq+2);
4472  // sface_lface is set later
4473  }
4474 
4475  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4476  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4477  I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4478  I_group_squad[group+1] = I_group_squad[group] + group_quads.Size();
4479 
4480  group_verts.CopyTo(J_group_svert + I_group_svert[group]);
4481  group_edges.CopyTo(J_group_sedge + I_group_sedge[group]);
4482  group_trias.CopyTo(J_group_stria + I_group_stria[group]);
4483  group_quads.CopyTo(J_group_squad + I_group_squad[group]);
4484  }
4485 
4486  FinalizeParTopo();
4487 
4488  group_svert.SetIJ(I_group_svert, J_group_svert);
4489  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4490  group_stria.SetIJ(I_group_stria, J_group_stria);
4491  group_squad.SetIJ(I_group_squad, J_group_squad);
4492 }
4493 
4495 {
4497 
4498  const int old_nv = NumOfVertices;
4499 
4500  // call Mesh::UniformRefinement2D so that it won't update the nodes
4501  {
4502  const bool update_nodes = false;
4503  Mesh::UniformRefinement2D_base(update_nodes);
4504  }
4505 
4506  // update the groups
4507  UniformRefineGroups2D(old_nv);
4508 
4509  UpdateNodes();
4510 
4511 #ifdef MFEM_DEBUG
4512  // If there are no Nodes, the orientation is checked in the call to
4513  // UniformRefinement2D_base() above.
4514  if (Nodes) { CheckElementOrientation(false); }
4515 #endif
4516 }
4517 
4519 {
4521 
4522  const int old_nv = NumOfVertices;
4523  const int old_nedges = NumOfEdges;
4524 
4525  DSTable v_to_v(NumOfVertices);
4526  GetVertexToVertexTable(v_to_v);
4527  STable3D *faces_tbl = GetFacesTable();
4528 
4529  // call Mesh::UniformRefinement3D_base so that it won't update the nodes
4530  Array<int> f2qf;
4531  {
4532  const bool update_nodes = false;
4533  UniformRefinement3D_base(&f2qf, &v_to_v, update_nodes);
4534  // Note: for meshes that have triangular faces, v_to_v is modified by the
4535  // above call to return different edge indices - this is used when
4536  // updating the groups. This is needed by ReorientTetMesh().
4537  }
4538 
4539  // update the groups
4540  UniformRefineGroups3D(old_nv, old_nedges, v_to_v, *faces_tbl,
4541  f2qf.Size() ? &f2qf : NULL);
4542  delete faces_tbl;
4543 
4544  UpdateNodes();
4545 }
4546 
4548 {
4549  if (MyRank == 0)
4550  {
4551  mfem::out << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
4552  }
4553 }
4554 
4555 void ParMesh::PrintXG(std::ostream &os) const
4556 {
4557  MFEM_ASSERT(Dim == spaceDim, "2D manifolds not supported");
4558  if (Dim == 3 && meshgen == 1)
4559  {
4560  int i, j, nv;
4561  const int *ind;
4562 
4563  os << "NETGEN_Neutral_Format\n";
4564  // print the vertices
4565  os << NumOfVertices << '\n';
4566  for (i = 0; i < NumOfVertices; i++)
4567  {
4568  for (j = 0; j < Dim; j++)
4569  {
4570  os << " " << vertices[i](j);
4571  }
4572  os << '\n';
4573  }
4574 
4575  // print the elements
4576  os << NumOfElements << '\n';
4577  for (i = 0; i < NumOfElements; i++)
4578  {
4579  nv = elements[i]->GetNVertices();
4580  ind = elements[i]->GetVertices();
4581  os << elements[i]->GetAttribute();
4582  for (j = 0; j < nv; j++)
4583  {
4584  os << " " << ind[j]+1;
4585  }
4586  os << '\n';
4587  }
4588 
4589  // print the boundary + shared faces information
4590  os << NumOfBdrElements + sface_lface.Size() << '\n';
4591  // boundary
4592  for (i = 0; i < NumOfBdrElements; i++)
4593  {
4594  nv = boundary[i]->GetNVertices();
4595  ind = boundary[i]->GetVertices();
4596  os << boundary[i]->GetAttribute();
4597  for (j = 0; j < nv; j++)
4598  {
4599  os << " " << ind[j]+1;
4600  }
4601  os << '\n';
4602  }
4603  // shared faces
4604  const int sf_attr =
4605  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4606  for (i = 0; i < shared_trias.Size(); i++)
4607  {
4608  ind = shared_trias[i].v;
4609  os << sf_attr;
4610  for (j = 0; j < 3; j++)
4611  {
4612  os << ' ' << ind[j]+1;
4613  }
4614  os << '\n';
4615  }
4616  // There are no quad shared faces
4617  }
4618 
4619  if (Dim == 3 && meshgen == 2)
4620  {
4621  int i, j, nv;
4622  const int *ind;
4623 
4624  os << "TrueGrid\n"
4625  << "1 " << NumOfVertices << " " << NumOfElements
4626  << " 0 0 0 0 0 0 0\n"
4627  << "0 0 0 1 0 0 0 0 0 0 0\n"
4628  << "0 0 " << NumOfBdrElements+sface_lface.Size()
4629  << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
4630  << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
4631  << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
4632 
4633  // print the vertices
4634  for (i = 0; i < NumOfVertices; i++)
4635  {
4636  os << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
4637  << " " << vertices[i](2) << " 0.0\n";
4638  }
4639 
4640  // print the elements
4641  for (i = 0; i < NumOfElements; i++)
4642  {
4643  nv = elements[i]->GetNVertices();
4644  ind = elements[i]->GetVertices();
4645  os << i+1 << " " << elements[i]->GetAttribute();
4646  for (j = 0; j < nv; j++)
4647  {
4648  os << " " << ind[j]+1;
4649  }
4650  os << '\n';
4651  }
4652 
4653  // print the boundary information
4654  for (i = 0; i < NumOfBdrElements; i++)
4655  {
4656  nv = boundary[i]->GetNVertices();
4657  ind = boundary[i]->GetVertices();
4658  os << boundary[i]->GetAttribute();
4659  for (j = 0; j < nv; j++)
4660  {
4661  os << " " << ind[j]+1;
4662  }
4663  os << " 1.0 1.0 1.0 1.0\n";
4664  }
4665 
4666  // print the shared faces information
4667  const int sf_attr =
4668  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4669  // There are no shared triangle faces
4670  for (i = 0; i < shared_quads.Size(); i++)
4671  {
4672  ind = shared_quads[i].v;
4673  os << sf_attr;
4674  for (j = 0; j < 4; j++)
4675  {
4676  os << ' ' << ind[j]+1;
4677  }
4678  os << " 1.0 1.0 1.0 1.0\n";
4679  }
4680  }
4681 
4682  if (Dim == 2)
4683  {
4684  int i, j, attr;
4685  Array<int> v;
4686 
4687  os << "areamesh2\n\n";
4688 
4689  // print the boundary + shared edges information
4690  os << NumOfBdrElements + shared_edges.Size() << '\n';
4691  // boundary
4692  for (i = 0; i < NumOfBdrElements; i++)
4693  {
4694  attr = boundary[i]->GetAttribute();
4695  boundary[i]->GetVertices(v);
4696  os << attr << " ";
4697  for (j = 0; j < v.Size(); j++)
4698  {
4699  os << v[j] + 1 << " ";
4700  }
4701  os << '\n';
4702  }
4703  // shared edges
4704  for (i = 0; i < shared_edges.Size(); i++)
4705  {
4706  attr = shared_edges[i]->GetAttribute();
4707  shared_edges[i]->GetVertices(v);
4708  os << attr << " ";
4709  for (j = 0; j < v.Size(); j++)
4710  {
4711  os << v[j] + 1 << " ";
4712  }
4713  os << '\n';
4714  }
4715 
4716  // print the elements
4717  os << NumOfElements << '\n';
4718  for (i = 0; i < NumOfElements; i++)
4719  {
4720  attr = elements[i]->GetAttribute();
4721  elements[i]->GetVertices(v);
4722 
4723  os << attr << " ";
4724  if ((j = GetElementType(i)) == Element::TRIANGLE)
4725  {
4726  os << 3 << " ";
4727  }
4728  else if (j == Element::QUADRILATERAL)
4729  {
4730  os << 4 << " ";
4731  }
4732  else if (j == Element::SEGMENT)
4733  {
4734  os << 2 << " ";
4735  }
4736  for (j = 0; j < v.Size(); j++)
4737  {
4738  os << v[j] + 1 << " ";
4739  }
4740  os << '\n';
4741  }
4742 
4743  // print the vertices
4744  os << NumOfVertices << '\n';
4745  for (i = 0; i < NumOfVertices; i++)
4746  {
4747  for (j = 0; j < Dim; j++)
4748  {
4749  os << vertices[i](j) << " ";
4750  }
4751  os << '\n';
4752  }
4753  }
4754  os.flush();
4755 }
4756 
4758 {
4759  // In 2D, this is a workaround for a CPU boundary rendering artifact. We need
4760  // to skip a shared master edge if one of its slaves has the same rank.
4761 
4762  const NCMesh::NCList &list = pncmesh->GetEdgeList();
4763  for (int i = master.slaves_begin; i < master.slaves_end; i++)
4764  {
4765  if (!pncmesh->IsGhost(1, list.slaves[i].index)) { return true; }
4766  }
4767  return false;
4768 }
4769 
4770 void ParMesh::Print(std::ostream &os) const
4771 {
4772  int shared_bdr_attr;
4773  Array<int> nc_shared_faces;
4774 
4775  if (NURBSext)
4776  {
4777  Printer(os); // does not print shared boundary
4778  return;
4779  }
4780 
4781  const Array<int>* s2l_face;
4782  if (!pncmesh)
4783  {
4784  s2l_face = ((Dim == 1) ? &svert_lvert :
4785  ((Dim == 2) ? &sedge_ledge : &sface_lface));
4786  }
4787  else
4788  {
4789  s2l_face = &nc_shared_faces;
4790  if (Dim >= 2)
4791  {
4792  // get a list of all shared non-ghost faces
4793  const NCMesh::NCList& sfaces =
4794  (Dim == 3) ? pncmesh->GetSharedFaces() : pncmesh->GetSharedEdges();
4795  const int nfaces = GetNumFaces();
4796  for (int i = 0; i < sfaces.conforming.Size(); i++)
4797  {
4798  int index = sfaces.conforming[i].index;
4799  if (index < nfaces) { nc_shared_faces.Append(index); }
4800  }
4801  for (int i = 0; i < sfaces.masters.Size(); i++)
4802  {
4803  if (Dim == 2 && WantSkipSharedMaster(sfaces.masters[i])) { continue; }
4804  int index = sfaces.masters[i].index;
4805  if (index < nfaces) { nc_shared_faces.Append(index); }
4806  }
4807  for (int i = 0; i < sfaces.slaves.Size(); i++)
4808  {
4809  int index = sfaces.slaves[i].index;
4810  if (index < nfaces) { nc_shared_faces.Append(index); }
4811  }
4812  }
4813  }
4814 
4815  os << "MFEM mesh v1.0\n";
4816 
4817  // optional
4818  os <<
4819  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4820  "# POINT = 0\n"
4821  "# SEGMENT = 1\n"
4822  "# TRIANGLE = 2\n"
4823  "# SQUARE = 3\n"
4824  "# TETRAHEDRON = 4\n"
4825  "# CUBE = 5\n"
4826  "# PRISM = 6\n"
4827  "#\n";
4828 
4829  os << "\ndimension\n" << Dim
4830  << "\n\nelements\n" << NumOfElements << '\n';
4831  for (int i = 0; i < NumOfElements; i++)
4832  {
4833  PrintElement(elements[i], os);
4834  }
4835 
4836  int num_bdr_elems = NumOfBdrElements;
4837  if (print_shared && Dim > 1)
4838  {
4839  num_bdr_elems += s2l_face->Size();
4840  }
4841  os << "\nboundary\n" << num_bdr_elems << '\n';
4842  for (int i = 0; i < NumOfBdrElements; i++)
4843  {
4844  PrintElement(boundary[i], os);
4845  }
4846 
4847  if (print_shared && Dim > 1)
4848  {
4849  if (bdr_attributes.Size())
4850  {
4851  shared_bdr_attr = bdr_attributes.Max() + MyRank + 1;
4852  }
4853  else
4854  {
4855  shared_bdr_attr = MyRank + 1;
4856  }
4857  for (int i = 0; i < s2l_face->Size(); i++)
4858  {
4859  // Modify the attributes of the faces (not used otherwise?)
4860  faces[(*s2l_face)[i]]->SetAttribute(shared_bdr_attr);
4861  PrintElement(faces[(*s2l_face)[i]], os);
4862  }
4863  }
4864  os << "\nvertices\n" << NumOfVertices << '\n';
4865  if (Nodes == NULL)
4866  {
4867  os << spaceDim << '\n';
4868  for (int i = 0; i < NumOfVertices; i++)
4869  {
4870  os << vertices[i](0);
4871  for (int j = 1; j < spaceDim; j++)
4872  {
4873  os << ' ' << vertices[i](j);
4874  }
4875  os << '\n';
4876  }
4877  os.flush();
4878  }
4879  else
4880  {
4881  os << "\nnodes\n";
4882  Nodes->Save(os);
4883  }
4884 }
4885 
4886 void ParMesh::Save(const char *fname, int precision) const
4887 {
4888  ostringstream fname_with_suffix;
4889  fname_with_suffix << fname << "." << setfill('0') << setw(6) << MyRank;
4890  ofstream ofs(fname_with_suffix.str().c_str());
4891  ofs.precision(precision);
4892  Print(ofs);
4893 }
4894 
4895 #ifdef MFEM_USE_ADIOS2
4897 {
4898  Mesh::Print(os);
4899 }
4900 #endif
4901 
4902 static void dump_element(const Element* elem, Array<int> &data)
4903 {
4904  data.Append(elem->GetGeometryType());
4905 
4906  int nv = elem->GetNVertices();
4907  const int *v = elem->GetVertices();
4908  for (int i = 0; i < nv; i++)
4909  {
4910  data.Append(v[i]);
4911  }
4912 }
4913 
4914 void ParMesh::PrintAsOne(std::ostream &os) const
4915 {
4916  int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
4917  const int *v;
4918  MPI_Status status;
4919  Array<double> vert;
4920  Array<int> ints;
4921 
4922  if (MyRank == 0)
4923  {
4924  os << "MFEM mesh v1.0\n";
4925 
4926  // optional
4927  os <<
4928  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4929  "# POINT = 0\n"
4930  "# SEGMENT = 1\n"
4931  "# TRIANGLE = 2\n"
4932  "# SQUARE = 3\n"
4933  "# TETRAHEDRON = 4\n"
4934  "# CUBE = 5\n"
4935  "# PRISM = 6\n"
4936  "#\n";
4937 
4938  os << "\ndimension\n" << Dim;
4939  }
4940 
4941  nv = NumOfElements;
4942  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
4943  if (MyRank == 0)
4944  {
4945  os << "\n\nelements\n" << ne << '\n';
4946  for (i = 0; i < NumOfElements; i++)
4947  {
4948  // processor number + 1 as attribute and geometry type
4949  os << 1 << ' ' << elements[i]->GetGeometryType();
4950  // vertices
4951  nv = elements[i]->GetNVertices();
4952  v = elements[i]->GetVertices();
4953  for (j = 0; j < nv; j++)
4954  {
4955  os << ' ' << v[j];
4956  }
4957  os << '\n';
4958  }
4959  vc = NumOfVertices;
4960  for (p = 1; p < NRanks; p++)
4961  {
4962  MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
4963  ints.SetSize(ne);
4964  if (ne)
4965  {
4966  MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
4967  }
4968  for (i = 0; i < ne; )
4969  {
4970  // processor number + 1 as attribute and geometry type
4971  os << p+1 << ' ' << ints[i];
4972  // vertices
4973  k = Geometries.GetVertices(ints[i++])->GetNPoints();
4974  for (j = 0; j < k; j++)
4975  {
4976  os << ' ' << vc + ints[i++];
4977  }
4978  os << '\n';
4979  }
4980  vc += nv;
4981  }
4982  }
4983  else
4984  {
4985  // for each element send its geometry type and its vertices
4986  ne = 0;
4987  for (i = 0; i < NumOfElements; i++)
4988  {
4989  ne += 1 + elements[i]->GetNVertices();
4990  }
4991  nv = NumOfVertices;
4992  MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
4993 
4994  ints.Reserve(ne);
4995  ints.SetSize(0);
4996  for (i = 0; i < NumOfElements; i++)
4997  {
4998  dump_element(elements[i], ints);
4999  }
5000  MFEM_ASSERT(ints.Size() == ne, "");
5001  if (ne)
5002  {
5003  MPI_Send(&ints[0], ne, MPI_INT, 0, 445, MyComm);
5004  }
5005  }
5006 
5007  // boundary + shared boundary
5008  ne = NumOfBdrElements;
5009  if (!pncmesh)
5010  {
5011  ne += GetNSharedFaces();
5012  }
5013  else if (Dim > 1)
5014  {
5015  const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5016  ne += list.conforming.Size() + list.masters.Size() + list.slaves.Size();
5017  // In addition to the number returned by GetNSharedFaces(), include the
5018  // the master shared faces as well.
5019  }
5020  ints.Reserve(ne * (1 + 2*(Dim-1))); // just an upper bound
5021  ints.SetSize(0);
5022 
5023  // for each boundary and shared boundary element send its geometry type
5024  // and its vertices
5025  ne = 0;
5026  for (i = j = 0; i < NumOfBdrElements; i++)
5027  {
5028  dump_element(boundary[i], ints); ne++;
5029  }
5030  if (!pncmesh)
5031  {
5032  switch (Dim)
5033  {
5034  case 1:
5035  for (i = 0; i < svert_lvert.Size(); i++)
5036  {
5037  ints.Append(Geometry::POINT);
5038  ints.Append(svert_lvert[i]);
5039  ne++;
5040  }
5041  break;
5042 
5043  case 2:
5044  for (i = 0; i < shared_edges.Size(); i++)
5045  {
5046  dump_element(shared_edges[i], ints); ne++;
5047  }
5048  break;
5049 
5050  case 3:
5051  for (i = 0; i < shared_trias.Size(); i++)
5052  {
5053  ints.Append(Geometry::TRIANGLE);
5054  ints.Append(shared_trias[i].v, 3);
5055  ne++;
5056  }
5057  for (i = 0; i < shared_quads.Size(); i++)
5058  {
5059  ints.Append(Geometry::SQUARE);
5060  ints.Append(shared_quads[i].v, 4);
5061  ne++;
5062  }
5063  break;
5064 
5065  default:
5066  MFEM_ABORT("invalid dimension: " << Dim);
5067  }
5068  }
5069  else if (Dim > 1)
5070  {
5071  const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5072  const int nfaces = GetNumFaces();
5073  for (i = 0; i < list.conforming.Size(); i++)
5074  {
5075  int index = list.conforming[i].index;
5076  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5077  }
5078  for (i = 0; i < list.masters.Size(); i++)
5079  {
5080  int index = list.masters[i].index;
5081  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5082  }
5083  for (i = 0; i < list.slaves.Size(); i++)
5084  {
5085  int index = list.slaves[i].index;
5086  if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5087  }
5088  }
5089 
5090  MPI_Reduce(&ne, &k, 1, MPI_INT, MPI_SUM, 0, MyComm);
5091  if (MyRank == 0)
5092  {
5093  os << "\nboundary\n" << k << '\n';
5094  vc = 0;
5095  for (p = 0; p < NRanks; p++)
5096  {
5097  if (p)
5098  {
5099  MPI_Recv(nv_ne, 2, MPI_INT, p, 446, MyComm, &status);
5100  ints.SetSize(ne);
5101  if (ne)
5102  {
5103  MPI_Recv(ints.GetData(), ne, MPI_INT, p, 447, MyComm, &status);
5104  }
5105  }
5106  else
5107  {
5108  ne = ints.Size();
5109  nv = NumOfVertices;
5110  }
5111  for (i = 0; i < ne; )
5112  {
5113  // processor number + 1 as bdr. attr. and bdr. geometry type
5114  os << p+1 << ' ' << ints[i];
5115  k = Geometries.NumVerts[ints[i++]];
5116  // vertices
5117  for (j = 0; j < k; j++)
5118  {
5119  os << ' ' << vc + ints[i++];
5120  }
5121  os << '\n';
5122  }
5123  vc += nv;
5124  }
5125  }
5126  else
5127  {
5128  nv = NumOfVertices;
5129  ne = ints.Size();
5130  MPI_Send(nv_ne, 2, MPI_INT, 0, 446, MyComm);
5131  if (ne)
5132  {
5133  MPI_Send(ints.GetData(), ne, MPI_INT, 0, 447, MyComm);
5134  }
5135  }
5136 
5137  // vertices / nodes
5138  MPI_Reduce(&NumOfVertices, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5139  if (MyRank == 0)
5140  {
5141  os << "\nvertices\n" << nv << '\n';
5142  }
5143  if (Nodes == NULL)
5144  {
5145  if (MyRank == 0)
5146  {
5147  os << spaceDim << '\n';
5148  for (i = 0; i < NumOfVertices; i++)
5149  {
5150  os << vertices[i](0);
5151  for (j = 1; j < spaceDim; j++)
5152  {
5153  os << ' ' << vertices[i](j);
5154  }
5155  os << '\n';
5156  }
5157  for (p = 1; p < NRanks; p++)
5158  {
5159  MPI_Recv(&nv, 1, MPI_INT, p, 448, MyComm, &status);
5160  vert.SetSize(nv*spaceDim);
5161  if (nv)
5162  {
5163  MPI_Recv(&vert[0], nv*spaceDim, MPI_DOUBLE, p, 449, MyComm, &status);
5164  }
5165  for (i = 0; i < nv; i++)
5166  {
5167  os << vert[i*spaceDim];
5168  for (j = 1; j < spaceDim; j++)
5169  {
5170  os << ' ' << vert[i*spaceDim+j];
5171  }
5172  os << '\n';
5173  }
5174  }
5175  os.flush();
5176  }
5177  else
5178  {
5179  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 448, MyComm);
5181  for (i = 0; i < NumOfVertices; i++)
5182  {
5183  for (j = 0; j < spaceDim; j++)
5184  {
5185  vert[i*spaceDim+j] = vertices[i](j);
5186  }
5187  }
5188  if (NumOfVertices)
5189  {
5190  MPI_Send(&vert[0], NumOfVertices*spaceDim, MPI_DOUBLE, 0, 449, MyComm);
5191  }
5192  }
5193  }
5194  else
5195  {
5196  if (MyRank == 0)
5197  {
5198  os << "\nnodes\n";
5199  }
5200  ParGridFunction *pnodes = dynamic_cast<ParGridFunction *>(Nodes);
5201  if (pnodes)
5202  {
5203  pnodes->SaveAsOne(os);
5204  }
5205  else
5206  {
5207  ParFiniteElementSpace *pfes =
5208  dynamic_cast<ParFiniteElementSpace *>(Nodes->FESpace());
5209  if (pfes)
5210  {
5211  // create a wrapper ParGridFunction
5212  ParGridFunction ParNodes(pfes, Nodes);
5213  ParNodes.SaveAsOne(os);
5214  }
5215  else
5216  {
5217  mfem_error("ParMesh::PrintAsOne : Nodes have no parallel info!");
5218  }
5219  }
5220  }
5221 }
5222 
5223 void ParMesh::PrintAsSerial(std::ostream &os) const
5224 {
5225  int save_rank = 0;
5226  Mesh serialmesh = GetSerialMesh(save_rank);
5227  if (MyRank == save_rank)
5228  {
5229  serialmesh.Printer(os);
5230  }
5231  MPI_Barrier(MyComm);
5232 }
5233 
5234 Mesh ParMesh::GetSerialMesh(int save_rank) const
5235 {
5236  if (pncmesh || NURBSext)
5237  {
5238  MFEM_ABORT("Nonconforming meshes and NURBS meshes are not yet supported.");
5239  }
5240 
5241  // Define linear H1 space for vertex numbering
5242  H1_FECollection fec_linear(1, Dim);
5243  ParMesh *pm = const_cast<ParMesh *>(this);
5244  ParFiniteElementSpace pfespace_linear(pm, &fec_linear);
5245 
5246  long long ne_glob_l = GetGlobalNE(); // needs to be called by all ranks
5247  MFEM_VERIFY(int(ne_glob_l) == ne_glob_l,
5248  "overflow in the number of elements!");
5249  int ne_glob = (save_rank == MyRank) ? int(ne_glob_l) : 0;
5250 
5251  long long nvertices = pfespace_linear.GetTrueVSize();
5252  long long nvertices_glob_l = 0;
5253  MPI_Reduce(&nvertices, &nvertices_glob_l, 1, MPI_LONG_LONG, MPI_SUM,
5254  save_rank, MyComm);
5255  int nvertices_glob = int(nvertices_glob_l);
5256  MFEM_VERIFY(nvertices_glob == nvertices_glob_l,
5257  "overflow in the number of vertices!");
5258 
5259  long long nbe = NumOfBdrElements;
5260  long long nbe_glob_l = 0;
5261  MPI_Reduce(&nbe, &nbe_glob_l, 1, MPI_LONG_LONG, MPI_SUM, save_rank, MyComm);
5262  int nbe_glob = int(nbe_glob_l);
5263  MFEM_VERIFY(nbe_glob == nbe_glob_l,
5264  "overflow in the number of boundary elements!");
5265 
5266  // On ranks other than save_rank, the *_glob variables are 0, so the serial
5267  // mesh is empty.
5268  Mesh serialmesh(Dim, nvertices_glob, ne_glob, nbe_glob, spaceDim);
5269 
5270  int n_send_recv;
5271  MPI_Status status;
5272  Array<double> vert;
5273  Array<int> ints, dofs;
5274 
5275  // First set the connectivity of serial mesh using the True Dofs from
5276  // the linear H1 space.
5277  if (MyRank == save_rank)
5278  {
5279  for (int e = 0; e < NumOfElements; e++)
5280  {
5281  const int attr = elements[e]->GetAttribute();
5282  const int geom_type = elements[e]->GetGeometryType();
5283  pfespace_linear.GetElementDofs(e, dofs);
5284  for (int j = 0; j < dofs.Size(); j++)
5285  {
5286  dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
5287  }
5288  Element *elem = serialmesh.NewElement(geom_type);
5289  elem->SetAttribute(attr);
5290  elem->SetVertices(dofs);
5291  serialmesh.AddElement(elem);
5292  }
5293 
5294  for (int p = 0; p < NRanks; p++)
5295  {
5296  if (p == save_rank) { continue; }
5297  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 444, MyComm, &status);
5298  ints.SetSize(n_send_recv);
5299  if (n_send_recv)
5300  {
5301  MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 445, MyComm, &status);
5302  }
5303  for (int i = 0; i < n_send_recv; )
5304  {
5305  int attr = ints[i++];
5306  int geom_type = ints[i++];
5307  Element *elem = serialmesh.NewElement(geom_type);
5308  elem->SetAttribute(attr);
5309  elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
5310  serialmesh.AddElement(elem);
5311  }
5312  }
5313  }
5314  else
5315  {
5316  n_send_recv = 0;
5317  for (int e = 0; e < NumOfElements; e++)
5318  {
5319  n_send_recv += 2 + elements[e]->GetNVertices();
5320  }
5321  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 444, MyComm);
5322  ints.Reserve(n_send_recv);
5323  ints.SetSize(0);
5324  for (int e = 0; e < NumOfElements; e++)
5325  {
5326  const int attr = elements[e]->GetAttribute();
5327  const int geom_type = elements[e]->GetGeometryType();;
5328  ints.Append(attr);
5329  ints.Append(geom_type);
5330  pfespace_linear.GetElementDofs(e, dofs);
5331  for (int j = 0; j < dofs.Size(); j++)
5332  {
5333  ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
5334  }
5335  }
5336  if (n_send_recv)
5337  {
5338  MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 445, MyComm);
5339  }
5340  }
5341 
5342  // Write out boundary elements
5343  if (MyRank == save_rank)
5344  {
5345  for (int e = 0; e < NumOfBdrElements; e++)
5346  {
5347  const int attr = boundary[e]->GetAttribute();
5348  const int geom_type = boundary[e]->GetGeometryType();
5349  pfespace_linear.GetBdrElementDofs(e, dofs);
5350  for (int j = 0; j < dofs.Size(); j++)
5351  {
5352  dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
5353  }
5354  Element *elem = serialmesh.NewElement(geom_type);
5355  elem->SetAttribute(attr);
5356  elem->SetVertices(dofs);
5357  serialmesh.AddBdrElement(elem);
5358  }
5359 
5360  for (int p = 0; p < NRanks; p++)
5361  {
5362  if (p == save_rank) { continue; }
5363  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 446, MyComm, &status);
5364  ints.SetSize(n_send_recv);
5365  if (n_send_recv)
5366  {
5367  MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 447, MyComm, &status);
5368  }
5369  for (int i = 0; i < n_send_recv; )
5370  {
5371  int attr = ints[i++];
5372  int geom_type = ints[i++];
5373  Element *elem = serialmesh.NewElement(geom_type);
5374  elem->SetAttribute(attr);
5375  elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
5376  serialmesh.AddBdrElement(elem);
5377  }
5378  }
5379  } // MyRank == save_rank
5380  else
5381  {
5382  n_send_recv = 0;
5383  for (int e = 0; e < NumOfBdrElements; e++)
5384  {
5385  n_send_recv += 2 + GetBdrElement(e)->GetNVertices();
5386  }
5387  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 446, MyComm);
5388  ints.Reserve(n_send_recv);
5389  ints.SetSize(0);
5390  for (int e = 0; e < NumOfBdrElements; e++)
5391  {
5392  const int attr = boundary[e]->GetAttribute();
5393  const int geom_type = boundary[e]->GetGeometryType();
5394  ints.Append(attr);
5395  ints.Append(geom_type);
5396  pfespace_linear.GetBdrElementDofs(e, dofs);
5397  for (int j = 0; j < dofs.Size(); j++)
5398  {
5399  ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
5400  }
5401  }
5402  if (n_send_recv)
5403  {
5404  MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 447, MyComm);
5405  }
5406  } // MyRank != save_rank
5407 
5408  if (MyRank == save_rank)
5409  {
5410  for (int v = 0; v < nvertices_glob; v++)
5411  {
5412  serialmesh.AddVertex(0.0); // all other coordinates are 0 by default
5413  }
5414  serialmesh.FinalizeTopology();
5415  }
5416 
5417  // From each processor, we send element-wise vertex/dof locations and
5418  // overwrite the vertex/dof locations of the serial mesh.
5419  if (MyRank == save_rank && Nodes)
5420  {
5421  FiniteElementSpace *fespace_serial = NULL;
5422  // Duplicate the FE collection to make sure the serial mesh is completely
5423  // independent of the parallel mesh:
5424  auto fec_serial = FiniteElementCollection::New(
5425  GetNodalFESpace()->FEColl()->Name());
5426  fespace_serial = new FiniteElementSpace(&serialmesh,
5427  fec_serial,
5428  spaceDim,
5429  GetNodalFESpace()->GetOrdering());
5430  serialmesh.SetNodalFESpace(fespace_serial);
5431  serialmesh.GetNodes()->MakeOwner(fec_serial);
5432  // The serial mesh owns its Nodes and they, in turn, own fec_serial and
5433  // fespace_serial.
5434  }
5435 
5436  int elem_count = 0; // To keep track of element count in serial mesh
5437  if (MyRank == save_rank)
5438  {
5439  Vector nodeloc;
5440  Array<int> ints_serial;
5441  for (int e = 0; e < NumOfElements; e++)
5442  {
5443  if (Nodes)
5444  {
5445  Nodes->GetElementDofValues(e, nodeloc);
5446  serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
5447  serialmesh.GetNodes()->SetSubVector(dofs, nodeloc);
5448  }
5449  else
5450  {
5451  GetElementVertices(e, ints);
5452  serialmesh.GetElementVertices(elem_count++, ints_serial);
5453  for (int i = 0; i < ints.Size(); i++)
5454  {
5455  const double *vdata = GetVertex(ints[i]);
5456  double *vdata_serial = serialmesh.GetVertex(ints_serial[i]);
5457  for (int d = 0; d < spaceDim; d++)
5458  {
5459  vdata_serial[d] = vdata[d];
5460  }
5461  }
5462  }
5463  }
5464 
5465  for (int p = 0; p < NRanks; p++)
5466  {
5467  if (p == save_rank) { continue; }
5468  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, MyComm, &status);
5469  vert.SetSize(n_send_recv);
5470  if (n_send_recv)
5471  {
5472  MPI_Recv(&vert[0], n_send_recv, MPI_DOUBLE, p, 449, MyComm, &status);
5473  }
5474  for (int i = 0; i < n_send_recv; )
5475  {
5476  if (Nodes)
5477  {
5478  serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
5479  serialmesh.GetNodes()->SetSubVector(dofs, &vert[i]);
5480  i += dofs.Size();
5481  }
5482  else
5483  {
5484  serialmesh.GetElementVertices(elem_count++, ints_serial);
5485  for (int j = 0; j < ints_serial.Size(); j++)
5486  {
5487  double *vdata_serial = serialmesh.GetVertex(ints_serial[j]);
5488  for (int d = 0; d < spaceDim; d++)
5489  {
5490  vdata_serial[d] = vert[i++];
5491  }
5492  }
5493  }
5494  }
5495  }
5496  } // MyRank == save_rank
5497  else
5498  {
5499  n_send_recv = 0;
5500  Vector nodeloc;
5501  for (int e = 0; e < NumOfElements; e++)
5502  {
5503  if (Nodes)
5504  {
5505  const FiniteElement *fe = Nodes->FESpace()->GetFE(e);
5506  n_send_recv += spaceDim*fe->GetDof();
5507  }
5508  else
5509  {
5510  n_send_recv += elements[e]->GetNVertices()*spaceDim;
5511  }
5512  }
5513  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, MyComm);
5514  vert.Reserve(n_send_recv);
5515  vert.SetSize(0);
5516  for (int e = 0; e < NumOfElements; e++)
5517  {
5518  if (Nodes)
5519  {
5520  Nodes->GetElementDofValues(e, nodeloc);
5521  for (int j = 0; j < nodeloc.Size(); j++)
5522  {
5523  vert.Append(nodeloc(j));
5524  }
5525  }
5526  else
5527  {
5528  GetElementVertices(e, ints);
5529  for (int i = 0; i < ints.Size(); i++)
5530  {
5531  const double *vdata = GetVertex(ints[i]);
5532  for (int d = 0; d < spaceDim; d++)
5533  {
5534  vert.Append(vdata[d]);
5535  }
5536  }
5537  }
5538  }
5539  if (n_send_recv)
5540  {
5541  MPI_Send(&vert[0], n_send_recv, MPI_DOUBLE, save_rank, 449, MyComm);
5542  }
5543  }
5544 
5545  MPI_Barrier(MyComm);
5546  return serialmesh;
5547 }
5548 
5549 void ParMesh::SaveAsOne(const char *fname, int precision) const
5550 {
5551  ofstream ofs;
5552  if (MyRank == 0)
5553  {
5554  ofs.open(fname);
5555  ofs.precision(precision);
5556  }
5557  PrintAsOne(ofs);
5558 }
5559 
5560 void ParMesh::PrintAsOneXG(std::ostream &os)
5561 {
5562  MFEM_ASSERT(Dim == spaceDim, "2D Manifolds not supported.");
5563  if (Dim == 3 && meshgen == 1)
5564  {
5565  int i, j, k, nv, ne, p;
5566  const int *ind, *v;
5567  MPI_Status status;
5568  Array<double> vert;
5569  Array<int> ints;
5570 
5571  if (MyRank == 0)
5572  {
5573  os << "NETGEN_Neutral_Format\n";
5574  // print the vertices
5575  ne = NumOfVertices;
5576  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5577  os << nv << '\n';
5578  for (i = 0; i < NumOfVertices; i++)
5579  {
5580  for (j = 0; j < Dim; j++)
5581  {
5582  os << " " << vertices[i](j);
5583  }
5584  os << '\n';
5585  }
5586  for (p = 1; p < NRanks; p++)
5587  {
5588  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5589  vert.SetSize(Dim*nv);
5590  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
5591  for (i = 0; i < nv; i++)
5592  {
5593  for (j = 0; j < Dim; j++)
5594  {
5595  os << " " << vert[Dim*i+j];
5596  }
5597  os << '\n';
5598  }
5599  }
5600 
5601  // print the elements
5602  nv = NumOfElements;
5603  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5604  os << ne << '\n';
5605  for (i = 0; i < NumOfElements; i++)
5606  {
5607  nv = elements[i]->GetNVertices();
5608  ind = elements[i]->GetVertices();
5609  os << 1;
5610  for (j = 0; j < nv; j++)
5611  {
5612  os << " " << ind[j]+1;
5613  }
5614  os << '\n';
5615  }
5616  k = NumOfVertices;
5617  for (p = 1; p < NRanks; p++)
5618  {
5619  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5620  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5621  ints.SetSize(4*ne);
5622  MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
5623  for (i = 0; i < ne; i++)
5624  {
5625  os << p+1;
5626  for (j = 0; j < 4; j++)
5627  {
5628  os << " " << k+ints[i*4+j]+1;
5629  }
5630  os << '\n';
5631  }
5632  k += nv;
5633  }
5634  // print the boundary + shared faces information
5636  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5637  os << ne << '\n';
5638  // boundary
5639  for (i = 0; i < NumOfBdrElements; i++)
5640  {
5641  nv = boundary[i]->GetNVertices();
5642  ind = boundary[i]->GetVertices();
5643  os << 1;
5644  for (j = 0; j < nv; j++)
5645  {
5646  os << " " << ind[j]+1;
5647  }
5648  os << '\n';
5649  }
5650  // shared faces
5651  const int sf_attr =
5652  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
5653  for (i = 0; i < shared_trias.Size(); i++)
5654  {
5655  ind = shared_trias[i].v;
5656  os << sf_attr;
5657  for (j = 0; j < 3; j++)
5658  {
5659  os << ' ' << ind[j]+1;
5660  }
5661  os << '\n';
5662  }
5663  // There are no quad shared faces
5664  k = NumOfVertices;
5665  for (p = 1; p < NRanks; p++)
5666  {
5667  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5668  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5669  ints.SetSize(3*ne);
5670  MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
5671  for (i = 0; i < ne; i++)
5672  {
5673  os << p+1;
5674  for (j = 0; j < 3; j++)
5675  {
5676  os << ' ' << k+ints[i*3+j]+1;
5677  }
5678  os << '\n';
5679  }
5680  k += nv;
5681  }
5682  }
5683  else
5684  {
5685  ne = NumOfVertices;
5686  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5687  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5688  vert.SetSize(Dim*NumOfVertices);
5689  for (i = 0; i < NumOfVertices; i++)
5690  for (j = 0; j < Dim; j++)
5691  {
5692  vert[Dim*i+j] = vertices[i](j);
5693  }
5694  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
5695  0, 445, MyComm);
5696  // elements
5697  ne = NumOfElements;
5698  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5699  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5700  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
5701  ints.SetSize(NumOfElements*4);
5702  for (i = 0; i < NumOfElements; i++)
5703  {
5704  v = elements[i]->GetVertices();
5705  for (j = 0; j < 4; j++)
5706  {
5707  ints[4*i+j] = v[j];
5708  }
5709  }
5710  MPI_Send(&ints[0], 4*NumOfElements, MPI_INT, 0, 447, MyComm);
5711  // boundary + shared faces
5713  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5714  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5716  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
5717  ints.SetSize(3*ne);
5718  for (i = 0; i < NumOfBdrElements; i++)
5719  {
5720  v = boundary[i]->GetVertices();
5721  for (j = 0; j < 3; j++)
5722  {
5723  ints[3*i+j] = v[j];
5724  }
5725  }
5726  for ( ; i < ne; i++)
5727  {
5728  v = shared_trias[i-NumOfBdrElements].v; // tet mesh
5729  for (j = 0; j < 3; j++)
5730  {
5731  ints[3*i+j] = v[j];
5732  }
5733  }
5734  MPI_Send(&ints[0], 3*ne, MPI_INT, 0, 447, MyComm);
5735  }
5736  }
5737 
5738  if (Dim == 3 && meshgen == 2)
5739  {
5740  int i, j, k, nv, ne, p;
5741  const int *ind, *v;
5742  MPI_Status status;
5743  Array<double> vert;
5744  Array<int> ints;
5745 
5746  int TG_nv, TG_ne, TG_nbe;
5747 
5748  if (MyRank == 0)
5749  {
5750  MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5751  MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5753  MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
5754 
5755  os << "TrueGrid\n"
5756  << "1 " << TG_nv << " " << TG_ne << " 0 0 0 0 0 0 0\n"
5757  << "0 0 0 1 0 0 0 0 0 0 0\n"
5758  << "0 0 " << TG_nbe << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
5759  << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
5760  << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
5761 
5762  // print the vertices
5763  nv = TG_nv;
5764  for (i = 0; i < NumOfVertices; i++)
5765  {
5766  os << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
5767  << " " << vertices[i](2) << " 0.0\n";
5768  }
5769  for (p = 1; p < NRanks; p++)
5770  {
5771  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5772  vert.SetSize(Dim*nv);
5773  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
5774  for (i = 0; i < nv; i++)
5775  {
5776  os << i+1 << " 0.0 " << vert[Dim*i] << " " << vert[Dim*i+1]
5777  << " " << vert[Dim*i+2] << " 0.0\n";
5778  }
5779  }
5780 
5781  // print the elements
5782  ne = TG_ne;
5783  for (i = 0; i < NumOfElements; i++)
5784  {
5785  nv = elements[i]->GetNVertices();
5786  ind = elements[i]->GetVertices();
5787  os << i+1 << " " << 1;
5788  for (j = 0; j < nv; j++)
5789  {
5790  os << " " << ind[j]+1;
5791  }
5792  os << '\n';
5793  }
5794  k = NumOfVertices;
5795  for (p = 1; p < NRanks; p++)
5796  {
5797  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5798  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5799  ints.SetSize(8*ne);
5800  MPI_Recv(&ints[0], 8*ne, MPI_INT, p, 447, MyComm, &status);
5801  for (i = 0; i < ne; i++)
5802  {
5803  os << i+1 << " " << p+1;
5804  for (j = 0; j < 8; j++)
5805  {
5806  os << " " << k+ints[i*8+j]+1;
5807  }
5808  os << '\n';
5809  }
5810  k += nv;
5811  }
5812 
5813  // print the boundary + shared faces information
5814  ne = TG_nbe;
5815  // boundary
5816  for (i = 0; i < NumOfBdrElements; i++)
5817  {
5818  nv = boundary[i]->GetNVertices();
5819  ind = boundary[i]->GetVertices();
5820  os << 1;
5821  for (j = 0; j < nv; j++)
5822  {
5823  os << " " << ind[j]+1;
5824  }
5825  os << " 1.0 1.0 1.0 1.0\n";
5826  }
5827  // shared faces
5828  const int sf_attr =
5829  MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
5830  // There are no shared triangle faces
5831  for (i = 0; i < shared_quads.Size(); i++)
5832  {
5833  ind = shared_quads[i].v;
5834  os << sf_attr;
5835  for (j = 0; j < 4; j++)
5836  {
5837  os << ' ' << ind[j]+1;
5838  }
5839  os << " 1.0 1.0 1.0 1.0\n";
5840  }
5841  k = NumOfVertices;
5842  for (p = 1; p < NRanks; p++)
5843  {
5844  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5845  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5846  ints.SetSize(4*ne);
5847  MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
5848  for (i = 0; i < ne; i++)
5849  {
5850  os << p+1;
5851  for (j = 0; j < 4; j++)
5852  {
5853  os << " " << k+ints[i*4+j]+1;
5854  }
5855  os << " 1.0 1.0 1.0 1.0\n";
5856  }
5857  k += nv;
5858  }
5859  }
5860  else
5861  {
5862  MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5863  MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5865  MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
5866 
5867  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5868  vert.SetSize(Dim*NumOfVertices);
5869  for (i = 0; i < NumOfVertices; i++)
5870  for (j = 0; j < Dim; j++)
5871  {
5872  vert[Dim*i+j] = vertices[i](j);
5873  }
5874  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE, 0, 445, MyComm);
5875  // elements
5876  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5877  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
5878  ints.SetSize(NumOfElements*8);
5879  for (i = 0; i < NumOfElements; i++)
5880  {
5881  v = elements[i]->GetVertices();
5882  for (j = 0; j < 8; j++)
5883  {
5884  ints[8*i+j] = v[j];
5885  }
5886  }
5887  MPI_Send(&ints[0], 8*NumOfElements, MPI_INT, 0, 447, MyComm);
5888  // boundary + shared faces
5889  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
5891  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
5892  ints.SetSize(4*ne);
5893  for (i = 0; i < NumOfBdrElements; i++)
5894  {
5895  v = boundary[i]->GetVertices();
5896  for (j = 0; j < 4; j++)
5897  {
5898  ints[4*i+j] = v[j];
5899  }
5900  }
5901  for ( ; i < ne; i++)
5902  {
5903  v = shared_quads[i-NumOfBdrElements].v; // hex mesh
5904  for (j = 0; j < 4; j++)
5905  {
5906  ints[4*i+j] = v[j];
5907  }
5908  }
5909  MPI_Send(&ints[0], 4*ne, MPI_INT, 0, 447, MyComm);
5910  }
5911  }
5912 
5913  if (Dim == 2)
5914  {
5915  int i, j, k, attr, nv, ne, p;
5916  Array<int> v;
5917  MPI_Status status;
5918  Array<double> vert;
5919  Array<int> ints;
5920 
5921  if (MyRank == 0)
5922  {
5923  os << "areamesh2\n\n";
5924 
5925  // print the boundary + shared edges information
5926  nv = NumOfBdrElements + shared_edges.Size();
5927  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5928  os << ne << '\n';
5929  // boundary
5930  for (i = 0; i < NumOfBdrElements; i++)
5931  {
5932  attr = boundary[i]->GetAttribute();
5933  boundary[i]->GetVertices(v);
5934  os << attr << " ";
5935  for (j = 0; j < v.Size(); j++)
5936  {
5937  os << v[j] + 1 << " ";
5938  }
5939  os << '\n';
5940  }
5941  // shared edges
5942  for (i = 0; i < shared_edges.Size(); i++)
5943  {
5944  attr = shared_edges[i]->GetAttribute();
5945  shared_edges[i]->GetVertices(v);
5946  os << attr << " ";
5947  for (j = 0; j < v.Size(); j++)
5948  {
5949  os << v[j] + 1 << " ";
5950  }
5951  os << '\n';
5952  }
5953  k = NumOfVertices;
5954  for (p = 1; p < NRanks; p++)
5955  {
5956  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5957  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5958  ints.SetSize(2*ne);
5959  MPI_Recv(&ints[0], 2*ne, MPI_INT, p, 447, MyComm, &status);
5960  for (i = 0; i < ne; i++)
5961  {
5962  os << p+1;
5963  for (j = 0; j < 2; j++)
5964  {
5965  os << " " << k+ints[i*2+j]+1;
5966  }
5967  os << '\n';
5968  }
5969  k += nv;
5970  }
5971 
5972  // print the elements
5973  nv = NumOfElements;
5974  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
5975  os << ne << '\n';
5976  for (i = 0; i < NumOfElements; i++)
5977  {
5978  // attr = elements[i]->GetAttribute(); // not used
5979  elements[i]->GetVertices(v);
5980  os << 1 << " " << 3 << " ";
5981  for (j = 0; j < v.Size(); j++)
5982  {
5983  os << v[j] + 1 << " ";
5984  }
5985  os << '\n';
5986  }
5987  k = NumOfVertices;
5988  for (p = 1; p < NRanks; p++)
5989  {
5990  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
5991  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
5992  ints.SetSize(3*ne);
5993  MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
5994  for (i = 0; i < ne; i++)
5995  {
5996  os << p+1 << " " << 3;
5997  for (j = 0; j < 3; j++)
5998  {
5999  os << " " << k+ints[i*3+j]+1;
6000  }
6001  os << '\n';
6002  }
6003  k += nv;
6004  }
6005 
6006  // print the vertices
6007  ne = NumOfVertices;
6008  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
6009  os << nv << '\n';
6010  for (i = 0; i < NumOfVertices; i++)
6011  {
6012  for (j = 0; j < Dim; j++)
6013  {
6014  os << vertices[i](j) << " ";
6015  }
6016  os << '\n';
6017  }
6018  for (p = 1; p < NRanks; p++)
6019  {
6020  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
6021  vert.SetSize(Dim*nv);
6022  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
6023  for (i = 0; i < nv; i++)
6024  {
6025  for (j = 0; j < Dim; j++)
6026  {
6027  os << " " << vert[Dim*i+j];
6028  }
6029  os << '\n';
6030  }
6031  }
6032  }
6033  else
6034  {
6035  // boundary + shared faces
6036  nv = NumOfBdrElements + shared_edges.Size();
6037  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
6038  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
6039  ne = NumOfBdrElements + shared_edges.Size();
6040  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
6041  ints.SetSize(2*ne);
6042  for (i = 0; i < NumOfBdrElements; i++)
6043  {
6044  boundary[i]->GetVertices(v);
6045  for (j = 0; j < 2; j++)
6046  {
6047  ints[2*i+j] = v[j];
6048  }
6049  }
6050  for ( ; i < ne; i++)
6051  {
6052  shared_edges[i-NumOfBdrElements]->GetVertices(v);
6053  for (j = 0; j < 2; j++)
6054  {
6055  ints[2*i+j] = v[j];
6056  }
6057  }
6058  MPI_Send(&ints[0], 2*ne, MPI_INT, 0, 447, MyComm);
6059  // elements
6060  ne = NumOfElements;
6061  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
6062  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
6063  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
6064  ints.SetSize(NumOfElements*3);
6065  for (i = 0; i < NumOfElements; i++)
6066  {
6067  elements[i]->GetVertices(v);
6068  for (j = 0; j < 3; j++)
6069  {
6070  ints[3*i+j] = v[j];
6071  }
6072  }
6073  MPI_Send(&ints[0], 3*NumOfElements, MPI_INT, 0, 447, MyComm);
6074  // vertices
6075  ne = NumOfVertices;
6076  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
6077  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
6078  vert.SetSize(Dim*NumOfVertices);
6079  for (i = 0; i < NumOfVertices; i++)
6080  for (j = 0; j < Dim; j++)
6081  {
6082  vert[Dim*i+j] = vertices[i](j);
6083  }
6084  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
6085  0, 445, MyComm);
6086  }
6087  }
6088 }
6089 
6090 void ParMesh::GetBoundingBox(Vector &gp_min, Vector &gp_max, int ref)
6091 {
6092  int sdim;
6093  Vector p_min, p_max;
6094 
6095  this->Mesh::GetBoundingBox(p_min, p_max, ref);
6096 
6097  sdim = SpaceDimension();
6098 
6099  gp_min.SetSize(sdim);
6100  gp_max.SetSize(sdim);
6101 
6102  MPI_Allreduce(p_min.GetData(), gp_min, sdim, MPI_DOUBLE, MPI_MIN, MyComm);
6103  MPI_Allreduce(p_max.GetData(), gp_max, sdim, MPI_DOUBLE, MPI_MAX, MyComm);
6104 }
6105 
6106 void ParMesh::GetCharacteristics(double &gh_min, double &gh_max,
6107  double &gk_min, double &gk_max)
6108 {
6109  double h_min, h_max, kappa_min, kappa_max;
6110 
6111  this->Mesh::GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
6112 
6113  MPI_Allreduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN, MyComm);
6114  MPI_Allreduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX, MyComm);
6115  MPI_Allreduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN, MyComm);
6116  MPI_Allreduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX, MyComm);
6117 }
6118 
6119 void ParMesh::PrintInfo(std::ostream &os)
6120 {
6121  int i;
6122  DenseMatrix J(Dim);
6123  double h_min, h_max, kappa_min, kappa_max, h, kappa;
6124 
6125  if (MyRank == 0)
6126  {
6127  os << "Parallel Mesh Stats:" << '\n';
6128  }
6129 
6130  for (i = 0; i < NumOfElements; i++)
6131  {
6132  GetElementJacobian(i, J);
6133  h = pow(fabs(J.Weight()), 1.0/double(Dim));
6134  kappa = (Dim == spaceDim) ?
6135  J.CalcSingularvalue(0) / J.CalcSingularvalue(Dim-1) : -1.0;
6136  if (i == 0)
6137  {
6138  h_min = h_max = h;
6139  kappa_min = kappa_max = kappa;
6140  }
6141  else
6142  {
6143  if (h < h_min) { h_min = h; }
6144  if (h > h_max) { h_max = h; }
6145  if (kappa < kappa_min) { kappa_min = kappa; }
6146  if (kappa > kappa_max) { kappa_max = kappa; }
6147  }
6148  }
6149 
6150  double gh_min, gh_max, gk_min, gk_max;
6151  MPI_Reduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
6152  MPI_Reduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
6153  MPI_Reduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
6154  MPI_Reduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
6155 
6156  // TODO: collect and print stats by geometry
6157 
6158  long long ldata[5]; // vert, edge, face, elem, neighbors;
6159  long long mindata[5], maxdata[5], sumdata[5];
6160 
6161  // count locally owned vertices, edges, and faces
6162  ldata[0] = GetNV();
6163  ldata[1] = GetNEdges();
6164  ldata[2] = GetNFaces();
6165  ldata[3] = GetNE();
6166  ldata[4] = gtopo.GetNumNeighbors()-1;
6167  for (int gr = 1; gr < GetNGroups(); gr++)
6168  {
6169  if (!gtopo.IAmMaster(gr)) // we are not the master
6170  {
6171  ldata[0] -= group_svert.RowSize(gr-1);
6172  ldata[1] -= group_sedge.RowSize(gr-1);
6173  ldata[2] -= group_stria.RowSize(gr-1);
6174  ldata[2] -= group_squad.RowSize(gr-1);
6175  }
6176  }
6177 
6178  MPI_Reduce(ldata, mindata, 5, MPI_LONG_LONG, MPI_MIN, 0, MyComm);
6179  MPI_Reduce(ldata, sumdata, 5, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
6180  MPI_Reduce(ldata, maxdata, 5, MPI_LONG_LONG, MPI_MAX, 0, MyComm);
6181 
6182  if (MyRank == 0)
6183  {
6184  os << '\n'
6185  << " "
6186  << setw(12) << "minimum"
6187  << setw(12) << "average"
6188  << setw(12) << "maximum"
6189  << setw(12) << "total" << '\n';
6190  os << " vertices "
6191  << setw(12) << mindata[0]
6192  << setw(12) << sumdata[0]/NRanks
6193  << setw(12) << maxdata[0]
6194  << setw(12) << sumdata[0] << '\n';
6195  os << " edges "
6196  << setw(12) << mindata[1]
6197  << setw(12) << sumdata[1]/NRanks
6198  << setw(12) << maxdata[1]
6199  << setw(12) << sumdata[1] << '\n';
6200  if (Dim == 3)
6201  {
6202  os << " faces "
6203  << setw(12) << mindata[2]
6204  << setw(12) << sumdata[2]/NRanks
6205  << setw(12) << maxdata[2]
6206  << setw(12) << sumdata[2] << '\n';
6207  }
6208  os << " elements "
6209  << setw(12) << mindata[3]
6210  << setw(12) << sumdata[3]/NRanks
6211  << setw(12) << maxdata[3]
6212  << setw(12) << sumdata[3] << '\n';
6213  os << " neighbors "
6214  << setw(12) << mindata[4]
6215  << setw(12) << sumdata[4]/NRanks
6216  << setw(12) << maxdata[4] << '\n';
6217  os << '\n'
6218  << " "
6219  << setw(12) << "minimum"
6220  << setw(12) << "maximum" << '\n';
6221  os << " h "
6222  << setw(12) << gh_min
6223  << setw(12) << gh_max << '\n';
6224  os << " kappa "
6225  << setw(12) << gk_min
6226  << setw(12) << gk_max << '\n';
6227  os << std::flush;
6228  }
6229 }
6230 
6231 long long ParMesh::ReduceInt(int value) const
6232 {
6233  long long local = value, global;
6234  MPI_Allreduce(&local, &global, 1, MPI_LONG_LONG, MPI_SUM, MyComm);
6235  return global;
6236 }
6237 
6238 void ParMesh::ParPrint(ostream &os) const
6239 {
6240  if (NURBSext)
6241  {
6242  // TODO: NURBS meshes.
6243  Print(os); // use the serial MFEM v1.0 format for now
6244  return;
6245  }
6246 
6247  if (Nonconforming())
6248  {
6249  // the NC mesh format works both in serial and in parallel
6250  Printer(os);
6251  return;
6252  }
6253 
6254  // Write out serial mesh. Tell serial mesh to deliniate the end of it's
6255  // output with 'mfem_serial_mesh_end' instead of 'mfem_mesh_end', as we will
6256  // be adding additional parallel mesh information.
6257  Printer(os, "mfem_serial_mesh_end");
6258 
6259  // write out group topology info.
6260  gtopo.Save(os);
6261 
6262  os << "\ntotal_shared_vertices " << svert_lvert.Size() << '\n';
6263  if (Dim >= 2)
6264  {
6265  os << "total_shared_edges " << shared_edges.Size() << '\n';
6266  }
6267  if (Dim >= 3)
6268  {
6269  os << "total_shared_faces " << sface_lface.Size() << '\n';
6270  }
6271  for (int gr = 1; gr < GetNGroups(); gr++)
6272  {
6273  {
6274  const int nv = group_svert.RowSize(gr-1);
6275  const int *sv = group_svert.GetRow(gr-1);
6276  os << "\n# group " << gr << "\nshared_vertices " << nv << '\n';
6277  for (int i = 0; i < nv; i++)
6278  {
6279  os << svert_lvert[sv[i]] << '\n';
6280  }
6281  }
6282  if (Dim >= 2)
6283  {
6284  const int ne = group_sedge.RowSize(gr-1);
6285  const int *se = group_sedge.GetRow(gr-1);
6286  os << "\nshared_edges " << ne << '\n';
6287  for (int i = 0; i < ne; i++)
6288  {
6289  const int *v = shared_edges[se[i]]->GetVertices();
6290  os << v[0] << ' ' << v[1] << '\n';
6291  }
6292  }
6293  if (Dim >= 3)
6294  {
6295  const int nt = group_stria.RowSize(gr-1);
6296  const int *st = group_stria.GetRow(gr-1);
6297  const int nq = group_squad.RowSize(gr-1);
6298  const int *sq = group_squad.GetRow(gr-1);
6299  os << "\nshared_faces " << nt+nq << '\n';
6300  for (int i = 0; i < nt; i++)
6301  {
6302  os << Geometry::TRIANGLE;
6303  const int *v = shared_trias[st[i]].v;
6304  for (int j = 0; j < 3; j++) { os << ' ' << v[j]; }
6305  os << '\n';
6306  }
6307  for (int i = 0; i < nq; i++)
6308  {
6309  os << Geometry::SQUARE;
6310  const int *v = shared_quads[sq[i]].v;
6311  for (int j = 0; j < 4; j++) { os << ' ' << v[j]; }
6312  os << '\n';
6313  }
6314  }
6315  }
6316 
6317  // Write out section end tag for mesh.
6318  os << "\nmfem_mesh_end" << endl;
6319 }
6320 
6321 void ParMesh::PrintVTU(std::string pathname,
6322  VTKFormat format,
6323  bool high_order_output,
6324  int compression_level,
6325  bool bdr)
6326 {
6327  int pad_digits_rank = 6;
6328  DataCollection::create_directory(pathname, this, MyRank);
6329 
6330  std::string::size_type pos = pathname.find_last_of('/');
6331  std::string fname
6332  = (pos == std::string::npos) ? pathname : pathname.substr(pos+1);
6333 
6334  if (MyRank == 0)
6335  {
6336  std::string pvtu_name = pathname + "/" + fname + ".pvtu";
6337  std::ofstream os(pvtu_name);
6338 
6339  std::string data_type = (format == VTKFormat::BINARY32) ? "Float32" : "Float64";
6340  std::string data_format = (format == VTKFormat::ASCII) ? "ascii" : "binary";
6341 
6342  os << "<?xml version=\"1.0\"?>\n";
6343  os << "<VTKFile type=\"PUnstructuredGrid\"";
6344  os << " version =\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
6345  os << "<PUnstructuredGrid GhostLevel=\"0\">\n";
6346 
6347  os << "<PPoints>\n";
6348  os << "\t<PDataArray type=\"" << data_type << "\" ";
6349  os << " Name=\"Points\" NumberOfComponents=\"3\""
6350  << " format=\"" << data_format << "\"/>\n";
6351  os << "</PPoints>\n";
6352 
6353  os << "<PCells>\n";
6354  os << "\t<PDataArray type=\"Int32\" ";
6355  os << " Name=\"connectivity\" NumberOfComponents=\"1\""
6356  << " format=\"" << data_format << "\"/>\n";
6357  os << "\t<PDataArray type=\"Int32\" ";
6358  os << " Name=\"offsets\" NumberOfComponents=\"1\""
6359  << " format=\"" << data_format << "\"/>\n";
6360  os << "\t<PDataArray type=\"UInt8\" ";
6361  os << " Name=\"types\" NumberOfComponents=\"1\""
6362  << " format=\"" << data_format << "\"/>\n";
6363  os << "</PCells>\n";
6364 
6365  os << "<PCellData>\n";
6366  os << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
6367  << "\" NumberOfComponents=\"1\""
6368  << " format=\"" << data_format << "\"/>\n";
6369  os << "</PCellData>\n";
6370 
6371  for (int ii=0; ii<NRanks; ii++)
6372  {
6373  std::string piece = fname + ".proc"
6374  + to_padded_string(ii, pad_digits_rank) + ".vtu";
6375  os << "<Piece Source=\"" << piece << "\"/>\n";
6376  }
6377 
6378  os << "</PUnstructuredGrid>\n";
6379  os << "</VTKFile>\n";
6380  os.close();
6381  }
6382 
6383  std::string vtu_fname = pathname + "/" + fname + ".proc"
6384  + to_padded_string(MyRank, pad_digits_rank);
6385  Mesh::PrintVTU(vtu_fname, format, high_order_output, compression_level, bdr);
6386 }
6387 
6388 int ParMesh::FindPoints(DenseMatrix& point_mat, Array<int>& elem_id,
6389  Array<IntegrationPoint>& ip, bool warn,
6390  InverseElementTransformation *inv_trans)
6391 {
6392  const int npts = point_mat.Width();
6393  if (npts == 0) { return 0; }
6394 
6395  const bool no_warn = false;
6396  Mesh::FindPoints(point_mat, elem_id, ip, no_warn, inv_trans);
6397 
6398  // If multiple processors find the same point, we need to choose only one of
6399  // the processors to mark that point as found.
6400  // Here, we choose the processor with the minimal rank.
6401 
6402  Array<int> my_point_rank(npts), glob_point_rank(npts);
6403  for (int k = 0; k < npts; k++)
6404  {
6405  my_point_rank[k] = (elem_id[k] == -1) ? NRanks : MyRank;
6406  }
6407 
6408  MPI_Allreduce(my_point_rank.GetData(), glob_point_rank.GetData(), npts,
6409  MPI_INT, MPI_MIN, MyComm);
6410 
6411  int pts_found = 0;
6412  for (int k = 0; k < npts; k++)
6413  {
6414  if (glob_point_rank[k] == NRanks) { elem_id[k] = -1; }
6415  else
6416  {
6417  pts_found++;
6418  if (glob_point_rank[k] != MyRank) { elem_id[k] = -2; }
6419  }
6420  }
6421  if (warn && pts_found != npts && MyRank == 0)
6422  {
6423  MFEM_WARNING((npts-pts_found) << " points were not found");
6424  }
6425  return pts_found;
6426 }
6427 
6428 static void PrintVertex(const Vertex &v, int space_dim, ostream &os)
6429 {
6430  os << v(0);
6431  for (int d = 1; d < space_dim; d++)
6432  {
6433  os << ' ' << v(d);
6434  }
6435 }
6436 
6437 void ParMesh::PrintSharedEntities(const char *fname_prefix) const
6438 {
6439  stringstream out_name;
6440  out_name << fname_prefix << '_' << setw(5) << setfill('0') << MyRank
6441  << ".shared_entities";
6442  ofstream os(out_name.str().c_str());
6443  os.precision(16);
6444 
6445  gtopo.Save(out);
6446 
6447  os << "\ntotal_shared_vertices " << svert_lvert.Size() << '\n';
6448  if (Dim >= 2)
6449  {
6450  os << "total_shared_edges " << shared_edges.Size() << '\n';
6451  }
6452  if (Dim >= 3)
6453  {
6454  os << "total_shared_faces " << sface_lface.Size() << '\n';
6455  }
6456  for (int gr = 1; gr < GetNGroups(); gr++)
6457  {
6458  {
6459  const int nv = group_svert.RowSize(gr-1);
6460  const int *sv = group_svert.GetRow(gr-1);
6461  os << "\n# group " << gr << "\n\nshared_vertices " << nv << '\n';
6462  for (int i = 0; i < nv; i++)
6463  {
6464  const int lvi = svert_lvert[sv[i]];
6465  // os << lvi << '\n';
6466  PrintVertex(vertices[lvi], spaceDim, os);
6467  os << '\n';
6468  }
6469  }
6470  if (Dim >= 2)
6471  {
6472  const int ne = group_sedge.RowSize(gr-1);
6473  const int *se = group_sedge.GetRow(gr-1);
6474  os << "\nshared_edges " << ne << '\n';
6475  for (int i = 0; i < ne; i++)
6476  {
6477  const int *v = shared_edges[se[i]]->GetVertices();
6478  // os << v[0] << ' ' << v[1] << '\n';
6479  PrintVertex(vertices[v[0]], spaceDim, os);
6480  os << " | ";
6481  PrintVertex(vertices[v[1]], spaceDim, os);
6482  os << '\n';
6483  }
6484  }
6485  if (Dim >= 3)
6486  {
6487  const int nt = group_stria.RowSize(gr-1);
6488  const int *st = group_stria.GetRow(gr-1);
6489  const int nq = group_squad.RowSize(gr-1);
6490  const int *sq = group_squad.GetRow(gr-1);
6491  os << "\nshared_faces " << nt+nq << '\n';
6492  for (int i = 0; i < nt; i++)
6493  {
6494  const int *v = shared_trias[st[i]].v;
6495 #if 0
6496  os << Geometry::TRIANGLE;
6497  for (int j = 0; j < 3; j++) { os << ' ' << v[j]; }
6498  os << '\n';
6499 #endif
6500  for (int j = 0; j < 3; j++)
6501  {
6502  PrintVertex(vertices[v[j]], spaceDim, os);
6503  (j < 2) ? os << " | " : os << '\n';
6504  }
6505  }
6506  for (int i = 0; i < nq; i++)
6507  {
6508  const int *v = shared_quads[sq[i]].v;
6509 #if 0
6510  os << Geometry::SQUARE;
6511  for (int j = 0; j < 4; j++) { os << ' ' << v[j]; }
6512  os << '\n';
6513 #endif
6514  for (int j = 0; j < 4; j++)
6515  {
6516  PrintVertex(vertices[v[j]], spaceDim, os);
6517  (j < 3) ? os << " | " : os << '\n';
6518  }
6519  }
6520  }
6521  }
6522 }
6523 
6525 {
6526  H1_FECollection fec(1, Dim); // Order 1, mesh dimension (not spatial dimension).
6527  ParMesh *pm = const_cast<ParMesh *>(this);
6528  ParFiniteElementSpace fespace(pm, &fec);
6529 
6530  gi.SetSize(GetNV());
6531 
6532  Array<int> dofs;
6533  for (int i=0; i<GetNV(); ++i)
6534  {
6535  fespace.GetVertexDofs(i, dofs);
6536  gi[i] = fespace.GetGlobalTDofNumber(dofs[0]);
6537  }
6538 }
6539 
6541 {
6542  if (Dim == 1)
6543  {
6545  return;
6546  }
6547 
6548  ND_FECollection fec(1, Dim); // Order 1, mesh dimension (not spatial dimension).
6549  ParMesh *pm = const_cast<ParMesh *>(this);
6550  ParFiniteElementSpace fespace(pm, &fec);
6551 
6552  gi.SetSize(GetNEdges());
6553 
6554  Array<int> dofs;
6555  for (int i=0; i<GetNEdges(); ++i)
6556  {
6557  fespace.GetEdgeDofs(i, dofs);
6558  const int ldof = (dofs[0] >= 0) ? dofs[0] : -1 - dofs[0];
6559  gi[i] = fespace.GetGlobalTDofNumber(ldof);
6560  }
6561 }
6562 
6564 {
6565  if (Dim == 2)
6566  {
6568  return;
6569  }
6570  else if (Dim == 1)
6571  {
6573  return;
6574  }
6575 
6576  RT_FECollection fec(0, Dim); // Order 0, mesh dimension (not spatial dimension).
6577  ParMesh *pm = const_cast<ParMesh *>(this);
6578  ParFiniteElementSpace fespace(pm, &fec);
6579 
6580  gi.SetSize(GetNFaces());
6581 
6582  Array<int> dofs;
6583  for (int i=0; i<GetNFaces(); ++i)
6584  {
6585  fespace.GetFaceDofs(i, dofs);
6586  const int ldof = (dofs[0] >= 0) ? dofs[0] : -1 - dofs[0];
6587  gi[i] = fespace.GetGlobalTDofNumber(ldof);
6588  }
6589 }
6590 
6592 {
6594 
6595  // Cast from long long to HYPRE_BigInt
6596  const HYPRE_BigInt offset = glob_elem_offset;
6597 
6598  gi.SetSize(GetNE());
6599  for (int i=0; i<GetNE(); ++i)
6600  {
6601  gi[i] = offset + i;
6602  }
6603 }
6604 
6606 {
6607  Mesh::Swap(other, true);
6608 
6609  mfem::Swap(MyComm, other.MyComm);
6610  mfem::Swap(NRanks, other.NRanks);
6611  mfem::Swap(MyRank, other.MyRank);
6612 
6615 
6616  gtopo.Swap(other.gtopo);
6617 
6618  group_svert.Swap(other.group_svert);
6619  group_sedge.Swap(other.group_sedge);
6620  group_stria.Swap(other.group_stria);
6621  group_squad.Swap(other.group_squad);
6622 
6629 
6630  // Swap face-neighbor data
6639 
6640  // Nodes, NCMesh, and NURBSExtension are taken care of by Mesh::Swap
6641  mfem::Swap(pncmesh, other.pncmesh);
6642 
6643  print_shared = other.print_shared;
6644 }
6645 
6647 {
6648  delete pncmesh;
6649  ncmesh = pncmesh = NULL;
6650 
6652 
6653  for (int i = 0; i < shared_edges.Size(); i++)
6654  {
6656  }
6657  shared_edges.DeleteAll();
6658 
6659  delete face_nbr_el_to_face;
6660  face_nbr_el_to_face = NULL;
6661 }
6662 
6664 {
6665  ParMesh::Destroy();
6666 
6667  // The Mesh destructor is called automatically
6668 }
6669 
6670 }
6671 
6672 #endif
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
void UniformRefinement3D() override
Refine a mixed 3D mesh uniformly.
Definition: pmesh.cpp:4518
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
int GetGroupMasterRank(int g) const
Return the rank of the group master for group &#39;g&#39;.
int GetNFaceNeighbors() const
Definition: pmesh.hpp:462
void PrintAsOneXG(std::ostream &out=mfem::out)
Old mesh format (Netgen/Truegrid) version of &#39;PrintAsOne&#39;.
Definition: pmesh.cpp:5560
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:3966
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
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:2742
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void GroupQuadrilateral(int group, int i, int &face, int &o) const
Definition: pmesh.cpp:1626
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
void GetConformingSharedStructures(class ParMesh &pmesh)
Definition: pncmesh.cpp:787
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:2963
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetGlobalFaceIndices(Array< HYPRE_BigInt > &gi) const
AMR meshes are not supported.
Definition: pmesh.cpp:6563
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:706
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1012
Element on side 1 is configured.
Definition: eltrans.hpp:513
virtual ~ParMesh()
Definition: pmesh.cpp:6663
void LoadSharedEntities(std::istream &input)
Definition: pmesh.cpp:978
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
long glob_offset_sequence
Definition: pmesh.hpp:87
void UniformRefineGroups2D(int old_nv)
Definition: pmesh.cpp:4294
void FreeElement(Element *E)
Definition: mesh.cpp:11873
FaceElementTransformations * GetSharedFaceTransformationsByLocalIndex(int FaceNo, bool fill2=true)
Definition: pmesh.cpp:3019
int CheckElementOrientation(bool fix_it=true)
Check (and optionally attempt to fix) the orientation of the elements.
Definition: mesh.cpp:5420
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 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
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:6016
int FindSharedVertices(const int *partition, Table *vertex_element, ListOfIntegerSets &groups)
Definition: pmesh.cpp:585
int GroupNVertices(int group) const
Definition: pmesh.hpp:395
void AddColumnsInRow(int r, int ncol)
Definition: table.hpp:78
void UniformRefinement3D_base(Array< int > *f2qf=NULL, DSTable *v_to_v_p=NULL, bool update_nodes=true)
Definition: mesh.cpp:8165
void PrintInfo(std::ostream &out=mfem::out) override
Print various parallel mesh stats.
Definition: pmesh.cpp:6119
int GroupNEdges(int group) const
Definition: pmesh.hpp:396
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:81
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:125
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
int GroupNTriangles(int group) const
Definition: pmesh.hpp:397
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
Array< Element * > boundary
Definition: mesh.hpp:91
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:3173
friend class ParNCMesh
Definition: pmesh.hpp:653
int * GeneratePartitioning(int nparts, int part_method=1)
Definition: mesh.cpp:7033
int BuildLocalVertices(const Mesh &global_mesh, const int *partitioning, Array< int > &vert_global_local)
Fills out partitioned Mesh::vertices.
Definition: pmesh.cpp:302
int GetNGroups() const
Definition: pmesh.hpp:392
void UniformRefinement2D() override
Refine a mixed 2D mesh uniformly.
Definition: pmesh.cpp:4494
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:6171
CoarseFineTransformations CoarseFineTr
Definition: mesh.hpp:235
int own_nodes
Definition: mesh.hpp:241
bool Conforming() const
Definition: mesh.hpp:1650
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:926
IsoparametricTransformation Transformation
Definition: mesh.hpp:229
void GetCharacteristics(double &h_min, double &h_max, double &kappa_min, double &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
Definition: mesh.cpp:193
Array< Element * > face_nbr_elements
Definition: pmesh.hpp:382
static FiniteElement * GetTransformationFEforElementType(Element::Type)
Definition: mesh.cpp:327
void GetSharedVertexCommunicator(int ordering, GroupCommunicator &svert_comm) const
Get the shared vertices GroupCommunicator.
Definition: pmesh.cpp:1661
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
Definition: mesh.cpp:7935
void GetFaceNeighbors(class ParMesh &pmesh)
Definition: pncmesh.cpp:886
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
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
int GetGroupSize(int g) const
Get the number of processors in a group.
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3127
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:6250
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
void ShiftRight(int &a, int &b, int &c)
Definition: mesh.hpp:1966
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:206
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
const int * GetGroup(int g) const
Return a pointer to a list of neighbors for a given group. Neighbor 0 is the local processor...
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
Definition: mesh.cpp:6255
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
T * GetData()
Returns the data.
Definition: array.hpp:112
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:246
virtual DofTransformation * GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: pfespace.cpp:495
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
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
GridFunction * Nodes
Definition: mesh.hpp:240
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int NumOfElements
Definition: mesh.hpp:69
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
int GetLocalElementNum(long long global_element_num) const
Definition: pmesh.cpp:1530
Point transformation for side 2 is configured.
Definition: eltrans.hpp:516
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:524
void ApplyLocalSlaveTransformation(FaceElementTransformations &FT, const FaceInfo &fi, bool is_ghost)
Definition: mesh.cpp:1056
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:6194
void GetCharacteristics(double &h_min, double &h_max, double &kappa_min, double &kappa_max)
Definition: pmesh.cpp:6106
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:7970
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
void SynchronizeDerefinementData(Array< Type > &elem_data, const Table &deref_table)
Definition: pncmesh.cpp:1545
const Element * GetFace(int i) const
Definition: mesh.hpp:1048
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:2985
void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0) override
This function is not public anymore. Use GeneralRefinement instead.
Definition: pmesh.cpp:3879
void EnsureParNodes()
If the mesh is curved, make sure &#39;Nodes&#39; is ParGridFunction.
Definition: pmesh.cpp:2095
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
Definition: element.cpp:17
void UniformRefineGroups3D(int old_nv, int old_nedges, const DSTable &old_v_to_v, const STable3D &old_faces, Array< int > *f2qf)
Definition: pmesh.cpp:4345
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:6388
virtual void Derefine(const Array< int > &derefs)
Definition: pncmesh.cpp:1373
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
const FiniteElement * GetFE() const
Get the current element used to compute the transformations.
Definition: eltrans.hpp:389
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.hpp:236
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:124
Data type for vertex.
Definition: vertex.hpp:22
void SetMeshGen()
Determine the mesh generator bitmask meshgen, see MeshGenerator().
Definition: mesh.cpp:3913
Array< int > RefEdges
Definition: geom.hpp:283
Array< int > face_nbr_group
Definition: pmesh.hpp:379
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:5616
void GetVertices(Vector &vert_coord) const
Definition: mesh.cpp:7837
bool NonconformingDerefinement(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1) override
NC version of GeneralDerefinement.
Definition: pmesh.cpp:3931
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:3999
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
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
void Swap(GroupTopology &other)
Swap the internal data with another GroupTopology object.
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:961
int Index(int r, int c, int f) const
Definition: stable3d.cpp:117
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
int Size_of_connections() const
Definition: table.hpp:98
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
STable3D * GetSharedFacesTable()
Definition: pmesh.cpp:2661
void GetVertexToVertexTable(DSTable &) const
Definition: mesh.cpp:6331
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:354
Operation last_operation
Definition: mesh.hpp:284
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:846
Table group_stria
Definition: pmesh.hpp:71
void AddConnections(int r, const int *c, int nc)
Definition: table.cpp:104
void InitRefinementTransforms()
Definition: mesh.cpp:9906
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:8006
Array< NCFaceInfo > nc_faces_info
Definition: mesh.hpp:218
bool IAmMaster(int g) const
Return true if I am master for group &#39;g&#39;.
bool WantSkipSharedMaster(const NCMesh::Master &master) const
Definition: pmesh.cpp:4757
void OnMeshUpdated(Mesh *mesh)
Definition: ncmesh.cpp:2398
A parallel extension of the NCMesh class.
Definition: pncmesh.hpp:64
void GetSharedEdgeDofs(int group, int ei, Array< int > &dofs) const
Definition: pfespace.cpp:584
Element * NewElement(int geom)
Definition: mesh.cpp:3841
Table * el_to_face
Definition: mesh.hpp:221
FaceElementTransformations FaceElemTr
Definition: mesh.hpp:232
ParNCMesh * pncmesh
Definition: pmesh.hpp:388
Element::Type GetFaceElementType(int Face) const
Definition: mesh.cpp:1432
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
Definition: mesh.cpp:5362
void CopyTo(U *dest)
STL-like copyTo dest from begin to end.
Definition: array.hpp:279
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:289
int FindSharedEdges(const Mesh &mesh, const int *partition, Table *&edge_element, ListOfIntegerSets &groups)
Definition: pmesh.cpp:531
int GetNRanks() const
Definition: pmesh.hpp:352
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5376
void MakeRefined_(Mesh &orig_mesh, const Array< int > ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
Definition: mesh.cpp:4333
void ExchangeFaceNbrData()
Definition: pmesh.cpp:2117
void BuildSharedVertMapping(int nvert, const Table *vert_element, const Array< int > &vert_global_local)
Definition: pmesh.cpp:834
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
virtual void GetElementDofValues(int el, Vector &dof_vals) const
Definition: gridfunc.cpp:1844
void ReduceMeshGen()
Definition: pmesh.cpp:881
double Weight() const
Definition: densemat.cpp:506
void NURBSUniformRefinement() override
Refine NURBS mesh.
Definition: pmesh.cpp:4547
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
Definition: segment.cpp:40
void GetElementJacobian(int i, DenseMatrix &J)
Definition: mesh.cpp:60
Geometry Geometries
Definition: fe.cpp:49
void Rebalance()
Definition: pmesh.cpp:3989
void GetVertexDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:3033
void MarkEdge(DenseMatrix &pmat)
Definition: triangle.cpp:53
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1613
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
Definition: fe_coll.cpp:1917
bool Nonconforming() const
Definition: mesh.hpp:1651
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 ...
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1517
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
Array< int > face_nbr_elements_offset
Definition: pmesh.hpp:380
void GetEdgeOrdering(DSTable &v_to_v, Array< int > &order)
Definition: mesh.cpp:2433
MPI_Comm MyComm
Definition: pmesh.hpp:37
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition: pmesh.cpp:1516
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
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:6090
const Table & GetDerefinementTable()
Definition: ncmesh.cpp:1807
Mesh GetSerialMesh(int save_rank) const
Definition: pmesh.cpp:5234
DenseTensor point_matrices[Geometry::NumGeom]
Definition: ncmesh.hpp:77
Array< Element * > shared_edges
Definition: pmesh.hpp:61
virtual void SetAttributes()
Definition: mesh.cpp:1559
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
void Destroy()
Definition: pmesh.cpp:6646
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:524
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3146
int GetNumNeighbors() const
Return the number of neighbors including the local processor.
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)
Returns the transformation defining the given face element in a user-defined variable.
Definition: mesh.cpp:497
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31) override
Definition: pmesh.cpp:2995
static const int NumVerts[NumGeom]
Definition: geom.hpp:49
void AddConnection(int r, int c)
Definition: table.hpp:80
STable3D * GetElementToFaceTable(int ret_ftbl=0)
Definition: mesh.cpp:6798
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3011
void LoseData()
NULL-ifies the data.
Definition: array.hpp:132
void Reset()
Force the reevaluation of the Jacobian in the next call.
Definition: eltrans.hpp:89
Table send_face_nbr_vertices
Definition: pmesh.hpp:386
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
int InitialPartition(int index) const
Helper to get the partitioning when the serial mesh gets split initially.
Definition: pncmesh.hpp:309
This structure stores the low level information necessary to interpret the configuration of elements ...
Definition: mesh.hpp:153
VTKFormat
Data array format for VTK and VTU files.
Definition: vtk.hpp:96
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
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
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:156
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:6321
void GetSharedQuadrilateralDofs(int group, int fi, Array< int > &dofs) const
Definition: pfespace.cpp:631
int AddElement(Element *elem)
The parameter elem should be allocated using the NewElement() method.
Definition: mesh.cpp:1809
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
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
int GetElementToEdgeTable(Table &, Array< int > &)
Definition: mesh.cpp:6356
virtual MFEM_DEPRECATED int GetNFaces(int &nFaceVertices) const =0
void MarkCoarseLevel()
Definition: ncmesh.cpp:4302
const Element * GetElement(int i) const
Definition: mesh.hpp:1040
void ResetLazyData()
Definition: mesh.cpp:1549
IntegrationRule RefPts
Definition: geom.hpp:282
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:3168
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
void GetGlobalVertexIndices(Array< HYPRE_BigInt > &gi) const
AMR meshes are not supported.
Definition: pmesh.cpp:6524
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:896
void GetSharedTriangleDofs(int group, int fi, Array< int > &dofs) const
Definition: pfespace.cpp:607
IsoparametricTransformation Transformation2
Definition: mesh.hpp:229
Table * face_nbr_el_to_face
Definition: pmesh.hpp:80
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
FiniteElementCollection * OwnFEC()
Definition: gridfunc.hpp:126
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
long long ReduceInt(int value) const override
Utility function: sum integers from all processors (Allreduce).
Definition: pmesh.cpp:6231
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
int Dimension() const
Definition: mesh.hpp:1006
int NumOfBdrElements
Definition: mesh.hpp:69
Table * el_to_edge
Definition: mesh.hpp:220
int Size()
Return the number of integer sets in the list.
Definition: sets.hpp:70
STable3D * GetFaceNbrElementToFaceTable(int ret_ftbl=0)
Definition: pmesh.cpp:2710
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:471
Array< Vert4 > shared_quads
Definition: pmesh.hpp:66
void Rebalance(const Array< int > *custom_partition=NULL)
Definition: pncmesh.cpp:1683
int slaves_end
slave faces
Definition: ncmesh.hpp:205
Data type tetrahedron element.
Definition: tetrahedron.hpp:22
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
void PrintXG(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4555
Point transformation for side 1 is configured.
Definition: eltrans.hpp:515
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:10503
Vector & FaceNbrData()
Definition: pgridfunc.hpp:208
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
Definition: mesh.cpp:5797
int SpaceDimension() const
Definition: mesh.hpp:1007
Table face_nbr_el_ori
Definition: pmesh.hpp:81
int GroupNQuadrilaterals(int group) const
Definition: pmesh.hpp:398
void Save(std::ostream &out) const
Save the data in a stream.
HYPRE_BigInt GetGlobalTDofNumber(int ldof) const
Returns the global tdof number of the given local degree of freedom.
Definition: pfespace.cpp:1083
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1394
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void FindSharedFaces(const Mesh &mesh, const int *partition, Array< int > &face_group, ListOfIntegerSets &groups)
Definition: pmesh.cpp:504
int GetMyRank() const
Definition: pmesh.hpp:353
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:3375
void RefineGroups(const DSTable &v_to_v, int *middle)
Update the groups after triangle refinement.
Definition: pmesh.cpp:4037
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
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:96
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
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 GetSharedTriCommunicator(int ordering, GroupCommunicator &stria_comm) const
Get the shared face triangles GroupCommunicator.
Definition: pmesh.cpp:1709
static void PrintElement(const Element *, std::ostream &)
Definition: mesh.cpp:3907
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
int GetNeighborRank(int i) const
Return the MPI rank of neighbor &#39;i&#39;.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void Swap(ParMesh &other)
Definition: pmesh.cpp:6605
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2890
void AddAColumnInRow(int r)
Definition: table.hpp:77
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
Definition: mesh.hpp:1541
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
Definition: mesh.cpp:2481
void PrintSharedEntities(const char *fname_prefix) const
Debugging method.
Definition: pmesh.cpp:6437
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
bool IsSlaveFace(const FaceInfo &fi) const
Definition: mesh.cpp:1051
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
void DeleteLast()
Delete the last entry of the array.
Definition: array.hpp:193
void GetFaceNbrElementFaces(int i, Array< int > &fcs, Array< int > &cor) const
Definition: pmesh.cpp:2875
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
Definition: pmesh.cpp:1580
double GetLength(int i, int j) const
Return the length of the segment from node i to node j.
Definition: mesh.cpp:6294
void ShiftUpI()
Definition: table.cpp:115
int meshgen
Definition: mesh.hpp:77
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:228
void GetSharedEdgeCommunicator(int ordering, GroupCommunicator &sedge_comm) const
Get the shared edges GroupCommunicator.
Definition: pmesh.cpp:1637
void BuildVertexGroup(int ngroups, const Table &vert_element)
Definition: pmesh.cpp:694
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:51
void GetFaceNbrElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: pmesh.cpp:1976
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:920
Array< int > be_to_edge
Definition: mesh.hpp:223
Face transformation is configured.
Definition: eltrans.hpp:517
virtual const char * Name() const
Definition: fe_coll.hpp:65
void DistributeAttributes(Array< int > &attr)
Ensure that bdr_attributes and attributes agree across processors.
Definition: pmesh.cpp:1544
long long GetGlobalElementNum(int local_element_num) const
Map a local element number to a global element number.
Definition: pmesh.cpp:1538
void GroupEdge(int group, int i, int &edge, int &o) const
Definition: pmesh.cpp:1607
virtual void Refine(const Array< Refinement > &refinements)
Definition: pncmesh.cpp:1246
void Swap(Mesh &other, bool non_geometry)
Definition: mesh.cpp:9285
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
void GetSharedQuadCommunicator(int ordering, GroupCommunicator &squad_comm) const
Get the shared face quadrilaterals GroupCommunicator.
Definition: pmesh.cpp:1685
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
Table group_sedge
Definition: pmesh.hpp:70
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:5338
void ComputeGlobalElementOffset() const
Definition: pmesh.cpp:868
void ExchangeFaceNbrNodes()
Definition: pmesh.cpp:2601
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:2983
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:193
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
virtual int GetFaceDofs(int i, Array< int > &dofs, int variant=0) const
Definition: pfespace.cpp:519
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 &quot;BigEndian&quot; or &quot;LittleEndian&quot;.
Definition: vtk.cpp:602
void MakeJ()
Definition: table.cpp:91
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:955
Table group_svert
Shared objects in each group.
Definition: pmesh.hpp:69
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1) override
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
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
void GroupTriangle(int group, int i, int &face, int &o) const
Definition: pmesh.cpp:1615
long sequence
Definition: mesh.hpp:83
void GetGlobalElementIndices(Array< HYPRE_BigInt > &gi) const
AMR meshes are supported.
Definition: pmesh.cpp:6591
IsoparametricTransformation Transf
Definition: eltrans.hpp:464
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:2781
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
bool IsGhost(int entity, int index) const
Return true if the specified vertex/edge/face is a ghost.
Definition: pncmesh.hpp:184
ElementTransformation * Elem1
Definition: eltrans.hpp:522
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1229
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
Definition: mesh.cpp:4569
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
void GetGlobalEdgeIndices(Array< HYPRE_BigInt > &gi) const
AMR meshes are not supported.
Definition: pmesh.cpp:6540
STable3D * GetFacesTable()
Definition: mesh.cpp:6735
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:273
int GroupVertex(int group, int i) const
Definition: pmesh.hpp:400
int AddBdrElement(Element *elem)
Definition: mesh.cpp:1816
T & Last()
Return the last element in the array.
Definition: array.hpp:784
int Dim
Definition: mesh.hpp:66
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
Definition: pncmesh.cpp:1625
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
Geometry::Type GetFaceGeometryType(int Face) const
Definition: mesh.cpp:1413
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:929
double AggregateError(const Array< double > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
Definition: mesh.cpp:9145
void GetMarkedFace(const int face, int *fv)
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:932
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:11898
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:415
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
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:5392
int NumberOfEntries() const
Definition: table.hpp:266
static ParMesh MakeSimplicial(ParMesh &orig_mesh)
Definition: pmesh.cpp:1374
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:5568
void PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4914
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:608
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
Array< int > svert_lvert
Shared to local index mapping.
Definition: pmesh.hpp:75
void GenerateFaces()
Definition: mesh.cpp:6573
const NCList & GetSharedFaces()
Definition: pncmesh.hpp:125
int MyRank
Definition: pmesh.hpp:38
void SaveAsOne(const char *fname, int precision=16) const
Definition: pmesh.cpp:5549
int spaceDim
Definition: mesh.hpp:67
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
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:2547
int RowSize(int i) const
Definition: table.hpp:108
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:5332
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
Table * GetFaceToAllElementTable() const
Definition: pmesh.cpp:2905
void FinalizeParTopo()
Definition: pmesh.cpp:887
void DeleteFaceNbrData()
Definition: pmesh.cpp:2035
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
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
Definition: mesh.cpp:9241
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:9571
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 DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
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:699
void ParPrint(std::ostream &out) const
Save the mesh in a parallel mesh format.
Definition: pmesh.cpp:6238
Class for parallel meshes.
Definition: pmesh.hpp:32
void PrintAsSerial(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:5223
Abstract data type element.
Definition: element.hpp:28
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1480
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:2857
const Table & ElementToEdgeTable() const
Definition: mesh.cpp:6440
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:9779
Table group_squad
Definition: pmesh.hpp:72
void GenerateNCFaceInfo()
Definition: mesh.cpp:6674
Array< int > RefGeoms
Definition: geom.hpp:283
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1109
void Save(const char *fname, int precision=16) const override
Definition: pmesh.cpp:4886
void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[], Array< HYPRE_BigInt > *offsets[]) const
Definition: pmesh.cpp:1933
virtual Type GetType() const =0
Returns element&#39;s type.
virtual void LimitNCLevel(int max_nc_level)
Parallel version of NCMesh::LimitNCLevel.
Definition: pncmesh.cpp:1324
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:268
class Mesh * mesh
The Mesh object containing the element.
Definition: eltrans.hpp:84
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
Definition: mesh.cpp:6079
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1044
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:172
void Load(std::istream &in)
Load the data from a stream.
void UpdateNodes()
Update the nodes of a curved mesh after refinement.
Definition: mesh.cpp:7991
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:49
MFEM_DEPRECATED void ReorientTetMesh() override
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:3200
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:6148
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1464
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
Definition: mesh.hpp:360
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
void Printer(std::ostream &out=mfem::out, std::string section_delimiter="") const
Definition: mesh.cpp:10155
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
int NumOfFaces
Definition: mesh.hpp:70