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