How to deal with AbstractGrid Classes

File galeri/example/Grid.cpp:

// @HEADER
// ************************************************************************
//
//           Galeri: Finite Element and Matrix Generation Package
//                 Copyright (2006) ETHZ/Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
//
// Questions about Galeri? Contact Marzio Sala (marzio.sala _AT_ gmail.com)
//
// ************************************************************************
// @HEADER

#include "Galeri_ConfigDefs.h"
#include "Galeri_Exception.h"
#include "Galeri_FiniteElements.h"
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif

using namespace Galeri;
using namespace Galeri::FiniteElements;

// =========== //
// main driver //
// =========== //

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  try {
    
    // Prepares the computational domain. For simplicity,           
    // the computation domain has (nx * NumProcs, ny, nz) elements, 
    // and it is partitioned into (NumProcs, 1, 1) subdomains. 
    
    int nx = 2 * Comm.NumProc();
    int ny = 2; 
    int nz = 2;
    int mx = Comm.NumProc(), my = 1, mz = 1;
    double lx = 1.0 * Comm.NumProc();
    double ly = 1.0;
    double lz = 1.0;

    // Uncomment one of the following:
    TriangleRectangleGrid Grid(Comm, nx, ny, mx, my, lx, ly);
    //QuadRectangleGrid Grid(Comm, nx, ny, mx, my, lx, ly);
    //TetCubeGrid Grid(Comm, nx, ny, nz, mx, my, mz, lx, ly, lz);
    //HexCubeGrid Grid(Comm, nx, ny, nz, mx, my, mz, lx, ly, lz);

    // just work on processor 0, only to reduce the output
    if (!Comm.MyPID()) 
    {
      // Extracts the information from the Grid using AbstractGrid 
      // methods. First, some general information. 

      cout << "Number of dimensions = " << Grid.NumDimensions() << endl;
      cout << "Number of vertices per element = " << Grid.NumVerticesPerElement() << endl;
      cout << "Number of faces per element = " << Grid.NumFacesPerElement() << endl;
      cout << "Number of vertices per face = " << Grid.NumVerticesPerFace() << endl;
      cout << "Element type = " << Grid.ElementType() << endl;

      cout << "Number of elements: global = " << Grid.NumGlobalElements();
      cout << ", on proc 0 = " << Grid.NumMyElements() << endl;
      cout << "Number of vertices: global = " << Grid.NumGlobalVertices();
      cout << ", on proc 0 = " << Grid.NumMyVertices() << endl;
      cout << "Number of boundary faces: "
        "on proc 0 = " << Grid.NumMyBoundaryFaces() << endl;

      // Now loop over all the local elements, and print  
      // the local ID of each vertex.

      vector<double> coord(3);
      vector<int> vertices(Grid.NumVerticesPerElement());

      for (int i = 0 ; i < Grid.NumMyElements() ; ++i) 
      {
        cout << "Element " << i << ": ";
        Grid.ElementVertices(i, &vertices[0]);
        for (int j = 0 ; j < Grid.NumVerticesPerElement() ; ++j)
          cout << vertices[j] << " ";
        cout << endl;
      }

      // Loop over all local vertices, and print out the coordinates

      for (int i = 0 ; i < Grid.NumMyVertices() ; ++i)
      {
        Grid.VertexCoord(i, &coord[0]);
        cout << "Vertex " << i << ": (" << coord[0];
        cout << ", " << coord[1] << ", " << coord[2] << ")" << endl;
      }

      // Loop over all local boundary faces, and print out the local ID of each
      // vertex belonging to the face

      for (int i = 0 ; i < Grid.NumMyBoundaryFaces() ; ++i)
      {
        cout << "Boundary face " << i << ": ";
        int tag;
        Grid.FaceVertices(i, tag, &vertices[0]);
        for (int j = 0 ; j < Grid.NumVerticesPerFace() ; ++j)
          cout << vertices[j] << " ";
        cout << endl;
      }
    }
    
    // The following instructions can be used to visualize with MEDIT
#if 0
    MEDITInterface MEDIT(Comm);

    // We need to define a vector to plot, in this case constant
    Epetra_Vector Vector(Grid.RowMap());
    MEDIT.Write(Grid, "grid", Vector);
#endif

  }
  catch (Exception& rhs)
  {
    if (Comm.MyPID() == 0)
    {
      cerr << "Caught exception: ";
      rhs.Print();
    }
  }

#ifdef HAVE_MPI
  MPI_Finalize();
#endif
}