MFEM  v3.0
 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.googlecode.com.
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 <iostream>
21 using namespace std;
22 
23 namespace mfem
24 {
25 
26 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
27  int part_method)
28  : gtopo(comm)
29 {
30  int i, j;
31  int *partitioning;
32  Array<bool> activeBdrElem;
33 
34  MyComm = comm;
35  MPI_Comm_size(MyComm, &NRanks);
36  MPI_Comm_rank(MyComm, &MyRank);
37 
38  Dim = mesh.Dim;
39  spaceDim = mesh.spaceDim;
40 
41  if (partitioning_)
42  partitioning = partitioning_;
43  else
44  partitioning = mesh.GeneratePartitioning(NRanks, part_method);
45 
46  // re-enumerate the partitions to better map to actual processor
47  // interconnect topology !?
48 
49  Array<int> vert;
50  Array<int> vert_global_local(mesh.GetNV());
51  int vert_counter, element_counter, bdrelem_counter;
52 
53  // build vert_global_local
54  vert_global_local = -1;
55 
56  element_counter = 0;
57  vert_counter = 0;
58  for (i = 0; i < mesh.GetNE(); i++)
59  if (partitioning[i] == MyRank)
60  {
61  mesh.GetElementVertices(i, vert);
62  element_counter++;
63  for (j = 0; j < vert.Size(); j++)
64  if (vert_global_local[vert[j]] < 0)
65  vert_global_local[vert[j]] = vert_counter++;
66  }
67 
68  NumOfVertices = vert_counter;
69  NumOfElements = element_counter;
70  vertices.SetSize(NumOfVertices);
71 
72  // re-enumerate the local vertices to preserve the global ordering
73  for (i = vert_counter = 0; i < vert_global_local.Size(); i++)
74  if (vert_global_local[i] >= 0)
75  vert_global_local[i] = vert_counter++;
76 
77  // determine vertices
78  for (i = 0; i < vert_global_local.Size(); i++)
79  if (vert_global_local[i] >= 0)
80  vertices[vert_global_local[i]].SetCoords(mesh.GetVertex(i));
81 
82  // determine elements
83  element_counter = 0;
84  elements.SetSize(NumOfElements);
85  for (i = 0; i < mesh.GetNE(); i++)
86  if (partitioning[i] == MyRank)
87  {
88  elements[element_counter] = mesh.GetElement(i)->Duplicate(this);
89  int *v = elements[element_counter]->GetVertices();
90  int nv = elements[element_counter]->GetNVertices();
91  for (j = 0; j < nv; j++)
92  v[j] = vert_global_local[v[j]];
93  element_counter++;
94  }
95 
96  Table *edge_element = NULL;
97  if (mesh.NURBSext)
98  {
99  activeBdrElem.SetSize(mesh.GetNBE());
100  activeBdrElem = false;
101  }
102  // build boundary elements
103  if (Dim == 3)
104  {
105  NumOfBdrElements = 0;
106  for (i = 0; i < mesh.GetNBE(); i++)
107  {
108  int face, o, el1, el2;
109  mesh.GetBdrElementFace(i, &face, &o);
110  mesh.GetFaceElements(face, &el1, &el2);
111  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
112  {
113  NumOfBdrElements++;
114  if (mesh.NURBSext)
115  activeBdrElem[i] = true;
116  }
117  }
118 
119  bdrelem_counter = 0;
120  boundary.SetSize(NumOfBdrElements);
121  for (i = 0; i < mesh.GetNBE(); i++)
122  {
123  int face, o, el1, el2;
124  mesh.GetBdrElementFace(i, &face, &o);
125  mesh.GetFaceElements(face, &el1, &el2);
126  if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
127  {
128  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
129  int *v = boundary[bdrelem_counter]->GetVertices();
130  int nv = boundary[bdrelem_counter]->GetNVertices();
131  for (j = 0; j < nv; j++)
132  v[j] = vert_global_local[v[j]];
133  bdrelem_counter++;
134  }
135  }
136  }
137  else if (Dim == 2)
138  {
139  edge_element = new Table;
140  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
141 
142  NumOfBdrElements = 0;
143  for (i = 0; i < mesh.GetNBE(); i++)
144  {
145  int edge = mesh.GetBdrElementEdgeIndex(i);
146  int el1 = edge_element->GetRow(edge)[0];
147  if (partitioning[el1] == MyRank)
148  {
149  NumOfBdrElements++;
150  if (mesh.NURBSext)
151  activeBdrElem[i] = true;
152  }
153  }
154 
155  bdrelem_counter = 0;
156  boundary.SetSize(NumOfBdrElements);
157  for (i = 0; i < mesh.GetNBE(); i++)
158  {
159  int edge = mesh.GetBdrElementEdgeIndex(i);
160  int el1 = edge_element->GetRow(edge)[0];
161  if (partitioning[el1] == MyRank)
162  {
163  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
164  int *v = boundary[bdrelem_counter]->GetVertices();
165  int nv = boundary[bdrelem_counter]->GetNVertices();
166  for (j = 0; j < nv; j++)
167  v[j] = vert_global_local[v[j]];
168  bdrelem_counter++;
169  }
170  }
171  }
172  else if (Dim == 1)
173  {
174  NumOfBdrElements = 0;
175  for (i = 0; i < mesh.GetNBE(); i++)
176  {
177  int vert = mesh.boundary[i]->GetVertices()[0];
178  int el1, el2;
179  mesh.GetFaceElements(vert, &el1, &el2);
180  if (partitioning[el1] == MyRank)
181  {
182  NumOfBdrElements++;
183  }
184  }
185 
186  bdrelem_counter = 0;
187  boundary.SetSize(NumOfBdrElements);
188  for (i = 0; i < mesh.GetNBE(); i++)
189  {
190  int vert = mesh.boundary[i]->GetVertices()[0];
191  int el1, el2;
192  mesh.GetFaceElements(vert, &el1, &el2);
193  if (partitioning[el1] == MyRank)
194  {
195  boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
196  int *v = boundary[bdrelem_counter]->GetVertices();
197  v[0] = vert_global_local[v[0]];
198  bdrelem_counter++;
199  }
200  }
201  }
202 
203  meshgen = mesh.MeshGenerator();
204 
205  mesh.attributes.Copy(attributes);
206  mesh.bdr_attributes.Copy(bdr_attributes);
207 
208  // this is called by the default Mesh constructor
209  // InitTables();
210 
211  if (Dim > 1)
212  {
213  el_to_edge = new Table;
214  NumOfEdges = Mesh::GetElementToEdgeTable(*el_to_edge, be_to_edge);
215  }
216  else
217  NumOfEdges = 0;
218 
219  STable3D *faces_tbl = NULL;
220  if (Dim == 3)
221  faces_tbl = GetElementToFaceTable(1);
222  else
223  NumOfFaces = 0;
224  GenerateFaces();
225 
226  c_el_to_edge = NULL;
227 
228  ListOfIntegerSets groups;
229  IntegerSet group;
230 
231  // the first group is the local one
232  group.Recreate(1, &MyRank);
233  groups.Insert(group);
234 
235 #ifdef MFEM_DEBUG
236  if (Dim < 3 && mesh.GetNFaces() != 0)
237  {
238  cerr << "ParMesh::ParMesh (proc " << MyRank << ") : "
239  "(Dim < 3 && mesh.GetNFaces() != 0) is true!" << endl;
240  mfem_error();
241  }
242 #endif
243  // determine shared faces
244  int sface_counter = 0;
245  Array<int> face_group(mesh.GetNFaces());
246  for (i = 0; i < face_group.Size(); i++)
247  {
248  int el[2];
249  face_group[i] = -1;
250  mesh.GetFaceElements(i, &el[0], &el[1]);
251  if (el[1] >= 0)
252  {
253  el[0] = partitioning[el[0]];
254  el[1] = partitioning[el[1]];
255  if ((el[0] == MyRank && el[1] != MyRank) ||
256  (el[0] != MyRank && el[1] == MyRank))
257  {
258  group.Recreate(2, el);
259  face_group[i] = groups.Insert(group) - 1;
260  sface_counter++;
261  }
262  }
263  }
264 
265  // determine shared edges
266  int sedge_counter = 0;
267  if (!edge_element)
268  {
269  edge_element = new Table;
270  if (Dim == 1)
271  edge_element->SetDims(0,0);
272  else
273  Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
274  }
275  for (i = 0; i < edge_element->Size(); i++)
276  {
277  int me = 0, others = 0;
278  for (j = edge_element->GetI()[i]; j < edge_element->GetI()[i+1]; j++)
279  {
280  edge_element->GetJ()[j] = partitioning[edge_element->GetJ()[j]];
281  if (edge_element->GetJ()[j] == MyRank)
282  me = 1;
283  else
284  others = 1;
285  }
286 
287  if (me && others)
288  {
289  sedge_counter++;
290  group.Recreate(edge_element->RowSize(i), edge_element->GetRow(i));
291  edge_element->GetRow(i)[0] = groups.Insert(group) - 1;
292  }
293  else
294  edge_element->GetRow(i)[0] = -1;
295  }
296 
297  // determine shared vertices
298  int svert_counter = 0;
299  Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
300 
301  for (i = 0; i < vert_element->Size(); i++)
302  {
303  int me = 0, others = 0;
304  for (j = vert_element->GetI()[i]; j < vert_element->GetI()[i+1]; j++)
305  {
306  vert_element->GetJ()[j] = partitioning[vert_element->GetJ()[j]];
307  if (vert_element->GetJ()[j] == MyRank)
308  me = 1;
309  else
310  others = 1;
311  }
312 
313  if (me && others)
314  {
315  svert_counter++;
316  group.Recreate(vert_element->RowSize(i), vert_element->GetRow(i));
317  vert_element->GetI()[i] = groups.Insert(group) - 1;
318  }
319  else
320  vert_element->GetI()[i] = -1;
321  }
322 
323  // build group_sface
324  group_sface.MakeI(groups.Size()-1);
325 
326  for (i = 0; i < face_group.Size(); i++)
327  if (face_group[i] >= 0)
328  group_sface.AddAColumnInRow(face_group[i]);
329 
330  group_sface.MakeJ();
331 
332  sface_counter = 0;
333  for (i = 0; i < face_group.Size(); i++)
334  if (face_group[i] >= 0)
335  group_sface.AddConnection(face_group[i], sface_counter++);
336 
337  group_sface.ShiftUpI();
338 
339  // build group_sedge
340  group_sedge.MakeI(groups.Size()-1);
341 
342  for (i = 0; i < edge_element->Size(); i++)
343  if (edge_element->GetRow(i)[0] >= 0)
344  group_sedge.AddAColumnInRow(edge_element->GetRow(i)[0]);
345 
346  group_sedge.MakeJ();
347 
348  sedge_counter = 0;
349  for (i = 0; i < edge_element->Size(); i++)
350  if (edge_element->GetRow(i)[0] >= 0)
351  group_sedge.AddConnection(edge_element->GetRow(i)[0],
352  sedge_counter++);
353 
354  group_sedge.ShiftUpI();
355 
356  // build group_svert
357  group_svert.MakeI(groups.Size()-1);
358 
359  for (i = 0; i < vert_element->Size(); i++)
360  if (vert_element->GetI()[i] >= 0)
361  group_svert.AddAColumnInRow(vert_element->GetI()[i]);
362 
363  group_svert.MakeJ();
364 
365  svert_counter = 0;
366  for (i = 0; i < vert_element->Size(); i++)
367  if (vert_element->GetI()[i] >= 0)
368  group_svert.AddConnection(vert_element->GetI()[i],
369  svert_counter++);
370 
371  group_svert.ShiftUpI();
372 
373  // build shared_faces and sface_lface
374  shared_faces.SetSize(sface_counter);
375  sface_lface. SetSize(sface_counter);
376 
377  if (Dim == 3)
378  {
379  sface_counter = 0;
380  for (i = 0; i < face_group.Size(); i++)
381  if (face_group[i] >= 0)
382  {
383  shared_faces[sface_counter] = mesh.GetFace(i)->Duplicate(this);
384  int *v = shared_faces[sface_counter]->GetVertices();
385  int nv = shared_faces[sface_counter]->GetNVertices();
386  for (j = 0; j < nv; j++)
387  v[j] = vert_global_local[v[j]];
388  switch (shared_faces[sface_counter]->GetType())
389  {
390  case Element::TRIANGLE:
391  sface_lface[sface_counter] = (*faces_tbl)(v[0], v[1], v[2]);
392  // mark the shared face for refinement by reorienting
393  // it according to the refinement flag in the tetradron
394  // to which this shared face belongs to.
395  {
396  int lface = sface_lface[sface_counter];
397  Tetrahedron *tet =
398  (Tetrahedron *)(elements[faces_info[lface].Elem1No]);
399  int re[2], type, flag, *tv;
400  tet->ParseRefinementFlag(re, type, flag);
401  tv = tet->GetVertices();
402  switch (faces_info[lface].Elem1Inf/64)
403  {
404  case 0:
405  switch (re[1])
406  {
407  case 1: v[0] = tv[1]; v[1] = tv[2]; v[2] = tv[3]; break;
408  case 4: v[0] = tv[3]; v[1] = tv[1]; v[2] = tv[2]; break;
409  case 5: v[0] = tv[2]; v[1] = tv[3]; v[2] = tv[1]; break;
410  }
411  break;
412  case 1:
413  switch (re[0])
414  {
415  case 2: v[0] = tv[2]; v[1] = tv[0]; v[2] = tv[3]; break;
416  case 3: v[0] = tv[0]; v[1] = tv[3]; v[2] = tv[2]; break;
417  case 5: v[0] = tv[3]; v[1] = tv[2]; v[2] = tv[0]; break;
418  }
419  break;
420  case 2:
421  v[0] = tv[0]; v[1] = tv[1]; v[2] = tv[3];
422  break;
423  case 3:
424  v[0] = tv[1]; v[1] = tv[0]; v[2] = tv[2];
425  break;
426  }
427  // flip the shared face in the processor that owns the
428  // second element (in 'mesh')
429  {
430  int gl_el1, gl_el2;
431  mesh.GetFaceElements(i, &gl_el1, &gl_el2);
432  if (MyRank == partitioning[gl_el2])
433  {
434  const int t = v[0]; v[0] = v[1]; v[1] = t;
435  }
436  }
437  }
438  break;
439  case Element::QUADRILATERAL:
440  sface_lface[sface_counter] =
441  (*faces_tbl)(v[0], v[1], v[2], v[3]);
442  break;
443  }
444  sface_counter++;
445  }
446 
447  delete faces_tbl;
448  }
449 
450  // build shared_edges and sedge_ledge
451  shared_edges.SetSize(sedge_counter);
452  sedge_ledge. SetSize(sedge_counter);
453 
454  {
455  DSTable v_to_v(NumOfVertices);
456  GetVertexToVertexTable(v_to_v);
457 
458  sedge_counter = 0;
459  for (i = 0; i < edge_element->Size(); i++)
460  if (edge_element->GetRow(i)[0] >= 0)
461  {
462  mesh.GetEdgeVertices(i, vert);
463 
464  shared_edges[sedge_counter] =
465  new Segment(vert_global_local[vert[0]],
466  vert_global_local[vert[1]], 1);
467 
468  if ((sedge_ledge[sedge_counter] =
469  v_to_v(vert_global_local[vert[0]],
470  vert_global_local[vert[1]])) < 0)
471  {
472  cerr << "\n\n\n" << MyRank << ": ParMesh::ParMesh: "
473  << "ERROR in v_to_v\n\n" << endl;
474  mfem_error();
475  }
476 
477  sedge_counter++;
478  }
479  }
480 
481  delete edge_element;
482 
483  // build svert_lvert
484  svert_lvert.SetSize(svert_counter);
485 
486  svert_counter = 0;
487  for (i = 0; i < vert_element->Size(); i++)
488  if (vert_element->GetI()[i] >= 0)
489  svert_lvert[svert_counter++] = vert_global_local[i];
490 
491  delete vert_element;
492 
493  // build the group communication topology
494  gtopo.Create(groups, 822);
495 
496  if (mesh.NURBSext)
497  {
498  NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
499  activeBdrElem);
500  }
501 
502  if (mesh.GetNodes()) // curved mesh
503  {
504  Nodes = new ParGridFunction(this, mesh.GetNodes());
505  own_nodes = 1;
506 
507  Array<int> gvdofs, lvdofs;
508  Vector lnodes;
509  element_counter = 0;
510  for (i = 0; i < mesh.GetNE(); i++)
511  if (partitioning[i] == MyRank)
512  {
513  Nodes->FESpace()->GetElementVDofs(element_counter, lvdofs);
514  mesh.GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
515  mesh.GetNodes()->GetSubVector(gvdofs, lnodes);
516  Nodes->SetSubVector(lvdofs, lnodes);
517  element_counter++;
518  }
519  }
520 
521  if (partitioning_ == NULL)
522  delete [] partitioning;
523 
524  have_face_nbr_data = false;
525 }
526 
527 void ParMesh::GroupEdge(int group, int i, int &edge, int &o)
528 {
529  int sedge = group_sedge.GetJ()[group_sedge.GetI()[group-1]+i];
530  edge = sedge_ledge[sedge];
531  int *v = shared_edges[sedge]->GetVertices();
532  o = (v[0] < v[1]) ? (+1) : (-1);
533 }
534 
535 void ParMesh::GroupFace(int group, int i, int &face, int &o)
536 {
537  int sface = group_sface.GetJ()[group_sface.GetI()[group-1]+i];
538  face = sface_lface[sface];
539  // face gives the base orientation
540  if (faces[face]->GetType() == Element::TRIANGLE)
541  o = GetTriOrientation(faces[face]->GetVertices(),
542  shared_faces[sface]->GetVertices());
543  if (faces[face]->GetType() == Element::QUADRILATERAL)
544  o = GetQuadOrientation(faces[face]->GetVertices(),
545  shared_faces[sface]->GetVertices());
546 }
547 
548 // For a line segment with vertices v[0] and v[1], return a number with
549 // the following meaning:
550 // 0 - the edge was not refined
551 // 1 - the edge e was refined once by splitting v[0],v[1]
552 int ParMesh::GetEdgeSplittings(Element *edge, const DSTable &v_to_v,
553  int *middle)
554 {
555  int m, *v = edge->GetVertices();
556 
557  if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
558  return 1;
559  else
560  return 0;
561 }
562 
563 // For a triangular face with (correctly ordered) vertices v[0], v[1], v[2]
564 // return a number with the following meaning:
565 // 0 - the face was not refined
566 // 1 - the face was refined once by splitting v[0],v[1]
567 // 2 - the face was refined twice by splitting v[0],v[1] and then v[1],v[2]
568 // 3 - the face was refined twice by splitting v[0],v[1] and then v[0],v[2]
569 // 4 - the face was refined three times (as in 2+3)
570 int ParMesh::GetFaceSplittings(Element *face, const DSTable &v_to_v,
571  int *middle)
572 {
573  int m, right = 0;
574  int number_of_splittings = 0;
575  int *v = face->GetVertices();
576 
577  if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
578  {
579  number_of_splittings++;
580  if ((m = v_to_v(v[1], v[2])) != -1 && middle[m] != -1)
581  {
582  right = 1;
583  number_of_splittings++;
584  }
585  if ((m = v_to_v(v[2], v[0])) != -1 && middle[m] != -1)
586  number_of_splittings++;
587 
588  switch (number_of_splittings)
589  {
590  case 2:
591  if (right == 0)
592  number_of_splittings++;
593  break;
594  case 3:
595  number_of_splittings++;
596  break;
597  }
598  }
599 
600  return number_of_splittings;
601 }
602 
603 void ParMesh::GetFaceNbrElementTransformation(
604  int i, IsoparametricTransformation *ElTr)
605 {
606  DenseMatrix &pointmat = ElTr->GetPointMat();
607  Element *elem = face_nbr_elements[i];
608 
609  ElTr->Attribute = elem->GetAttribute();
610  ElTr->ElementNo = NumOfElements + i;
611 
612  if (Nodes == NULL)
613  {
614  const int nv = elem->GetNVertices();
615  const int *v = elem->GetVertices();
616 
617  pointmat.SetSize(spaceDim, nv);
618  for (int k = 0; k < spaceDim; k++)
619  for (int j = 0; j < nv; j++)
620  pointmat(k, j) = face_nbr_vertices[v[j]](k);
621 
622  ElTr->SetFE(GetTransformationFEforElementType(elem->GetType()));
623  }
624  else
625  {
626  Array<int> vdofs;
627  ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
628  if (pNodes)
629  {
630  pNodes->ParFESpace()->GetFaceNbrElementVDofs(i, vdofs);
631  int n = vdofs.Size()/spaceDim;
632  pointmat.SetSize(spaceDim, n);
633  for (int k = 0; k < spaceDim; k++)
634  for (int j = 0; j < n; j++)
635  pointmat(k,j) = (pNodes->FaceNbrData())(vdofs[n*k+j]);
636 
637  ElTr->SetFE(pNodes->ParFESpace()->GetFaceNbrFE(i));
638  }
639  else
640  mfem_error("ParMesh::GetFaceNbrElementTransformation : "
641  "Nodes are not ParGridFunction!");
642  }
643 }
644 
645 void ParMesh::DeleteFaceNbrData()
646 {
647  if (!have_face_nbr_data)
648  return;
649 
650  have_face_nbr_data = false;
651  face_nbr_group.DeleteAll();
652  face_nbr_elements_offset.DeleteAll();
653  face_nbr_vertices_offset.DeleteAll();
654  for (int i = 0; i < face_nbr_elements.Size(); i++)
655  FreeElement(face_nbr_elements[i]);
656  face_nbr_elements.DeleteAll();
657  face_nbr_vertices.DeleteAll();
658  send_face_nbr_elements.Clear();
659  send_face_nbr_vertices.Clear();
660 }
661 
662 void ParMesh::ExchangeFaceNbrData()
663 {
664  if (have_face_nbr_data)
665  return;
666 
667  Table *gr_sface;
668  int *s2l_face;
669  if (Dim == 1)
670  {
671  gr_sface = &group_svert;
672  s2l_face = svert_lvert;
673  }
674  else if (Dim == 2)
675  {
676  gr_sface = &group_sedge;
677  s2l_face = sedge_ledge;
678  }
679  else
680  {
681  gr_sface = &group_sface;
682  s2l_face = sface_lface;
683  }
684 
685  int num_face_nbrs = 0;
686  for (int g = 1; g < GetNGroups(); g++)
687  if (gr_sface->RowSize(g-1) > 0)
688  num_face_nbrs++;
689 
690  face_nbr_group.SetSize(num_face_nbrs);
691 
692  if (num_face_nbrs == 0)
693  {
694  have_face_nbr_data = true;
695  return;
696  }
697 
698  {
699  // sort face-neighbors by processor rank
700  Array<Pair<int, int> > rank_group(num_face_nbrs);
701 
702  for (int g = 1, counter = 0; g < GetNGroups(); g++)
703  if (gr_sface->RowSize(g-1) > 0)
704  {
705 #ifdef MFEM_DEBUG
706  if (gtopo.GetGroupSize(g) != 2)
707  mfem_error("ParMesh::ExchangeFaceNbrData() : "
708  "group size is not 2!");
709 #endif
710  const int *nbs = gtopo.GetGroup(g);
711  int lproc = (nbs[0]) ? nbs[0] : nbs[1];
712  rank_group[counter].one = gtopo.GetNeighborRank(lproc);
713  rank_group[counter].two = g;
714  counter++;
715  }
716 
717  SortPairs<int, int>(rank_group, rank_group.Size());
718 
719  for (int fn = 0; fn < num_face_nbrs; fn++)
720  face_nbr_group[fn] = rank_group[fn].two;
721  }
722 
723  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
724  MPI_Request *send_requests = requests;
725  MPI_Request *recv_requests = requests + num_face_nbrs;
726  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
727 
728  int *nbr_data = new int[6*num_face_nbrs];
729  int *nbr_send_data = nbr_data;
730  int *nbr_recv_data = nbr_data + 3*num_face_nbrs;
731 
732  Array<int> el_marker(GetNE());
733  Array<int> vertex_marker(GetNV());
734  el_marker = -1;
735  vertex_marker = -1;
736 
737  Table send_face_nbr_elemdata, send_face_nbr_facedata;
738 
739  send_face_nbr_elements.MakeI(num_face_nbrs);
740  send_face_nbr_vertices.MakeI(num_face_nbrs);
741  send_face_nbr_elemdata.MakeI(num_face_nbrs);
742  send_face_nbr_facedata.MakeI(num_face_nbrs);
743  for (int fn = 0; fn < num_face_nbrs; fn++)
744  {
745  int nbr_group = face_nbr_group[fn];
746  int num_sfaces = gr_sface->RowSize(nbr_group-1);
747  int *sface = gr_sface->GetRow(nbr_group-1);
748  for (int i = 0; i < num_sfaces; i++)
749  {
750  int lface = s2l_face[sface[i]];
751  int el = faces_info[lface].Elem1No;
752  if (el_marker[el] != fn)
753  {
754  el_marker[el] = fn;
755  send_face_nbr_elements.AddAColumnInRow(fn);
756 
757  const int nv = elements[el]->GetNVertices();
758  const int *v = elements[el]->GetVertices();
759  for (int j = 0; j < nv; j++)
760  if (vertex_marker[v[j]] != fn)
761  {
762  vertex_marker[v[j]] = fn;
763  send_face_nbr_vertices.AddAColumnInRow(fn);
764  }
765 
766  send_face_nbr_elemdata.AddColumnsInRow(fn, nv + 2);
767  }
768  }
769  send_face_nbr_facedata.AddColumnsInRow(fn, 2*num_sfaces);
770 
771  nbr_send_data[3*fn ] = send_face_nbr_elements.GetI()[fn];
772  nbr_send_data[3*fn+1] = send_face_nbr_vertices.GetI()[fn];
773  nbr_send_data[3*fn+2] = send_face_nbr_elemdata.GetI()[fn];
774 
775  int nbr_rank = GetFaceNbrRank(fn);
776  int tag = 0;
777 
778  MPI_Isend(&nbr_send_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
779  &send_requests[fn]);
780  MPI_Irecv(&nbr_recv_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
781  &recv_requests[fn]);
782  }
783  send_face_nbr_elements.MakeJ();
784  send_face_nbr_vertices.MakeJ();
785  send_face_nbr_elemdata.MakeJ();
786  send_face_nbr_facedata.MakeJ();
787  el_marker = -1;
788  vertex_marker = -1;
789  for (int fn = 0; fn < num_face_nbrs; fn++)
790  {
791  int nbr_group = face_nbr_group[fn];
792  int num_sfaces = gr_sface->RowSize(nbr_group-1);
793  int *sface = gr_sface->GetRow(nbr_group-1);
794  for (int i = 0; i < num_sfaces; i++)
795  {
796  int lface = s2l_face[sface[i]];
797  int el = faces_info[lface].Elem1No;
798  if (el_marker[el] != fn)
799  {
800  el_marker[el] = fn;
801  send_face_nbr_elements.AddConnection(fn, el);
802 
803  const int nv = elements[el]->GetNVertices();
804  const int *v = elements[el]->GetVertices();
805  for (int j = 0; j < nv; j++)
806  if (vertex_marker[v[j]] != fn)
807  {
808  vertex_marker[v[j]] = fn;
809  send_face_nbr_vertices.AddConnection(fn, v[j]);
810  }
811 
812  send_face_nbr_elemdata.AddConnection(fn, GetAttribute(el));
813  send_face_nbr_elemdata.AddConnection(
814  fn, GetElementBaseGeometry(el));
815  send_face_nbr_elemdata.AddConnections(fn, v, nv);
816  }
817  send_face_nbr_facedata.AddConnection(fn, el);
818  int info = faces_info[lface].Elem1Inf;
819  // change the orientation in info to be relative to the shared face
820  // in 1D and 2D keep the orientation equal to 0
821  if (Dim == 3)
822  {
823  Element *lf = faces[lface];
824  const int *sf_v = shared_faces[sface[i]]->GetVertices();
825 
826  if (lf->GetGeometryType() == Geometry::TRIANGLE)
827  info += GetTriOrientation(sf_v, lf->GetVertices());
828  else
829  info += GetQuadOrientation(sf_v, lf->GetVertices());
830  }
831  send_face_nbr_facedata.AddConnection(fn, info);
832  }
833  }
834  send_face_nbr_elements.ShiftUpI();
835  send_face_nbr_vertices.ShiftUpI();
836  send_face_nbr_elemdata.ShiftUpI();
837  send_face_nbr_facedata.ShiftUpI();
838 
839  // convert the vertex indices in send_face_nbr_elemdata
840  // convert the element indices in send_face_nbr_facedata
841  for (int fn = 0; fn < num_face_nbrs; fn++)
842  {
843  int num_elems = send_face_nbr_elements.RowSize(fn);
844  int *elems = send_face_nbr_elements.GetRow(fn);
845  int num_verts = send_face_nbr_vertices.RowSize(fn);
846  int *verts = send_face_nbr_vertices.GetRow(fn);
847  int *elemdata = send_face_nbr_elemdata.GetRow(fn);
848  int num_sfaces = send_face_nbr_facedata.RowSize(fn)/2;
849  int *facedata = send_face_nbr_facedata.GetRow(fn);
850 
851  for (int i = 0; i < num_verts; i++)
852  vertex_marker[verts[i]] = i;
853 
854  for (int el = 0; el < num_elems; el++)
855  {
856  const int nv = elements[el]->GetNVertices();
857  elemdata += 2; // skip the attribute and the geometry type
858  for (int j = 0; j < nv; j++)
859  elemdata[j] = vertex_marker[elemdata[j]];
860  elemdata += nv;
861 
862  el_marker[elems[el]] = el;
863  }
864 
865  for (int i = 0; i < num_sfaces; i++)
866  facedata[2*i] = el_marker[facedata[2*i]];
867  }
868 
869  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
870 
871  Array<int> recv_face_nbr_facedata;
872  Table recv_face_nbr_elemdata;
873 
874  // fill-in face_nbr_elements_offset, face_nbr_vertices_offset
875  face_nbr_elements_offset.SetSize(num_face_nbrs + 1);
876  face_nbr_vertices_offset.SetSize(num_face_nbrs + 1);
877  recv_face_nbr_elemdata.MakeI(num_face_nbrs);
878  face_nbr_elements_offset[0] = 0;
879  face_nbr_vertices_offset[0] = 0;
880  for (int fn = 0; fn < num_face_nbrs; fn++)
881  {
882  face_nbr_elements_offset[fn+1] =
883  face_nbr_elements_offset[fn] + nbr_recv_data[3*fn];
884  face_nbr_vertices_offset[fn+1] =
885  face_nbr_vertices_offset[fn] + nbr_recv_data[3*fn+1];
886  recv_face_nbr_elemdata.AddColumnsInRow(fn, nbr_recv_data[3*fn+2]);
887  }
888  recv_face_nbr_elemdata.MakeJ();
889 
890  MPI_Waitall(num_face_nbrs, send_requests, statuses);
891 
892  // send and receive the element data
893  for (int fn = 0; fn < num_face_nbrs; fn++)
894  {
895  int nbr_rank = GetFaceNbrRank(fn);
896  int tag = 0;
897 
898  MPI_Isend(send_face_nbr_elemdata.GetRow(fn),
899  send_face_nbr_elemdata.RowSize(fn),
900  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
901 
902  MPI_Irecv(recv_face_nbr_elemdata.GetRow(fn),
903  recv_face_nbr_elemdata.RowSize(fn),
904  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
905  }
906 
907  // convert the element data into face_nbr_elements
908  face_nbr_elements.SetSize(face_nbr_elements_offset[num_face_nbrs]);
909  while (true)
910  {
911  int fn;
912  MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
913 
914  if (fn == MPI_UNDEFINED)
915  break;
916 
917  int vert_off = face_nbr_vertices_offset[fn];
918  int elem_off = face_nbr_elements_offset[fn];
919  int num_elems = face_nbr_elements_offset[fn+1] - elem_off;
920  int *recv_elemdata = recv_face_nbr_elemdata.GetRow(fn);
921 
922  for (int i = 0; i < num_elems; i++)
923  {
924  Element *el = NewElement(recv_elemdata[1]);
925  el->SetAttribute(recv_elemdata[0]);
926  recv_elemdata += 2;
927  int nv = el->GetNVertices();
928  for (int j = 0; j < nv; j++)
929  recv_elemdata[j] += vert_off;
930  el->SetVertices(recv_elemdata);
931  recv_elemdata += nv;
932  face_nbr_elements[elem_off++] = el;
933  }
934  }
935 
936  MPI_Waitall(num_face_nbrs, send_requests, statuses);
937 
938  // send and receive the face data
939  recv_face_nbr_facedata.SetSize(
940  send_face_nbr_facedata.Size_of_connections());
941  for (int fn = 0; fn < num_face_nbrs; fn++)
942  {
943  int nbr_rank = GetFaceNbrRank(fn);
944  int tag = 0;
945 
946  MPI_Isend(send_face_nbr_facedata.GetRow(fn),
947  send_face_nbr_facedata.RowSize(fn),
948  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
949 
950  // the size of the send and receive face data is the same
951  MPI_Irecv(&recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]],
952  send_face_nbr_facedata.RowSize(fn),
953  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
954  }
955 
956  // transfer the received face data into faces_info
957  while (true)
958  {
959  int fn;
960  MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
961 
962  if (fn == MPI_UNDEFINED)
963  break;
964 
965  int elem_off = face_nbr_elements_offset[fn];
966  int nbr_group = face_nbr_group[fn];
967  int num_sfaces = gr_sface->RowSize(nbr_group-1);
968  int *sface = gr_sface->GetRow(nbr_group-1);
969  int *facedata =
970  &recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]];
971 
972  for (int i = 0; i < num_sfaces; i++)
973  {
974  int lface = s2l_face[sface[i]];
975  FaceInfo &face_info = faces_info[lface];
976  face_info.Elem2No = -1 - (facedata[2*i] + elem_off);
977  int info = facedata[2*i+1];
978  // change the orientation in info to be relative to the local face
979  if (Dim < 3)
980  {
981  info++; // orientation 0 --> orientation 1
982  }
983  else
984  {
985  int nbr_ori = info%64, nbr_v[4];
986  Element *lf = faces[lface];
987  const int *sf_v = shared_faces[sface[i]]->GetVertices();
988 
989  if (lf->GetGeometryType() == Geometry::TRIANGLE)
990  {
991  // apply the nbr_ori to sf_v to get nbr_v
992  const int *perm = tri_orientations[nbr_ori];
993  for (int j = 0; j < 3; j++)
994  nbr_v[perm[j]] = sf_v[j];
995  // get the orientation of nbr_v w.r.t. the local face
996  nbr_ori = GetTriOrientation(lf->GetVertices(), nbr_v);
997  }
998  else
999  {
1000  // apply the nbr_ori to sf_v to get nbr_v
1001  const int *perm = quad_orientations[nbr_ori];
1002  for (int j = 0; j < 4; j++)
1003  nbr_v[perm[j]] = sf_v[j];
1004  // get the orientation of nbr_v w.r.t. the local face
1005  nbr_ori = GetQuadOrientation(lf->GetVertices(), nbr_v);
1006  }
1007 
1008  info = 64*(info/64) + nbr_ori;
1009  }
1010  face_info.Elem2Inf = info;
1011  }
1012  }
1013 
1014  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1015 
1016  // allocate the face_nbr_vertices
1017  face_nbr_vertices.SetSize(face_nbr_vertices_offset[num_face_nbrs]);
1018 
1019  delete [] nbr_data;
1020 
1021  delete [] statuses;
1022  delete [] requests;
1023 
1024  have_face_nbr_data = true;
1025 
1026  ExchangeFaceNbrNodes();
1027 }
1028 
1029 void ParMesh::ExchangeFaceNbrNodes()
1030 {
1031  if (!have_face_nbr_data)
1032  {
1033  ExchangeFaceNbrData(); // calls this method at the end
1034  }
1035  else if (Nodes == NULL)
1036  {
1037  int num_face_nbrs = GetNFaceNeighbors();
1038 
1039  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
1040  MPI_Request *send_requests = requests;
1041  MPI_Request *recv_requests = requests + num_face_nbrs;
1042  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
1043 
1044  // allocate buffer and copy the vertices to be sent
1045  Array<Vertex> send_vertices(send_face_nbr_vertices.Size_of_connections());
1046  for (int i = 0; i < send_vertices.Size(); i++)
1047  send_vertices[i] = vertices[send_face_nbr_vertices.GetJ()[i]];
1048 
1049  // send and receive the vertices
1050  for (int fn = 0; fn < num_face_nbrs; fn++)
1051  {
1052  int nbr_rank = GetFaceNbrRank(fn);
1053  int tag = 0;
1054 
1055  MPI_Isend(send_vertices[send_face_nbr_vertices.GetI()[fn]](),
1056  3*send_face_nbr_vertices.RowSize(fn),
1057  MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
1058 
1059  MPI_Irecv(face_nbr_vertices[face_nbr_vertices_offset[fn]](),
1060  3*(face_nbr_vertices_offset[fn+1] -
1061  face_nbr_vertices_offset[fn]),
1062  MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
1063  }
1064 
1065  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1066  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1067 
1068  delete [] statuses;
1069  delete [] requests;
1070  }
1071  else
1072  {
1073  ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
1074  if (pNodes)
1075  pNodes->ExchangeFaceNbrData();
1076  else
1077  mfem_error("ParMesh::ExchangeFaceNbrNodes() : "
1078  "Nodes are not ParGridFunction!");
1079  }
1080 }
1081 
1082 int ParMesh::GetFaceNbrRank(int fn) const
1083 {
1084  int nbr_group = face_nbr_group[fn];
1085  const int *nbs = gtopo.GetGroup(nbr_group);
1086  int nbr_lproc = (nbs[0]) ? nbs[0] : nbs[1];
1087  int nbr_rank = gtopo.GetNeighborRank(nbr_lproc);
1088 
1089  return nbr_rank;
1090 }
1091 
1092 Table *ParMesh::GetFaceToAllElementTable() const
1093 {
1094  const Array<int> *s2l_face;
1095  if (Dim == 1)
1096  s2l_face = &svert_lvert;
1097  else if (Dim == 2)
1098  s2l_face = &sedge_ledge;
1099  else
1100  s2l_face = &sface_lface;
1101 
1102  Table *face_elem = new Table;
1103 
1104  face_elem->MakeI(faces_info.Size());
1105 
1106  for (int i = 0; i < faces_info.Size(); i++)
1107  {
1108  if (faces_info[i].Elem2No >= 0)
1109  face_elem->AddColumnsInRow(i, 2);
1110  else
1111  face_elem->AddAColumnInRow(i);
1112  }
1113  for (int i = 0; i < s2l_face->Size(); i++)
1114  face_elem->AddAColumnInRow((*s2l_face)[i]);
1115 
1116  face_elem->MakeJ();
1117 
1118  for (int i = 0; i < faces_info.Size(); i++)
1119  {
1120  face_elem->AddConnection(i, faces_info[i].Elem1No);
1121  if (faces_info[i].Elem2No >= 0)
1122  face_elem->AddConnection(i, faces_info[i].Elem2No);
1123  }
1124  for (int i = 0; i < s2l_face->Size(); i++)
1125  {
1126  int lface = (*s2l_face)[i];
1127  int nbr_elem_idx = -1 - faces_info[lface].Elem2No;
1128  face_elem->AddConnection(lface, NumOfElements + nbr_elem_idx);
1129  }
1130 
1131  face_elem->ShiftUpI();
1132 
1133  return face_elem;
1134 }
1135 
1136 FaceElementTransformations *ParMesh::GetSharedFaceTransformations(int sf)
1137 {
1138  int FaceNo;
1139 
1140  if (Dim == 1)
1141  FaceNo = svert_lvert[sf];
1142  else if (Dim == 2)
1143  FaceNo = sedge_ledge[sf];
1144  else
1145  FaceNo = sface_lface[sf];
1146 
1147  // fill-in the face and the first (local) element data into FaceElemTr
1148  GetFaceElementTransformations(FaceNo);
1149 
1150  FaceInfo &face_info = faces_info[FaceNo];
1151 
1152  // setup the transformation for the second (neighbor) element
1153  FaceElemTr.Elem2No = -1 - face_info.Elem2No;
1154  GetFaceNbrElementTransformation(FaceElemTr.Elem2No, &Transformation2);
1155  FaceElemTr.Elem2 = &Transformation2;
1156 
1157  // setup Loc2
1158  int face_type = (Dim == 1) ? Element::POINT : faces[FaceNo]->GetType();
1159  switch (face_type)
1160  {
1161  case Element::POINT:
1162  GetLocalPtToSegTransformation(FaceElemTr.Loc2.Transf,
1163  face_info.Elem2Inf);
1164  break;
1165 
1166  case Element::SEGMENT:
1167  if (face_nbr_elements[FaceElemTr.Elem2No]->GetType() == Element::TRIANGLE)
1168  GetLocalSegToTriTransformation(FaceElemTr.Loc2.Transf,
1169  face_info.Elem2Inf);
1170  else // assume the element is a quad
1171  GetLocalSegToQuadTransformation(FaceElemTr.Loc2.Transf,
1172  face_info.Elem2Inf);
1173  break;
1174 
1175  case Element::TRIANGLE:
1176  // --------- assumes the face is a triangle -- face of a tetrahedron
1177  GetLocalTriToTetTransformation(FaceElemTr.Loc2.Transf,
1178  face_info.Elem2Inf);
1179  break;
1180 
1181  case Element::QUADRILATERAL:
1182  // --------- assumes the face is a quad -- face of a hexahedron
1183  GetLocalQuadToHexTransformation(FaceElemTr.Loc2.Transf,
1184  face_info.Elem2Inf);
1185  break;
1186  }
1187 
1188  return &FaceElemTr;
1189 }
1190 
1191 int ParMesh::GetNSharedFaces() const
1192 {
1193  if (Dim == 1)
1194  return svert_lvert.Size();
1195  if (Dim == 2)
1196  return sedge_ledge.Size();
1197  return sface_lface.Size();
1198 }
1199 
1200 void ParMesh::ReorientTetMesh()
1201 {
1202  if (Dim != 3 || !(meshgen & 1))
1203  return;
1204 
1205  Mesh::ReorientTetMesh();
1206 
1207  int *v;
1208 
1209  // The local edge and face numbering is changed therefore we need to
1210  // update sedge_ledge and sface_lface.
1211  {
1212  DSTable v_to_v(NumOfVertices);
1213  GetVertexToVertexTable(v_to_v);
1214  for (int i = 0; i < shared_edges.Size(); i++)
1215  {
1216  v = shared_edges[i]->GetVertices();
1217  sedge_ledge[i] = v_to_v(v[0], v[1]);
1218  }
1219  }
1220 
1221  // Rotate shared faces and update sface_lface.
1222  // Note that no communication is needed to ensure that the shared
1223  // faces are rotated in the same way in both processors. This is
1224  // automatic due to various things, e.g. the global to local vertex
1225  // mapping preserves the global order; also the way new vertices
1226  // are introduced during refinement is essential.
1227  {
1228  STable3D *faces_tbl = GetFacesTable();
1229  for (int i = 0; i < shared_faces.Size(); i++)
1230  if (shared_faces[i]->GetType() == Element::TRIANGLE)
1231  {
1232  v = shared_faces[i]->GetVertices();
1233 
1234  Rotate3(v[0], v[1], v[2]);
1235 
1236  sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
1237  }
1238  delete faces_tbl;
1239  }
1240 }
1241 
1242 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
1243 {
1244  int i, j, wtls = WantTwoLevelState;
1245 
1246  if (Nodes) // curved mesh
1247  {
1248  UseTwoLevelState(1);
1249  }
1250 
1251  SetState(Mesh::NORMAL);
1252  DeleteCoarseTables();
1253  DeleteFaceNbrData();
1254 
1255  if (Dim == 3)
1256  {
1257  if (WantTwoLevelState)
1258  {
1259  c_NumOfVertices = NumOfVertices;
1260  c_NumOfEdges = NumOfEdges;
1261  c_NumOfFaces = NumOfFaces;
1262  c_NumOfElements = NumOfElements;
1263  c_NumOfBdrElements = NumOfBdrElements;
1264  }
1265 
1266  int uniform_refinement = 0;
1267  if (type < 0)
1268  {
1269  type = -type;
1270  uniform_refinement = 1;
1271  }
1272 
1273  // 1. Get table of vertex to vertex connections.
1274  DSTable v_to_v(NumOfVertices);
1275  GetVertexToVertexTable(v_to_v);
1276 
1277  // 2. Get edge to element connections in arrays edge1 and edge2
1278  Array<int> middle(v_to_v.NumberOfEntries());
1279  middle = -1;
1280 
1281  // 3. Do the red refinement.
1282  switch (type)
1283  {
1284  case 1:
1285  for (i = 0; i < marked_el.Size(); i++)
1286  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1287  break;
1288  case 2:
1289  for (i = 0; i < marked_el.Size(); i++)
1290  {
1291  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1292 
1293  Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
1294  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1295  }
1296  break;
1297  case 3:
1298  for (i = 0; i < marked_el.Size(); i++)
1299  {
1300  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1301 
1302  j = NumOfElements - 1;
1303  Bisection(j, v_to_v, NULL, NULL, middle);
1304  Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
1305  Bisection(j, v_to_v, NULL, NULL, middle);
1306 
1307  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1308  Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
1309  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
1310  }
1311  break;
1312  }
1313 
1314  if (WantTwoLevelState)
1315  {
1316  RefinedElement::State = RefinedElement::FINE;
1317  State = Mesh::TWO_LEVEL_FINE;
1318  }
1319 
1320  // 4. Do the green refinement (to get conforming mesh).
1321  int need_refinement;
1322  int refined_edge[5][3] = {{0, 0, 0},
1323  {1, 0, 0},
1324  {1, 1, 0},
1325  {1, 0, 1},
1326  {1, 1, 1}};
1327  int faces_in_group, max_faces_in_group = 0;
1328  // face_splittings identify how the shared faces have been split
1329  int **face_splittings = new int*[GetNGroups()-1];
1330  for (i = 0; i < GetNGroups()-1; i++)
1331  {
1332  faces_in_group = GroupNFaces(i+1);
1333  face_splittings[i] = new int[faces_in_group];
1334  if (faces_in_group > max_faces_in_group)
1335  max_faces_in_group = faces_in_group;
1336  }
1337  int neighbor, *iBuf = new int[max_faces_in_group];
1338 
1339  Array<int> group_faces;
1340  Vertex V;
1341 
1342  MPI_Request request;
1343  MPI_Status status;
1344 
1345 #ifdef MFEM_DEBUG
1346  int ref_loops_all = 0, ref_loops_par = 0;
1347 #endif
1348  do
1349  {
1350  need_refinement = 0;
1351  for (i = 0; i < NumOfElements; i++)
1352  {
1353  if (elements[i]->NeedRefinement(v_to_v, middle))
1354  {
1355  need_refinement = 1;
1356  Bisection(i, v_to_v, NULL, NULL, middle);
1357  }
1358  }
1359 #ifdef MFEM_DEBUG
1360  ref_loops_all++;
1361 #endif
1362 
1363  if (uniform_refinement)
1364  continue;
1365 
1366  // if the mesh is locally conforming start making it globally
1367  // conforming
1368  if (need_refinement == 0)
1369  {
1370 #ifdef MFEM_DEBUG
1371  ref_loops_par++;
1372 #endif
1373  // MPI_Barrier(MyComm);
1374 
1375  // (a) send the type of interface splitting
1376  for (i = 0; i < GetNGroups()-1; i++)
1377  {
1378  group_sface.GetRow(i, group_faces);
1379  faces_in_group = group_faces.Size();
1380  // it is enough to communicate through the faces
1381  if (faces_in_group != 0)
1382  {
1383  for (j = 0; j < faces_in_group; j++)
1384  face_splittings[i][j] =
1385  GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
1386  middle);
1387  const int *nbs = gtopo.GetGroup(i+1);
1388  if (nbs[0] == 0)
1389  neighbor = gtopo.GetNeighborRank(nbs[1]);
1390  else
1391  neighbor = gtopo.GetNeighborRank(nbs[0]);
1392  MPI_Isend(face_splittings[i], faces_in_group, MPI_INT,
1393  neighbor, 0, MyComm, &request);
1394  }
1395  }
1396 
1397  // (b) receive the type of interface splitting
1398  for (i = 0; i < GetNGroups()-1; i++)
1399  {
1400  group_sface.GetRow(i, group_faces);
1401  faces_in_group = group_faces.Size();
1402  if (faces_in_group != 0)
1403  {
1404  const int *nbs = gtopo.GetGroup(i+1);
1405  if (nbs[0] == 0)
1406  neighbor = gtopo.GetNeighborRank(nbs[1]);
1407  else
1408  neighbor = gtopo.GetNeighborRank(nbs[0]);
1409  MPI_Recv(iBuf, faces_in_group, MPI_INT, neighbor,
1410  MPI_ANY_TAG, MyComm, &status);
1411 
1412  for (j = 0; j < faces_in_group; j++)
1413  if (iBuf[j] != face_splittings[i][j])
1414  {
1415  int *v = shared_faces[group_faces[j]]->GetVertices();
1416  for (int k = 0; k < 3; k++)
1417  if (refined_edge[iBuf[j]][k] == 1 &&
1418  refined_edge[face_splittings[i][j]][k] == 0)
1419  {
1420  int ii = v_to_v(v[k], v[(k+1)%3]);
1421  if (middle[ii] == -1)
1422  {
1423  need_refinement = 1;
1424  middle[ii] = NumOfVertices++;
1425  for (int c = 0; c < 3; c++)
1426  V(c) = 0.5 * (vertices[v[k]](c) +
1427  vertices[v[(k+1)%3]](c));
1428  vertices.Append(V);
1429  }
1430  }
1431  }
1432  }
1433  }
1434 
1435  i = need_refinement;
1436  MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
1437  }
1438  }
1439  while (need_refinement == 1);
1440 
1441 #ifdef MFEM_DEBUG
1442  i = ref_loops_all;
1443  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
1444  if (MyRank == 0)
1445  {
1446  cout << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
1447  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
1448  << '\n' << endl;
1449  }
1450 #endif
1451 
1452  delete [] iBuf;
1453  for (i = 0; i < GetNGroups()-1; i++)
1454  delete [] face_splittings[i];
1455  delete [] face_splittings;
1456 
1457 
1458  // 5. Update the boundary elements.
1459  do
1460  {
1461  need_refinement = 0;
1462  for (i = 0; i < NumOfBdrElements; i++)
1463  if (boundary[i]->NeedRefinement(v_to_v, middle))
1464  {
1465  need_refinement = 1;
1466  Bisection(i, v_to_v, middle);
1467  }
1468  }
1469  while (need_refinement == 1);
1470 
1471  if (NumOfBdrElements != boundary.Size())
1472  mfem_error("ParMesh::LocalRefinement :"
1473  " (NumOfBdrElements != boundary.Size())");
1474 
1475  // 5a. Update the groups after refinement.
1476  if (el_to_face != NULL)
1477  {
1478  if (WantTwoLevelState)
1479  {
1480  c_el_to_face = el_to_face;
1481  el_to_face = NULL;
1482  mfem::Swap(faces_info, fc_faces_info);
1483  }
1484  RefineGroups(v_to_v, middle);
1485  // GetElementToFaceTable(); // Called by RefineGroups
1486  GenerateFaces();
1487  if (WantTwoLevelState)
1488  {
1489  f_el_to_face = el_to_face;
1490  }
1491  }
1492 
1493  // 6. Un-mark the Pf elements.
1494  int refinement_edges[2], type, flag;
1495  for (i = 0; i < NumOfElements; i++)
1496  {
1497  Element *El = elements[i];
1498  while (El->GetType() == Element::BISECTED)
1499  El = ((BisectedElement *) El)->FirstChild;
1500  ((Tetrahedron *) El)->ParseRefinementFlag(refinement_edges,
1501  type, flag);
1502  if (type == Tetrahedron::TYPE_PF)
1503  ((Tetrahedron *) El)->CreateRefinementFlag(refinement_edges,
1504  Tetrahedron::TYPE_PU,
1505  flag);
1506  }
1507 
1508  // 7. Free the allocated memory.
1509  middle.DeleteAll();
1510 
1511  if (el_to_edge != NULL)
1512  {
1513  if (WantTwoLevelState)
1514  {
1515  c_el_to_edge = el_to_edge;
1516  f_el_to_edge = new Table;
1517  c_bel_to_edge = bel_to_edge;
1518  bel_to_edge = NULL;
1519  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
1520  el_to_edge = f_el_to_edge;
1521  f_bel_to_edge = bel_to_edge;
1522  }
1523  else
1524  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1525  }
1526 
1527  if (WantTwoLevelState)
1528  {
1529  f_NumOfVertices = NumOfVertices;
1530  f_NumOfEdges = NumOfEdges;
1531  f_NumOfFaces = NumOfFaces;
1532  f_NumOfElements = NumOfElements;
1533  f_NumOfBdrElements = NumOfBdrElements;
1534  }
1535  } // 'if (Dim == 3)'
1536 
1537 
1538  if (Dim == 2)
1539  {
1540  if (WantTwoLevelState)
1541  {
1542  c_NumOfVertices = NumOfVertices;
1543  c_NumOfEdges = NumOfEdges;
1544  c_NumOfElements = NumOfElements;
1545  c_NumOfBdrElements = NumOfBdrElements;
1546  }
1547 
1548  int uniform_refinement = 0;
1549  if (type < 0)
1550  {
1551  type = -type;
1552  uniform_refinement = 1;
1553  }
1554 
1555  // 1. Get table of vertex to vertex connections.
1556  DSTable v_to_v(NumOfVertices);
1557  GetVertexToVertexTable(v_to_v);
1558 
1559  // 2. Get edge to element connections in arrays edge1 and edge2
1560  int nedges = v_to_v.NumberOfEntries();
1561  int *edge1 = new int[nedges];
1562  int *edge2 = new int[nedges];
1563  int *middle = new int[nedges];
1564 
1565  for (i = 0; i < nedges; i++)
1566  edge1[i] = edge2[i] = middle[i] = -1;
1567 
1568  for (i = 0; i < NumOfElements; i++)
1569  {
1570  int *v = elements[i]->GetVertices();
1571  for (j = 0; j < 3; j++)
1572  {
1573  int ind = v_to_v(v[j], v[(j+1)%3]);
1574  (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
1575  }
1576  }
1577 
1578  // 3. Do the red refinement.
1579  for (i = 0; i < marked_el.Size(); i++)
1580  RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
1581 
1582  if (WantTwoLevelState)
1583  {
1584  RefinedElement::State = RefinedElement::FINE;
1585  State = Mesh::TWO_LEVEL_FINE;
1586  }
1587 
1588  // 4. Do the green refinement (to get conforming mesh).
1589  int need_refinement;
1590  int edges_in_group, max_edges_in_group = 0;
1591  // edge_splittings identify how the shared edges have been split
1592  int **edge_splittings = new int*[GetNGroups()-1];
1593  for (i = 0; i < GetNGroups()-1; i++)
1594  {
1595  edges_in_group = GroupNEdges(i+1);
1596  edge_splittings[i] = new int[edges_in_group];
1597  if (edges_in_group > max_edges_in_group)
1598  max_edges_in_group = edges_in_group;
1599  }
1600  int neighbor, *iBuf = new int[max_edges_in_group];
1601 
1602  Array<int> group_edges;
1603 
1604  MPI_Request request;
1605  MPI_Status status;
1606  Vertex V;
1607  V(2) = 0.0;
1608 
1609 #ifdef MFEM_DEBUG
1610  int ref_loops_all = 0, ref_loops_par = 0;
1611 #endif
1612  do
1613  {
1614  need_refinement = 0;
1615  for (i = 0; i < nedges; i++)
1616  if (middle[i] != -1 && edge1[i] != -1)
1617  {
1618  need_refinement = 1;
1619  GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
1620  }
1621 #ifdef MFEM_DEBUG
1622  ref_loops_all++;
1623 #endif
1624 
1625  if (uniform_refinement)
1626  continue;
1627 
1628  // if the mesh is locally conforming start making it globally
1629  // conforming
1630  if (need_refinement == 0)
1631  {
1632 #ifdef MFEM_DEBUG
1633  ref_loops_par++;
1634 #endif
1635  // MPI_Barrier(MyComm);
1636 
1637  // (a) send the type of interface splitting
1638  for (i = 0; i < GetNGroups()-1; i++)
1639  {
1640  group_sedge.GetRow(i, group_edges);
1641  edges_in_group = group_edges.Size();
1642  // it is enough to communicate through the edges
1643  if (edges_in_group != 0)
1644  {
1645  for (j = 0; j < edges_in_group; j++)
1646  edge_splittings[i][j] =
1647  GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
1648  middle);
1649  const int *nbs = gtopo.GetGroup(i+1);
1650  if (nbs[0] == 0)
1651  neighbor = gtopo.GetNeighborRank(nbs[1]);
1652  else
1653  neighbor = gtopo.GetNeighborRank(nbs[0]);
1654  MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
1655  neighbor, 0, MyComm, &request);
1656  }
1657  }
1658 
1659  // (b) receive the type of interface splitting
1660  for (i = 0; i < GetNGroups()-1; i++)
1661  {
1662  group_sedge.GetRow(i, group_edges);
1663  edges_in_group = group_edges.Size();
1664  if (edges_in_group != 0)
1665  {
1666  const int *nbs = gtopo.GetGroup(i+1);
1667  if (nbs[0] == 0)
1668  neighbor = gtopo.GetNeighborRank(nbs[1]);
1669  else
1670  neighbor = gtopo.GetNeighborRank(nbs[0]);
1671  MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
1672  MPI_ANY_TAG, MyComm, &status);
1673 
1674  for (j = 0; j < edges_in_group; j++)
1675  if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
1676  {
1677  int *v = shared_edges[group_edges[j]]->GetVertices();
1678  int ii = v_to_v(v[0], v[1]);
1679 #ifdef MFEM_DEBUG
1680  if (middle[ii] != -1)
1681  mfem_error("ParMesh::LocalRefinement (triangles) : "
1682  "Oops!");
1683 #endif
1684  need_refinement = 1;
1685  middle[ii] = NumOfVertices++;
1686  for (int c = 0; c < 2; c++)
1687  V(c) = 0.5 * (vertices[v[0]](c) + vertices[v[1]](c));
1688  vertices.Append(V);
1689  }
1690  }
1691  }
1692 
1693  i = need_refinement;
1694  MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
1695  }
1696  }
1697  while (need_refinement == 1);
1698 
1699 #ifdef MFEM_DEBUG
1700  i = ref_loops_all;
1701  MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
1702  if (MyRank == 0)
1703  {
1704  cout << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
1705  << ref_loops_all << ", ref_loops_par = " << ref_loops_par
1706  << '\n' << endl;
1707  }
1708 #endif
1709 
1710  for (i = 0; i < GetNGroups()-1; i++)
1711  delete [] edge_splittings[i];
1712  delete [] edge_splittings;
1713 
1714  delete [] iBuf;
1715 
1716  // 5. Update the boundary elements.
1717  int v1[2], v2[2], bisect, temp;
1718  temp = NumOfBdrElements;
1719  for (i = 0; i < temp; i++)
1720  {
1721  int *v = boundary[i]->GetVertices();
1722  bisect = v_to_v(v[0], v[1]);
1723  if (middle[bisect] != -1)
1724  { // the element was refined (needs updating)
1725  if (boundary[i]->GetType() == Element::SEGMENT)
1726  {
1727  v1[0] = v[0]; v1[1] = middle[bisect];
1728  v2[0] = middle[bisect]; v2[1] = v[1];
1729 
1730  if (WantTwoLevelState)
1731  {
1732  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
1733 #ifdef MFEM_USE_MEMALLOC
1734  BisectedElement *aux = BEMemory.Alloc();
1735  aux->SetCoarseElem(boundary[i]);
1736 #else
1737  BisectedElement *aux = new BisectedElement(boundary[i]);
1738 #endif
1739  aux->FirstChild =
1740  new Segment(v1, boundary[i]->GetAttribute());
1741  aux->SecondChild = NumOfBdrElements;
1742  boundary[i] = aux;
1743  NumOfBdrElements++;
1744  }
1745  else
1746  {
1747  boundary[i]->SetVertices(v1);
1748  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
1749  }
1750  }
1751  else
1752  mfem_error("Only bisection of segment is implemented for bdr"
1753  " elem.");
1754  }
1755  }
1756  NumOfBdrElements = boundary.Size();
1757 
1758  // 5a. Update the groups after refinement.
1759  RefineGroups(v_to_v, middle);
1760 
1761  // 6. Free the allocated memory.
1762  delete [] edge1;
1763  delete [] edge2;
1764  delete [] middle;
1765 
1766  if (WantTwoLevelState)
1767  {
1768  f_NumOfVertices = NumOfVertices;
1769  f_NumOfElements = NumOfElements;
1770  f_NumOfBdrElements = NumOfBdrElements;
1771  RefinedElement::State = RefinedElement::FINE;
1772  State = Mesh::TWO_LEVEL_FINE;
1773  }
1774 
1775  if (el_to_edge != NULL)
1776  {
1777  if (WantTwoLevelState)
1778  {
1779  c_el_to_edge = el_to_edge;
1780  mfem::Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
1781  f_el_to_edge = new Table;
1782  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
1783  el_to_edge = f_el_to_edge;
1784  f_NumOfEdges = NumOfEdges;
1785  }
1786  else
1787  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1788  GenerateFaces();
1789  }
1790  } // 'if (Dim == 2)'
1791 
1792  if (Dim == 1) // --------------------------------------------------------
1793  {
1794  if (WantTwoLevelState)
1795  {
1796  c_NumOfVertices = NumOfVertices;
1797  c_NumOfElements = NumOfElements;
1798  c_NumOfBdrElements = NumOfBdrElements;
1799  c_NumOfEdges = 0;
1800  }
1801  int cne = NumOfElements, cnv = NumOfVertices;
1802  NumOfVertices += marked_el.Size();
1803  NumOfElements += marked_el.Size();
1804  vertices.SetSize(NumOfVertices);
1805  elements.SetSize(NumOfElements);
1806  for (j = 0; j < marked_el.Size(); j++)
1807  {
1808  i = marked_el[j];
1809  Segment *c_seg = (Segment *)elements[i];
1810  int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
1811  int new_v = cnv + j, new_e = cne + j;
1812  AverageVertices(vert, 2, new_v);
1813  elements[new_e] = new Segment(new_v, vert[1], attr);
1814  if (WantTwoLevelState)
1815  {
1816 #ifdef MFEM_USE_MEMALLOC
1817  BisectedElement *aux = BEMemory.Alloc();
1818  aux->SetCoarseElem(c_seg);
1819 #else
1820  BisectedElement *aux = new BisectedElement(c_seg);
1821 #endif
1822  aux->FirstChild = new Segment(vert[0], new_v, attr);
1823  aux->SecondChild = new_e;
1824  elements[i] = aux;
1825  }
1826  else
1827  {
1828  vert[1] = new_v;
1829  }
1830  }
1831  if (WantTwoLevelState)
1832  {
1833  f_NumOfVertices = NumOfVertices;
1834  f_NumOfElements = NumOfElements;
1835  f_NumOfBdrElements = NumOfBdrElements;
1836  f_NumOfEdges = 0;
1837 
1838  RefinedElement::State = RefinedElement::FINE;
1839  State = Mesh::TWO_LEVEL_FINE;
1840  }
1841  GenerateFaces();
1842  } // end of 'if (Dim == 1)'
1843 
1844  if (Nodes) // curved mesh
1845  {
1846  UpdateNodes();
1847  UseTwoLevelState(wtls);
1848  }
1849 
1850 #ifdef MFEM_DEBUG
1851  CheckElementOrientation(false);
1852  CheckBdrElementOrientation(false);
1853 #endif
1854 }
1855 
1856 void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
1857 {
1858  int i, attr, newv[3], ind, f_ind, *v;
1859 
1860  int group;
1861  Array<int> group_verts, group_edges, group_faces;
1862 
1863  // To update the groups after a refinement, we observe that:
1864  // - every (new and old) vertex, edge and face belongs to exactly one group
1865  // - the refinement does not create new groups
1866  // - a new vertex appears only as the middle of a refined edge
1867  // - a face can be refined 2, 3 or 4 times producing new edges and faces
1868 
1869  int *I_group_svert, *J_group_svert;
1870  int *I_group_sedge, *J_group_sedge;
1871  int *I_group_sface, *J_group_sface;
1872 
1873  I_group_svert = new int[GetNGroups()+1];
1874  I_group_sedge = new int[GetNGroups()+1];
1875  if (Dim == 3)
1876  I_group_sface = new int[GetNGroups()+1];
1877  else
1878  I_group_sface = NULL;
1879 
1880  I_group_svert[0] = I_group_svert[1] = 0;
1881  I_group_sedge[0] = I_group_sedge[1] = 0;
1882  if (Dim == 3)
1883  I_group_sface[0] = I_group_sface[1] = 0;
1884 
1885  // overestimate the size of the J arrays
1886  if (Dim == 3)
1887  {
1888  J_group_svert = new int[group_svert.Size_of_connections()
1889  + group_sedge.Size_of_connections()];
1890  J_group_sedge = new int[2*group_sedge.Size_of_connections()
1891  + 3*group_sface.Size_of_connections()];
1892  J_group_sface = new int[4*group_sface.Size_of_connections()];
1893  }
1894  else if (Dim == 2)
1895  {
1896  J_group_svert = new int[group_svert.Size_of_connections()
1897  + group_sedge.Size_of_connections()];
1898  J_group_sedge = new int[2*group_sedge.Size_of_connections()];
1899  J_group_sface = NULL;
1900  }
1901  else
1902  {
1903  J_group_svert = J_group_sedge = J_group_sface = NULL;
1904  }
1905 
1906  for (group = 0; group < GetNGroups()-1; group++)
1907  {
1908  // Get the group shared objects
1909  group_svert.GetRow(group, group_verts);
1910  group_sedge.GetRow(group, group_edges);
1911  group_sface.GetRow(group, group_faces);
1912 
1913  // Check which edges have been refined
1914  for (i = 0; i < group_sedge.RowSize(group); i++)
1915  {
1916  v = shared_edges[group_edges[i]]->GetVertices();
1917  ind = middle[v_to_v(v[0], v[1])];
1918  if (ind != -1)
1919  {
1920  // add a vertex
1921  group_verts.Append(svert_lvert.Append(ind)-1);
1922  // update the edges
1923  attr = shared_edges[group_edges[i]]->GetAttribute();
1924  shared_edges.Append(new Segment(v[1], ind, attr));
1925  group_edges.Append(sedge_ledge.Append(-1)-1);
1926  v[1] = ind;
1927  }
1928  }
1929 
1930  // Check which faces have been refined
1931  for (i = 0; i < group_sface.RowSize(group); i++)
1932  {
1933  v = shared_faces[group_faces[i]]->GetVertices();
1934  ind = middle[v_to_v(v[0], v[1])];
1935  if (ind != -1)
1936  {
1937  attr = shared_faces[group_faces[i]]->GetAttribute();
1938  // add the refinement edge
1939  shared_edges.Append(new Segment(v[2], ind, attr));
1940  group_edges.Append(sedge_ledge.Append(-1)-1);
1941  // add a face
1942  f_ind = group_faces.Size();
1943  shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
1944  group_faces.Append(sface_lface.Append(-1)-1);
1945  newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
1946  shared_faces[group_faces[i]]->SetVertices(newv);
1947 
1948  // check if the left face has also been refined
1949  // v = shared_faces[group_faces[i]]->GetVertices();
1950  ind = middle[v_to_v(v[0], v[1])];
1951  if (ind != -1)
1952  {
1953  // add the refinement edge
1954  shared_edges.Append(new Segment(v[2], ind, attr));
1955  group_edges.Append(sedge_ledge.Append(-1)-1);
1956  // add a face
1957  shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
1958  group_faces.Append(sface_lface.Append(-1)-1);
1959  newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
1960  shared_faces[group_faces[i]]->SetVertices(newv);
1961  }
1962 
1963  // check if the right face has also been refined
1964  v = shared_faces[group_faces[f_ind]]->GetVertices();
1965  ind = middle[v_to_v(v[0], v[1])];
1966  if (ind != -1)
1967  {
1968  // add the refinement edge
1969  shared_edges.Append(new Segment(v[2], ind, attr));
1970  group_edges.Append(sedge_ledge.Append(-1)-1);
1971  // add a face
1972  shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
1973  group_faces.Append(sface_lface.Append(-1)-1);
1974  newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
1975  shared_faces[group_faces[f_ind]]->SetVertices(newv);
1976  }
1977  }
1978  }
1979 
1980  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
1981  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
1982  if (Dim == 3)
1983  I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
1984 
1985  int *J;
1986  J = J_group_svert+I_group_svert[group];
1987  for (i = 0; i < group_verts.Size(); i++)
1988  J[i] = group_verts[i];
1989  J = J_group_sedge+I_group_sedge[group];
1990  for (i = 0; i < group_edges.Size(); i++)
1991  J[i] = group_edges[i];
1992  if (Dim == 3)
1993  {
1994  J = J_group_sface+I_group_sface[group];
1995  for (i = 0; i < group_faces.Size(); i++)
1996  J[i] = group_faces[i];
1997  }
1998  }
1999 
2000  // Fix the local numbers of shared edges and faces
2001  {
2002  DSTable new_v_to_v(NumOfVertices);
2003  GetVertexToVertexTable(new_v_to_v);
2004  for (i = 0; i < shared_edges.Size(); i++)
2005  {
2006  v = shared_edges[i]->GetVertices();
2007  sedge_ledge[i] = new_v_to_v(v[0], v[1]);
2008  }
2009  }
2010  if (Dim == 3)
2011  {
2012  STable3D *faces_tbl = GetElementToFaceTable(1);
2013  for (i = 0; i < shared_faces.Size(); i++)
2014  {
2015  v = shared_faces[i]->GetVertices();
2016  sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
2017  }
2018  delete faces_tbl;
2019  }
2020 
2021  group_svert.SetIJ(I_group_svert, J_group_svert);
2022  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
2023  if (Dim == 3)
2024  group_sface.SetIJ(I_group_sface, J_group_sface);
2025 }
2026 
2027 void ParMesh::QuadUniformRefinement()
2028 {
2029  SetState(Mesh::NORMAL);
2030  DeleteFaceNbrData();
2031 
2032  int oedge = NumOfVertices, wtls = WantTwoLevelState;
2033 
2034  if (Nodes) // curved mesh
2035  UseTwoLevelState(1);
2036 
2037  // call Mesh::QuadUniformRefinement so that it won't update the nodes
2038  {
2039  GridFunction *nodes = Nodes;
2040  Nodes = NULL;
2041  Mesh::QuadUniformRefinement();
2042  Nodes = nodes;
2043  }
2044 
2045  // update the groups
2046  {
2047  int i, attr, ind, *v;
2048 
2049  int group;
2050  Array<int> sverts, sedges;
2051 
2052  int *I_group_svert, *J_group_svert;
2053  int *I_group_sedge, *J_group_sedge;
2054 
2055  I_group_svert = new int[GetNGroups()+1];
2056  I_group_sedge = new int[GetNGroups()+1];
2057 
2058  I_group_svert[0] = I_group_svert[1] = 0;
2059  I_group_sedge[0] = I_group_sedge[1] = 0;
2060 
2061  // compute the size of the J arrays
2062  J_group_svert = new int[group_svert.Size_of_connections()
2063  + group_sedge.Size_of_connections()];
2064  J_group_sedge = new int[2*group_sedge.Size_of_connections()];
2065 
2066  for (group = 0; group < GetNGroups()-1; group++)
2067  {
2068  // Get the group shared objects
2069  group_svert.GetRow(group, sverts);
2070  group_sedge.GetRow(group, sedges);
2071 
2072  // Process all the edges
2073  for (i = 0; i < group_sedge.RowSize(group); i++)
2074  {
2075  v = shared_edges[sedges[i]]->GetVertices();
2076  ind = oedge + sedge_ledge[sedges[i]];
2077  // add a vertex
2078  sverts.Append(svert_lvert.Append(ind)-1);
2079  // update the edges
2080  attr = shared_edges[sedges[i]]->GetAttribute();
2081  shared_edges.Append(new Segment(v[1], ind, attr));
2082  sedges.Append(sedge_ledge.Append(-1)-1);
2083  v[1] = ind;
2084  }
2085 
2086  I_group_svert[group+1] = I_group_svert[group] + sverts.Size();
2087  I_group_sedge[group+1] = I_group_sedge[group] + sedges.Size();
2088 
2089  int *J;
2090  J = J_group_svert+I_group_svert[group];
2091  for (i = 0; i < sverts.Size(); i++)
2092  J[i] = sverts[i];
2093  J = J_group_sedge+I_group_sedge[group];
2094  for (i = 0; i < sedges.Size(); i++)
2095  J[i] = sedges[i];
2096  }
2097 
2098  // Fix the local numbers of shared edges
2099  DSTable v_to_v(NumOfVertices);
2100  GetVertexToVertexTable(v_to_v);
2101  for (i = 0; i < shared_edges.Size(); i++)
2102  {
2103  v = shared_edges[i]->GetVertices();
2104  sedge_ledge[i] = v_to_v(v[0], v[1]);
2105  }
2106 
2107  group_svert.SetIJ(I_group_svert, J_group_svert);
2108  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
2109  }
2110 
2111  if (Nodes) // curved mesh
2112  {
2113  UpdateNodes();
2114  UseTwoLevelState(wtls);
2115  }
2116 }
2117 
2118 void ParMesh::HexUniformRefinement()
2119 {
2120  SetState(Mesh::NORMAL);
2121  DeleteFaceNbrData();
2122 
2123  int wtls = WantTwoLevelState;
2124  int oedge = NumOfVertices;
2125  int oface = oedge + NumOfEdges;
2126 
2127  DSTable v_to_v(NumOfVertices);
2128  GetVertexToVertexTable(v_to_v);
2129  STable3D *faces_tbl = GetFacesTable();
2130 
2131  if (Nodes) // curved mesh
2132  UseTwoLevelState(1);
2133 
2134  // call Mesh::HexUniformRefinement so that it won't update the nodes
2135  {
2136  GridFunction *nodes = Nodes;
2137  Nodes = NULL;
2138  Mesh::HexUniformRefinement();
2139  Nodes = nodes;
2140  }
2141 
2142  // update the groups
2143  {
2144  int i, attr, newv[4], ind, m[5];
2145  Array<int> v;
2146 
2147  int group;
2148  Array<int> group_verts, group_edges, group_faces;
2149 
2150  int *I_group_svert, *J_group_svert;
2151  int *I_group_sedge, *J_group_sedge;
2152  int *I_group_sface, *J_group_sface;
2153 
2154  I_group_svert = new int[GetNGroups()+1];
2155  I_group_sedge = new int[GetNGroups()+1];
2156  I_group_sface = new int[GetNGroups()+1];
2157 
2158  I_group_svert[0] = I_group_svert[1] = 0;
2159  I_group_sedge[0] = I_group_sedge[1] = 0;
2160  I_group_sface[0] = I_group_sface[1] = 0;
2161 
2162  // compute the size of the J arrays
2163  J_group_svert = new int[group_svert.Size_of_connections()
2164  + group_sedge.Size_of_connections()
2165  + group_sface.Size_of_connections()];
2166  J_group_sedge = new int[2*group_sedge.Size_of_connections()
2167  + 4*group_sface.Size_of_connections()];
2168  J_group_sface = new int[4*group_sface.Size_of_connections()];
2169 
2170  for (group = 0; group < GetNGroups()-1; group++)
2171  {
2172  // Get the group shared objects
2173  group_svert.GetRow(group, group_verts);
2174  group_sedge.GetRow(group, group_edges);
2175  group_sface.GetRow(group, group_faces);
2176 
2177  // Process the edges that have been refined
2178  for (i = 0; i < group_sedge.RowSize(group); i++)
2179  {
2180  shared_edges[group_edges[i]]->GetVertices(v);
2181  ind = oedge + v_to_v(v[0], v[1]);
2182  // add a vertex
2183  group_verts.Append(svert_lvert.Append(ind)-1);
2184  // update the edges
2185  attr = shared_edges[group_edges[i]]->GetAttribute();
2186  shared_edges.Append(new Segment(v[1], ind, attr));
2187  group_edges.Append(sedge_ledge.Append(-1)-1);
2188  newv[0] = v[0]; newv[1] = ind;
2189  shared_edges[group_edges[i]]->SetVertices(newv);
2190  }
2191 
2192  // Process the faces that have been refined
2193  for (i = 0; i < group_sface.RowSize(group); i++)
2194  {
2195  shared_faces[group_faces[i]]->GetVertices(v);
2196  m[0] = oface+(*faces_tbl)(v[0], v[1], v[2], v[3]);
2197  // add a vertex
2198  group_verts.Append(svert_lvert.Append(m[0])-1);
2199  // add the refinement edges
2200  attr = shared_faces[group_faces[i]]->GetAttribute();
2201  m[1] = oedge + v_to_v(v[0], v[1]);
2202  m[2] = oedge + v_to_v(v[1], v[2]);
2203  m[3] = oedge + v_to_v(v[2], v[3]);
2204  m[4] = oedge + v_to_v(v[3], v[0]);
2205  shared_edges.Append(new Segment(m[1], m[0], attr));
2206  group_edges.Append(sedge_ledge.Append(-1)-1);
2207  shared_edges.Append(new Segment(m[2], m[0], attr));
2208  group_edges.Append(sedge_ledge.Append(-1)-1);
2209  shared_edges.Append(new Segment(m[3], m[0], attr));
2210  group_edges.Append(sedge_ledge.Append(-1)-1);
2211  shared_edges.Append(new Segment(m[4], m[0], attr));
2212  group_edges.Append(sedge_ledge.Append(-1)-1);
2213  // update faces
2214  newv[0] = v[0]; newv[1] = m[1]; newv[2] = m[0]; newv[3] = m[4];
2215  shared_faces[group_faces[i]]->SetVertices(newv);
2216  shared_faces.Append(new Quadrilateral(m[1],v[1],m[2],m[0],attr));
2217  group_faces.Append(sface_lface.Append(-1)-1);
2218  shared_faces.Append(new Quadrilateral(m[0],m[2],v[2],m[3],attr));
2219  group_faces.Append(sface_lface.Append(-1)-1);
2220  shared_faces.Append(new Quadrilateral(m[4],m[0],m[3],v[3],attr));
2221  group_faces.Append(sface_lface.Append(-1)-1);
2222  }
2223 
2224  I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
2225  I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
2226  I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
2227 
2228  int *J;
2229  J = J_group_svert+I_group_svert[group];
2230  for (i = 0; i < group_verts.Size(); i++)
2231  J[i] = group_verts[i];
2232  J = J_group_sedge+I_group_sedge[group];
2233  for (i = 0; i < group_edges.Size(); i++)
2234  J[i] = group_edges[i];
2235  J = J_group_sface+I_group_sface[group];
2236  for (i = 0; i < group_faces.Size(); i++)
2237  J[i] = group_faces[i];
2238  }
2239 
2240  // Fix the local numbers of shared edges and faces
2241  DSTable new_v_to_v(NumOfVertices);
2242  GetVertexToVertexTable(new_v_to_v);
2243  for (i = 0; i < shared_edges.Size(); i++)
2244  {
2245  shared_edges[i]->GetVertices(v);
2246  sedge_ledge[i] = new_v_to_v(v[0], v[1]);
2247  }
2248 
2249  delete faces_tbl;
2250  faces_tbl = GetFacesTable();
2251  for (i = 0; i < shared_faces.Size(); i++)
2252  {
2253  shared_faces[i]->GetVertices(v);
2254  sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
2255  }
2256  delete faces_tbl;
2257 
2258  group_svert.SetIJ(I_group_svert, J_group_svert);
2259  group_sedge.SetIJ(I_group_sedge, J_group_sedge);
2260  group_sface.SetIJ(I_group_sface, J_group_sface);
2261  }
2262 
2263  if (Nodes) // curved mesh
2264  {
2265  UpdateNodes();
2266  UseTwoLevelState(wtls);
2267  }
2268 }
2269 
2270 void ParMesh::NURBSUniformRefinement()
2271 {
2272  if (MyRank == 0)
2273  cout << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
2274 }
2275 
2276 void ParMesh::PrintXG(std::ostream &out) const
2277 {
2278  MFEM_ASSERT(Dim == spaceDim, "2D manifolds not supported");
2279  if (Dim == 3 && meshgen == 1)
2280  {
2281  int i, j, nv;
2282  const int *ind;
2283 
2284  out << "NETGEN_Neutral_Format\n";
2285  // print the vertices
2286  out << NumOfVertices << '\n';
2287  for (i = 0; i < NumOfVertices; i++)
2288  {
2289  for (j = 0; j < Dim; j++)
2290  out << " " << vertices[i](j);
2291  out << '\n';
2292  }
2293 
2294  // print the elements
2295  out << NumOfElements << '\n';
2296  for (i = 0; i < NumOfElements; i++)
2297  {
2298  nv = elements[i]->GetNVertices();
2299  ind = elements[i]->GetVertices();
2300  out << elements[i]->GetAttribute();
2301  for (j = 0; j < nv; j++)
2302  out << " " << ind[j]+1;
2303  out << '\n';
2304  }
2305 
2306  // print the boundary + shared faces information
2307  out << NumOfBdrElements + shared_faces.Size() << '\n';
2308  // boundary
2309  for (i = 0; i < NumOfBdrElements; i++)
2310  {
2311  nv = boundary[i]->GetNVertices();
2312  ind = boundary[i]->GetVertices();
2313  out << boundary[i]->GetAttribute();
2314  for (j = 0; j < nv; j++)
2315  out << " " << ind[j]+1;
2316  out << '\n';
2317  }
2318  // shared faces
2319  for (i = 0; i < shared_faces.Size(); i++)
2320  {
2321  nv = shared_faces[i]->GetNVertices();
2322  ind = shared_faces[i]->GetVertices();
2323  out << shared_faces[i]->GetAttribute();
2324  for (j = 0; j < nv; j++)
2325  out << " " << ind[j]+1;
2326  out << '\n';
2327  }
2328  }
2329 
2330  if (Dim == 3 && meshgen == 2)
2331  {
2332  int i, j, nv;
2333  const int *ind;
2334 
2335  out << "TrueGrid\n"
2336  << "1 " << NumOfVertices << " " << NumOfElements << " 0 0 0 0 0 0 0\n"
2337  << "0 0 0 1 0 0 0 0 0 0 0\n"
2338  << "0 0 " << NumOfBdrElements+shared_faces.Size()
2339  << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
2340  << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
2341  << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
2342 
2343  // print the vertices
2344  for (i = 0; i < NumOfVertices; i++)
2345  out << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
2346  << " " << vertices[i](2) << " 0.0\n";
2347 
2348  // print the elements
2349  for (i = 0; i < NumOfElements; i++)
2350  {
2351  nv = elements[i]->GetNVertices();
2352  ind = elements[i]->GetVertices();
2353  out << i+1 << " " << elements[i]->GetAttribute();
2354  for (j = 0; j < nv; j++)
2355  out << " " << ind[j]+1;
2356  out << '\n';
2357  }
2358 
2359  // print the boundary information
2360  for (i = 0; i < NumOfBdrElements; i++)
2361  {
2362  nv = boundary[i]->GetNVertices();
2363  ind = boundary[i]->GetVertices();
2364  out << boundary[i]->GetAttribute();
2365  for (j = 0; j < nv; j++)
2366  out << " " << ind[j]+1;
2367  out << " 1.0 1.0 1.0 1.0\n";
2368  }
2369 
2370  // print the shared faces information
2371  for (i = 0; i < shared_faces.Size(); i++)
2372  {
2373  nv = shared_faces[i]->GetNVertices();
2374  ind = shared_faces[i]->GetVertices();
2375  out << shared_faces[i]->GetAttribute();
2376  for (j = 0; j < nv; j++)
2377  out << " " << ind[j]+1;
2378  out << " 1.0 1.0 1.0 1.0\n";
2379  }
2380  }
2381 
2382  if (Dim == 2)
2383  {
2384  int i, j, attr;
2385  Array<int> v;
2386 
2387  out << "areamesh2\n\n";
2388 
2389  // print the boundary + shared edges information
2390  out << NumOfBdrElements + shared_edges.Size() << '\n';
2391  // boundary
2392  for (i = 0; i < NumOfBdrElements; i++)
2393  {
2394  attr = boundary[i]->GetAttribute();
2395  boundary[i]->GetVertices(v);
2396  out << attr << " ";
2397  for (j = 0; j < v.Size(); j++)
2398  out << v[j] + 1 << " ";
2399  out << '\n';
2400  }
2401  // shared edges
2402  for (i = 0; i < shared_edges.Size(); i++)
2403  {
2404  attr = shared_edges[i]->GetAttribute();
2405  shared_edges[i]->GetVertices(v);
2406  out << attr << " ";
2407  for (j = 0; j < v.Size(); j++)
2408  out << v[j] + 1 << " ";
2409  out << '\n';
2410  }
2411 
2412  // print the elements
2413  out << NumOfElements << '\n';
2414  for (i = 0; i < NumOfElements; i++)
2415  {
2416  attr = elements[i]->GetAttribute();
2417  elements[i]->GetVertices(v);
2418 
2419  out << attr << " ";
2420  if ((j = GetElementType(i)) == Element::TRIANGLE)
2421  out << 3 << " ";
2422  else
2423  if (j == Element::QUADRILATERAL)
2424  out << 4 << " ";
2425  else
2426  if (j == Element::SEGMENT)
2427  out << 2 << " ";
2428  for (j = 0; j < v.Size(); j++)
2429  out << v[j] + 1 << " ";
2430  out << '\n';
2431  }
2432 
2433  // print the vertices
2434  out << NumOfVertices << '\n';
2435  for (i = 0; i < NumOfVertices; i++)
2436  {
2437  for (j = 0; j < Dim; j++)
2438  out << vertices[i](j) << " ";
2439  out << '\n';
2440  }
2441  }
2442 }
2443 
2444 void ParMesh::Print(std::ostream &out) const
2445 {
2446  bool print_shared = true;
2447  int i, j, shared_bdr_attr;
2448  const Array<int> &s2l_face = ((Dim == 1) ? svert_lvert :
2449  ((Dim == 2) ? sedge_ledge : sface_lface));
2450 
2451  if (NURBSext)
2452  {
2453  Mesh::Print(out); // does not print shared boundary
2454  return;
2455  }
2456 
2457  out << "MFEM mesh v1.0\n";
2458 
2459  // optional
2460  out <<
2461  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
2462  "# POINT = 0\n"
2463  "# SEGMENT = 1\n"
2464  "# TRIANGLE = 2\n"
2465  "# SQUARE = 3\n"
2466  "# TETRAHEDRON = 4\n"
2467  "# CUBE = 5\n"
2468  "#\n";
2469 
2470  out << "\ndimension\n" << Dim
2471  << "\n\nelements\n" << NumOfElements << '\n';
2472  for (i = 0; i < NumOfElements; i++)
2473  PrintElement(elements[i], out);
2474 
2475  int num_bdr_elems = NumOfBdrElements;
2476  if (print_shared && Dim > 1)
2477  num_bdr_elems += s2l_face.Size();
2478  out << "\nboundary\n" << num_bdr_elems << '\n';
2479  for (i = 0; i < NumOfBdrElements; i++)
2480  PrintElement(boundary[i], out);
2481 
2482  if (print_shared && Dim > 1)
2483  {
2484  if (bdr_attributes.Size())
2485  shared_bdr_attr = bdr_attributes.Max() + MyRank + 1;
2486  else
2487  shared_bdr_attr = MyRank + 1;
2488  for (i = 0; i < s2l_face.Size(); i++)
2489  {
2490  // Modify the attrributes of the faces (not used otherwise?)
2491  faces[s2l_face[i]]->SetAttribute(shared_bdr_attr);
2492  PrintElement(faces[s2l_face[i]], out);
2493  }
2494  }
2495  out << "\nvertices\n" << NumOfVertices << '\n';
2496  if (Nodes == NULL)
2497  {
2498  out << spaceDim << '\n';
2499  for (i = 0; i < NumOfVertices; i++)
2500  {
2501  out << vertices[i](0);
2502  for (j = 1; j < spaceDim; j++)
2503  out << ' ' << vertices[i](j);
2504  out << '\n';
2505  }
2506  }
2507  else
2508  {
2509  out << "\nnodes\n";
2510  Nodes->Save(out);
2511  }
2512 }
2513 
2514 void ParMesh::PrintAsOne(std::ostream &out)
2515 {
2516  int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
2517  const int *v;
2518  MPI_Status status;
2519  Array<double> vert;
2520  Array<int> ints;
2521 
2522  if (MyRank == 0)
2523  {
2524  out << "MFEM mesh v1.0\n";
2525 
2526  // optional
2527  out <<
2528  "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
2529  "# POINT = 0\n"
2530  "# SEGMENT = 1\n"
2531  "# TRIANGLE = 2\n"
2532  "# SQUARE = 3\n"
2533  "# TETRAHEDRON = 4\n"
2534  "# CUBE = 5\n"
2535  "#\n";
2536 
2537  out << "\ndimension\n" << Dim;
2538  }
2539 
2540  nv = NumOfElements;
2541  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2542  if (MyRank == 0)
2543  {
2544  out << "\n\nelements\n" << ne << '\n';
2545  for (i = 0; i < NumOfElements; i++)
2546  {
2547  // processor number + 1 as attribute and geometry type
2548  out << 1 << ' ' << elements[i]->GetGeometryType();
2549  // vertices
2550  nv = elements[i]->GetNVertices();
2551  v = elements[i]->GetVertices();
2552  for (j = 0; j < nv; j++)
2553  out << ' ' << v[j];
2554  out << '\n';
2555  }
2556  vc = NumOfVertices;
2557  for (p = 1; p < NRanks; p++)
2558  {
2559  MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
2560  ints.SetSize(ne);
2561  MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
2562  for (i = 0; i < ne; )
2563  {
2564  // processor number + 1 as attribute and geometry type
2565  out << p+1 << ' ' << ints[i];
2566  // vertices
2567  k = Geometries.GetVertices(ints[i++])->GetNPoints();
2568  for (j = 0; j < k; j++)
2569  out << ' ' << vc + ints[i++];
2570  out << '\n';
2571  }
2572  vc += nv;
2573  }
2574  }
2575  else
2576  {
2577  // for each element send its geometry type and its vertices
2578  ne = 0;
2579  for (i = 0; i < NumOfElements; i++)
2580  ne += 1 + elements[i]->GetNVertices();
2581  nv = NumOfVertices;
2582  MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
2583  ints.SetSize(ne);
2584  for (i = j = 0; i < NumOfElements; i++)
2585  {
2586  ints[j++] = elements[i]->GetGeometryType();
2587  nv = elements[i]->GetNVertices();
2588  v = elements[i]->GetVertices();
2589  for (k = 0; k < nv; k++)
2590  ints[j++] = v[k];
2591  }
2592  MPI_Send(&ints[0], ne, MPI_INT, 0, 445, MyComm);
2593  }
2594 
2595  // boundary + shared boundary
2596  Array<Element *> &shared_boundary =
2597  (Dim == 2) ? shared_edges : shared_faces;
2598  nv = NumOfBdrElements + shared_boundary.Size();
2599  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2600  if (MyRank == 0)
2601  {
2602  out << "\nboundary\n" << ne << '\n';
2603  // actual boundary
2604  for (i = 0; i < NumOfBdrElements; i++)
2605  {
2606  // processor number + 1 as bdr. attr. and bdr. geometry type
2607  out << 1 << ' ' << boundary[i]->GetGeometryType();
2608  // vertices
2609  nv = boundary[i]->GetNVertices();
2610  v = boundary[i]->GetVertices();
2611  for (j = 0; j < nv; j++)
2612  out << ' ' << v[j];
2613  out << '\n';
2614  }
2615  // shared boundary (interface)
2616  for (i = 0; i < shared_boundary.Size(); i++)
2617  {
2618  // processor number + 1 as bdr. attr. and bdr. geometry type
2619  out << 1 << ' ' << shared_boundary[i]->GetGeometryType();
2620  // vertices
2621  nv = shared_boundary[i]->GetNVertices();
2622  v = shared_boundary[i]->GetVertices();
2623  for (j = 0; j < nv; j++)
2624  out << ' ' << v[j];
2625  out << '\n';
2626  }
2627  vc = NumOfVertices;
2628  for (p = 1; p < NRanks; p++)
2629  {
2630  MPI_Recv(nv_ne, 2, MPI_INT, p, 446, MyComm, &status);
2631  ints.SetSize(ne);
2632  MPI_Recv(&ints[0], ne, MPI_INT, p, 447, MyComm, &status);
2633  for (i = 0; i < ne; )
2634  {
2635  // processor number + 1 as bdr. attr. and bdr. geometry type
2636  out << p+1 << ' ' << ints[i];
2637  k = Geometries.GetVertices(ints[i++])->GetNPoints();
2638  // vertices
2639  for (j = 0; j < k; j++)
2640  out << ' ' << vc + ints[i++];
2641  out << '\n';
2642  }
2643  vc += nv;
2644  }
2645  }
2646  else
2647  {
2648  // for each boundary and shared boundary element send its
2649  // geometry type and its vertices
2650  ne = 0;
2651  for (i = 0; i < NumOfBdrElements; i++)
2652  ne += 1 + boundary[i]->GetNVertices();
2653  for (i = 0; i < shared_boundary.Size(); i++)
2654  ne += 1 + shared_boundary[i]->GetNVertices();
2655  nv = NumOfVertices;
2656  MPI_Send(nv_ne, 2, MPI_INT, 0, 446, MyComm);
2657  ints.SetSize(ne);
2658  // boundary
2659  for (i = j = 0; i < NumOfBdrElements; i++)
2660  {
2661  ints[j++] = boundary[i]->GetGeometryType();
2662  nv = boundary[i]->GetNVertices();
2663  v = boundary[i]->GetVertices();
2664  for (k = 0; k < nv; k++)
2665  ints[j++] = v[k];
2666  }
2667  // shared boundary
2668  for (i = 0; i < shared_boundary.Size(); i++)
2669  {
2670  ints[j++] = shared_boundary[i]->GetGeometryType();
2671  nv = shared_boundary[i]->GetNVertices();
2672  v = shared_boundary[i]->GetVertices();
2673  for (k = 0; k < nv; k++)
2674  ints[j++] = v[k];
2675  }
2676  MPI_Send(&ints[0], ne, MPI_INT, 0, 447, MyComm);
2677  }
2678 
2679  // vertices / nodes
2680  MPI_Reduce(&NumOfVertices, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2681  if (MyRank == 0)
2682  out << "\nvertices\n" << nv << '\n';
2683  if (Nodes == NULL)
2684  {
2685  if (MyRank == 0)
2686  {
2687  out << spaceDim << '\n';
2688  for (i = 0; i < NumOfVertices; i++)
2689  {
2690  out << vertices[i](0);
2691  for (j = 1; j < spaceDim; j++)
2692  out << ' ' << vertices[i](j);
2693  out << '\n';
2694  }
2695  for (p = 1; p < NRanks; p++)
2696  {
2697  MPI_Recv(&nv, 1, MPI_INT, p, 448, MyComm, &status);
2698  vert.SetSize(nv*spaceDim);
2699  MPI_Recv(&vert[0], nv*spaceDim, MPI_DOUBLE, p, 449, MyComm, &status);
2700  for (i = 0; i < nv; i++)
2701  {
2702  out << vert[i*spaceDim];
2703  for (j = 1; j < spaceDim; j++)
2704  out << ' ' << vert[i*spaceDim+j];
2705  out << '\n';
2706  }
2707  }
2708  }
2709  else
2710  {
2711  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 448, MyComm);
2712  vert.SetSize(NumOfVertices*spaceDim);
2713  for (i = 0; i < NumOfVertices; i++)
2714  for (j = 0; j < spaceDim; j++)
2715  vert[i*spaceDim+j] = vertices[i](j);
2716  MPI_Send(&vert[0], NumOfVertices*spaceDim, MPI_DOUBLE, 0, 449, MyComm);
2717  }
2718  }
2719  else
2720  {
2721  if (MyRank == 0)
2722  out << "\nnodes\n";
2723  ParGridFunction *pnodes = dynamic_cast<ParGridFunction *>(Nodes);
2724  if (pnodes)
2725  {
2726  pnodes->SaveAsOne(out);
2727  }
2728  else
2729  {
2730  ParFiniteElementSpace *pfes =
2731  dynamic_cast<ParFiniteElementSpace *>(Nodes->FESpace());
2732  if (pfes)
2733  {
2734  // create a wrapper ParGridFunction
2735  ParGridFunction ParNodes(pfes, Nodes);
2736  ParNodes.SaveAsOne(out);
2737  }
2738  else
2739  mfem_error("ParMesh::PrintAsOne : Nodes have no parallel info!");
2740  }
2741  }
2742 }
2743 
2744 void ParMesh::PrintAsOneXG(std::ostream &out)
2745 {
2746  MFEM_ASSERT(Dim == spaceDim, "2D Manifolds not supported.");
2747  if (Dim == 3 && meshgen == 1)
2748  {
2749  int i, j, k, nv, ne, p;
2750  const int *ind, *v;
2751  MPI_Status status;
2752  Array<double> vert;
2753  Array<int> ints;
2754 
2755  if (MyRank == 0)
2756  {
2757  out << "NETGEN_Neutral_Format\n";
2758  // print the vertices
2759  ne = NumOfVertices;
2760  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2761  out << nv << '\n';
2762  for (i = 0; i < NumOfVertices; i++)
2763  {
2764  for (j = 0; j < Dim; j++)
2765  out << " " << vertices[i](j);
2766  out << '\n';
2767  }
2768  for (p = 1; p < NRanks; p++)
2769  {
2770  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2771  vert.SetSize(Dim*nv);
2772  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
2773  for (i = 0; i < nv; i++)
2774  {
2775  for (j = 0; j < Dim; j++)
2776  out << " " << vert[Dim*i+j];
2777  out << '\n';
2778  }
2779  }
2780 
2781  // print the elements
2782  nv = NumOfElements;
2783  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2784  out << ne << '\n';
2785  for (i = 0; i < NumOfElements; i++)
2786  {
2787  nv = elements[i]->GetNVertices();
2788  ind = elements[i]->GetVertices();
2789  out << 1;
2790  for (j = 0; j < nv; j++)
2791  out << " " << ind[j]+1;
2792  out << '\n';
2793  }
2794  k = NumOfVertices;
2795  for (p = 1; p < NRanks; p++)
2796  {
2797  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2798  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
2799  ints.SetSize(4*ne);
2800  MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
2801  for (i = 0; i < ne; i++)
2802  {
2803  out << p+1;
2804  for (j = 0; j < 4; j++)
2805  out << " " << k+ints[i*4+j]+1;
2806  out << '\n';
2807  }
2808  k += nv;
2809  }
2810  // print the boundary + shared faces information
2811  nv = NumOfBdrElements + shared_faces.Size();
2812  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2813  out << ne << '\n';
2814  // boundary
2815  for (i = 0; i < NumOfBdrElements; i++)
2816  {
2817  nv = boundary[i]->GetNVertices();
2818  ind = boundary[i]->GetVertices();
2819  out << 1;
2820  for (j = 0; j < nv; j++)
2821  out << " " << ind[j]+1;
2822  out << '\n';
2823  }
2824  // shared faces
2825  for (i = 0; i < shared_faces.Size(); i++)
2826  {
2827  nv = shared_faces[i]->GetNVertices();
2828  ind = shared_faces[i]->GetVertices();
2829  out << 1;
2830  for (j = 0; j < nv; j++)
2831  out << " " << ind[j]+1;
2832  out << '\n';
2833  }
2834  k = NumOfVertices;
2835  for (p = 1; p < NRanks; p++)
2836  {
2837  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2838  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
2839  ints.SetSize(3*ne);
2840  MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
2841  for (i = 0; i < ne; i++)
2842  {
2843  out << p+1;
2844  for (j = 0; j < 3; j++)
2845  out << " " << k+ints[i*3+j]+1;
2846  out << '\n';
2847  }
2848  k += nv;
2849  }
2850  }
2851  else
2852  {
2853  ne = NumOfVertices;
2854  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2855  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
2856  vert.SetSize(Dim*NumOfVertices);
2857  for (i = 0; i < NumOfVertices; i++)
2858  for (j = 0; j < Dim; j++)
2859  vert[Dim*i+j] = vertices[i](j);
2860  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
2861  0, 445, MyComm);
2862  // elements
2863  ne = NumOfElements;
2864  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2865  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
2866  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
2867  ints.SetSize(NumOfElements*4);
2868  for (i = 0; i < NumOfElements; i++)
2869  {
2870  v = elements[i]->GetVertices();
2871  for (j = 0; j < 4; j++)
2872  ints[4*i+j] = v[j];
2873  }
2874  MPI_Send(&ints[0], 4*NumOfElements, MPI_INT, 0, 447, MyComm);
2875  // boundary + shared faces
2876  nv = NumOfBdrElements + shared_faces.Size();
2877  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2878  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
2879  ne = NumOfBdrElements + shared_faces.Size();
2880  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
2881  ints.SetSize(3*ne);
2882  for (i = 0; i < NumOfBdrElements; i++)
2883  {
2884  v = boundary[i]->GetVertices();
2885  for (j = 0; j < 3; j++)
2886  ints[3*i+j] = v[j];
2887  }
2888  for ( ; i < ne; i++)
2889  {
2890  v = shared_faces[i-NumOfBdrElements]->GetVertices();
2891  for (j = 0; j < 3; j++)
2892  ints[3*i+j] = v[j];
2893  }
2894  MPI_Send(&ints[0], 3*ne, MPI_INT, 0, 447, MyComm);
2895  }
2896  }
2897 
2898  if (Dim == 3 && meshgen == 2)
2899  {
2900  int i, j, k, nv, ne, p;
2901  const int *ind, *v;
2902  MPI_Status status;
2903  Array<double> vert;
2904  Array<int> ints;
2905 
2906  int TG_nv, TG_ne, TG_nbe;
2907 
2908  if (MyRank == 0)
2909  {
2910  MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2911  MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
2912  nv = NumOfBdrElements + shared_faces.Size();
2913  MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
2914 
2915  out << "TrueGrid\n"
2916  << "1 " << TG_nv << " " << TG_ne << " 0 0 0 0 0 0 0\n"
2917  << "0 0 0 1 0 0 0 0 0 0 0\n"
2918  << "0 0 " << TG_nbe << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
2919  << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
2920  << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
2921 
2922  // print the vertices
2923  nv = TG_nv;
2924  for (i = 0; i < NumOfVertices; i++)
2925  out << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
2926  << " " << vertices[i](2) << " 0.0\n";
2927  for (p = 1; p < NRanks; p++)
2928  {
2929  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2930  vert.SetSize(Dim*nv);
2931  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
2932  for (i = 0; i < nv; i++)
2933  out << i+1 << " 0.0 " << vert[Dim*i] << " " << vert[Dim*i+1]
2934  << " " << vert[Dim*i+2] << " 0.0\n";
2935  }
2936 
2937  // print the elements
2938  ne = TG_ne;
2939  for (i = 0; i < NumOfElements; i++)
2940  {
2941  nv = elements[i]->GetNVertices();
2942  ind = elements[i]->GetVertices();
2943  out << i+1 << " " << 1;
2944  for (j = 0; j < nv; j++)
2945  out << " " << ind[j]+1;
2946  out << '\n';
2947  }
2948  k = NumOfVertices;
2949  for (p = 1; p < NRanks; p++)
2950  {
2951  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2952  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
2953  ints.SetSize(8*ne);
2954  MPI_Recv(&ints[0], 8*ne, MPI_INT, p, 447, MyComm, &status);
2955  for (i = 0; i < ne; i++)
2956  {
2957  out << i+1 << " " << p+1;
2958  for (j = 0; j < 8; j++)
2959  out << " " << k+ints[i*8+j]+1;
2960  out << '\n';
2961  }
2962  k += nv;
2963  }
2964 
2965  // print the boundary + shared faces information
2966  ne = TG_nbe;
2967  // boundary
2968  for (i = 0; i < NumOfBdrElements; i++)
2969  {
2970  nv = boundary[i]->GetNVertices();
2971  ind = boundary[i]->GetVertices();
2972  out << 1;
2973  for (j = 0; j < nv; j++)
2974  out << " " << ind[j]+1;
2975  out << " 1.0 1.0 1.0 1.0\n";
2976  }
2977  // shared faces
2978  for (i = 0; i < shared_faces.Size(); i++)
2979  {
2980  nv = shared_faces[i]->GetNVertices();
2981  ind = shared_faces[i]->GetVertices();
2982  out << 1;
2983  for (j = 0; j < nv; j++)
2984  out << " " << ind[j]+1;
2985  out << " 1.0 1.0 1.0 1.0\n";
2986  }
2987  k = NumOfVertices;
2988  for (p = 1; p < NRanks; p++)
2989  {
2990  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
2991  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
2992  ints.SetSize(4*ne);
2993  MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
2994  for (i = 0; i < ne; i++)
2995  {
2996  out << p+1;
2997  for (j = 0; j < 4; j++)
2998  out << " " << k+ints[i*4+j]+1;
2999  out << " 1.0 1.0 1.0 1.0\n";
3000  }
3001  k += nv;
3002  }
3003  }
3004  else
3005  {
3006  MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
3007  MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
3008  nv = NumOfBdrElements + shared_faces.Size();
3009  MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
3010 
3011  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3012  vert.SetSize(Dim*NumOfVertices);
3013  for (i = 0; i < NumOfVertices; i++)
3014  for (j = 0; j < Dim; j++)
3015  vert[Dim*i+j] = vertices[i](j);
3016  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE, 0, 445, MyComm);
3017  // elements
3018  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3019  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
3020  ints.SetSize(NumOfElements*8);
3021  for (i = 0; i < NumOfElements; i++)
3022  {
3023  v = elements[i]->GetVertices();
3024  for (j = 0; j < 8; j++)
3025  ints[8*i+j] = v[j];
3026  }
3027  MPI_Send(&ints[0], 8*NumOfElements, MPI_INT, 0, 447, MyComm);
3028  // boundary + shared faces
3029  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3030  ne = NumOfBdrElements + shared_faces.Size();
3031  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
3032  ints.SetSize(4*ne);
3033  for (i = 0; i < NumOfBdrElements; i++)
3034  {
3035  v = boundary[i]->GetVertices();
3036  for (j = 0; j < 4; j++)
3037  ints[4*i+j] = v[j];
3038  }
3039  for ( ; i < ne; i++)
3040  {
3041  v = shared_faces[i-NumOfBdrElements]->GetVertices();
3042  for (j = 0; j < 4; j++)
3043  ints[4*i+j] = v[j];
3044  }
3045  MPI_Send(&ints[0], 4*ne, MPI_INT, 0, 447, MyComm);
3046  }
3047  }
3048 
3049  if (Dim == 2)
3050  {
3051  int i, j, k, attr, nv, ne, p;
3052  Array<int> v;
3053  MPI_Status status;
3054  Array<double> vert;
3055  Array<int> ints;
3056 
3057 
3058  if (MyRank == 0)
3059  {
3060  out << "areamesh2\n\n";
3061 
3062  // print the boundary + shared edges information
3063  nv = NumOfBdrElements + shared_edges.Size();
3064  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
3065  out << ne << '\n';
3066  // boundary
3067  for (i = 0; i < NumOfBdrElements; i++)
3068  {
3069  attr = boundary[i]->GetAttribute();
3070  boundary[i]->GetVertices(v);
3071  out << attr << " ";
3072  for (j = 0; j < v.Size(); j++)
3073  out << v[j] + 1 << " ";
3074  out << '\n';
3075  }
3076  // shared edges
3077  for (i = 0; i < shared_edges.Size(); i++)
3078  {
3079  attr = shared_edges[i]->GetAttribute();
3080  shared_edges[i]->GetVertices(v);
3081  out << attr << " ";
3082  for (j = 0; j < v.Size(); j++)
3083  out << v[j] + 1 << " ";
3084  out << '\n';
3085  }
3086  k = NumOfVertices;
3087  for (p = 1; p < NRanks; p++)
3088  {
3089  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
3090  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
3091  ints.SetSize(2*ne);
3092  MPI_Recv(&ints[0], 2*ne, MPI_INT, p, 447, MyComm, &status);
3093  for (i = 0; i < ne; i++)
3094  {
3095  out << p+1;
3096  for (j = 0; j < 2; j++)
3097  out << " " << k+ints[i*2+j]+1;
3098  out << '\n';
3099  }
3100  k += nv;
3101  }
3102 
3103  // print the elements
3104  nv = NumOfElements;
3105  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
3106  out << ne << '\n';
3107  for (i = 0; i < NumOfElements; i++)
3108  {
3109  attr = elements[i]->GetAttribute();
3110  elements[i]->GetVertices(v);
3111  out << 1 << " " << 3 << " ";
3112  for (j = 0; j < v.Size(); j++)
3113  out << v[j] + 1 << " ";
3114  out << '\n';
3115  }
3116  k = NumOfVertices;
3117  for (p = 1; p < NRanks; p++)
3118  {
3119  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
3120  MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
3121  ints.SetSize(3*ne);
3122  MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
3123  for (i = 0; i < ne; i++)
3124  {
3125  out << p+1 << " " << 3;
3126  for (j = 0; j < 3; j++)
3127  out << " " << k+ints[i*3+j]+1;
3128  out << '\n';
3129  }
3130  k += nv;
3131  }
3132 
3133  // print the vertices
3134  ne = NumOfVertices;
3135  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
3136  out << nv << '\n';
3137  for (i = 0; i < NumOfVertices; i++)
3138  {
3139  for (j = 0; j < Dim; j++)
3140  out << vertices[i](j) << " ";
3141  out << '\n';
3142  }
3143  for (p = 1; p < NRanks; p++)
3144  {
3145  MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
3146  vert.SetSize(Dim*nv);
3147  MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
3148  for (i = 0; i < nv; i++)
3149  {
3150  for (j = 0; j < Dim; j++)
3151  out << " " << vert[Dim*i+j];
3152  out << '\n';
3153  }
3154  }
3155  }
3156  else
3157  {
3158  // boundary + shared faces
3159  nv = NumOfBdrElements + shared_edges.Size();
3160  MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
3161  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3162  ne = NumOfBdrElements + shared_edges.Size();
3163  MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
3164  ints.SetSize(2*ne);
3165  for (i = 0; i < NumOfBdrElements; i++)
3166  {
3167  boundary[i]->GetVertices(v);
3168  for (j = 0; j < 2; j++)
3169  ints[2*i+j] = v[j];
3170  }
3171  for ( ; i < ne; i++)
3172  {
3173  shared_edges[i-NumOfBdrElements]->GetVertices(v);
3174  for (j = 0; j < 2; j++)
3175  ints[2*i+j] = v[j];
3176  }
3177  MPI_Send(&ints[0], 2*ne, MPI_INT, 0, 447, MyComm);
3178  // elements
3179  ne = NumOfElements;
3180  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
3181  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3182  MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
3183  ints.SetSize(NumOfElements*3);
3184  for (i = 0; i < NumOfElements; i++)
3185  {
3186  elements[i]->GetVertices(v);
3187  for (j = 0; j < 3; j++)
3188  ints[3*i+j] = v[j];
3189  }
3190  MPI_Send(&ints[0], 3*NumOfElements, MPI_INT, 0, 447, MyComm);
3191  // vertices
3192  ne = NumOfVertices;
3193  MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
3194  MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
3195  vert.SetSize(Dim*NumOfVertices);
3196  for (i = 0; i < NumOfVertices; i++)
3197  for (j = 0; j < Dim; j++)
3198  vert[Dim*i+j] = vertices[i](j);
3199  MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
3200  0, 445, MyComm);
3201  }
3202  }
3203 }
3204 
3205 void ParMesh::PrintInfo(std::ostream &out)
3206 {
3207  int i;
3208  DenseMatrix J(Dim);
3209  double h_min, h_max, kappa_min, kappa_max, h, kappa;
3210 
3211  if (MyRank == 0)
3212  out << "Parallel Mesh Stats:" << endl;
3213 
3214  for (i = 0; i < NumOfElements; i++)
3215  {
3216  GetElementJacobian(i, J);
3217  h = pow(fabs(J.Det()), 1.0/double(Dim));
3218  kappa = J.CalcSingularvalue(0) / J.CalcSingularvalue(Dim-1);
3219  if (i == 0)
3220  {
3221  h_min = h_max = h;
3222  kappa_min = kappa_max = kappa;
3223  }
3224  else
3225  {
3226  if (h < h_min) h_min = h;
3227  if (h > h_max) h_max = h;
3228  if (kappa < kappa_min) kappa_min = kappa;
3229  if (kappa > kappa_max) kappa_max = kappa;
3230  }
3231  }
3232 
3233  double gh_min, gh_max, gk_min, gk_max;
3234  MPI_Reduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
3235  MPI_Reduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
3236  MPI_Reduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
3237  MPI_Reduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
3238 
3239  int ldata[5]; // vert, edge, face, elem, neighbors;
3240  int mindata[5], maxdata[5], sumdata[5];
3241 
3242  // count locally owned vertices, edges, and faces
3243  ldata[0] = GetNV();
3244  ldata[1] = GetNEdges();
3245  ldata[2] = GetNFaces();
3246  ldata[3] = GetNE();
3247  ldata[4] = gtopo.GetNumNeighbors()-1;
3248  for (int gr = 1; gr < GetNGroups(); gr++)
3249  if (!gtopo.IAmMaster(gr)) // we are not the master
3250  {
3251  ldata[0] -= group_svert.RowSize(gr-1);
3252  ldata[1] -= group_sedge.RowSize(gr-1);
3253  ldata[2] -= group_sface.RowSize(gr-1);
3254  }
3255 
3256  MPI_Reduce(ldata, mindata, 5, MPI_INT, MPI_MIN, 0, MyComm);
3257  MPI_Reduce(ldata, sumdata, 5, MPI_INT, MPI_SUM, 0, MyComm); // overflow?
3258  MPI_Reduce(ldata, maxdata, 5, MPI_INT, MPI_MAX, 0, MyComm);
3259 
3260  if (MyRank == 0)
3261  {
3262  out << '\n'
3263  << " "
3264  << setw(12) << "minimum"
3265  << setw(12) << "average"
3266  << setw(12) << "maximum"
3267  << setw(12) << "total" << '\n';
3268  out << " vertices "
3269  << setw(12) << mindata[0]
3270  << setw(12) << sumdata[0]/NRanks
3271  << setw(12) << maxdata[0]
3272  << setw(12) << sumdata[0] << '\n';
3273  out << " edges "
3274  << setw(12) << mindata[1]
3275  << setw(12) << sumdata[1]/NRanks
3276  << setw(12) << maxdata[1]
3277  << setw(12) << sumdata[1] << '\n';
3278  if (Dim == 3)
3279  out << " faces "
3280  << setw(12) << mindata[2]
3281  << setw(12) << sumdata[2]/NRanks
3282  << setw(12) << maxdata[2]
3283  << setw(12) << sumdata[2] << '\n';
3284  out << " elements "
3285  << setw(12) << mindata[3]
3286  << setw(12) << sumdata[3]/NRanks
3287  << setw(12) << maxdata[3]
3288  << setw(12) << sumdata[3] << '\n';
3289  out << " neighbors "
3290  << setw(12) << mindata[4]
3291  << setw(12) << sumdata[4]/NRanks
3292  << setw(12) << maxdata[4] << '\n';
3293  out << '\n'
3294  << " "
3295  << setw(12) << "minimum"
3296  << setw(12) << "maximum" << '\n';
3297  out << " h "
3298  << setw(12) << gh_min
3299  << setw(12) << gh_max << '\n';
3300  out << " kappa "
3301  << setw(12) << gk_min
3302  << setw(12) << gk_max << '\n';
3303  }
3304 }
3305 
3306 ParMesh::~ParMesh()
3307 {
3308  int i;
3309 
3310  DeleteFaceNbrData();
3311 
3312  for (i = 0; i < shared_faces.Size(); i++)
3313  FreeElement(shared_faces[i]);
3314  for (i = 0; i < shared_edges.Size(); i++)
3315  FreeElement(shared_edges[i]);
3316 
3317  // The Mesh destructor is called automatically
3318 }
3319 
3320 }
3321 
3322 #endif
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:202
mfem::ParMesh::GetNFaceNeighbors
int GetNFaceNeighbors() const
Definition: pmesh.hpp:99
mfem::GroupTopology::Create
void Create(ListOfIntegerSets &groups, int mpitag)
Definition: communication.cpp:54
mfem::Mesh::SetState
void SetState(int s)
Change the mesh state to NORMAL, TWO_LEVEL_COARSE, TWO_LEVEL_FINE.
Definition: mesh.cpp:6298
mfem::Array::Size
int Size() const
Logical size of the array.
Definition: array.hpp:108
mfem::Mesh::f_el_to_face
Table * f_el_to_face
Definition: mesh.hpp:85
mfem::IntegerSet::Recreate
void Recreate(const int n, const int *p)
Definition: sets.cpp:51
mfem::Table::GetJ
int * GetJ()
Definition: table.hpp:88
mfem::Mesh::GetVertex
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:421
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
mfem::ParMesh::~ParMesh
virtual ~ParMesh()
Definition: pmesh.cpp:3306
mfem::ParGridFunction::SaveAsOne
void SaveAsOne(std::ostream &out=std::cout)
Merge the local grid functions.
Definition: pgridfunc.cpp:263
mfem::Mesh::FreeElement
void FreeElement(Element *E)
Definition: mesh.cpp:8206
mfem::Element::Duplicate
virtual Element * Duplicate(Mesh *m) const =0
mfem::Mesh::GetEdgeVertices
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:3421
mfem::Mesh::f_NumOfEdges
int f_NumOfEdges
Definition: mesh.hpp:62
mfem::Table::AddColumnsInRow
void AddColumnsInRow(int r, int ncol)
Definition: table.hpp:52
mfem::Mesh::f_NumOfElements
int f_NumOfElements
Definition: mesh.hpp:60
mfem::Table::MakeI
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:57
mfem::Element::GetVertices
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
mfem::Mesh::c_bel_to_edge
Table * c_bel_to_edge
Definition: mesh.hpp:83
mfem::Mesh::GetFaceElements
void GetFaceElements(int Face, int *Elem1, int *Elem2)
Definition: mesh.cpp:615
mfem::Mesh::boundary
Array< Element * > boundary
Definition: mesh.hpp:68
mfem::Mesh::GeneratePartitioning
int * GeneratePartitioning(int nparts, int part_method=1)
Definition: mesh.cpp:4252
mfem::ParMesh::ReorientTetMesh
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:1200
mfem::Mesh::own_nodes
int own_nodes
Definition: mesh.hpp:96
mfem::Mesh::GetNBE
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:399
mfem::Mesh::GetLocalPtToSegTransformation
void GetLocalPtToSegTransformation(IsoparametricTransformation &, int)
Used in GetFaceElementTransformations (...)
Definition: mesh.cpp:374
mfem::ParMesh::face_nbr_elements
Array< Element * > face_nbr_elements
Definition: pmesh.hpp:79
mfem::DenseMatrix::Det
double Det() const
Calculates the determinant of the matrix (for 2x2 or 3x3 matrices)
Definition: densemat.cpp:321
mfem::Mesh::NumOfEdges
int NumOfEdges
Definition: mesh.hpp:52
mfem::GroupTopology::GetGroupSize
int GetGroupSize(int g) const
Definition: communication.hpp:69
mfem::ParMesh::GetNSharedFaces
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:1191
mfem::Table::SetDims
void SetDims(int rows, int nnz)
Definition: table.cpp:107
mfem::Mesh::GetElementBaseGeometry
int GetElementBaseGeometry(int i) const
Definition: mesh.hpp:436
mfem::GroupTopology::GetGroup
const int * GetGroup(int g) const
Definition: communication.hpp:72
mfem::Array::Copy
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:155
mfem::Mesh::WantTwoLevelState
int WantTwoLevelState
Definition: mesh.hpp:54
mfem::ParMesh::LocalRefinement
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
Refine the marked elements.
Definition: pmesh.cpp:1242
mfem::DenseMatrix
Data type dense matrix.
Definition: densemat.hpp:22
mfem::Mesh::Nodes
GridFunction * Nodes
Definition: mesh.hpp:95
mfem::Mesh::NumOfElements
int NumOfElements
Definition: mesh.hpp:51
mfem::Table::GetRow
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:148
mfem::FaceElementTransformations::Loc2
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:102
mfem::Mesh::GetBdrElementEdgeIndex
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:3625
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
mfem::Mesh::GetFace
const Element * GetFace(int i) const
Definition: mesh.hpp:432
mfem::Element::SetVertices
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
Definition: element.cpp:17
mfem::Mesh::GetElementVertices
void GetElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of element i.
Definition: mesh.hpp:443
mfem::ParFiniteElementSpace
Abstract parallel finite element space.
Definition: pfespace.hpp:28
mfem::ParMesh::have_face_nbr_data
bool have_face_nbr_data
Definition: pmesh.hpp:75
mfem::ParMesh::face_nbr_vertices_offset
Array< int > face_nbr_vertices_offset
Definition: pmesh.hpp:78
mfem::Mesh::DeleteCoarseTables
void DeleteCoarseTables()
Definition: mesh.cpp:671
mfem::Vertex
Data type for vertex.
Definition: vertex.hpp:21
mfem::ParMesh::face_nbr_group
Array< int > face_nbr_group
Definition: pmesh.hpp:76
mfem::Mesh::GetQuadOrientation
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:3230
mfem::Mesh::GetVertices
void GetVertices(Vector &vert_coord) const
Definition: mesh.cpp:4893
mfem::Mesh::c_el_to_face
Table * c_el_to_face
Definition: mesh.hpp:85
mfem::FaceElementTransformations::Elem2
ElementTransformation * Elem2
Definition: eltrans.hpp:101
mfem::Mesh::CheckBdrElementOrientation
void CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
Definition: mesh.cpp:3266
mfem::Mesh::GetFaceElementTransformations
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:499
mfem::Table::Size_of_connections
int Size_of_connections() const
Definition: table.hpp:72
mfem::Mesh::faces
Array< Element * > faces
Definition: mesh.hpp:69
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Definition: geom.cpp:172
mfem::Mesh::GetVertexToVertexTable
void GetVertexToVertexTable(DSTable &) const
Definition: mesh.cpp:3723
mfem::Mesh::RedRefinement
void RedRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
Definition: mesh.hpp:153
mfem::Element::GetGeometryType
int GetGeometryType() const
Definition: element.hpp:49
mfem::Array::DeleteAll
void DeleteAll()
Delete whole array.
Definition: array.hpp:407
mfem::Table::AddConnections
void AddConnections(int r, const int *c, int nc)
Definition: table.cpp:75
mfem::Mesh::fc_faces_info
Array< FaceInfo > fc_faces_info
Definition: mesh.hpp:86
mfem::GroupTopology::IAmMaster
bool IAmMaster(int g) const
Definition: communication.hpp:59
mfem::Mesh::NewElement
Element * NewElement(int geom)
Definition: mesh.cpp:1670
mfem::Mesh::el_to_face
Table * el_to_face
Definition: mesh.hpp:75
mfem::Mesh::FaceElemTr
FaceElementTransformations FaceElemTr
Definition: mesh.hpp:90
mfem::SortPairs< int, int >
template void SortPairs< int, int >(Pair< int, int > *, int)
mfem::ParMesh::ExchangeFaceNbrData
void ExchangeFaceNbrData()
Definition: pmesh.cpp:662
mfem::GridFunction::Save
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:1797
mfem::ParMesh::PrintAsOne
void PrintAsOne(std::ostream &out=std::cout)
Definition: pmesh.cpp:2514
mfem::Segment::GetVertices
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
Definition: segment.cpp:35
mfem::Mesh::GetElementJacobian
void GetElementJacobian(int i, DenseMatrix &J)
Definition: mesh.cpp:31
mfem::Geometries
Geometry Geometries
Definition: geom.cpp:267
mfem::ParMesh::face_nbr_elements_offset
Array< int > face_nbr_elements_offset
Definition: pmesh.hpp:77
mfem::STable3D
Symmetric 3D Table.
Definition: stable3d.hpp:28
mfem::Mesh::quad_orientations
static const int quad_orientations[8][4]
Definition: mesh.hpp:106
mfem::Array::Append
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:330
mfem::Mesh::State
int State
Definition: mesh.hpp:54
mfem::Mesh::f_el_to_edge
Table * f_el_to_edge
Definition: mesh.hpp:83
mfem::GroupTopology::GetNumNeighbors
int GetNumNeighbors() const
Definition: communication.hpp:56
mfem::Table::Clear
void Clear()
Definition: table.cpp:265
mfem::Mesh::GetLocalSegToTriTransformation
void GetLocalSegToTriTransformation(IsoparametricTransformation &loc, int i)
Definition: mesh.cpp:388
mfem::Mesh::GetTransformationFEforElementType
static FiniteElement * GetTransformationFEforElementType(int)
Definition: mesh.cpp:141
mfem::Table::AddConnection
void AddConnection(int r, int c)
Definition: table.hpp:54
mfem::Mesh::GetElementToFaceTable
STable3D * GetElementToFaceTable(int ret_ftbl=0)
Definition: mesh.cpp:4053
mfem::Mesh::Rotate3
static void Rotate3(int &, int &, int &)
Definition: mesh.hpp:801
mfem::Mesh::MeshGenerator
int MeshGenerator()
Truegrid or NetGen?
Definition: mesh.hpp:389
mfem::ParMesh::send_face_nbr_vertices
Table send_face_nbr_vertices
Definition: pmesh.hpp:83
mfem::Array::Max
T Max() const
Definition: array.cpp:78
mfem::Mesh::CheckElementOrientation
void CheckElementOrientation(bool fix_it=true)
Check the orientation of the elements.
Definition: mesh.cpp:3098
mesh_headers.hpp
mfem::ListOfIntegerSets::Insert
int Insert(IntegerSet &s)
Definition: sets.cpp:70
mfem::Mesh::f_bel_to_edge
Table * f_bel_to_edge
Definition: mesh.hpp:83
mfem::Transpose
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:305
mfem::Element::GetAttribute
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:79
mfem::Mesh::GetElementToEdgeTable
int GetElementToEdgeTable(Table &, Array< int > &)
Definition: mesh.cpp:3748
mfem::Mesh::c_NumOfEdges
int c_NumOfEdges
Definition: mesh.hpp:61
mfem::Mesh::GetElementType
int GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:3632
mfem::Triangle
Data type triangle element.
Definition: triangle.hpp:22
mfem::RefinedElement::FirstChild
Element * FirstChild
Definition: element.hpp:101
mfem::Mesh::GetElement
const Element * GetElement(int i) const
Definition: mesh.hpp:424
mfem::Mesh::GetLocalQuadToHexTransformation
void GetLocalQuadToHexTransformation(IsoparametricTransformation &loc, int i)
Used in GetFaceElementTransformations (...)
Definition: mesh.cpp:478
mfem::RefinedElement::State
static int State
Definition: element.hpp:99
mfem::Mesh::Transformation2
IsoparametricTransformation Transformation2
Definition: mesh.hpp:88
mfem::Mesh::f_NumOfVertices
int f_NumOfVertices
Definition: mesh.hpp:60
mfem::Table::Size
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:66
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:246
mfem::Mesh::ReorientTetMesh
virtual void ReorientTetMesh()
Definition: mesh.cpp:4118
mfem::Mesh::NumOfBdrElements
int NumOfBdrElements
Definition: mesh.hpp:51
mfem::Mesh::el_to_edge
Table * el_to_edge
Definition: mesh.hpp:74
mfem::Mesh::c_NumOfBdrElements
int c_NumOfBdrElements
Definition: mesh.hpp:59
mfem::Tetrahedron
Data type tetrahedron element.
Definition: tetrahedron.hpp:22
mfem::Mesh::c_NumOfVertices
int c_NumOfVertices
Definition: mesh.hpp:59
mfem::Mesh::GetBdrElementFace
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:3560
kappa
const double kappa
Definition: ex3.cpp:186
mfem::Mesh::Print
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
Definition: mesh.cpp:7136
mfem::Swap
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:284
mfem::Mesh::bdr_attributes
Array< int > bdr_attributes
Definition: mesh.hpp:305
mfem::Mesh::UseTwoLevelState
void UseTwoLevelState(int use)
Definition: mesh.hpp:684
mfem::RefinedElement::SetCoarseElem
void SetCoarseElem(Element *ce)
Definition: element.hpp:110
mfem::ParMesh::RefineGroups
void RefineGroups(const DSTable &v_to_v, int *middle)
Update the groups after tet refinement.
Definition: pmesh.cpp:1856
mfem::ParMesh::face_nbr_vertices
Array< Vertex > face_nbr_vertices
Definition: pmesh.hpp:80
mfem::Mesh::PrintElement
static void PrintElement(const Element *, std::ostream &)
Definition: mesh.cpp:1727
mfem::GroupTopology::GetNeighborRank
int GetNeighborRank(int i) const
Definition: communication.hpp:57
mfem::Table::AddAColumnInRow
void AddAColumnInRow(int r)
Definition: table.hpp:51
mfem::mfem_error
void mfem_error(const char *msg)
Definition: error.cpp:23
mfem::Array::SetSize
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:293
mfem::Vector::SetSubVector
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Definition: vector.cpp:427
mfem::Mesh::vertices
Array< Vertex > vertices
Definition: mesh.hpp:67
mfem::FiniteElementSpace::GetElementVDofs
void GetElementVDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:117
mfem::Mesh::BEMemory
MemAlloc< BisectedElement, 1024 > BEMemory
Definition: mesh.hpp:112
mfem::ParMesh::PrintInfo
void PrintInfo(std::ostream &out=std::cout)
Print various parallel mesh stats.
Definition: pmesh.cpp:3205
mfem::Mesh::QuadUniformRefinement
virtual void QuadUniformRefinement()
Refine quadrilateral mesh.
Definition: mesh.cpp:5009
mfem::Mesh::HexUniformRefinement
virtual void HexUniformRefinement()
Refine hexahedral mesh.
Definition: mesh.cpp:5156
mfem::Mesh::elements
Array< Element * > elements
Definition: mesh.hpp:64
mfem::Table::ShiftUpI
void ShiftUpI()
Definition: table.cpp:84
mfem::Mesh::meshgen
int meshgen
Definition: mesh.hpp:57
mfem::Mesh::bel_to_edge
Table * bel_to_edge
Definition: mesh.hpp:78
mfem::ParMesh::PrintXG
virtual void PrintXG(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2276
mfem::Mesh::fc_be_to_edge
Array< int > fc_be_to_edge
Definition: mesh.hpp:84
mfem::Tetrahedron::GetVertices
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
Definition: tetrahedron.cpp:243
mfem::Mesh::GetLocalTriToTetTransformation
void GetLocalTriToTetTransformation(IsoparametricTransformation &loc, int i)
Used in GetFaceElementTransformations (...)
Definition: mesh.cpp:455
mfem::Mesh::NURBSext
NURBSExtension * NURBSext
Definition: mesh.hpp:307
mfem::Mesh::GetNV
int GetNV() const
Definition: mesh.hpp:393
mfem::Mesh::be_to_edge
Array< int > be_to_edge
Definition: mesh.hpp:77
mfem::Mesh::c_el_to_edge
Table * c_el_to_edge
Definition: mesh.hpp:83
mfem::ParMesh::GroupFace
void GroupFace(int group, int i, int &face, int &o)
Definition: pmesh.cpp:535
mfem::Mesh::f_NumOfFaces
int f_NumOfFaces
Definition: mesh.hpp:62
mfem::Mesh::GetLocalSegToQuadTransformation
void GetLocalSegToQuadTransformation(IsoparametricTransformation &loc, int i)
Definition: mesh.cpp:411
mfem::Mesh::c_NumOfElements
int c_NumOfElements
Definition: mesh.hpp:59
mfem::ParMesh::ExchangeFaceNbrNodes
void ExchangeFaceNbrNodes()
Definition: pmesh.cpp:1029
mfem::Mesh::c_NumOfFaces
int c_NumOfFaces
Definition: mesh.hpp:61
mfem::IntegerSet
A set of integers.
Definition: sets.hpp:23
mfem::Table::SetIJ
void SetIJ(int *newI, int *newJ, int newsize=-1)
Definition: table.cpp:157
mfem::Table::MakeJ
void MakeJ()
Definition: table.cpp:65
mfem::ParMesh::GroupNFaces
int GroupNFaces(int group)
Definition: pmesh.hpp:90
mfem::IntegrationPointTransformation::Transf
IsoparametricTransformation Transf
Definition: eltrans.hpp:92
mfem::DenseMatrix::CalcSingularvalue
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1451
mfem::Mesh::faces_info
Array< FaceInfo > faces_info
Definition: mesh.hpp:72
mfem::Element::GetNVertices
virtual int GetNVertices() const =0
mfem::Element::SetAttribute
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:82
mfem::Mesh::GetFacesTable
STable3D * GetFacesTable()
Definition: mesh.cpp:4020
mfem::Mesh::Dim
int Dim
Definition: mesh.hpp:48
mfem::Mesh::GetNEdges
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:402
mfem::ParMesh::GroupEdge
void GroupEdge(int group, int i, int &edge, int &o)
Definition: pmesh.cpp:527
mfem::Mesh::GetNFaces
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:405
mfem::Vector
Vector data type.
Definition: vector.hpp:29
mfem::Mesh::GetNodes
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:4948
mfem::Mesh::AverageVertices
void AverageVertices(int *indexes, int n, int result)
Definition: mesh.cpp:4988
mfem::DSTable::NumberOfEntries
int NumberOfEntries() const
Definition: table.hpp:189
mfem::Mesh::f_NumOfBdrElements
int f_NumOfBdrElements
Definition: mesh.hpp:60
mfem::Mesh::GetTriOrientation
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of &quot;test&quot; relative to &quot;base&quot;.
Definition: mesh.cpp:3200
mfem::Table::GetI
int * GetI()
Definition: table.hpp:87
mfem::ParMesh::PrintAsOneXG
void PrintAsOneXG(std::ostream &out=std::cout)
Old mesh format (Netgen/Truegrid) version of &#39;PrintAsOne&#39;.
Definition: pmesh.cpp:2744
mfem::Mesh::GenerateFaces
void GenerateFaces()
Definition: mesh.cpp:3954
mfem::Mesh::spaceDim
int spaceDim
Definition: mesh.hpp:49
mfem::Tetrahedron::ParseRefinementFlag
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Definition: tetrahedron.cpp:22
mfem::ParGridFunction
Class for parallel grid function.
Definition: pgridfunc.hpp:31
mfem::Table::RowSize
int RowSize(int i) const
Definition: table.hpp:82
mfem::ParMesh::send_face_nbr_elements
Table send_face_nbr_elements
Definition: pmesh.hpp:82
mfem::ListOfIntegerSets
List of integer sets.
Definition: sets.hpp:49
mfem::Mesh::NumOfVertices
int NumOfVertices
Definition: mesh.hpp:51
mfem::ParMesh::GetFaceToAllElementTable
Table * GetFaceToAllElementTable() const
Definition: pmesh.cpp:1092
mfem::ParMesh::gtopo
GroupTopology gtopo
Definition: pmesh.hpp:72
mfem::ParMesh::GetSharedFaceTransformations
FaceElementTransformations * GetSharedFaceTransformations(int)
Definition: pmesh.cpp:1136
mfem::Mesh::Bisection
void Bisection(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:5942
mfem::ParMesh::GetNGroups
int GetNGroups()
Definition: pmesh.hpp:85
mfem::Element
Abstract data type element.
Definition: element.hpp:27
mfem::Mesh::GetAttribute
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:583
mfem::ParMesh::GetFaceNbrRank
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:1082
mfem::Mesh::ElementToEdgeTable
const Table & ElementToEdgeTable() const
Definition: mesh.cpp:3855
mfem::Segment
Data type line segment element.
Definition: segment.hpp:22
mfem::Mesh::attributes
Array< int > attributes
Definition: mesh.hpp:304
mfem::Mesh::GetVertexToElementTable
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
Definition: mesh.cpp:3475
mfem::Mesh::GetBdrElement
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:428
mfem::Element::GetType
virtual int GetType() const =0
Returns element&#39;s type.
mfem::Mesh::UpdateNodes
void UpdateNodes()
Update the nodes of a curved mesh after refinement.
Definition: mesh.cpp:5003
mfem::ParMesh::GroupNEdges
int GroupNEdges(int group)
Definition: pmesh.hpp:89
mfem::Mesh::GreenRefinement
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
Definition: mesh.hpp:159
mfem::Mesh::tri_orientations
static const int tri_orientations[6][3]
Definition: mesh.hpp:105
mfem::ParMesh::Print
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2444
mfem::Mesh::NumOfFaces
int NumOfFaces
Definition: mesh.hpp:52
Generated on Fri Feb 18 2022 18:57:40 for MFEM by   doxygen 1.8.5