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