MFEM v2.0
pmesh.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 #ifdef MFEM_USE_MPI
00013 
00014 #include "mesh_headers.hpp"
00015 #include "../fem/fem.hpp"
00016 #include "../general/sets.hpp"
00017 
00018 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
00019                  int part_method)
00020    : gtopo(comm)
00021 {
00022    int i, j;
00023    int *partitioning;
00024    Array<bool> activeBdrElem;
00025 
00026    MyComm = comm;
00027    MPI_Comm_size(MyComm, &NRanks);
00028    MPI_Comm_rank(MyComm, &MyRank);
00029 
00030    Dim = mesh.Dim;
00031 
00032    if (partitioning_)
00033       partitioning = partitioning_;
00034    else
00035       partitioning = mesh.GeneratePartitioning(NRanks, part_method);
00036 
00037    // re-enumerate the partitions to better map to actual processor
00038    // interconnect topology !?
00039 
00040    Array<int> vert;
00041    Array<int> vert_global_local(mesh.GetNV());
00042    int vert_counter, element_counter, bdrelem_counter;
00043 
00044    // build vert_global_local
00045    vert_global_local = -1;
00046 
00047    element_counter = 0;
00048    vert_counter = 0;
00049    for (i = 0; i < mesh.GetNE(); i++)
00050       if (partitioning[i] == MyRank)
00051       {
00052          mesh.GetElementVertices(i, vert);
00053          element_counter++;
00054          for (j = 0; j < vert.Size(); j++)
00055             if (vert_global_local[vert[j]] < 0)
00056                vert_global_local[vert[j]] = vert_counter++;
00057       }
00058 
00059    NumOfVertices = vert_counter;
00060    NumOfElements = element_counter;
00061    vertices.SetSize(NumOfVertices);
00062 
00063    // re-enumerate the local vertices to preserve the global ordering
00064    for (i = vert_counter = 0; i < vert_global_local.Size(); i++)
00065       if (vert_global_local[i] >= 0)
00066          vert_global_local[i] = vert_counter++;
00067 
00068    // determine vertices
00069    for (i = 0; i < vert_global_local.Size(); i++)
00070       if (vert_global_local[i] >= 0)
00071          vertices[vert_global_local[i]].SetCoords(mesh.GetVertex(i));
00072 
00073    // determine elements
00074    element_counter = 0;
00075    elements.SetSize(NumOfElements);
00076    for (i = 0; i < mesh.GetNE(); i++)
00077       if (partitioning[i] == MyRank)
00078       {
00079          elements[element_counter] = mesh.GetElement(i)->Duplicate(this);
00080          int *v = elements[element_counter]->GetVertices();
00081          int nv = elements[element_counter]->GetNVertices();
00082          for (j = 0; j < nv; j++)
00083             v[j] = vert_global_local[v[j]];
00084          element_counter++;
00085       }
00086 
00087    Table *edge_element = NULL;
00088    if (mesh.NURBSext)
00089    {
00090       activeBdrElem.SetSize(mesh.GetNBE());
00091       activeBdrElem = false;
00092    }
00093    // build boundary elements
00094    if (Dim == 3)
00095    {
00096       NumOfBdrElements = 0;
00097       for (i = 0; i < mesh.GetNBE(); i++)
00098       {
00099          int face = mesh.GetBdrElementEdgeIndex(i);
00100          int el1, el2;
00101          mesh.GetFaceElements(face, &el1, &el2);
00102          if (partitioning[el1] == MyRank)
00103          {
00104             NumOfBdrElements++;
00105             if (mesh.NURBSext)
00106                activeBdrElem[i] = true;
00107          }
00108       }
00109 
00110       bdrelem_counter = 0;
00111       boundary.SetSize(NumOfBdrElements);
00112       for (i = 0; i < mesh.GetNBE(); i++)
00113       {
00114          int face = mesh.GetBdrElementEdgeIndex(i);
00115          int el1, el2;
00116          mesh.GetFaceElements(face, &el1, &el2);
00117          if (partitioning[el1] == MyRank)
00118          {
00119             boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
00120             int *v = boundary[bdrelem_counter]->GetVertices();
00121             int nv = boundary[bdrelem_counter]->GetNVertices();
00122             for (j = 0; j < nv; j++)
00123                v[j] = vert_global_local[v[j]];
00124             bdrelem_counter++;
00125          }
00126       }
00127 
00128    }
00129    else if (Dim == 2)
00130    {
00131       edge_element = new Table;
00132       Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
00133 
00134       NumOfBdrElements = 0;
00135       for (i = 0; i < mesh.GetNBE(); i++)
00136       {
00137          int edge = mesh.GetBdrElementEdgeIndex(i);
00138          int el1 = edge_element->GetRow(edge)[0];
00139          if (partitioning[el1] == MyRank)
00140          {
00141             NumOfBdrElements++;
00142             if (mesh.NURBSext)
00143                activeBdrElem[i] = true;
00144          }
00145       }
00146 
00147       bdrelem_counter = 0;
00148       boundary.SetSize(NumOfBdrElements);
00149       for (i = 0; i < mesh.GetNBE(); i++)
00150       {
00151          int edge = mesh.GetBdrElementEdgeIndex(i);
00152          int el1 = edge_element->GetRow(edge)[0];
00153          if (partitioning[el1] == MyRank)
00154          {
00155             boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
00156             int *v = boundary[bdrelem_counter]->GetVertices();
00157             int nv = boundary[bdrelem_counter]->GetNVertices();
00158             for (j = 0; j < nv; j++)
00159                v[j] = vert_global_local[v[j]];
00160             bdrelem_counter++;
00161          }
00162       }
00163    }
00164 
00165    meshgen = mesh.MeshGenerator();
00166 
00167    mesh.attributes.Copy(attributes);
00168    mesh.bdr_attributes.Copy(bdr_attributes);
00169 
00170    // this is called by the default Mesh constructor
00171    // InitTables();
00172 
00173    el_to_edge = new Table;
00174    NumOfEdges = Mesh::GetElementToEdgeTable(*el_to_edge, be_to_edge);
00175 
00176    STable3D *faces_tbl = NULL;
00177    if (Dim == 3)
00178       faces_tbl = GetElementToFaceTable(1);
00179    else
00180       NumOfFaces = 0;
00181    GenerateFaces();
00182 
00183    c_el_to_edge = NULL;
00184 
00185    ListOfIntegerSets  groups;
00186    IntegerSet         group;
00187 
00188    // the first group is the local one
00189    group.Recreate(1, &MyRank);
00190    groups.Insert(group);
00191 
00192 #ifdef MFEM_DEBUG
00193    if (Dim < 3 && mesh.GetNFaces() != 0)
00194    {
00195       cerr << "ParMesh::ParMesh (proc " << MyRank << ") : "
00196          "(Dim < 3 && mesh.GetNFaces() != 0) is true!" << endl;
00197       mfem_error();
00198    }
00199 #endif
00200    // determine shared faces
00201    int sface_counter = 0;
00202    Array<int> face_group(mesh.GetNFaces());
00203    for (i = 0; i < face_group.Size(); i++)
00204    {
00205       int el[2];
00206       face_group[i] = -1;
00207       mesh.GetFaceElements(i, &el[0], &el[1]);
00208       if (el[1] >= 0)
00209       {
00210          el[0] = partitioning[el[0]];
00211          el[1] = partitioning[el[1]];
00212          if ((el[0] == MyRank && el[1] != MyRank) ||
00213              (el[0] != MyRank && el[1] == MyRank))
00214          {
00215             group.Recreate(2, el);
00216             face_group[i] = groups.Insert(group) - 1;
00217             sface_counter++;
00218          }
00219       }
00220    }
00221 
00222    // determine shared edges
00223    int sedge_counter = 0;
00224    if (!edge_element)
00225    {
00226       edge_element = new Table;
00227       Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
00228    }
00229    for (i = 0; i < edge_element->Size(); i++)
00230    {
00231       int me = 0, others = 0;
00232       for (j = edge_element->GetI()[i]; j < edge_element->GetI()[i+1]; j++)
00233       {
00234          edge_element->GetJ()[j] = partitioning[edge_element->GetJ()[j]];
00235          if (edge_element->GetJ()[j] == MyRank)
00236             me = 1;
00237          else
00238             others = 1;
00239       }
00240 
00241       if (me && others)
00242       {
00243          sedge_counter++;
00244          group.Recreate(edge_element->RowSize(i), edge_element->GetRow(i));
00245          edge_element->GetRow(i)[0] = groups.Insert(group) - 1;
00246       }
00247       else
00248          edge_element->GetRow(i)[0] = -1;
00249    }
00250 
00251    // determine shared vertices
00252    int svert_counter = 0;
00253    Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
00254 
00255    for (i = 0; i < vert_element->Size(); i++)
00256    {
00257       int me = 0, others = 0;
00258       for (j = vert_element->GetI()[i]; j < vert_element->GetI()[i+1]; j++)
00259       {
00260          vert_element->GetJ()[j] = partitioning[vert_element->GetJ()[j]];
00261          if (vert_element->GetJ()[j] == MyRank)
00262             me = 1;
00263          else
00264             others = 1;
00265       }
00266 
00267       if (me && others)
00268       {
00269          svert_counter++;
00270          group.Recreate(vert_element->RowSize(i), vert_element->GetRow(i));
00271          vert_element->GetRow(i)[0] = groups.Insert(group) - 1;
00272       }
00273       else
00274          vert_element->GetRow(i)[0] = -1;
00275    }
00276 
00277    // build group_sface
00278    group_sface.MakeI(groups.Size()-1);
00279 
00280    for (i = 0; i < face_group.Size(); i++)
00281       if (face_group[i] >= 0)
00282          group_sface.AddAColumnInRow(face_group[i]);
00283 
00284    group_sface.MakeJ();
00285 
00286    sface_counter = 0;
00287    for (i = 0; i < face_group.Size(); i++)
00288       if (face_group[i] >= 0)
00289          group_sface.AddConnection(face_group[i], sface_counter++);
00290 
00291    group_sface.ShiftUpI();
00292 
00293    // build group_sedge
00294    group_sedge.MakeI(groups.Size()-1);
00295 
00296    for (i = 0; i < edge_element->Size(); i++)
00297       if (edge_element->GetRow(i)[0] >= 0)
00298          group_sedge.AddAColumnInRow(edge_element->GetRow(i)[0]);
00299 
00300    group_sedge.MakeJ();
00301 
00302    sedge_counter = 0;
00303    for (i = 0; i < edge_element->Size(); i++)
00304       if (edge_element->GetRow(i)[0] >= 0)
00305          group_sedge.AddConnection(edge_element->GetRow(i)[0],
00306                                    sedge_counter++);
00307 
00308    group_sedge.ShiftUpI();
00309 
00310    // build group_svert
00311    group_svert.MakeI(groups.Size()-1);
00312 
00313    for (i = 0; i < vert_element->Size(); i++)
00314       if (vert_element->GetRow(i)[0] >= 0)
00315          group_svert.AddAColumnInRow(vert_element->GetRow(i)[0]);
00316 
00317    group_svert.MakeJ();
00318 
00319    svert_counter = 0;
00320    for (i = 0; i < vert_element->Size(); i++)
00321       if (vert_element->GetRow(i)[0] >= 0)
00322          group_svert.AddConnection(vert_element->GetRow(i)[0],
00323                                    svert_counter++);
00324 
00325    group_svert.ShiftUpI();
00326 
00327    // build shared_faces and sface_lface
00328    shared_faces.SetSize(sface_counter);
00329    sface_lface. SetSize(sface_counter);
00330 
00331    if (Dim == 3)
00332    {
00333       sface_counter = 0;
00334       for (i = 0; i < face_group.Size(); i++)
00335          if (face_group[i] >= 0)
00336          {
00337             shared_faces[sface_counter] = mesh.GetFace(i)->Duplicate(this);
00338             int *v = shared_faces[sface_counter]->GetVertices();
00339             int nv = shared_faces[sface_counter]->GetNVertices();
00340             for (j = 0; j < nv; j++)
00341                v[j] = vert_global_local[v[j]];
00342             switch (shared_faces[sface_counter]->GetType())
00343             {
00344             case Element::TRIANGLE:
00345                sface_lface[sface_counter] = (*faces_tbl)(v[0], v[1], v[2]);
00346                // mark the shared face for refinement by reorienting
00347                // it according to the refinement flag in the tetradron
00348                // to which this shared face belongs to.
00349                {
00350                   int lface = sface_lface[sface_counter];
00351                   Tetrahedron *tet =
00352                      (Tetrahedron *)(elements[faces_info[lface].Elem1No]);
00353                   int re[2], type, flag, *tv;
00354                   tet->ParseRefinementFlag(re, type, flag);
00355                   tv = tet->GetVertices();
00356                   switch (faces_info[lface].Elem1Inf/64)
00357                   {
00358                   case 0:
00359                      switch (re[1])
00360                      {
00361                      case 1: v[0] = tv[1]; v[1] = tv[2]; v[2] = tv[3]; break;
00362                      case 4: v[0] = tv[3]; v[1] = tv[1]; v[2] = tv[2]; break;
00363                      case 5: v[0] = tv[2]; v[1] = tv[3]; v[2] = tv[1]; break;
00364                      }
00365                      break;
00366                   case 1:
00367                      switch (re[0])
00368                      {
00369                      case 2: v[0] = tv[2]; v[1] = tv[0]; v[2] = tv[3]; break;
00370                      case 3: v[0] = tv[0]; v[1] = tv[3]; v[2] = tv[2]; break;
00371                      case 5: v[0] = tv[3]; v[1] = tv[2]; v[2] = tv[0]; break;
00372                      }
00373                      break;
00374                   case 2:
00375                      v[0] = tv[0]; v[1] = tv[1]; v[2] = tv[3];
00376                      break;
00377                   case 3:
00378                      v[0] = tv[1]; v[1] = tv[0]; v[2] = tv[2];
00379                      break;
00380                   }
00381                   // flip the shared face in the processor that owns the
00382                   // second element (in 'mesh')
00383                   {
00384                      int gl_el1, gl_el2;
00385                      mesh.GetFaceElements(i, &gl_el1, &gl_el2);
00386                      if (MyRank == partitioning[gl_el2])
00387                      {
00388                         const int t = v[0]; v[0] = v[1]; v[1] = t;
00389                      }
00390                   }
00391                }
00392                break;
00393             case Element::QUADRILATERAL:
00394                sface_lface[sface_counter] =
00395                   (*faces_tbl)(v[0], v[1], v[2], v[3]);
00396                break;
00397             }
00398             sface_counter++;
00399          }
00400 
00401       delete faces_tbl;
00402    }
00403 
00404    // build shared_edges and sedge_ledge
00405    shared_edges.SetSize(sedge_counter);
00406    sedge_ledge. SetSize(sedge_counter);
00407 
00408    {
00409       DSTable v_to_v(NumOfVertices);
00410       GetVertexToVertexTable(v_to_v);
00411 
00412       sedge_counter = 0;
00413       for (i = 0; i < edge_element->Size(); i++)
00414          if (edge_element->GetRow(i)[0] >= 0)
00415          {
00416             mesh.GetEdgeVertices(i, vert);
00417 
00418             shared_edges[sedge_counter] =
00419                new Segment(vert_global_local[vert[0]],
00420                            vert_global_local[vert[1]], 1);
00421 
00422             if ((sedge_ledge[sedge_counter] =
00423                  v_to_v(vert_global_local[vert[0]],
00424                         vert_global_local[vert[1]])) < 0)
00425             {
00426                cerr << "\n\n\n" << MyRank << ": ParMesh::ParMesh: "
00427                     << "ERROR in v_to_v\n\n" << endl;
00428                mfem_error();
00429             }
00430 
00431             sedge_counter++;
00432          }
00433    }
00434 
00435    delete edge_element;
00436 
00437    // build svert_lvert
00438    svert_lvert.SetSize(svert_counter);
00439 
00440    svert_counter = 0;
00441    for (i = 0; i < vert_element->Size(); i++)
00442       if (vert_element->GetRow(i)[0] >= 0)
00443          svert_lvert[svert_counter++] = vert_global_local[i];
00444 
00445    delete vert_element;
00446 
00447    // build the group communication topology
00448    gtopo.Create(groups, 822);
00449 
00450    if (mesh.NURBSext)
00451    {
00452       NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
00453                                        activeBdrElem);
00454    }
00455 
00456    if (mesh.GetNodes()) // curved mesh
00457    {
00458       Nodes = new ParGridFunction(this, mesh.GetNodes());
00459       own_nodes = 1;
00460 
00461       Array<int> gvdofs, lvdofs;
00462       Vector lnodes;
00463       element_counter = 0;
00464       for (i = 0; i < mesh.GetNE(); i++)
00465          if (partitioning[i] == MyRank)
00466          {
00467             Nodes->FESpace()->GetElementVDofs(element_counter, lvdofs);
00468             mesh.GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
00469             mesh.GetNodes()->GetSubVector(gvdofs, lnodes);
00470             Nodes->SetSubVector(lvdofs, lnodes);
00471             element_counter++;
00472          }
00473    }
00474 
00475    if (NURBSext == NULL)
00476       delete [] partitioning;
00477 }
00478 
00479 void ParMesh::GroupEdge(int group, int i, int &edge, int &o)
00480 {
00481    int sedge = group_sedge.GetJ()[group_sedge.GetI()[group-1]+i];
00482    edge = sedge_ledge[sedge];
00483    int *v = shared_edges[sedge]->GetVertices();
00484    o = (v[0] < v[1]) ? (+1) : (-1);
00485 }
00486 
00487 void ParMesh::GroupFace(int group, int i, int &face, int &o)
00488 {
00489    int sface = group_sface.GetJ()[group_sface.GetI()[group-1]+i];
00490    face = sface_lface[sface];
00491    // face gives the base orientation
00492    if (faces[face]->GetType() == Element::TRIANGLE)
00493       o = GetTriOrientation(faces[face]->GetVertices(),
00494                             shared_faces[sface]->GetVertices());
00495    if (faces[face]->GetType() == Element::QUADRILATERAL)
00496       o = GetQuadOrientation(faces[face]->GetVertices(),
00497                              shared_faces[sface]->GetVertices());
00498 }
00499 
00500 // For a line segment with vertices v[0] and v[1], return a number with
00501 // the following meaning:
00502 // 0 - the edge was not refined
00503 // 1 - the edge e was refined once by splitting v[0],v[1]
00504 int ParMesh::GetEdgeSplittings(Element *edge, const DSTable &v_to_v,
00505                                int *middle)
00506 {
00507    int m, *v = edge->GetVertices();
00508 
00509    if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
00510       return 1;
00511    else
00512       return 0;
00513 }
00514 
00515 // For a triangular face with (correctly ordered) vertices v[0], v[1], v[2]
00516 // return a number with the following meaning:
00517 // 0 - the face was not refined
00518 // 1 - the face was refined once by splitting v[0],v[1]
00519 // 2 - the face was refined twice by splitting v[0],v[1] and then v[1],v[2]
00520 // 3 - the face was refined twice by splitting v[0],v[1] and then v[0],v[2]
00521 // 4 - the face was refined three times (as in 2+3)
00522 int ParMesh::GetFaceSplittings(Element *face, const DSTable &v_to_v,
00523                                int *middle)
00524 {
00525    int m, right = 0;
00526    int number_of_splittings = 0;
00527    int *v = face->GetVertices();
00528 
00529    if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
00530    {
00531       number_of_splittings++;
00532       if ((m = v_to_v(v[1], v[2])) != -1 && middle[m] != -1)
00533       {
00534          right = 1;
00535          number_of_splittings++;
00536       }
00537       if ((m = v_to_v(v[2], v[0])) != -1 && middle[m] != -1)
00538          number_of_splittings++;
00539 
00540       switch (number_of_splittings)
00541       {
00542       case 2:
00543          if (right == 0)
00544             number_of_splittings++;
00545          break;
00546       case 3:
00547          number_of_splittings++;
00548          break;
00549       }
00550    }
00551 
00552    return number_of_splittings;
00553 }
00554 
00555 void ParMesh::ReorientTetMesh()
00556 {
00557    if (Dim != 3 || !(meshgen & 1))
00558       return;
00559 
00560    Mesh::ReorientTetMesh();
00561 
00562    int *v;
00563 
00564    // The local edge and face numbering is changed therefore we need to
00565    // update sedge_ledge and sface_lface.
00566    {
00567       DSTable v_to_v(NumOfVertices);
00568       GetVertexToVertexTable(v_to_v);
00569       for (int i = 0; i < shared_edges.Size(); i++)
00570       {
00571          v = shared_edges[i]->GetVertices();
00572          sedge_ledge[i] = v_to_v(v[0], v[1]);
00573       }
00574    }
00575 
00576    // Rotate shared faces and update sface_lface.
00577    // Note that no communication is needed to ensure that the shared
00578    // faces are rotated in the same way in both processors. This is
00579    // automatic due to various things, e.g. the global to local vertex
00580    // mapping preserves the global order; also the way new vertices
00581    // are introduced during refinement is essential.
00582    {
00583       STable3D *faces_tbl = GetFacesTable();
00584       for (int i = 0; i < shared_faces.Size(); i++)
00585          if (shared_faces[i]->GetType() == Element::TRIANGLE)
00586          {
00587             v = shared_faces[i]->GetVertices();
00588 
00589             Rotate3(v[0], v[1], v[2]);
00590 
00591             sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
00592          }
00593       delete faces_tbl;
00594    }
00595 }
00596 
00597 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
00598 {
00599    int i, j, wtls = WantTwoLevelState;
00600 
00601    if (Nodes)  // curved mesh
00602    {
00603       UseTwoLevelState(1);
00604    }
00605 
00606    SetState(Mesh::NORMAL);
00607    DeleteCoarseTables();
00608 
00609    if (Dim == 3)
00610    {
00611       if (WantTwoLevelState)
00612       {
00613          c_NumOfVertices    = NumOfVertices;
00614          c_NumOfEdges       = NumOfEdges;
00615          c_NumOfFaces       = NumOfFaces;
00616          c_NumOfElements    = NumOfElements;
00617          c_NumOfBdrElements = NumOfBdrElements;
00618       }
00619 
00620       int uniform_refinement = 0;
00621       if (type < 0)
00622       {
00623          type = -type;
00624          uniform_refinement = 1;
00625       }
00626 
00627       // 1. Get table of vertex to vertex connections.
00628       DSTable v_to_v(NumOfVertices);
00629       GetVertexToVertexTable(v_to_v);
00630 
00631       // 2. Get edge to element connections in arrays edge1 and edge2
00632       Array<int> middle(v_to_v.NumberOfEntries());
00633       middle = -1;
00634 
00635       // 3. Do the red refinement.
00636       switch (type)
00637       {
00638       case 1:
00639          for (i = 0; i < marked_el.Size(); i++)
00640             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00641          break;
00642       case 2:
00643          for (i = 0; i < marked_el.Size(); i++)
00644          {
00645             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00646 
00647             Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
00648             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00649          }
00650          break;
00651       case 3:
00652          for (i = 0; i < marked_el.Size(); i++)
00653          {
00654             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00655 
00656             j = NumOfElements - 1;
00657             Bisection(j, v_to_v, NULL, NULL, middle);
00658             Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
00659             Bisection(j, v_to_v, NULL, NULL, middle);
00660 
00661             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00662             Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
00663             Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
00664          }
00665          break;
00666       }
00667 
00668       if (WantTwoLevelState)
00669       {
00670          RefinedElement::State = RefinedElement::FINE;
00671          State = Mesh::TWO_LEVEL_FINE;
00672       }
00673 
00674       // 4. Do the green refinement (to get conforming mesh).
00675       int need_refinement;
00676       int refined_edge[5][3] = {{0, 0, 0},
00677                                 {1, 0, 0},
00678                                 {1, 1, 0},
00679                                 {1, 0, 1},
00680                                 {1, 1, 1}};
00681       int faces_in_group, max_faces_in_group = 0;
00682       // face_splittings identify how the shared faces have been split
00683       int **face_splittings = new int*[GetNGroups()-1];
00684       for (i = 0; i < GetNGroups()-1; i++)
00685       {
00686          faces_in_group = GroupNFaces(i+1);
00687          face_splittings[i] = new int[faces_in_group];
00688          if (faces_in_group > max_faces_in_group)
00689             max_faces_in_group = faces_in_group;
00690       }
00691       int neighbor, *iBuf = new int[max_faces_in_group];
00692 
00693       Array<int> group_faces;
00694       Vertex V;
00695 
00696       MPI_Request request;
00697       MPI_Status  status;
00698 
00699 #ifdef MFEM_DEBUG
00700       int ref_loops_all = 0, ref_loops_par = 0;
00701 #endif
00702       do
00703       {
00704          need_refinement = 0;
00705          for (i = 0; i < NumOfElements; i++)
00706          {
00707             if (elements[i]->NeedRefinement(v_to_v, middle))
00708             {
00709                need_refinement = 1;
00710                Bisection(i, v_to_v, NULL, NULL, middle);
00711             }
00712          }
00713 #ifdef MFEM_DEBUG
00714          ref_loops_all++;
00715 #endif
00716 
00717          if (uniform_refinement)
00718             continue;
00719 
00720          // if the mesh is locally conforming start making it globally
00721          // conforming
00722          if (need_refinement == 0)
00723          {
00724 #ifdef MFEM_DEBUG
00725             ref_loops_par++;
00726 #endif
00727             // MPI_Barrier(MyComm);
00728 
00729             // (a) send the type of interface splitting
00730             for (i = 0; i < GetNGroups()-1; i++)
00731             {
00732                group_sface.GetRow(i, group_faces);
00733                faces_in_group = group_faces.Size();
00734                // it is enough to communicate through the faces
00735                if (faces_in_group != 0)
00736                {
00737                   for (j = 0; j < faces_in_group; j++)
00738                      face_splittings[i][j] =
00739                         GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
00740                                           middle);
00741                   const int *nbs = gtopo.GetGroup(i+1);
00742                   if (nbs[0] == 0)
00743                      neighbor = gtopo.GetNeighborRank(nbs[1]);
00744                   else
00745                      neighbor = gtopo.GetNeighborRank(nbs[0]);
00746                   MPI_Isend(face_splittings[i], faces_in_group, MPI_INT,
00747                             neighbor, 0, MyComm, &request);
00748                }
00749             }
00750 
00751             // (b) receive the type of interface splitting
00752             for (i = 0; i < GetNGroups()-1; i++)
00753             {
00754                group_sface.GetRow(i, group_faces);
00755                faces_in_group = group_faces.Size();
00756                if (faces_in_group != 0)
00757                {
00758                   const int *nbs = gtopo.GetGroup(i+1);
00759                   if (nbs[0] == 0)
00760                      neighbor = gtopo.GetNeighborRank(nbs[1]);
00761                   else
00762                      neighbor = gtopo.GetNeighborRank(nbs[0]);
00763                   MPI_Recv(iBuf, faces_in_group, MPI_INT, neighbor,
00764                            MPI_ANY_TAG, MyComm, &status);
00765 
00766                   for (j = 0; j < faces_in_group; j++)
00767                      if (iBuf[j] != face_splittings[i][j])
00768                      {
00769                         int *v = shared_faces[group_faces[j]]->GetVertices();
00770                         for (int k = 0; k < 3; k++)
00771                            if (refined_edge[iBuf[j]][k] == 1 &&
00772                                refined_edge[face_splittings[i][j]][k] == 0)
00773                            {
00774                               int ii = v_to_v(v[k], v[(k+1)%3]);
00775                               if (middle[ii] == -1)
00776                               {
00777                                  need_refinement = 1;
00778                                  middle[ii] = NumOfVertices++;
00779                                  for (int c = 0; c < 3; c++)
00780                                     V(c) = 0.5 * (vertices[v[k]](c) +
00781                                                   vertices[v[(k+1)%3]](c));
00782                                  vertices.Append(V);
00783                               }
00784                            }
00785                      }
00786                }
00787             }
00788 
00789             i = need_refinement;
00790             MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
00791          }
00792       }
00793       while (need_refinement == 1);
00794 
00795 #ifdef MFEM_DEBUG
00796       i = ref_loops_all;
00797       MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
00798       if (MyRank == 0)
00799       {
00800          cout << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
00801               << ref_loops_all << ", ref_loops_par = " << ref_loops_par
00802               << '\n' << endl;
00803       }
00804 #endif
00805 
00806       delete [] iBuf;
00807       for (i = 0; i < GetNGroups()-1; i++)
00808          delete [] face_splittings[i];
00809       delete [] face_splittings;
00810 
00811 
00812       // 5. Update the boundary elements.
00813       do
00814       {
00815          need_refinement = 0;
00816          for (i = 0; i < NumOfBdrElements; i++)
00817             if (boundary[i]->NeedRefinement(v_to_v, middle))
00818             {
00819                need_refinement = 1;
00820                Bisection(i, v_to_v, middle);
00821             }
00822       }
00823       while (need_refinement == 1);
00824 
00825       if (NumOfBdrElements != boundary.Size())
00826          mfem_error("ParMesh::LocalRefinement :"
00827                     " (NumOfBdrElements != boundary.Size())");
00828 
00829       // 5a. Update the groups after refinement.
00830       if (el_to_face != NULL)
00831       {
00832          if (WantTwoLevelState)
00833          {
00834             c_el_to_face = el_to_face;
00835             el_to_face = NULL;
00836             Swap(faces_info, fc_faces_info);
00837          }
00838          RefineGroups(v_to_v, middle);
00839          // GetElementToFaceTable(); // Called by RefineGroups
00840          GenerateFaces();
00841          if (WantTwoLevelState)
00842          {
00843             f_el_to_face = el_to_face;
00844          }
00845       }
00846 
00847       // 6. Un-mark the Pf elements.
00848       int refinement_edges[2], type, flag;
00849       for (i = 0; i < NumOfElements; i++)
00850       {
00851          Element *El = elements[i];
00852          while (El->GetType() == Element::BISECTED)
00853             El = ((BisectedElement *) El)->FirstChild;
00854          ((Tetrahedron *) El)->ParseRefinementFlag(refinement_edges,
00855                                                    type, flag);
00856          if (type == Tetrahedron::TYPE_PF)
00857             ((Tetrahedron *) El)->CreateRefinementFlag(refinement_edges,
00858                                                        Tetrahedron::TYPE_PU,
00859                                                        flag);
00860       }
00861 
00862       // 7. Free the allocated memory.
00863       middle.DeleteAll();
00864 
00865 #ifdef MFEM_DEBUG
00866       CheckElementOrientation();
00867 #endif
00868 
00869       if (el_to_edge != NULL)
00870       {
00871          if (WantTwoLevelState)
00872          {
00873             c_el_to_edge = el_to_edge;
00874             f_el_to_edge = new Table;
00875             c_bel_to_edge = bel_to_edge;
00876             bel_to_edge = NULL;
00877             NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
00878             el_to_edge = f_el_to_edge;
00879             f_bel_to_edge = bel_to_edge;
00880          }
00881          else
00882             NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
00883       }
00884 
00885       if (WantTwoLevelState)
00886       {
00887          f_NumOfVertices    = NumOfVertices;
00888          f_NumOfEdges       = NumOfEdges;
00889          f_NumOfFaces       = NumOfFaces;
00890          f_NumOfElements    = NumOfElements;
00891          f_NumOfBdrElements = NumOfBdrElements;
00892       }
00893    } //  'if (Dim == 3)'
00894 
00895 
00896    if (Dim == 2)
00897    {
00898       if (WantTwoLevelState)
00899       {
00900          c_NumOfVertices    = NumOfVertices;
00901          c_NumOfEdges       = NumOfEdges;
00902          c_NumOfElements    = NumOfElements;
00903          c_NumOfBdrElements = NumOfBdrElements;
00904       }
00905 
00906       int uniform_refinement = 0;
00907       if (type < 0)
00908       {
00909          type = -type;
00910          uniform_refinement = 1;
00911       }
00912 
00913       // 1. Get table of vertex to vertex connections.
00914       DSTable v_to_v(NumOfVertices);
00915       GetVertexToVertexTable(v_to_v);
00916 
00917       // 2. Get edge to element connections in arrays edge1 and edge2
00918       int nedges  = v_to_v.NumberOfEntries();
00919       int *edge1  = new int[nedges];
00920       int *edge2  = new int[nedges];
00921       int *middle = new int[nedges];
00922 
00923       for (i = 0; i < nedges; i++)
00924          edge1[i] = edge2[i] = middle[i] = -1;
00925 
00926       for (i = 0; i < NumOfElements; i++)
00927       {
00928          int *v = elements[i]->GetVertices();
00929          for (j = 0; j < 3; j++)
00930          {
00931             int ind = v_to_v(v[j], v[(j+1)%3]);
00932             (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
00933          }
00934       }
00935 
00936       // 3. Do the red refinement.
00937       for (i = 0; i < marked_el.Size(); i++)
00938          RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
00939 
00940       if (WantTwoLevelState)
00941       {
00942          RefinedElement::State = RefinedElement::FINE;
00943          State = Mesh::TWO_LEVEL_FINE;
00944       }
00945 
00946       // 4. Do the green refinement (to get conforming mesh).
00947       int need_refinement;
00948       int edges_in_group, max_edges_in_group = 0;
00949       // edge_splittings identify how the shared edges have been split
00950       int **edge_splittings = new int*[GetNGroups()-1];
00951       for (i = 0; i < GetNGroups()-1; i++)
00952       {
00953          edges_in_group = GroupNEdges(i+1);
00954          edge_splittings[i] = new int[edges_in_group];
00955          if (edges_in_group > max_edges_in_group)
00956             max_edges_in_group = edges_in_group;
00957       }
00958       int neighbor, *iBuf = new int[max_edges_in_group];
00959 
00960       Array<int> group_edges;
00961 
00962       MPI_Request request;
00963       MPI_Status  status;
00964       Vertex V;
00965       V(2) = 0.0;
00966 
00967 #ifdef MFEM_DEBUG
00968       int ref_loops_all = 0, ref_loops_par = 0;
00969 #endif
00970       do
00971       {
00972          need_refinement = 0;
00973          for (i = 0; i < nedges; i++)
00974             if (middle[i] != -1 && edge1[i] != -1)
00975             {
00976                need_refinement = 1;
00977                GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
00978             }
00979 #ifdef MFEM_DEBUG
00980          ref_loops_all++;
00981 #endif
00982 
00983          if (uniform_refinement)
00984             continue;
00985 
00986          // if the mesh is locally conforming start making it globally
00987          // conforming
00988          if (need_refinement == 0)
00989          {
00990 #ifdef MFEM_DEBUG
00991             ref_loops_par++;
00992 #endif
00993             // MPI_Barrier(MyComm);
00994 
00995             // (a) send the type of interface splitting
00996             for (i = 0; i < GetNGroups()-1; i++)
00997             {
00998                group_sedge.GetRow(i, group_edges);
00999                edges_in_group = group_edges.Size();
01000                // it is enough to communicate through the edges
01001                if (edges_in_group != 0)
01002                {
01003                   for (j = 0; j < edges_in_group; j++)
01004                      edge_splittings[i][j] =
01005                         GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
01006                                           middle);
01007                   const int *nbs = gtopo.GetGroup(i+1);
01008                   if (nbs[0] == 0)
01009                      neighbor = gtopo.GetNeighborRank(nbs[1]);
01010                   else
01011                      neighbor = gtopo.GetNeighborRank(nbs[0]);
01012                   MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
01013                             neighbor, 0, MyComm, &request);
01014                }
01015             }
01016 
01017             // (b) receive the type of interface splitting
01018             for (i = 0; i < GetNGroups()-1; i++)
01019             {
01020                group_sedge.GetRow(i, group_edges);
01021                edges_in_group = group_edges.Size();
01022                if (edges_in_group != 0)
01023                {
01024                   const int *nbs = gtopo.GetGroup(i+1);
01025                   if (nbs[0] == 0)
01026                      neighbor = gtopo.GetNeighborRank(nbs[1]);
01027                   else
01028                      neighbor = gtopo.GetNeighborRank(nbs[0]);
01029                   MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
01030                            MPI_ANY_TAG, MyComm, &status);
01031 
01032                   for (j = 0; j < edges_in_group; j++)
01033                      if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
01034                      {
01035                         int *v = shared_edges[group_edges[j]]->GetVertices();
01036                         int ii = v_to_v(v[0], v[1]);
01037 #ifdef MFEM_DEBUG
01038                         if (middle[ii] != -1)
01039                            mfem_error("ParMesh::LocalRefinement (triangles) : "
01040                                       "Oops!");
01041 #endif
01042                         need_refinement = 1;
01043                         middle[ii] = NumOfVertices++;
01044                         for (int c = 0; c < 2; c++)
01045                            V(c) = 0.5 * (vertices[v[0]](c) + vertices[v[1]](c));
01046                         vertices.Append(V);
01047                      }
01048                }
01049             }
01050 
01051             i = need_refinement;
01052             MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
01053          }
01054       }
01055       while (need_refinement == 1);
01056 
01057 #ifdef MFEM_DEBUG
01058       i = ref_loops_all;
01059       MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
01060       if (MyRank == 0)
01061       {
01062          cout << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
01063               << ref_loops_all << ", ref_loops_par = " << ref_loops_par
01064               << '\n' << endl;
01065       }
01066 #endif
01067 
01068       for (i = 0; i < GetNGroups()-1; i++)
01069          delete [] edge_splittings[i];
01070       delete [] edge_splittings;
01071 
01072       delete [] iBuf;
01073 
01074       // 5. Update the boundary elements.
01075       int v1[2], v2[2], bisect, temp;
01076       temp = NumOfBdrElements;
01077       for (i = 0; i < temp; i++)
01078       {
01079          int *v = boundary[i]->GetVertices();
01080          bisect = v_to_v(v[0], v[1]);
01081          if (middle[bisect] != -1)
01082          {  // the element was refined (needs updating)
01083             if (boundary[i]->GetType() == Element::SEGMENT)
01084             {
01085                v1[0] =           v[0]; v1[1] = middle[bisect];
01086                v2[0] = middle[bisect]; v2[1] =           v[1];
01087 
01088                if (WantTwoLevelState)
01089                {
01090                   boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
01091 #ifdef MFEM_USE_MEMALLOC
01092                   BisectedElement *aux = BEMemory.Alloc();
01093                   aux->SetCoarseElem(boundary[i]);
01094 #else
01095                   BisectedElement *aux = new BisectedElement(boundary[i]);
01096 #endif
01097                   aux->FirstChild =
01098                      new Segment(v1, boundary[i]->GetAttribute());
01099                   aux->SecondChild = NumOfBdrElements;
01100                   boundary[i] = aux;
01101                   NumOfBdrElements++;
01102                }
01103                else
01104                {
01105                   boundary[i]->SetVertices(v1);
01106                   boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
01107                }
01108             }
01109             else
01110                mfem_error("Only bisection of segment is implemented for bdr"
01111                           " elem.");
01112          }
01113       }
01114       NumOfBdrElements = boundary.Size();
01115 
01116       // 5a. Update the groups after refinement.
01117       RefineGroups(v_to_v, middle);
01118 
01119       // 6. Free the allocated memory.
01120       delete [] edge1;
01121       delete [] edge2;
01122       delete [] middle;
01123 
01124 #ifdef MFEM_DEBUG
01125       CheckElementOrientation();
01126 #endif
01127 
01128       if (WantTwoLevelState)
01129       {
01130          f_NumOfVertices    = NumOfVertices;
01131          f_NumOfElements    = NumOfElements;
01132          f_NumOfBdrElements = NumOfBdrElements;
01133          RefinedElement::State = RefinedElement::FINE;
01134          State = Mesh::TWO_LEVEL_FINE;
01135       }
01136 
01137       if (el_to_edge != NULL)
01138       {
01139          if (WantTwoLevelState)
01140          {
01141             c_el_to_edge = el_to_edge;
01142             Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
01143             f_el_to_edge = new Table;
01144             NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
01145             el_to_edge = f_el_to_edge;
01146             f_NumOfEdges = NumOfEdges;
01147          }
01148          else
01149             NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
01150          GenerateFaces();
01151       }
01152    } //  'if (Dim == 2)'
01153 
01154    if (Nodes)  // curved mesh
01155    {
01156       UpdateNodes();
01157       UseTwoLevelState(wtls);
01158    }
01159 }
01160 
01161 void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
01162 {
01163    int i, attr, newv[3], ind, f_ind, *v;
01164 
01165    int group;
01166    Array<int> group_verts, group_edges, group_faces;
01167 
01168    // To update the groups after a refinement, we observe that:
01169    // - every (new and old) vertex, edge and face belongs to exactly one group
01170    // - the refinement does not create new groups
01171    // - a new vertex appears only as the middle of a refined edge
01172    // - a face can be refined 2, 3 or 4 times producing new edges and faces
01173 
01174    int *I_group_svert, *J_group_svert;
01175    int *I_group_sedge, *J_group_sedge;
01176    int *I_group_sface, *J_group_sface;
01177 
01178    I_group_svert = new int[GetNGroups()+1];
01179    I_group_sedge = new int[GetNGroups()+1];
01180    if (Dim == 3)
01181       I_group_sface = new int[GetNGroups()+1];
01182 
01183    I_group_svert[0] = I_group_svert[1] = 0;
01184    I_group_sedge[0] = I_group_sedge[1] = 0;
01185    if (Dim == 3)
01186       I_group_sface[0] = I_group_sface[1] = 0;
01187 
01188    // overestimate the size of the J arrays
01189    if (Dim == 3)
01190    {
01191       J_group_svert = new int[group_svert.Size_of_connections()
01192                               + group_sedge.Size_of_connections()];
01193       J_group_sedge = new int[2*group_sedge.Size_of_connections()
01194                               + 3*group_sface.Size_of_connections()];
01195       J_group_sface = new int[4*group_sface.Size_of_connections()];
01196    }
01197    else if (Dim == 2)
01198    {
01199       J_group_svert = new int[group_svert.Size_of_connections()
01200                               + group_sedge.Size_of_connections()];
01201       J_group_sedge = new int[2*group_sedge.Size_of_connections()];
01202    }
01203 
01204    for (group = 0; group < GetNGroups()-1; group++)
01205    {
01206       // Get the group shared objects
01207       group_svert.GetRow(group, group_verts);
01208       group_sedge.GetRow(group, group_edges);
01209       group_sface.GetRow(group, group_faces);
01210 
01211       // Check which edges have been refined
01212       for (i = 0; i < group_sedge.RowSize(group); i++)
01213       {
01214          v = shared_edges[group_edges[i]]->GetVertices();
01215          ind = middle[v_to_v(v[0], v[1])];
01216          if (ind != -1)
01217          {
01218             // add a vertex
01219             group_verts.Append(svert_lvert.Append(ind)-1);
01220             // update the edges
01221             attr = shared_edges[group_edges[i]]->GetAttribute();
01222             shared_edges.Append(new Segment(v[1], ind, attr));
01223             group_edges.Append(sedge_ledge.Append(-1)-1);
01224             v[1] = ind;
01225          }
01226       }
01227 
01228       // Check which faces have been refined
01229       for (i = 0; i < group_sface.RowSize(group); i++)
01230       {
01231          v = shared_faces[group_faces[i]]->GetVertices();
01232          ind = middle[v_to_v(v[0], v[1])];
01233          if (ind != -1)
01234          {
01235             attr = shared_faces[group_faces[i]]->GetAttribute();
01236             // add the refinement edge
01237             shared_edges.Append(new Segment(v[2], ind, attr));
01238             group_edges.Append(sedge_ledge.Append(-1)-1);
01239             // add a face
01240             f_ind = group_faces.Size();
01241             shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
01242             group_faces.Append(sface_lface.Append(-1)-1);
01243             newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
01244             shared_faces[group_faces[i]]->SetVertices(newv);
01245 
01246             // check if the left face has also been refined
01247             // v = shared_faces[group_faces[i]]->GetVertices();
01248             ind = middle[v_to_v(v[0], v[1])];
01249             if (ind != -1)
01250             {
01251                // add the refinement edge
01252                shared_edges.Append(new Segment(v[2], ind, attr));
01253                group_edges.Append(sedge_ledge.Append(-1)-1);
01254                // add a face
01255                shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
01256                group_faces.Append(sface_lface.Append(-1)-1);
01257                newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
01258                shared_faces[group_faces[i]]->SetVertices(newv);
01259             }
01260 
01261             // check if the right face has also been refined
01262             v = shared_faces[group_faces[f_ind]]->GetVertices();
01263             ind = middle[v_to_v(v[0], v[1])];
01264             if (ind != -1)
01265             {
01266                // add the refinement edge
01267                shared_edges.Append(new Segment(v[2], ind, attr));
01268                group_edges.Append(sedge_ledge.Append(-1)-1);
01269                // add a face
01270                shared_faces.Append(new Triangle(v[1], v[2], ind, attr));
01271                group_faces.Append(sface_lface.Append(-1)-1);
01272                newv[0] = v[2]; newv[1] = v[0]; newv[2] = ind;
01273                shared_faces[group_faces[f_ind]]->SetVertices(newv);
01274             }
01275          }
01276       }
01277 
01278       I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
01279       I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
01280       if (Dim == 3)
01281          I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
01282 
01283       int *J;
01284       J = J_group_svert+I_group_svert[group];
01285       for (i = 0; i < group_verts.Size(); i++)
01286          J[i] = group_verts[i];
01287       J = J_group_sedge+I_group_sedge[group];
01288       for (i = 0; i < group_edges.Size(); i++)
01289          J[i] = group_edges[i];
01290       if (Dim == 3)
01291       {
01292          J = J_group_sface+I_group_sface[group];
01293          for (i = 0; i < group_faces.Size(); i++)
01294             J[i] = group_faces[i];
01295       }
01296    }
01297 
01298    // Fix the local numbers of shared edges and faces
01299    {
01300       DSTable new_v_to_v(NumOfVertices);
01301       GetVertexToVertexTable(new_v_to_v);
01302       for (i = 0; i < shared_edges.Size(); i++)
01303       {
01304          v = shared_edges[i]->GetVertices();
01305          sedge_ledge[i] = new_v_to_v(v[0], v[1]);
01306       }
01307    }
01308    if (Dim == 3)
01309    {
01310       STable3D *faces_tbl = GetElementToFaceTable(1);
01311       for (i = 0; i < shared_faces.Size(); i++)
01312       {
01313          v = shared_faces[i]->GetVertices();
01314          sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
01315       }
01316       delete faces_tbl;
01317    }
01318 
01319    group_svert.SetIJ(I_group_svert, J_group_svert);
01320    group_sedge.SetIJ(I_group_sedge, J_group_sedge);
01321    if (Dim == 3)
01322       group_sface.SetIJ(I_group_sface, J_group_sface);
01323 }
01324 
01325 void ParMesh::QuadUniformRefinement()
01326 {
01327    SetState(Mesh::NORMAL);
01328 
01329    int oedge = NumOfVertices, wtls = WantTwoLevelState;
01330 
01331    if (Nodes)  // curved mesh
01332       UseTwoLevelState(1);
01333 
01334    // call Mesh::QuadUniformRefinement so that it won't update the nodes
01335    {
01336       GridFunction *nodes = Nodes;
01337       Nodes = NULL;
01338       Mesh::QuadUniformRefinement();
01339       Nodes = nodes;
01340    }
01341 
01342    // update the groups
01343    {
01344       int i, attr, ind, *v;
01345 
01346       int group;
01347       Array<int> sverts, sedges;
01348 
01349       int *I_group_svert, *J_group_svert;
01350       int *I_group_sedge, *J_group_sedge;
01351 
01352       I_group_svert = new int[GetNGroups()+1];
01353       I_group_sedge = new int[GetNGroups()+1];
01354 
01355       I_group_svert[0] = I_group_svert[1] = 0;
01356       I_group_sedge[0] = I_group_sedge[1] = 0;
01357 
01358       // compute the size of the J arrays
01359       J_group_svert = new int[group_svert.Size_of_connections()
01360                               + group_sedge.Size_of_connections()];
01361       J_group_sedge = new int[2*group_sedge.Size_of_connections()];
01362 
01363       for (group = 0; group < GetNGroups()-1; group++)
01364       {
01365          // Get the group shared objects
01366          group_svert.GetRow(group, sverts);
01367          group_sedge.GetRow(group, sedges);
01368 
01369          // Process all the edges
01370          for (i = 0; i < group_sedge.RowSize(group); i++)
01371          {
01372             v = shared_edges[sedges[i]]->GetVertices();
01373             ind = oedge + sedge_ledge[sedges[i]];
01374             // add a vertex
01375             sverts.Append(svert_lvert.Append(ind)-1);
01376             // update the edges
01377             attr = shared_edges[sedges[i]]->GetAttribute();
01378             shared_edges.Append(new Segment(v[1], ind, attr));
01379             sedges.Append(sedge_ledge.Append(-1)-1);
01380             v[1] = ind;
01381          }
01382 
01383          I_group_svert[group+1] = I_group_svert[group] + sverts.Size();
01384          I_group_sedge[group+1] = I_group_sedge[group] + sedges.Size();
01385 
01386          int *J;
01387          J = J_group_svert+I_group_svert[group];
01388          for (i = 0; i < sverts.Size(); i++)
01389             J[i] = sverts[i];
01390          J = J_group_sedge+I_group_sedge[group];
01391          for (i = 0; i < sedges.Size(); i++)
01392             J[i] = sedges[i];
01393       }
01394 
01395       // Fix the local numbers of shared edges
01396       DSTable v_to_v(NumOfVertices);
01397       GetVertexToVertexTable(v_to_v);
01398       for (i = 0; i < shared_edges.Size(); i++)
01399       {
01400          v = shared_edges[i]->GetVertices();
01401          sedge_ledge[i] = v_to_v(v[0], v[1]);
01402       }
01403 
01404       group_svert.SetIJ(I_group_svert, J_group_svert);
01405       group_sedge.SetIJ(I_group_sedge, J_group_sedge);
01406    }
01407 
01408    if (Nodes)  // curved mesh
01409    {
01410       UpdateNodes();
01411       UseTwoLevelState(wtls);
01412    }
01413 }
01414 
01415 void ParMesh::HexUniformRefinement()
01416 {
01417    SetState(Mesh::NORMAL);
01418 
01419    int wtls = WantTwoLevelState;
01420    int oedge = NumOfVertices;
01421    int oface = oedge + NumOfEdges;
01422 
01423    DSTable v_to_v(NumOfVertices);
01424    GetVertexToVertexTable(v_to_v);
01425    STable3D *faces_tbl = GetFacesTable();
01426 
01427    if (Nodes)  // curved mesh
01428       UseTwoLevelState(1);
01429 
01430    // call Mesh::HexUniformRefinement so that it won't update the nodes
01431    {
01432       GridFunction *nodes = Nodes;
01433       Nodes = NULL;
01434       Mesh::HexUniformRefinement();
01435       Nodes = nodes;
01436    }
01437 
01438    // update the groups
01439    {
01440       int i, attr, newv[4], ind, m[5];
01441       Array<int> v;
01442 
01443       int group;
01444       Array<int> group_verts, group_edges, group_faces;
01445 
01446       int *I_group_svert, *J_group_svert;
01447       int *I_group_sedge, *J_group_sedge;
01448       int *I_group_sface, *J_group_sface;
01449 
01450       I_group_svert = new int[GetNGroups()+1];
01451       I_group_sedge = new int[GetNGroups()+1];
01452       I_group_sface = new int[GetNGroups()+1];
01453 
01454       I_group_svert[0] = I_group_svert[1] = 0;
01455       I_group_sedge[0] = I_group_sedge[1] = 0;
01456       I_group_sface[0] = I_group_sface[1] = 0;
01457 
01458       // compute the size of the J arrays
01459       J_group_svert = new int[group_svert.Size_of_connections()
01460                               + group_sedge.Size_of_connections()
01461                               + group_sface.Size_of_connections()];
01462       J_group_sedge = new int[2*group_sedge.Size_of_connections()
01463                               + 4*group_sface.Size_of_connections()];
01464       J_group_sface = new int[4*group_sface.Size_of_connections()];
01465 
01466       for (group = 0; group < GetNGroups()-1; group++)
01467       {
01468          // Get the group shared objects
01469          group_svert.GetRow(group, group_verts);
01470          group_sedge.GetRow(group, group_edges);
01471          group_sface.GetRow(group, group_faces);
01472 
01473          // Process the edges that have been refined
01474          for (i = 0; i < group_sedge.RowSize(group); i++)
01475          {
01476             shared_edges[group_edges[i]]->GetVertices(v);
01477             ind = oedge + v_to_v(v[0], v[1]);
01478             // add a vertex
01479             group_verts.Append(svert_lvert.Append(ind)-1);
01480             // update the edges
01481             attr = shared_edges[group_edges[i]]->GetAttribute();
01482             shared_edges.Append(new Segment(v[1], ind, attr));
01483             group_edges.Append(sedge_ledge.Append(-1)-1);
01484             newv[0] = v[0]; newv[1] = ind;
01485             shared_edges[group_edges[i]]->SetVertices(newv);
01486          }
01487 
01488          // Process the faces that have been refined
01489          for (i = 0; i < group_sface.RowSize(group); i++)
01490          {
01491             shared_faces[group_faces[i]]->GetVertices(v);
01492             m[0] = oface+(*faces_tbl)(v[0], v[1], v[2], v[3]);
01493             // add a vertex
01494             group_verts.Append(svert_lvert.Append(m[0])-1);
01495             // add the refinement edges
01496             attr = shared_faces[group_faces[i]]->GetAttribute();
01497             m[1] = oedge + v_to_v(v[0], v[1]);
01498             m[2] = oedge + v_to_v(v[1], v[2]);
01499             m[3] = oedge + v_to_v(v[2], v[3]);
01500             m[4] = oedge + v_to_v(v[3], v[0]);
01501             shared_edges.Append(new Segment(m[1], m[0], attr));
01502             group_edges.Append(sedge_ledge.Append(-1)-1);
01503             shared_edges.Append(new Segment(m[2], m[0], attr));
01504             group_edges.Append(sedge_ledge.Append(-1)-1);
01505             shared_edges.Append(new Segment(m[3], m[0], attr));
01506             group_edges.Append(sedge_ledge.Append(-1)-1);
01507             shared_edges.Append(new Segment(m[4], m[0], attr));
01508             group_edges.Append(sedge_ledge.Append(-1)-1);
01509             // update faces
01510             newv[0] = v[0]; newv[1] = m[1]; newv[2] = m[0]; newv[3] = m[4];
01511             shared_faces[group_faces[i]]->SetVertices(newv);
01512             shared_faces.Append(new Quadrilateral(m[1],v[1],m[2],m[0],attr));
01513             group_faces.Append(sface_lface.Append(-1)-1);
01514             shared_faces.Append(new Quadrilateral(m[0],m[2],v[2],m[3],attr));
01515             group_faces.Append(sface_lface.Append(-1)-1);
01516             shared_faces.Append(new Quadrilateral(m[4],m[0],m[3],v[3],attr));
01517             group_faces.Append(sface_lface.Append(-1)-1);
01518          }
01519 
01520          I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
01521          I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
01522          I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
01523 
01524          int *J;
01525          J = J_group_svert+I_group_svert[group];
01526          for (i = 0; i < group_verts.Size(); i++)
01527             J[i] = group_verts[i];
01528          J = J_group_sedge+I_group_sedge[group];
01529          for (i = 0; i < group_edges.Size(); i++)
01530             J[i] = group_edges[i];
01531          J = J_group_sface+I_group_sface[group];
01532          for (i = 0; i < group_faces.Size(); i++)
01533             J[i] = group_faces[i];
01534       }
01535 
01536       // Fix the local numbers of shared edges and faces
01537       DSTable new_v_to_v(NumOfVertices);
01538       GetVertexToVertexTable(new_v_to_v);
01539       for (i = 0; i < shared_edges.Size(); i++)
01540       {
01541          shared_edges[i]->GetVertices(v);
01542          sedge_ledge[i] = new_v_to_v(v[0], v[1]);
01543       }
01544 
01545       delete faces_tbl;
01546       faces_tbl = GetFacesTable();
01547       for (i = 0; i < shared_faces.Size(); i++)
01548       {
01549          shared_faces[i]->GetVertices(v);
01550          sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
01551       }
01552       delete faces_tbl;
01553 
01554       group_svert.SetIJ(I_group_svert, J_group_svert);
01555       group_sedge.SetIJ(I_group_sedge, J_group_sedge);
01556       group_sface.SetIJ(I_group_sface, J_group_sface);
01557    }
01558 
01559    if (Nodes)  // curved mesh
01560    {
01561       UpdateNodes();
01562       UseTwoLevelState(wtls);
01563    }
01564 }
01565 
01566 void ParMesh::NURBSUniformRefinement()
01567 {
01568    if (MyRank == 0)
01569       cout << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
01570 }
01571 
01572 void ParMesh::PrintXG(ostream &out) const
01573 {
01574    if (Dim == 3 && meshgen == 1)
01575    {
01576       int i, j, nv;
01577       const int *ind;
01578 
01579       out << "NETGEN_Neutral_Format\n";
01580       // print the vertices
01581       out << NumOfVertices << '\n';
01582       for (i = 0; i < NumOfVertices; i++)
01583       {
01584          for (j = 0; j < Dim; j++)
01585             out << " " << vertices[i](j);
01586          out << '\n';
01587       }
01588 
01589       // print the elements
01590       out << NumOfElements << '\n';
01591       for (i = 0; i < NumOfElements; i++)
01592       {
01593          nv = elements[i]->GetNVertices();
01594          ind = elements[i]->GetVertices();
01595          out << elements[i]->GetAttribute();
01596          for (j = 0; j < nv; j++)
01597             out << " " << ind[j]+1;
01598          out << '\n';
01599       }
01600 
01601       // print the boundary + shared faces information
01602       out << NumOfBdrElements + shared_faces.Size() << '\n';
01603       // boundary
01604       for (i = 0; i < NumOfBdrElements; i++)
01605       {
01606          nv = boundary[i]->GetNVertices();
01607          ind = boundary[i]->GetVertices();
01608          out << boundary[i]->GetAttribute();
01609          for (j = 0; j < nv; j++)
01610             out << " " << ind[j]+1;
01611          out << '\n';
01612       }
01613       // shared faces
01614       for (i = 0; i < shared_faces.Size(); i++)
01615       {
01616          nv = shared_faces[i]->GetNVertices();
01617          ind = shared_faces[i]->GetVertices();
01618          out << shared_faces[i]->GetAttribute();
01619          for (j = 0; j < nv; j++)
01620             out << " " << ind[j]+1;
01621          out << '\n';
01622       }
01623    }
01624 
01625    if (Dim == 3 && meshgen == 2)
01626    {
01627       int i, j, nv;
01628       const int *ind;
01629 
01630       out << "TrueGrid\n"
01631           << "1 " << NumOfVertices << " " << NumOfElements << " 0 0 0 0 0 0 0\n"
01632           << "0 0 0 1 0 0 0 0 0 0 0\n"
01633           << "0 0 " << NumOfBdrElements+shared_faces.Size()
01634           << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
01635           << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
01636           << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
01637 
01638       // print the vertices
01639       for (i = 0; i < NumOfVertices; i++)
01640          out << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
01641              << " " << vertices[i](2) << " 0.0\n";
01642 
01643       // print the elements
01644       for (i = 0; i < NumOfElements; i++)
01645       {
01646          nv = elements[i]->GetNVertices();
01647          ind = elements[i]->GetVertices();
01648          out << i+1 << " " << elements[i]->GetAttribute();
01649          for (j = 0; j < nv; j++)
01650             out << " " << ind[j]+1;
01651          out << '\n';
01652       }
01653 
01654       // print the boundary information
01655       for (i = 0; i < NumOfBdrElements; i++)
01656       {
01657          nv = boundary[i]->GetNVertices();
01658          ind = boundary[i]->GetVertices();
01659          out << boundary[i]->GetAttribute();
01660          for (j = 0; j < nv; j++)
01661             out << " " << ind[j]+1;
01662          out << " 1.0 1.0 1.0 1.0\n";
01663       }
01664 
01665       // print the shared faces information
01666       for (i = 0; i < shared_faces.Size(); i++)
01667       {
01668          nv = shared_faces[i]->GetNVertices();
01669          ind = shared_faces[i]->GetVertices();
01670          out << shared_faces[i]->GetAttribute();
01671          for (j = 0; j < nv; j++)
01672             out << " " << ind[j]+1;
01673          out << " 1.0 1.0 1.0 1.0\n";
01674       }
01675    }
01676 
01677    if (Dim == 2)
01678    {
01679       int i, j, attr;
01680       Array<int> v;
01681 
01682       out << "areamesh2\n\n";
01683 
01684       // print the boundary + shared edges information
01685       out << NumOfBdrElements + shared_edges.Size() << '\n';
01686       // boundary
01687       for (i = 0; i < NumOfBdrElements; i++)
01688       {
01689          attr = boundary[i]->GetAttribute();
01690          boundary[i]->GetVertices(v);
01691          out << attr << "     ";
01692          for (j = 0; j < v.Size(); j++)
01693             out << v[j] + 1 << "   ";
01694          out << '\n';
01695       }
01696       // shared edges
01697       for (i = 0; i < shared_edges.Size(); i++)
01698       {
01699          attr = shared_edges[i]->GetAttribute();
01700          shared_edges[i]->GetVertices(v);
01701          out << attr << "     ";
01702          for (j = 0; j < v.Size(); j++)
01703             out << v[j] + 1 << "   ";
01704          out << '\n';
01705       }
01706 
01707       // print the elements
01708       out << NumOfElements << '\n';
01709       for (i = 0; i < NumOfElements; i++)
01710       {
01711          attr = elements[i]->GetAttribute();
01712          elements[i]->GetVertices(v);
01713 
01714          out << attr << "   ";
01715          if ((j = GetElementType(i)) == Element::TRIANGLE)
01716             out << 3 << "   ";
01717          else
01718             if (j == Element::QUADRILATERAL)
01719                out << 4 << "   ";
01720             else
01721                if (j == Element::SEGMENT)
01722                   out << 2 << "   ";
01723          for (j = 0; j < v.Size(); j++)
01724             out << v[j] + 1 << "  ";
01725          out << '\n';
01726       }
01727 
01728       // print the vertices
01729       out << NumOfVertices << '\n';
01730       for (i = 0; i < NumOfVertices; i++)
01731       {
01732          for (j = 0; j < Dim; j++)
01733             out << vertices[i](j) << " ";
01734          out << '\n';
01735       }
01736    }
01737 }
01738 
01739 void ParMesh::Print(ostream &out) const
01740 {
01741    int i, j, shared_bdr_attr;
01742    const Array<Element *> &shared_bdr =
01743       (Dim == 3) ? shared_faces : shared_edges;
01744 
01745    if (NURBSext)
01746    {
01747       Mesh::Print(out); // does not print shared boundary
01748       return;
01749    }
01750 
01751    out << "MFEM mesh v1.0\n";
01752 
01753    // optional
01754    out <<
01755       "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
01756       "# POINT       = 0\n"
01757       "# SEGMENT     = 1\n"
01758       "# TRIANGLE    = 2\n"
01759       "# SQUARE      = 3\n"
01760       "# TETRAHEDRON = 4\n"
01761       "# CUBE        = 5\n"
01762       "#\n";
01763 
01764    out << "\ndimension\n" << Dim
01765        << "\n\nelements\n" << NumOfElements << '\n';
01766    for (i = 0; i < NumOfElements; i++)
01767       PrintElement(elements[i], out);
01768 
01769    out << "\nboundary\n" << NumOfBdrElements + shared_bdr.Size() << '\n';
01770    for (i = 0; i < NumOfBdrElements; i++)
01771       PrintElement(boundary[i], out);
01772 
01773    shared_bdr_attr = bdr_attributes.Max() + MyRank + 1;
01774    for (i = 0; i < shared_bdr.Size(); i++)
01775    {
01776       shared_bdr[i]->SetAttribute(shared_bdr_attr);
01777       PrintElement(shared_bdr[i], out);
01778    }
01779    out << "\nvertices\n" << NumOfVertices << '\n';
01780    if (Nodes == NULL)
01781    {
01782       out << Dim << '\n';
01783       for (i = 0; i < NumOfVertices; i++)
01784       {
01785          out << vertices[i](0);
01786          for (j = 1; j < Dim; j++)
01787             out << ' ' << vertices[i](j);
01788          out << '\n';
01789       }
01790    }
01791    else
01792    {
01793       out << "\nnodes\n";
01794       Nodes->Save(out);
01795    }
01796 }
01797 
01798 void ParMesh::PrintAsOne(ostream &out)
01799 {
01800    int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
01801    const int *v;
01802    MPI_Status status;
01803    Array<double> vert;
01804    Array<int> ints;
01805 
01806    if (MyRank == 0)
01807    {
01808       out << "MFEM mesh v1.0\n";
01809 
01810       // optional
01811       out <<
01812          "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
01813          "# POINT       = 0\n"
01814          "# SEGMENT     = 1\n"
01815          "# TRIANGLE    = 2\n"
01816          "# SQUARE      = 3\n"
01817          "# TETRAHEDRON = 4\n"
01818          "# CUBE        = 5\n"
01819          "#\n";
01820 
01821       out << "\ndimension\n" << Dim;
01822    }
01823 
01824    nv = NumOfElements;
01825    MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
01826    if (MyRank == 0)
01827    {
01828       out << "\n\nelements\n" << ne << '\n';
01829       for (i = 0; i < NumOfElements; i++)
01830       {
01831          // processor number + 1 as attribute and geometry type
01832          out << 1 << ' ' << elements[i]->GetGeometryType();
01833          // vertices
01834          nv = elements[i]->GetNVertices();
01835          v  = elements[i]->GetVertices();
01836          for (j = 0; j < nv; j++)
01837             out << ' ' << v[j];
01838          out << '\n';
01839       }
01840       vc = NumOfVertices;
01841       for (p = 1; p < NRanks; p++)
01842       {
01843          MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
01844          ints.SetSize(ne);
01845          MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
01846          for (i = 0; i < ne; )
01847          {
01848             // processor number + 1 as attribute and geometry type
01849             out << p+1 << ' ' << ints[i];
01850             // vertices
01851             k = Geometries.GetVertices(ints[i++])->GetNPoints();
01852             for (j = 0; j < k; j++)
01853                out << ' ' << vc + ints[i++];
01854             out << '\n';
01855          }
01856          vc += nv;
01857       }
01858    }
01859    else
01860    {
01861       // for each element send its geometry type and its vertices
01862       ne = 0;
01863       for (i = 0; i < NumOfElements; i++)
01864          ne += 1 + elements[i]->GetNVertices();
01865       nv = NumOfVertices;
01866       MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
01867       ints.SetSize(ne);
01868       for (i = j = 0; i < NumOfElements; i++)
01869       {
01870          ints[j++] = elements[i]->GetGeometryType();
01871          nv = elements[i]->GetNVertices();
01872          v  = elements[i]->GetVertices();
01873          for (k = 0; k < nv; k++)
01874             ints[j++] = v[k];
01875       }
01876       MPI_Send(&ints[0], ne, MPI_INT, 0, 445, MyComm);
01877    }
01878 
01879    // boundary + shared boundary
01880    Array<Element *> &shared_boundary =
01881       (Dim == 2) ? shared_edges : shared_faces;
01882    nv = NumOfBdrElements + shared_boundary.Size();
01883    MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
01884    if (MyRank == 0)
01885    {
01886       out << "\nboundary\n" << ne << '\n';
01887       // actual boundary
01888       for (i = 0; i < NumOfBdrElements; i++)
01889       {
01890          // processor number + 1 as bdr. attr. and bdr. geometry type
01891          out << 1 << ' ' << boundary[i]->GetGeometryType();
01892          // vertices
01893          nv = boundary[i]->GetNVertices();
01894          v  = boundary[i]->GetVertices();
01895          for (j = 0; j < nv; j++)
01896             out << ' ' << v[j];
01897          out << '\n';
01898       }
01899       // shared boundary (interface)
01900       for (i = 0; i < shared_boundary.Size(); i++)
01901       {
01902          // processor number + 1 as bdr. attr. and bdr. geometry type
01903          out << 1 << ' ' << shared_boundary[i]->GetGeometryType();
01904          // vertices
01905          nv = shared_boundary[i]->GetNVertices();
01906          v  = shared_boundary[i]->GetVertices();
01907          for (j = 0; j < nv; j++)
01908             out << ' ' << v[j];
01909          out << '\n';
01910       }
01911       vc = NumOfVertices;
01912       for (p = 1; p < NRanks; p++)
01913       {
01914          MPI_Recv(nv_ne, 2, MPI_INT, p, 446, MyComm, &status);
01915          ints.SetSize(ne);
01916          MPI_Recv(&ints[0], ne, MPI_INT, p, 447, MyComm, &status);
01917          for (i = 0; i < ne; )
01918          {
01919             // processor number + 1 as bdr. attr. and bdr. geometry type
01920             out << p+1 << ' ' << ints[i];
01921             k = Geometries.GetVertices(ints[i++])->GetNPoints();
01922             // vertices
01923             for (j = 0; j < k; j++)
01924                out << ' ' << vc + ints[i++];
01925             out << '\n';
01926          }
01927          vc += nv;
01928       }
01929    }
01930    else
01931    {
01932       // for each boundary and shared boundary element send its
01933       // geometry type and its vertices
01934       ne = 0;
01935       for (i = 0; i < NumOfBdrElements; i++)
01936          ne += 1 + boundary[i]->GetNVertices();
01937       for (i = 0; i < shared_boundary.Size(); i++)
01938          ne += 1 + shared_boundary[i]->GetNVertices();
01939       nv = NumOfVertices;
01940       MPI_Send(nv_ne, 2, MPI_INT, 0, 446, MyComm);
01941       ints.SetSize(ne);
01942       // boundary
01943       for (i = j = 0; i < NumOfBdrElements; i++)
01944       {
01945          ints[j++] = boundary[i]->GetGeometryType();
01946          nv = boundary[i]->GetNVertices();
01947          v  = boundary[i]->GetVertices();
01948          for (k = 0; k < nv; k++)
01949             ints[j++] = v[k];
01950       }
01951       // shared boundary
01952       for (i = 0; i < shared_boundary.Size(); i++)
01953       {
01954          ints[j++] = shared_boundary[i]->GetGeometryType();
01955          nv = shared_boundary[i]->GetNVertices();
01956          v  = shared_boundary[i]->GetVertices();
01957          for (k = 0; k < nv; k++)
01958             ints[j++] = v[k];
01959       }
01960       MPI_Send(&ints[0], ne, MPI_INT, 0, 447, MyComm);
01961    }
01962 
01963    // vertices / nodes
01964    MPI_Reduce(&NumOfVertices, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
01965    if (MyRank == 0)
01966       out << "\nvertices\n" << nv << '\n';
01967    if (Nodes == NULL)
01968    {
01969       if (MyRank == 0)
01970       {
01971          out << Dim << '\n';
01972          for (i = 0; i < NumOfVertices; i++)
01973          {
01974             out << vertices[i](0);
01975             for (j = 1; j < Dim; j++)
01976                out << ' ' << vertices[i](j);
01977             out << '\n';
01978          }
01979          for (p = 1; p < NRanks; p++)
01980          {
01981             MPI_Recv(&nv, 1, MPI_INT, p, 448, MyComm, &status);
01982             vert.SetSize(nv*Dim);
01983             MPI_Recv(&vert[0], nv*Dim, MPI_DOUBLE, p, 449, MyComm, &status);
01984             for (i = 0; i < nv; i++)
01985             {
01986                out << vert[i*Dim];
01987                for (j = 1; j < Dim; j++)
01988                   out << ' ' << vert[i*Dim+j];
01989                out << '\n';
01990             }
01991          }
01992       }
01993       else
01994       {
01995          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 448, MyComm);
01996          vert.SetSize(NumOfVertices*Dim);
01997          for (i = 0; i < NumOfVertices; i++)
01998             for (j = 0; j < Dim; j++)
01999                vert[i*Dim+j] = vertices[i](j);
02000          MPI_Send(&vert[0], NumOfVertices*Dim, MPI_DOUBLE, 0, 449, MyComm);
02001       }
02002    }
02003    else
02004    {
02005       if (MyRank == 0)
02006          out << "\nnodes\n";
02007       ParGridFunction *pnodes = dynamic_cast<ParGridFunction *>(Nodes);
02008       if (pnodes)
02009       {
02010          pnodes->SaveAsOne(out);
02011       }
02012       else
02013       {
02014          ParFiniteElementSpace *pfes =
02015             dynamic_cast<ParFiniteElementSpace *>(Nodes->FESpace());
02016          if (pfes)
02017          {
02018             // create a wrapper ParGridFunction
02019             ParGridFunction ParNodes(pfes, Nodes);
02020             ParNodes.SaveAsOne(out);
02021          }
02022          else
02023             mfem_error("ParMesh::PrintAsOne : Nodes have no parallel info!");
02024       }
02025    }
02026 }
02027 
02028 void ParMesh::PrintAsOneXG(ostream &out)
02029 {
02030    if (Dim == 3 && meshgen == 1)
02031    {
02032       int i, j, k, nv, ne, p;
02033       const int *ind, *v;
02034       MPI_Status status;
02035       Array<double> vert;
02036       Array<int> ints;
02037 
02038       if (MyRank == 0)
02039       {
02040          out << "NETGEN_Neutral_Format\n";
02041          // print the vertices
02042          ne = NumOfVertices;
02043          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02044          out << nv << '\n';
02045          for (i = 0; i < NumOfVertices; i++)
02046          {
02047             for (j = 0; j < Dim; j++)
02048                out << " " << vertices[i](j);
02049             out << '\n';
02050          }
02051          for (p = 1; p < NRanks; p++)
02052          {
02053             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02054             vert.SetSize(Dim*nv);
02055             MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
02056             for (i = 0; i < nv; i++)
02057             {
02058                for (j = 0; j < Dim; j++)
02059                   out << " " << vert[Dim*i+j];
02060                out << '\n';
02061             }
02062          }
02063 
02064          // print the elements
02065          nv = NumOfElements;
02066          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02067          out << ne << '\n';
02068          for (i = 0; i < NumOfElements; i++)
02069          {
02070             nv = elements[i]->GetNVertices();
02071             ind = elements[i]->GetVertices();
02072             out << 1;
02073             for (j = 0; j < nv; j++)
02074                out << " " << ind[j]+1;
02075             out << '\n';
02076          }
02077          k = NumOfVertices;
02078          for (p = 1; p < NRanks; p++)
02079          {
02080             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02081             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02082             ints.SetSize(4*ne);
02083             MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
02084             for (i = 0; i < ne; i++)
02085             {
02086                out << p+1;
02087                for (j = 0; j < 4; j++)
02088                   out << " " << k+ints[i*4+j]+1;
02089                out << '\n';
02090             }
02091             k += nv;
02092          }
02093          // print the boundary + shared faces information
02094          nv = NumOfBdrElements + shared_faces.Size();
02095          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02096          out << ne << '\n';
02097          // boundary
02098          for (i = 0; i < NumOfBdrElements; i++)
02099          {
02100             nv = boundary[i]->GetNVertices();
02101             ind = boundary[i]->GetVertices();
02102             out << 1;
02103             for (j = 0; j < nv; j++)
02104                out << " " << ind[j]+1;
02105             out << '\n';
02106          }
02107          // shared faces
02108          for (i = 0; i < shared_faces.Size(); i++)
02109          {
02110             nv = shared_faces[i]->GetNVertices();
02111             ind = shared_faces[i]->GetVertices();
02112             out << 1;
02113             for (j = 0; j < nv; j++)
02114                out << " " << ind[j]+1;
02115             out << '\n';
02116          }
02117          k = NumOfVertices;
02118          for (p = 1; p < NRanks; p++)
02119          {
02120             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02121             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02122             ints.SetSize(3*ne);
02123             MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
02124             for (i = 0; i < ne; i++)
02125             {
02126                out << p+1;
02127                for (j = 0; j < 3; j++)
02128                   out << " " << k+ints[i*3+j]+1;
02129                out << '\n';
02130             }
02131             k += nv;
02132          }
02133       }
02134       else
02135       {
02136          ne = NumOfVertices;
02137          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02138          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02139          vert.SetSize(Dim*NumOfVertices);
02140          for (i = 0; i < NumOfVertices; i++)
02141             for (j = 0; j < Dim; j++)
02142                vert[Dim*i+j] = vertices[i](j);
02143          MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
02144                   0, 445, MyComm);
02145          // elements
02146          ne = NumOfElements;
02147          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02148          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02149          MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
02150          ints.SetSize(NumOfElements*4);
02151          for (i = 0; i < NumOfElements; i++)
02152          {
02153             v = elements[i]->GetVertices();
02154             for (j = 0; j < 4; j++)
02155                ints[4*i+j] = v[j];
02156          }
02157          MPI_Send(&ints[0], 4*NumOfElements, MPI_INT, 0, 447, MyComm);
02158          // boundary + shared faces
02159          nv = NumOfBdrElements + shared_faces.Size();
02160          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02161          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02162          ne = NumOfBdrElements + shared_faces.Size();
02163          MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
02164          ints.SetSize(3*ne);
02165          for (i = 0; i < NumOfBdrElements; i++)
02166          {
02167             v = boundary[i]->GetVertices();
02168             for (j = 0; j < 3; j++)
02169                ints[3*i+j] = v[j];
02170          }
02171          for ( ; i < ne; i++)
02172          {
02173             v = shared_faces[i-NumOfBdrElements]->GetVertices();
02174             for (j = 0; j < 3; j++)
02175                ints[3*i+j] = v[j];
02176          }
02177          MPI_Send(&ints[0], 3*ne, MPI_INT, 0, 447, MyComm);
02178       }
02179    }
02180 
02181    if (Dim == 3 && meshgen == 2)
02182    {
02183       int i, j, k, nv, ne, p;
02184       const int *ind, *v;
02185       MPI_Status status;
02186       Array<double> vert;
02187       Array<int> ints;
02188 
02189       int TG_nv, TG_ne, TG_nbe;
02190 
02191       if (MyRank == 0)
02192       {
02193          MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02194          MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02195          nv = NumOfBdrElements + shared_faces.Size();
02196          MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
02197 
02198          out << "TrueGrid\n"
02199              << "1 " << TG_nv << " " << TG_ne << " 0 0 0 0 0 0 0\n"
02200              << "0 0 0 1 0 0 0 0 0 0 0\n"
02201              << "0 0 " << TG_nbe << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
02202              << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
02203              << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
02204 
02205          // print the vertices
02206          nv = TG_nv;
02207          for (i = 0; i < NumOfVertices; i++)
02208             out << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
02209                 << " " << vertices[i](2) << " 0.0\n";
02210          for (p = 1; p < NRanks; p++)
02211          {
02212             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02213             vert.SetSize(Dim*nv);
02214             MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
02215             for (i = 0; i < nv; i++)
02216                out << i+1 << " 0.0 " << vert[Dim*i] << " " << vert[Dim*i+1]
02217                    << " " << vert[Dim*i+2] << " 0.0\n";
02218          }
02219 
02220          // print the elements
02221          ne = TG_ne;
02222          for (i = 0; i < NumOfElements; i++)
02223          {
02224             nv = elements[i]->GetNVertices();
02225             ind = elements[i]->GetVertices();
02226             out << i+1 << " " << 1;
02227             for (j = 0; j < nv; j++)
02228                out << " " << ind[j]+1;
02229             out << '\n';
02230          }
02231          k = NumOfVertices;
02232          for (p = 1; p < NRanks; p++)
02233          {
02234             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02235             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02236             ints.SetSize(8*ne);
02237             MPI_Recv(&ints[0], 8*ne, MPI_INT, p, 447, MyComm, &status);
02238             for (i = 0; i < ne; i++)
02239             {
02240                out << i+1 << " " << p+1;
02241                for (j = 0; j < 8; j++)
02242                   out << " " << k+ints[i*8+j]+1;
02243                out << '\n';
02244             }
02245             k += nv;
02246          }
02247 
02248          // print the boundary + shared faces information
02249          ne = TG_nbe;
02250          // boundary
02251          for (i = 0; i < NumOfBdrElements; i++)
02252          {
02253             nv = boundary[i]->GetNVertices();
02254             ind = boundary[i]->GetVertices();
02255             out << 1;
02256             for (j = 0; j < nv; j++)
02257                out << " " << ind[j]+1;
02258             out << " 1.0 1.0 1.0 1.0\n";
02259          }
02260          // shared faces
02261          for (i = 0; i < shared_faces.Size(); i++)
02262          {
02263             nv = shared_faces[i]->GetNVertices();
02264             ind = shared_faces[i]->GetVertices();
02265             out << 1;
02266             for (j = 0; j < nv; j++)
02267                out << " " << ind[j]+1;
02268             out << " 1.0 1.0 1.0 1.0\n";
02269          }
02270          k = NumOfVertices;
02271          for (p = 1; p < NRanks; p++)
02272          {
02273             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02274             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02275             ints.SetSize(4*ne);
02276             MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447, MyComm, &status);
02277             for (i = 0; i < ne; i++)
02278             {
02279                out << p+1;
02280                for (j = 0; j < 4; j++)
02281                   out << " " << k+ints[i*4+j]+1;
02282                out << " 1.0 1.0 1.0 1.0\n";
02283             }
02284             k += nv;
02285          }
02286       }
02287       else
02288       {
02289          MPI_Reduce(&NumOfVertices, &TG_nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02290          MPI_Reduce(&NumOfElements, &TG_ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02291          nv = NumOfBdrElements + shared_faces.Size();
02292          MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0, MyComm);
02293 
02294          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02295          vert.SetSize(Dim*NumOfVertices);
02296          for (i = 0; i < NumOfVertices; i++)
02297             for (j = 0; j < Dim; j++)
02298                vert[Dim*i+j] = vertices[i](j);
02299          MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE, 0, 445, MyComm);
02300          // elements
02301          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02302          MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
02303          ints.SetSize(NumOfElements*8);
02304          for (i = 0; i < NumOfElements; i++)
02305          {
02306             v = elements[i]->GetVertices();
02307             for (j = 0; j < 8; j++)
02308                ints[8*i+j] = v[j];
02309          }
02310          MPI_Send(&ints[0], 8*NumOfElements, MPI_INT, 0, 447, MyComm);
02311          // boundary + shared faces
02312          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02313          ne = NumOfBdrElements + shared_faces.Size();
02314          MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
02315          ints.SetSize(4*ne);
02316          for (i = 0; i < NumOfBdrElements; i++)
02317          {
02318             v = boundary[i]->GetVertices();
02319             for (j = 0; j < 4; j++)
02320                ints[4*i+j] = v[j];
02321          }
02322          for ( ; i < ne; i++)
02323          {
02324             v = shared_faces[i-NumOfBdrElements]->GetVertices();
02325             for (j = 0; j < 4; j++)
02326                ints[4*i+j] = v[j];
02327          }
02328          MPI_Send(&ints[0], 4*ne, MPI_INT, 0, 447, MyComm);
02329       }
02330    }
02331 
02332    if (Dim == 2)
02333    {
02334       int i, j, k, attr, nv, ne, p;
02335       Array<int> v;
02336       MPI_Status status;
02337       Array<double> vert;
02338       Array<int> ints;
02339 
02340 
02341       if (MyRank == 0)
02342       {
02343          out << "areamesh2\n\n";
02344 
02345          // print the boundary + shared edges information
02346          nv = NumOfBdrElements + shared_edges.Size();
02347          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02348          out << ne << '\n';
02349          // boundary
02350          for (i = 0; i < NumOfBdrElements; i++)
02351          {
02352             attr = boundary[i]->GetAttribute();
02353             boundary[i]->GetVertices(v);
02354             out << attr << "     ";
02355             for (j = 0; j < v.Size(); j++)
02356                out << v[j] + 1 << "   ";
02357             out << '\n';
02358          }
02359          // shared edges
02360          for (i = 0; i < shared_edges.Size(); i++)
02361          {
02362             attr = shared_edges[i]->GetAttribute();
02363             shared_edges[i]->GetVertices(v);
02364             out << attr << "     ";
02365             for (j = 0; j < v.Size(); j++)
02366                out << v[j] + 1 << "   ";
02367             out << '\n';
02368          }
02369          k = NumOfVertices;
02370          for (p = 1; p < NRanks; p++)
02371          {
02372             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02373             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02374             ints.SetSize(2*ne);
02375             MPI_Recv(&ints[0], 2*ne, MPI_INT, p, 447, MyComm, &status);
02376             for (i = 0; i < ne; i++)
02377             {
02378                out << p+1;
02379                for (j = 0; j < 2; j++)
02380                   out << " " << k+ints[i*2+j]+1;
02381                out << '\n';
02382             }
02383             k += nv;
02384          }
02385 
02386          // print the elements
02387          nv = NumOfElements;
02388          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02389          out << ne << '\n';
02390          for (i = 0; i < NumOfElements; i++)
02391          {
02392             attr = elements[i]->GetAttribute();
02393             elements[i]->GetVertices(v);
02394             out << 1 << "   " << 3 << "   ";
02395             for (j = 0; j < v.Size(); j++)
02396                out << v[j] + 1 << "  ";
02397             out << '\n';
02398          }
02399          k = NumOfVertices;
02400          for (p = 1; p < NRanks; p++)
02401          {
02402             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02403             MPI_Recv(&ne, 1, MPI_INT, p, 446, MyComm, &status);
02404             ints.SetSize(3*ne);
02405             MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447, MyComm, &status);
02406             for (i = 0; i < ne; i++)
02407             {
02408                out << p+1 << " " << 3;
02409                for (j = 0; j < 3; j++)
02410                   out << " " << k+ints[i*3+j]+1;
02411                out << '\n';
02412             }
02413             k += nv;
02414          }
02415 
02416          // print the vertices
02417          ne = NumOfVertices;
02418          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02419          out << nv << '\n';
02420          for (i = 0; i < NumOfVertices; i++)
02421          {
02422             for (j = 0; j < Dim; j++)
02423                out << vertices[i](j) << " ";
02424             out << '\n';
02425          }
02426          for (p = 1; p < NRanks; p++)
02427          {
02428             MPI_Recv(&nv, 1, MPI_INT, p, 444, MyComm, &status);
02429             vert.SetSize(Dim*nv);
02430             MPI_Recv(&vert[0], Dim*nv, MPI_DOUBLE, p, 445, MyComm, &status);
02431             for (i = 0; i < nv; i++)
02432             {
02433                for (j = 0; j < Dim; j++)
02434                   out << " " << vert[Dim*i+j];
02435                out << '\n';
02436             }
02437          }
02438       }
02439       else
02440       {
02441          // boundary + shared faces
02442          nv = NumOfBdrElements + shared_edges.Size();
02443          MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
02444          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02445          ne = NumOfBdrElements + shared_edges.Size();
02446          MPI_Send(&ne, 1, MPI_INT, 0, 446, MyComm);
02447          ints.SetSize(2*ne);
02448          for (i = 0; i < NumOfBdrElements; i++)
02449          {
02450             boundary[i]->GetVertices(v);
02451             for (j = 0; j < 2; j++)
02452                ints[2*i+j] = v[j];
02453          }
02454          for ( ; i < ne; i++)
02455          {
02456             shared_edges[i-NumOfBdrElements]->GetVertices(v);
02457             for (j = 0; j < 2; j++)
02458                ints[2*i+j] = v[j];
02459          }
02460          MPI_Send(&ints[0], 2*ne, MPI_INT, 0, 447, MyComm);
02461          // elements
02462          ne = NumOfElements;
02463          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02464          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02465          MPI_Send(&NumOfElements, 1, MPI_INT, 0, 446, MyComm);
02466          ints.SetSize(NumOfElements*3);
02467          for (i = 0; i < NumOfElements; i++)
02468          {
02469             elements[i]->GetVertices(v);
02470             for (j = 0; j < 3; j++)
02471                ints[3*i+j] = v[j];
02472          }
02473          MPI_Send(&ints[0], 3*NumOfElements, MPI_INT, 0, 447, MyComm);
02474          // vertices
02475          ne = NumOfVertices;
02476          MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
02477          MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444, MyComm);
02478          vert.SetSize(Dim*NumOfVertices);
02479          for (i = 0; i < NumOfVertices; i++)
02480             for (j = 0; j < Dim; j++)
02481                vert[Dim*i+j] = vertices[i](j);
02482          MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
02483                   0, 445, MyComm);
02484       }
02485    }
02486 }
02487 
02488 void ParMesh::PrintInfo(ostream &out)
02489 {
02490    int i;
02491    DenseMatrix J(Dim);
02492    double h_min, h_max, kappa_min, kappa_max, h, kappa;
02493 
02494    if (MyRank == 0)
02495       out << "Parallel Mesh Stats:" << endl;
02496 
02497    for (i = 0; i < NumOfElements; i++)
02498    {
02499       GetElementJacobian(i, J);
02500       h = pow(fabs(J.Det()), 1.0/double(Dim));
02501       kappa = J.CalcSingularvalue(0) / J.CalcSingularvalue(Dim-1);
02502       if (i == 0)
02503       {
02504          h_min = h_max = h;
02505          kappa_min = kappa_max = kappa;
02506       }
02507       else
02508       {
02509          if (h < h_min)  h_min = h;
02510          if (h > h_max)  h_max = h;
02511          if (kappa < kappa_min)  kappa_min = kappa;
02512          if (kappa > kappa_max)  kappa_max = kappa;
02513       }
02514    }
02515 
02516    double gh_min, gh_max, gk_min, gk_max;
02517    MPI_Reduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
02518    MPI_Reduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
02519    MPI_Reduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN, 0, MyComm);
02520    MPI_Reduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX, 0, MyComm);
02521 
02522    int ldata[5]; // vert, edge, face, elem, neighbors;
02523    int mindata[5], maxdata[5], sumdata[5];
02524 
02525    // count locally owned vertices, edges, and faces
02526    ldata[0] = GetNV();
02527    ldata[1] = GetNEdges();
02528    ldata[2] = GetNFaces();
02529    ldata[3] = GetNE();
02530    ldata[4] = gtopo.GetNumNeighbors()-1;
02531    for (int gr = 1; gr < GetNGroups(); gr++)
02532       if (!gtopo.IAmMaster(gr)) // we are not the master
02533       {
02534          ldata[0] -= group_svert.RowSize(gr-1);
02535          ldata[1] -= group_sedge.RowSize(gr-1);
02536          ldata[2] -= group_sface.RowSize(gr-1);
02537       }
02538 
02539    MPI_Reduce(ldata, mindata, 5, MPI_INT, MPI_MIN, 0, MyComm);
02540    MPI_Reduce(ldata, sumdata, 5, MPI_INT, MPI_SUM, 0, MyComm); // overflow?
02541    MPI_Reduce(ldata, maxdata, 5, MPI_INT, MPI_MAX, 0, MyComm);
02542 
02543    if (MyRank == 0)
02544    {
02545       out << '\n'
02546           << "           "
02547           << setw(12) << "minimum"
02548           << setw(12) << "average"
02549           << setw(12) << "maximum"
02550           << setw(12) << "total" << '\n';
02551       out << " vertices  "
02552           << setw(12) << mindata[0]
02553           << setw(12) << sumdata[0]/NRanks
02554           << setw(12) << maxdata[0]
02555           << setw(12) << sumdata[0] << '\n';
02556       out << " edges     "
02557           << setw(12) << mindata[1]
02558           << setw(12) << sumdata[1]/NRanks
02559           << setw(12) << maxdata[1]
02560           << setw(12) << sumdata[1] << '\n';
02561       if (Dim == 3)
02562          out << " faces     "
02563              << setw(12) << mindata[2]
02564              << setw(12) << sumdata[2]/NRanks
02565              << setw(12) << maxdata[2]
02566              << setw(12) << sumdata[2] << '\n';
02567       out << " elements  "
02568           << setw(12) << mindata[3]
02569           << setw(12) << sumdata[3]/NRanks
02570           << setw(12) << maxdata[3]
02571           << setw(12) << sumdata[3] << '\n';
02572       out << " neighbors "
02573           << setw(12) << mindata[4]
02574           << setw(12) << sumdata[4]/NRanks
02575           << setw(12) << maxdata[4] << '\n';
02576       out << '\n'
02577           << "       "
02578           << setw(12) << "minimum"
02579           << setw(12) << "maximum" << '\n';
02580       out << " h     "
02581           << setw(12) << gh_min
02582           << setw(12) << gh_max << '\n';
02583       out << " kappa "
02584           << setw(12) << gk_min
02585           << setw(12) << gk_max << '\n';
02586    }
02587 }
02588 
02589 ParMesh::~ParMesh()
02590 {
02591    int i;
02592 
02593    for (i = 0; i < shared_faces.Size(); i++)
02594       FreeElement(shared_faces[i]);
02595    for (i = 0; i < shared_edges.Size(); i++)
02596       FreeElement(shared_edges[i]);
02597 
02598    // The Mesh destructor is called automatically
02599 }
02600 
02601 #endif
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines