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