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