Branch to develop further the element reordering capabilities in MFEM

fem/fespace.cpp

@@ -924,6 +924,35 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
    return R;
 }
 
+PermutationOperator *FiniteElementSpace::GetElementReorderingOperator(
+   const Table *old_elem_dof)
+{
+   PermutationOperator *pT = new PermutationOperator(GetVSize());
+   Array<int> perm;
+   perm.MakeRef(pT->GetPermutation());
+   const int *el_perm = mesh->GetElementPermutation();
+   for (int old_el = 0; old_el < GetNE(); ++old_el)
+   {
+      const int new_el = el_perm[old_el];
+      const int nd = old_elem_dof->RowSize(old_el);
+      const int *old_dofs = old_elem_dof->GetRow(old_el);
+      const int *new_dofs = elem_dof->GetRow(new_el);
+      MFEM_ASSERT(nd == elem_dof->RowSize(new_el), "internal error");
+      for (int i = 0; i < nd; ++i)
+      {
+         for (int vd = 0; vd < vdim; ++vd)
+         {
+            const int ovd = DofToVDof(old_dofs[i],vd);
+            const int nvd = DofToVDof(new_dofs[i],vd);
+            (ovd >= 0) ? (perm[ovd] = nvd) : (perm[-1-ovd] = -1-nvd);
+         }
+      }
+   }
+   MFEM_DEBUG_DO(int err = pT->CheckPermutation());
+   MFEM_ASSERT(!err, "internal error: " << err);
+   return pT;
+}
+
 FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
                                        const FiniteElementCollection *fec,
                                        int vdim, int ordering)
@@ -1544,8 +1573,19 @@ void FiniteElementSpace::Update(bool want_transform)
             break;
          }
 
-         default:
+         case Mesh::NONE:
+         case Mesh::REBALANCE:
             break; // T stays NULL
+
+         case Mesh::REORDER:
+         {
+            T = GetElementReorderingOperator(old_elem_dof);
+            break;
+         }
+
+         default:
+            MFEM_ABORT("Mesh::Operation not supported!");
+            break;
       }
    }
 

fem/fespace.hpp

@@ -127,6 +127,7 @@ protected:
    /// Calculate GridFunction restriction matrix after mesh derefinement.
    SparseMatrix* DerefinementMatrix(int old_ndofs, const Table* old_elem_dof);
 
+   PermutationOperator *GetElementReorderingOperator(const Table *old_elem_dof);
 
 public:
    FiniteElementSpace(Mesh *mesh, const FiniteElementCollection *fec,

fem/gridfunc.cpp

@@ -160,7 +160,7 @@ void GridFunction::Update()
    {
       Vector tmp(T->Height());
       T->Mult(*this, tmp);
-      *this = tmp;
+      this->Swap(tmp);
    }
    else
    {

fem/pfespace.cpp

@@ -2287,8 +2287,18 @@ void ParFiniteElementSpace::Update(bool want_transform)
             break;
          }
 
-         default:
+         case Mesh::REORDER:
+         {
+            T = GetElementReorderingOperator(old_elem_dof);
+            break;
+         }
+
+         case Mesh::NONE:
             break; // T stays NULL
+
+         default:
+            MFEM_ABORT("Mesh::Operation not supported!");
+            break;
       }
       delete old_elem_dof;
    }

fem/pgridfunc.cpp

@@ -263,8 +263,9 @@ const
       fes->GetElementDofs(i, dofs);
       fes->DofsToVDofs(vdim-1, dofs);
       DofVal.SetSize(dofs.Size());
-      const FiniteElement *fe = fes->GetFE(i);
-      MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
+      MFEM_DEBUG_DO(const FiniteElement *fe = fes->GetFE(i));
+      MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
+                  "invalid FE map type");
       fes->GetFE(i)->CalcShape(ip, DofVal);
       GetSubVector(dofs, LocVec);
    }

general/communication.cpp

@@ -33,6 +33,26 @@ using namespace std;
 namespace mfem
 {
 
+void ParTimer::GetParStats()
+{
+   double my_rt = RealTime();
+   MPI_Reduce(&my_rt, &min_rt, 1, MPI_DOUBLE, MPI_MIN, 0, comm);
+   MPI_Reduce(&my_rt, &avg_rt, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
+   MPI_Reduce(&my_rt, &max_rt, 1, MPI_DOUBLE, MPI_MAX, 0, comm);
+   avg_rt /= comm_size;
+}
+
+std::ostream &operator<<(std::ostream &out, ParTimer &pt)
+{
+   if (pt.CommRank() == 0)
+   {
+      out << "max: " << pt.RealTimeMax()
+          << "s, avg: " << pt.RealTimeAvg()
+          << "s, min: " << pt.RealTimeMin() << "s";
+   }
+   return out;
+}
+
 GroupTopology::GroupTopology(const GroupTopology &gt)
    : MyComm(gt.MyComm),
      group_lproc(gt.group_lproc)

general/communication.hpp

@@ -16,6 +16,7 @@
 
 #ifdef MFEM_USE_MPI
 
+#include "tic_toc.hpp"
 #include "array.hpp"
 #include "table.hpp"
 #include "sets.hpp"
@@ -49,6 +50,42 @@ public:
    bool Root() const { return world_rank == 0; }
 };
 
+
+/** @brief This class extends class StopWatch with an MPI communicator,
+    providing parallel timing statistics. */
+/** The statistics are: the minimal, average, and maximal real times from all
+    instances of the StopWatch, across the MPI communicator. The statistics are
+    accessible from rank 0 only. */
+class ParTimer : public StopWatch
+{
+protected:
+   MPI_Comm comm;
+   int comm_rank, comm_size;
+   double max_rt, avg_rt, min_rt;
+
+public:
+   ParTimer(MPI_Comm comm_) : comm(comm_)
+   {
+      MPI_Comm_rank(comm, &comm_rank);
+      MPI_Comm_size(comm, &comm_size);
+      max_rt = avg_rt = min_rt = 0.0;
+   }
+   int CommRank() const { return comm_rank; }
+   int CommSize() const { return comm_size; }
+   /** @brief This method computes the parallel statistics. It must be called on
+       all ranks. */
+   void GetParStats();
+   /// Shortcut for Stop() plus GetParStats().
+   void ParStop() { Stop(); GetParStats(); }
+   double RealTimeMax() const { return max_rt; }
+   double RealTimeAvg() const { return avg_rt; }
+   double RealTimeMin() const { return min_rt; }
+};
+
+/// Overload operator<< for std::ostream and ParTimer.
+std::ostream &operator<<(std::ostream &out, ParTimer &pt);
+
+
 class GroupTopology
 {
 private:

general/sort_pairs.hpp

@@ -25,6 +25,7 @@ class Pair
 public:
    A one;
    B two;
+   Pair(A a, B b) : one(a), two(b) { }
 };
 
 /// @brief Comparison operator for class Pair, based on the first element only.

general/table.cpp

@@ -11,13 +11,22 @@
 
 // Implementation of data types Table.
 
-#include <iostream>
-#include <iomanip>
-
 #include "array.hpp"
 #include "table.hpp"
+#include "sort_pairs.hpp"
 #include "error.hpp"
 
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <algorithm>
+
+// Include the METIS header, if using version 5. If using METIS 4, the needed
+// declarations are inlined below, i.e. no header is needed.
+#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
+#include "metis.h"
+#endif
+
 namespace mfem
 {
 
@@ -292,6 +301,250 @@ int Table::Width() const
    return width + 1;
 }
 
+void Table::GetCMReordering(Array<int> &ordering, bool reverse) const
+{
+   if (size <= 0)
+   {
+      ordering.SetSize(0);
+      return;
+   }
+
+   int num_el = size, stack_p, stack_top_p;
+   ordering.SetSize(num_el);
+   Array<Pair<int,int> > el_stack(num_el);
+   Array<int> el_layer;
+   el_layer.MakeRef(ordering);
+
+   // Assuming that either all diagonal entries are present or none are present.
+
+   // Choose starting element (for one connected component only).
+   int el0 = 0, min_nbrs = RowSize(el0);
+   for (int el = 1; el < num_el; el++)
+   {
+      const int num_nbrs = RowSize(el);
+      if (num_nbrs < min_nbrs)
+      {
+         el0 = el;
+         min_nbrs = num_nbrs;
+      }
+   }
+
+   el_layer = -1;
+   stack_p = stack_top_p = 0;
+   for (int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
+   {
+      if (el_layer[el] != -1) { continue; }
+
+      // FIXME: choose starting element for this connected component.
+
+      el_layer[el] = 0;
+      el_stack[stack_top_p++] = Pair<int,int>(RowSize(el),el);
+
+      int layer = 0, layer_start = stack_p;
+      for ( ; stack_p < stack_top_p; stack_p++)
+      {
+         const int i = el_stack[stack_p].two;
+         for (int j = I[i]; j < I[i+1]; j++)
+         {
+            int k = J[j];
+            if (el_layer[k] == -1)
+            {
+               el_layer[k] = el_layer[i] + 1;
+               el_stack[stack_top_p++] = Pair<int,int>(RowSize(k),k);
+            }
+         }
+         if (stack_p+1 == stack_top_p ||
+             layer < el_layer[el_stack[stack_p+1].two])
+         {
+            std::sort(&el_stack[layer_start], &el_stack[stack_p] + 1);
+            layer++;
+            layer_start = stack_p+1;
+         }
+      }
+   }
+
+   if (!reverse)
+   {
+      for (int i = 0; i < num_el; i++)
+      {
+         ordering[el_stack[i].two] = i;
+      }
+   }
+   else
+   {
+      for (int i = 0; i < num_el; i++)
+      {
+         ordering[el_stack[num_el-1-i].two] = i;
+      }
+   }
+}
+
+#ifdef MFEM_USE_GECKO
+void Table::GetGeckoReordering(const GeckoParameters &g_params,
+                               Array<int> &ordering) const
+{
+   Gecko::Graph graph;
+
+   // Run through all the elements and insert the nodes in the graph for them
+   for (int elemid = 0; elemid < size; ++elemid)
+   {
+      graph.insert();
+   }
+
+   // Run through all the elems and insert arcs to the graph for each element
+   // face Indices in Gecko are 1 based hence the +1 on the insertion
+   for (int elemid = 0; elemid < size; ++elemid)
+   {
+      const int num_neigh = RowSize(elemid);
+      const int *neighid = GetRow(elemid);
+      for (int i = 0; i < num_neigh; ++i)
+      {
+         if (elemid != neighid[i])
+         {
+            graph.insert(elemid + 1,  neighid[i] + 1);
+         }
+      }
+   }
+
+   // Get the reordering from Gecko and copy it into the ordering Array<int>
+   graph.order(g_params.functional,
+               g_params.iterations,
+               g_params.window,
+               g_params.period,
+               g_params.seed);
+   ordering.DeleteAll();
+   ordering.SetSize(size);
+   Gecko::Node::Index NE = size;
+   for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
+   {
+      ordering[gnodeid - 1] = graph.rank(gnodeid);
+   }
+}
+#endif // #ifdef MFEM_USE_GECKO
+
+#ifdef MFEM_USE_METIS
+#ifndef MFEM_USE_METIS_5
+// METIS 4 prototypes
+typedef int idxtype;
+extern "C" {
+   void METIS_EdgeND(int *, idxtype *, idxtype *, int *, int *, idxtype *,
+                     idxtype *);
+   void METIS_NodeND(int *, idxtype *, idxtype *, int *, int *, idxtype *,
+                     idxtype *);
+}
+#endif
+
+void Table::GetMetisReordering(Array<int> &ordering, int type,
+                               bool check_for_diag) const
+{
+   if (size <= 0)
+   {
+      ordering.SetSize(0);
+      return;
+   }
+
+#ifndef MFEM_USE_METIS_5
+   int numflag = 0;
+   int options[8];
+#else
+   int err;
+   int options[METIS_NOPTIONS];
+#endif
+
+#ifndef MFEM_USE_METIS_5
+   options[0] = 0; // use the default options
+#else
+   METIS_SetDefaultOptions(options);
+   // METIS_OPTION_CTYPE, METIS_OPTION_RTYPE, METIS_OPTION_NO2HOP,
+   // METIS_OPTION_NSEPS, METIS_OPTION_NITER, METIS_OPTION_UFACTOR,
+   // METIS_OPTION_COMPRESS, METIS_OPTION_CCORDER, METIS_OPTION_SEED,
+   // METIS_OPTION_PFACTOR, METIS_OPTION_NUMBERING, METIS_OPTION_DBGLVL
+#endif
+   int n = size, *mI = I, *mJ = J;
+
+   // Check if we need to remove any diagonal entries.
+   int num_diag = 0;
+   if (check_for_diag)
+   {
+      for (int row = 0; row < n; row++)
+      {
+         for (int j = mI[row]; j < mI[row+1]; j++)
+         {
+            if (row == mJ[j]) { num_diag++; }
+         }
+      }
+      if (num_diag)
+      {
+         // Remove the diagonal entries.
+         mI = new int[n+1];
+         mJ = new int[I[n]-num_diag];
+         mI[0] = 0;
+         for (int row = 0, counter = 0; row < n; row++)
+         {
+            for (int j = I[row]; j < I[row+1]; j++)
+            {
+               if (row != J[j]) { mJ[counter++] = J[j]; }
+            }
+            mI[row+1] = counter;
+         }
+      }
+   }
+
+   ordering.SetSize(n);
+   Array<int> inv_ordering(n);
+
+   if (type == 0 || type == 1)
+   {
+#ifndef MFEM_USE_METIS_5
+      // Metis 4
+      if (type == 0)
+      {
+         // From the manual: "This function computes fill reducing orderings of
+         // sparse matrices using the multilevel nested dissection algorithm".
+         // We create the reordering based on the element-to-element matrix as
+         // defined by the method ElementToElementTable().
+         METIS_NodeND(&n,
+                      (idxtype *) mI,
+                      (idxtype *) mJ,
+                      &numflag,
+                      options,
+                      (idxtype *) inv_ordering.GetData(),
+                      (idxtype *) ordering.GetData());
+      }
+      else
+      {
+         METIS_EdgeND(&n,
+                      (idxtype *) mI,
+                      (idxtype *) mJ,
+                      &numflag,
+                      options,
+                      (idxtype *) inv_ordering.GetData(),
+                      (idxtype *) ordering.GetData());
+      }
+#else // #ifndef MFEM_USE_METIS_5
+      // Metis 5
+      err = METIS_NodeND((idx_t *) &n,
+                         (idx_t *) mI,
+                         (idx_t *) mJ,
+                         NULL,         // vwgt, NULL - equal weights
+                         options,
+                         (idx_t *) inv_ordering.GetData(),
+                         (idx_t *) ordering.GetData());
+      MFEM_VERIFY(err == METIS_OK, "error in METIS_NodeND");
+#endif // #ifndef MFEM_USE_METIS_5
+   }
+   else
+   {
+      MFEM_ABORT("invalid parameter value: type = " << type);
+   }
+   if (num_diag)
+   {
+      delete [] mJ;
+      delete [] mI;
+   }
+}
+#endif // #ifdef MFEM_USE_METIS
+
 void Table::Print(std::ostream & out, int width) const
 {
    int i, j;
@@ -327,6 +580,100 @@ void Table::PrintMatlab(std::ostream & out) const
    out << flush;
 }
 
+void Table::PrintOrderingStats(std::ostream &out) const
+{
+   out << "Table ordering statistics:\n";
+   if (size <= 0)
+   {
+      out << "   (the Table is empty)\n";
+      return;
+   }
+   const int num_conn = I[size];
+   const int width = Width();
+   out << "   number of rows        = " << size << '\n'
+       << "   number of columns     = " << width << '\n'
+       << "   number of connections = " << num_conn << '\n';
+
+   const int bin_factor = 4, num_bins = 10;
+
+   int max_jump = 0, min_jump = width;
+   long sum_dist = 0, sum_jump = 0;
+
+   int bins[2*num_bins+1];
+   std::fill(bins, bins+2*num_bins+1, 0);
+
+   for (int j = 1; j < num_conn; j++)
+   {
+      const int jump = J[j] - J[j-1];
+      const int dist = std::abs(jump);
+      max_jump = std::max(max_jump, jump);
+      min_jump = std::min(min_jump, jump);
+      sum_jump += jump;
+      sum_dist += dist;
+      // Put 'jump' in the appropriate bin.
+      if (jump == 0)
+      {
+         bins[num_bins]++;
+         continue;
+      }
+      for (int bin_id = 0, bin_max = bin_factor; true;
+           bin_id++, bin_max *= bin_factor)
+      {
+         if (bin_id < num_bins-1)
+         {
+            if (dist < bin_max)
+            {
+               if (jump > 0) { bins[num_bins+1+bin_id]++; }
+               else { bins[num_bins-1-bin_id]++; }
+               break;
+            }
+         }
+         else
+         {
+            if (jump > 0) { bins[2*num_bins]++; }
+            else { bins[0]++; }
+            break;
+         }
+      }
+   }
+   // Save precision and flags.
+   streamsize old_prec = out.precision(4);
+   ios_base::fmtflags old_flags = out.flags();
+   out << fixed;
+   out << "  jumps between consecutive column indices:"
+       << "\n   minimal  = " << min_jump
+       << "\n   maximal  = " << max_jump
+       << "\n   average  = " << 1.*sum_jump/num_conn
+       << "\n   avg dist = " << 1.*sum_dist/num_conn
+       << "\n  distribution of the jumps, positive (+) and negative (-):";
+   out << "\n                {0} : " << right << setw(8)
+       << 100.*bins[num_bins]/num_conn << "% (" << bins[num_bins] << ")";
+   for (int bin_id = 0, bin_min = 1; bin_id < num_bins;
+        bin_id++, bin_min *= bin_factor)
+   {
+      out << "\n   [" << setw(6) << bin_min << ", ";
+      if (bin_id < num_bins-1)
+      {
+         out << setw(6) << bin_min*bin_factor;
+      }
+      else
+      {
+         out << "     ∞";
+      }
+      const int n_neg = bins[num_bins-1-bin_id];
+      const int n_pos = bins[num_bins+1+bin_id];
+      out << ") : " << right << setw(8) << 100.*(n_neg+n_pos)/num_conn
+          << "% = (+) " << setw(8)
+          << 100.*n_pos/num_conn << "% + (-) " << setw(8)
+          << 100.*n_neg/num_conn << "% ("
+          << (n_neg+n_pos) << " = " << n_pos << " + " << n_neg << ")";
+   }
+   out << endl;
+   // Restore precision and flags.
+   out.precision(old_prec);
+   out.flags(old_flags);
+}
+
 void Table::Save(std::ostream &out) const
 {
    out << size << '\n';

general/table.hpp

@@ -16,6 +16,9 @@
 
 #include "mem_alloc.hpp"
 #include "array.hpp"
+#ifdef MFEM_USE_GECKO
+#include <graph.h>
+#endif
 
 namespace mfem
 {
@@ -33,6 +36,27 @@ struct Connection
 };
 
 
+#ifdef MFEM_USE_GECKO
+class GeckoParameters
+{
+public:
+   Gecko::Functional *functional;
+   unsigned int iterations;  ///< number of V cycles
+   unsigned int window;      ///< initial window size
+   unsigned int period;      ///< iterations between window increment
+   unsigned int seed;        ///< random number seed
+
+   /// Constructor. Assumes ownership of @a f.
+   /** If @a f is NULL (default), an instance of Gecko::FunctionalGeometric is
+       used. */
+   GeckoParameters(Gecko::Functional *f = NULL)
+      : functional(f ? f : new Gecko::FunctionalGeometric()),
+        iterations(1), window(2), period(1), seed(0) { }
+   ~GeckoParameters() { delete functional; }
+};
+#endif // #ifdef MFEM_USE_GECKO
+
+
 /** Data type Table. Table stores the connectivity of elements of TYPE I
     to elements of TYPE II, for example, it may be Element-To-Face
     connectivity table, etc. */
@@ -140,10 +164,31 @@ public:
    /// Call this if data has been stolen.
    void LoseData() { size = -1; I = J = NULL; }
 
+   /** @brief Assuming a symmetric Table, compute a row (and column) reordering
+       using the Cuthill-McKee (CM) algorithm. */
+   void GetCMReordering(Array<int> &ordering, bool reverse = false) const;
+
+#ifdef MFEM_USE_GECKO
+   /** @brief Assuming a symmetric Table, compute a row (and column) reordering
+       using the Gecko library. */
+   void GetGeckoReordering(const GeckoParameters &g_params,
+                           Array<int> &ordering) const;
+#endif
+
+#ifdef MFEM_USE_METIS
+   /** @brief Assuming a symmetric Table, compute a row (and column) reordering
+       using the Metis library. */
+   void GetMetisReordering(Array<int> &ordering, int type = 0,
+                           bool check_for_diag = true) const;
+#endif
+
    /// Prints the table to stream out.
    void Print(std::ostream & out = std::cout, int width = 4) const;
    void PrintMatlab(std::ostream & out) const;
 
+   /// Print statistics about the ordering of the J array.
+   void PrintOrderingStats(std::ostream &out = std::cout) const;
+
    void Save(std::ostream &out) const;
    void Load(std::istream &in);
 

linalg/operator.cpp

@@ -146,4 +146,24 @@ void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
    }
 }
 
+
+int PermutationOperator::CheckPermutation() const
+{
+   // Make sure 'perm' is a permutation of [0,height)
+   const int n = Height();
+   if (Width() != n) { return 1; }
+   if (perm.Size() != n) { return 2; }
+   Array<int> inv_perm(n);
+   inv_perm = n;
+   for (int oi = 0; oi < n; oi++)
+   {
+      const int sni = perm[oi];
+      const int ni = sni >= 0 ? sni : -1-sni;
+      if (ni >= n) { return 3; }
+      if (inv_perm[ni] != n) { return 4; }
+      inv_perm[ni] = oi;
+   }
+   return 0;
+}
+
 }

linalg/operator.hpp

@@ -424,6 +424,67 @@ public:
    virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
 };
 
+
+/// Local permutation operator.
+class PermutationOperator : public Operator
+{
+protected:
+   Array<int> perm;
+
+public:
+   /// Construct a (square) permutation Operator.
+   /** This constructor allocates a new permutation array without initializing
+       it. Use the method GetPermutation() to set the permutation. */
+   explicit PermutationOperator(int n) : Operator(n), perm(n) { }
+
+   /// Construct a (square) permutation Operator.
+   /** This constructor wraps already existing permutation array, @a p. If the
+       parameter @a own_p is true, the PermutationOperator assumes ownership of
+       the array @a p. */
+   PermutationOperator(int n, int *p, bool own_p = true)
+      : Operator(n), perm(p, n) { if (own_p) { perm.MakeDataOwner(); } }
+
+   // Default destructor.
+
+   /// Read + write access to the permutation array.
+   /** The permutation array, `perm`, defines the operator as follows:
+
+           x_new(   perm[oi]) = +x_old(oi),   perm[oi] >= 0,
+           x_new(-1-perm[oi]) = -x_old(oi),   perm[oi] <  0.
+
+       This method should be used to initialize the permutation array. */
+   Array<int> &GetPermutation() { return perm; }
+
+   /// Read-only access to the permutation array.
+   /** The permutation array, `perm`, defines the operator as follows:
+
+           x_new(   perm[oi]) = +x_old(oi),   perm[oi] >= 0,
+           x_new(-1-perm[oi]) = -x_old(oi),   perm[oi] <  0.
+   */
+   const Array<int> &GetPermutation() const { return perm; }
+
+   /// Return 0 if the permutation array defines a valid permutation.
+   int CheckPermutation() const;
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      for (int oi = 0; oi < Width(); oi++)
+      {
+         const int ni = perm[oi];
+         ni >= 0 ? y(ni) = x(oi) : y(-1-ni) = -x(oi);
+      }
+   }
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   {
+      for (int oi = 0; oi < Width(); oi++)
+      {
+         const int ni = perm[oi];
+         ni >= 0 ? y(oi) = x(ni) : y(oi) = -x(-1-ni);
+      }
+   }
+};
+
 }
 
 #endif

linalg/vector.cpp

@@ -516,15 +516,18 @@ void Vector::SetSubVector(const Array<int> &dofs, const double value)
 void Vector::SetSubVector(const Array<int> &dofs, const Vector &elemvect)
 {
    int i, j, n = dofs.Size();
+   MFEM_ASSERT(n == elemvect.Size(), "");
 
    for (i = 0; i < n; i++)
    {
       if ((j=dofs[i]) >= 0)
       {
+         MFEM_ASSERT(j < size, "");
          data[j] = elemvect(i);
       }
       else
       {
+         MFEM_ASSERT(-1-j < size, "");
          data[-1-j] = -elemvect(i);
       }
    }

makefile

@@ -338,6 +338,12 @@ test:
 	   if [ 0 -ne $${ERR} ]; then echo "Some tests failed."; exit 1; \
 	   else echo "All tests passed."; fi
 
+define subdir_rule
+$(1)/%: lib
+	$$(MAKE) -C $$(BLD)$$(@D) $$(@F)
+endef
+$(foreach dir,$(EM_DIRS),$(eval $(call subdir_rule,$(dir))))
+
 ALL_CLEAN_SUBDIRS = $(addsuffix /clean,config $(EM_DIRS) doc)
 .PHONY: $(ALL_CLEAN_SUBDIRS) miniapps/clean
 miniapps/clean: $(addsuffix /clean,$(MINIAPP_DIRS))

mesh/mesh.cpp

@@ -31,9 +31,6 @@
 #include "metis.h"
 #endif
 
-#ifdef MFEM_USE_GECKO
-#include "graph.h"
-#endif
 
 using namespace std;
 
@@ -764,6 +761,7 @@ void Mesh::Init()
    NURBSext = NULL;
    ncmesh = NULL;
    last_operation = Mesh::NONE;
+   el_perm = NULL;
 }
 
 void Mesh::InitTables()
@@ -1131,54 +1129,10 @@ void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
    meshgen = 2;
 }
 
-
-#ifdef MFEM_USE_GECKO
-void Mesh::GetGeckoElementReordering(Array<int> &ordering)
-{
-   Gecko::Graph graph;
-
-   // We will put some accesors in for these later
-   Gecko::Functional *functional =
-      new Gecko::FunctionalGeometric(); // ordering functional
-   unsigned int iterations = 1;         // number of V cycles
-   unsigned int window = 2;             // initial window size
-   unsigned int period = 1;             // iterations between window increment
-   unsigned int seed = 0;               // random number seed
-
-   // Run through all the elements and insert the nodes in the graph for them
-   for (int elemid = 0; elemid < GetNE(); ++elemid)
-   {
-      graph.insert();
-   }
-
-   // Run through all the elems and insert arcs to the graph for each element
-   // face Indices in Gecko are 1 based hence the +1 on the insertion
-   const Table &my_el_to_el = ElementToElementTable();
-   for (int elemid = 0; elemid < GetNE(); ++elemid)
-   {
-      const int *neighid = my_el_to_el.GetRow(elemid);
-      for (int i = 0; i < my_el_to_el.RowSize(elemid); ++i)
-      {
-         graph.insert(elemid + 1,  neighid[i] + 1);
-      }
-   }
-
-   // Get the reordering from Gecko and copy it into the ordering Array<int>
-   graph.order(functional, iterations, window, period, seed);
-   ordering.DeleteAll();
-   ordering.SetSize(GetNE());
-   Gecko::Node::Index NE = GetNE();
-   for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
-   {
-      ordering[gnodeid - 1] = graph.rank(gnodeid);
-   }
-
-   delete functional;
-}
-#endif
-
-
-void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
+void Mesh::ReorderElements_internal(const Array<int> &ordering,
+                                    bool reorder_vertices,
+                                    Array<int> &vertex_ordering,
+                                    bool update_nodes)
 {
    if (NURBSext)
    {
@@ -1191,7 +1145,16 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
                    " supported.");
       return;
    }
-   MFEM_VERIFY(ordering.Size() == GetNE(), "invalid reordering array.")
+   MFEM_VERIFY(ordering.Size() == GetNE(), "invalid reordering array.");
+#ifdef MFEM_DEBUG
+   {
+      // Make sure 'ordering' is a permutation of [0,NumOfElements)
+      PermutationOperator po(GetNE(), const_cast<int*>(ordering.GetData()),
+                             false);
+      int err = po.CheckPermutation();
+      MFEM_VERIFY(!err, "invalid reordering array.");
+   }
+#endif
 
    // Data members that need to be updated:
 
@@ -1215,22 +1178,8 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
 
    // - Nodes
 
-   // Save the locations of the Nodes so we can rebuild them later
-   Array<Vector*> old_elem_node_vals;
-   FiniteElementSpace *nodes_fes = NULL;
-   if (Nodes)
-   {
-      old_elem_node_vals.SetSize(GetNE());
-      nodes_fes = Nodes->FESpace();
-      Array<int> old_dofs;
-      Vector vals;
-      for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
-      {
-         nodes_fes->GetElementVDofs(old_elid, old_dofs);
-         Nodes->GetSubVector(old_dofs, vals);
-         old_elem_node_vals[old_elid] = new Vector(vals);
-      }
-   }
+   // Destroy tables that need to be rebuild
+   DeleteTables();
 
    // Get the newly ordered elements
    Array<Element *> new_elements(GetNE());
@@ -1246,7 +1195,7 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
    {
       // Get the new vertex ordering permutation vectors and fill the new
       // vertices
-      Array<int> vertex_ordering(GetNV());
+      vertex_ordering.SetSize(GetNV());
       vertex_ordering = -1;
       Array<Vertex> new_vertices(GetNV());
       int new_vertex_ind = 0;
@@ -1292,9 +1241,6 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
       }
    }
 
-   // Destroy tables that need to be rebuild
-   DeleteTables();
-
    if (Dim > 1)
    {
       // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
@@ -1309,18 +1255,14 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
    // Update faces and faces_info
    GenerateFaces();
 
-   // Build the nodes from the saved locations if they were around before
-   if (Nodes)
+   // Element reordering is a `Mesh::Operation`.
+   last_operation = Mesh::REORDER;
+   sequence++;
+   el_perm = ordering.GetData();
+
+   if (Nodes && update_nodes)
    {
-      nodes_fes->Update();
-      Array<int> new_dofs;
-      for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
-      {
-         int new_elid = ordering[old_elid];
-         nodes_fes->GetElementVDofs(new_elid, new_dofs);
-         Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
-         delete old_elem_node_vals[old_elid];
-      }
+      Nodes->Update();
    }
 }
 
@@ -2296,6 +2238,7 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    // Create the new Mesh instance without a record of its refinement history
    sequence = 0;
    last_operation = Mesh::NONE;
+   el_perm = NULL;
 
    // Duplicate the elements
    elements.SetSize(NumOfElements);
@@ -3043,6 +2986,9 @@ void Mesh::KnotInsert(Array<KnotVector *> &kv)
 
    NURBSext->KnotInsert(kv);
 
+   last_operation = Mesh::REFINE;
+   sequence++;
+
    UpdateNURBS();
 }
 
@@ -3335,6 +3281,10 @@ static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
 
 int Mesh::CheckElementOrientation(bool fix_it)
 {
+   // Note: this public operation may change the mesh topology (rotate elements)
+   //       if 'fix_it' is true, however, it is not a 'Mesh::Operation', i.e. it
+   //       does not change 'last_operation' and 'sequence'.
+
    int i, j, k, wo = 0, fo = 0, *vi = 0;
    double *v[4];
 
@@ -3542,6 +3492,11 @@ int Mesh::GetQuadOrientation(const int *base, const int *test)
 
 int Mesh::CheckBdrElementOrientation(bool fix_it)
 {
+   // Note: this public operation may change the mesh topology (rotate boundary
+   //       elements) if 'fix_it' is true, however, it is not a
+   //       'Mesh::Operation', i.e. it does not change 'last_operation' and
+   //       'sequence'.
+
    int i, wo = 0;
 
    if (Dim == 2)
@@ -4120,10 +4075,13 @@ const Table & Mesh::ElementToElementTable()
 {
    if (el_to_el)
    {
+      MFEM_ASSERT(el_to_el->Size() == GetNE(),
+                  "internal error, el_to_el->Size() = "
+                  << el_to_el->Size() << ", GetNE() = " << GetNE());
       return *el_to_el;
    }
 
-   int num_faces = GetNumFaces();
+   MFEM_DEBUG_DO(int num_faces = GetNumFaces());
    // Note that, for ParNCMeshes, faces_info will contain also the ghost faces
    MFEM_ASSERT(faces_info.Size() >= num_faces, "faces were not generated!");
 
@@ -4159,6 +4117,7 @@ const Table & Mesh::ElementToFaceTable() const
    {
       mfem_error("Mesh::ElementToFaceTable()");
    }
+   MFEM_ASSERT(el_to_face->Size() == GetNE(), "internal error");
    return *el_to_face;
 }
 
@@ -4168,6 +4127,7 @@ const Table & Mesh::ElementToEdgeTable() const
    {
       mfem_error("Mesh::ElementToEdgeTable()");
    }
+   MFEM_ASSERT(el_to_edge->Size() == GetNE(), "internal error");
    return *el_to_edge;
 }
 
@@ -4416,10 +4376,7 @@ STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
    int i, *v;
    STable3D *faces_tbl;
 
-   if (el_to_face != NULL)
-   {
-      delete el_to_face;
-   }
+   delete el_to_face;
    el_to_face = new Table(NumOfElements, 6);  // must be 6 for hexahedra
    faces_tbl = new STable3D(NumOfVertices);
    for (i = 0; i < NumOfElements; i++)
@@ -4486,7 +4443,7 @@ STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
 
 void Mesh::ReorientTetMesh()
 {
-   int *v;
+   // This operation could be a 'Mesh::Operation'.
 
    if (Dim != 3 || !(meshgen & 1))
    {
@@ -4505,7 +4462,7 @@ void Mesh::ReorientTetMesh()
    {
       if (GetElementType(i) == Element::TETRAHEDRON)
       {
-         v = elements[i]->GetVertices();
+         int *v = elements[i]->GetVertices();
 
          Rotate3(v[0], v[1], v[2]);
          if (v[0] < v[3])
@@ -4523,7 +4480,7 @@ void Mesh::ReorientTetMesh()
    {
       if (GetBdrElementType(i) == Element::TRIANGLE)
       {
-         v = boundary[i]->GetVertices();
+         int *v = boundary[i]->GetVertices();
 
          Rotate3(v[0], v[1], v[2]);
       }
@@ -4630,7 +4587,7 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
 #else
       int ncon = 1;
       int err;
-      int options[40];
+      int options[METIS_NOPTIONS];
 #endif
       int edgecut;
 
@@ -4687,9 +4644,11 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
                                         options,
                                         &edgecut,
                                         partitioning);
-         if (err != 1)
+         if (err != METIS_OK)
+         {
             mfem_error("Mesh::GeneratePartitioning: "
                        " error in METIS_PartGraphRecursive!");
+         }
 #endif
       }
 
@@ -4723,9 +4682,11 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
                                    options,
                                    &edgecut,
                                    partitioning);
-         if (err != 1)
+         if (err != METIS_OK)
+         {
             mfem_error("Mesh::GeneratePartitioning: "
                        " error in METIS_PartGraphKway!");
+         }
 #endif
       }
 
@@ -4760,9 +4721,11 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
                                    options,
                                    &edgecut,
                                    partitioning);
-         if (err != 1)
+         if (err != METIS_OK)
+         {
             mfem_error("Mesh::GeneratePartitioning: "
                        " error in METIS_PartGraphKway!");
+         }
 #endif
       }
 
@@ -5469,6 +5432,14 @@ void Mesh::QuadUniformRefinement()
    int i, j, *v, vv[2], attr;
    const int *e;
 
+   // Call DeleteTables() but preserve el_to_edge, if built.
+   {
+      Table *elem_to_edge = el_to_edge;
+      el_to_edge = NULL;
+      DeleteTables();
+      el_to_edge = elem_to_edge;
+   }
+
    if (el_to_edge == NULL)
    {
       el_to_edge = new Table;
@@ -5573,10 +5544,24 @@ void Mesh::HexUniformRefinement()
    const int *e, *f;
    int vv[4];
 
+   // Call DeleteTables() but preserve el_to_edge, el_to_face, and bel_to_edge,
+   // if built.
+   {
+      Table *elem_to_edge = el_to_edge;
+      Table *elem_to_face = el_to_face;
+      Table *belem_to_edge = bel_to_edge;
+      el_to_edge = el_to_face = bel_to_edge = NULL;
+      DeleteTables();
+      el_to_edge = elem_to_edge;
+      el_to_face = elem_to_face;
+      bel_to_edge = belem_to_edge;
+   }
+
    if (el_to_edge == NULL)
    {
       el_to_edge = new Table;
       NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      // the above call creates 'bel_to_edge' too
    }
    if (el_to_face == NULL)
    {
@@ -5774,10 +5759,14 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
    } // end of 'if (Dim == 1)'
    else if (Dim == 2) // ---------------------------------------------------
    {
+      const bool have_el_to_edge = (el_to_edge != NULL);
+
       // 1. Get table of vertex to vertex connections.
       DSTable v_to_v(NumOfVertices);
       GetVertexToVertexTable(v_to_v);
 
+      DeleteTables();
+
       // 2. Get edge to element connections in arrays edge1 and edge2
       nedges = v_to_v.NumberOfEntries();
       int *edge1  = new int[nedges];
@@ -5841,8 +5830,10 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
                boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
             }
             else
+            {
                mfem_error("Only bisection of segment is implemented"
                           " for bdr elem.");
+            }
          }
       }
       NumOfBdrElements = boundary.Size();
@@ -5852,8 +5843,9 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
       delete [] edge2;
       delete [] middle;
 
-      if (el_to_edge != NULL)
+      if (have_el_to_edge)
       {
+         el_to_edge = new Table; // el_to_edge was deleted by DeleteTables()
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
          GenerateFaces();
       }
@@ -5861,10 +5853,15 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
    }
    else if (Dim == 3) // ---------------------------------------------------
    {
+      const bool have_el_to_edge = (el_to_edge != NULL);
+      const bool have_el_to_face = (el_to_face != NULL);
+
       // 1. Get table of vertex to vertex connections.
       DSTable v_to_v(NumOfVertices);
       GetVertexToVertexTable(v_to_v);
 
+      DeleteTables();
+
       // 2. Get edge to element connections in arrays edge1 and edge2
       nedges = v_to_v.NumberOfEntries();
       int *middle = new int[nedges];
@@ -5967,11 +5964,12 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 7. Free the allocated memory.
       delete [] middle;
 
-      if (el_to_edge != NULL)
+      if (have_el_to_edge)
       {
+         el_to_edge = new Table; // was deleted by DeleteTables()
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
-      if (el_to_face != NULL)
+      if (have_el_to_face)
       {
          GetElementToFaceTable();
          GenerateFaces();
@@ -8307,6 +8305,10 @@ void Mesh::Transform(VectorCoefficient &deformation)
 
 void Mesh::RemoveUnusedVertices()
 {
+   // Note: this public operation may change the mesh topology, however, it is
+   //       not a 'Mesh::Operation', i.e. it does not change 'last_operation'
+   //       and 'sequence'.
+
    if (NURBSext || ncmesh) { return; }
 
    Array<int> v2v(GetNV());
@@ -8414,6 +8416,10 @@ void Mesh::RemoveUnusedVertices()
 
 void Mesh::RemoveInternalBoundaries()
 {
+   // Note: this public operation may change the mesh topology, however, it is
+   //       not a 'Mesh::Operation', i.e. it does not change 'last_operation'
+   //       and 'sequence'.
+
    if (NURBSext || ncmesh) { return; }
 
    int num_bdr_elem = 0;

mesh/mesh.hpp

@@ -164,7 +164,9 @@ public:
    typedef Geometry::Constants<Geometry::TETRAHEDRON> tet_t;
    typedef Geometry::Constants<Geometry::CUBE>        hex_t;
 
-   enum Operation { NONE, REFINE, DEREFINE, REBALANCE };
+   /// Mesh operations that, generally, change the topology of the mesh.
+   /** This enumeration defines the values returned by GetLastOperation(). */
+   enum Operation { NONE, REFINE, DEREFINE, REBALANCE, REORDER };
 
    /// A list of all unique element attributes used by the Mesh.
    Array<int> attributes;
@@ -176,6 +178,7 @@ public:
 
 protected:
    Operation last_operation;
+   const int *el_perm;  ///< Pointer to the last element reordering array.
 
    void Init();
    void InitTables();
@@ -217,6 +220,11 @@ protected:
        reference element at the center of the element. */
    void GetElementJacobian(int i, DenseMatrix &J);
 
+   void ReorderElements_internal(const Array<int> &ordering,
+                                 bool reorder_vertices,
+                                 Array<int> &vertex_ordering,
+                                 bool update_nodes);
+
    void MarkForRefinement();
    void MarkTriMeshForRefinement();
    void GetEdgeOrdering(DSTable &v_to_v, Array<int> &order);
@@ -515,18 +523,53 @@ public:
 
    void SetAttributes();
 
+   /// Generate element reordering using the Cuthill-McKee (CM) algorithm.
+   /** The generated re-ordering can be applied to the Mesh using the method
+       ReorderElements(). */
+   void GetCMElementReordering(Array<int> &ordering, bool reverse = false)
+   {
+      ElementToElementTable().GetCMReordering(ordering, reverse);
+   }
+
 #ifdef MFEM_USE_GECKO
    /** This is our integration with the Gecko library.  This will call the
        Gecko library to find an element ordering that will increase memory
        coherency by putting elements that are in physical proximity closer in
        memory. */
-   void GetGeckoElementReordering(Array<int> &ordering);
+   void GetGeckoElementReordering(Array<int> &ordering)
+   {
+      ElementToElementTable().GetGeckoReordering(GeckoParameters(), ordering);
+   }
+#endif
+
+#ifdef MFEM_USE_METIS
+   /// Generate element reordering using the Metis library.
+   /** The generated re-ordering can be applied to the Mesh using the method
+       ReorderElements(). */
+   void GetMetisElementReordering(Array<int> &ordering, int type = 0)
+   {
+      const bool check_diag = false;
+      ElementToElementTable().GetMetisReordering(ordering, type, check_diag);
+   }
 #endif
 
    /** Rebuilds the mesh with a different order of elements.  The ordering
-       vector maps the old element number to the new element number.  This also
-       reorders the vertices and nodes edges and faces along with the elements.  */
-   void ReorderElements(const Array<int> &ordering, bool reorder_vertices = true);
+       vector maps the old element number to the new element number, i.e.
+       `new_element_id = ordering[old_element_id]`.  This also reorders the
+       vertices (if @a reorder_vertices is true), as well as the edges and faces
+       along with the elements. */
+   virtual void ReorderElements(const Array<int> &ordering,
+                                bool reorder_vertices = true)
+   {
+      const bool update_nodes = true;
+      Array<int> vertex_ordering;
+      ReorderElements_internal(ordering, reorder_vertices, vertex_ordering,
+                               update_nodes);
+   }
+
+   /** @brief Returns the permutation/ordering used in the last call to
+       ReorderElements(). */
+   const int *GetElementPermutation() const { return el_perm; }
 
    /** Creates mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into
        nx*ny*nz hexahedrals if type=HEXAHEDRON or into 6*nx*ny*nz tetrahedrons

mesh/pmesh.cpp

@@ -986,9 +986,15 @@ ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
    group_sedge.ShiftUpI();
    group_sface.ShiftUpI();
 
+   SetSharedToLocalMaps();
+}
+
+void ParMesh::SetSharedToLocalMaps()
+{
    // determine sedge_ledge
    if (shared_edges.Size() > 0)
    {
+      sedge_ledge.SetSize(shared_edges.Size());
       DSTable v_to_v(NumOfVertices);
       GetVertexToVertexTable(v_to_v);
       for (int se = 0; se < shared_edges.Size(); se++)
@@ -1003,6 +1009,7 @@ ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
    // determine sface_lface
    if (shared_faces.Size() > 0)
    {
+      sface_lface.SetSize(shared_faces.Size());
       STable3D *faces_tbl = GetFacesTable();
       for (int sf = 0; sf < shared_faces.Size(); sf++)
       {
@@ -1389,10 +1396,12 @@ void ParMesh::ExchangeFaceNbrData()
 
    int num_face_nbrs = 0;
    for (int g = 1; g < GetNGroups(); g++)
+   {
       if (gr_sface->RowSize(g-1) > 0)
       {
          num_face_nbrs++;
       }
+   }
 
    face_nbr_group.SetSize(num_face_nbrs);
 
@@ -1407,12 +1416,15 @@ void ParMesh::ExchangeFaceNbrData()
       Array<Pair<int, int> > rank_group(num_face_nbrs);
 
       for (int g = 1, counter = 0; g < GetNGroups(); g++)
+      {
          if (gr_sface->RowSize(g-1) > 0)
          {
 #ifdef MFEM_DEBUG
             if (gtopo.GetGroupSize(g) != 2)
+            {
                mfem_error("ParMesh::ExchangeFaceNbrData() : "
                           "group size is not 2!");
+            }
 #endif
             const int *nbs = gtopo.GetGroup(g);
             int lproc = (nbs[0]) ? nbs[0] : nbs[1];
@@ -1420,6 +1432,7 @@ void ParMesh::ExchangeFaceNbrData()
             rank_group[counter].two = g;
             counter++;
          }
+      }
 
       SortPairs<int, int>(rank_group, rank_group.Size());
 
@@ -1451,9 +1464,9 @@ void ParMesh::ExchangeFaceNbrData()
    send_face_nbr_facedata.MakeI(num_face_nbrs);
    for (int fn = 0; fn < num_face_nbrs; fn++)
    {
-      int nbr_group = face_nbr_group[fn];
-      int  num_sfaces = gr_sface->RowSize(nbr_group-1);
-      int *sface = gr_sface->GetRow(nbr_group-1);
+      int nbr_group    = face_nbr_group[fn];
+      int num_sfaces   = gr_sface->RowSize(nbr_group-1);
+      const int *sface = gr_sface->GetRow(nbr_group-1);
       for (int i = 0; i < num_sfaces; i++)
       {
          int lface = s2l_face[sface[i]];
@@ -1466,11 +1479,13 @@ void ParMesh::ExchangeFaceNbrData()
             const int nv = elements[el]->GetNVertices();
             const int *v = elements[el]->GetVertices();
             for (int j = 0; j < nv; j++)
+            {
                if (vertex_marker[v[j]] != fn)
                {
                   vertex_marker[v[j]] = fn;
                   send_face_nbr_vertices.AddAColumnInRow(fn);
                }
+            }
 
             send_face_nbr_elemdata.AddColumnsInRow(fn, nv + 2);
          }
@@ -1497,9 +1512,9 @@ void ParMesh::ExchangeFaceNbrData()
    vertex_marker = -1;
    for (int fn = 0; fn < num_face_nbrs; fn++)
    {
-      int nbr_group = face_nbr_group[fn];
-      int  num_sfaces = gr_sface->RowSize(nbr_group-1);
-      int *sface = gr_sface->GetRow(nbr_group-1);
+      int nbr_group    = face_nbr_group[fn];
+      int num_sfaces   = gr_sface->RowSize(nbr_group-1);
+      const int *sface = gr_sface->GetRow(nbr_group-1);
       for (int i = 0; i < num_sfaces; i++)
       {
          int lface = s2l_face[sface[i]];
@@ -1512,11 +1527,13 @@ void ParMesh::ExchangeFaceNbrData()
             const int nv = elements[el]->GetNVertices();
             const int *v = elements[el]->GetVertices();
             for (int j = 0; j < nv; j++)
+            {
                if (vertex_marker[v[j]] != fn)
                {
                   vertex_marker[v[j]] = fn;
                   send_face_nbr_vertices.AddConnection(fn, v[j]);
                }
+            }
 
             send_face_nbr_elemdata.AddConnection(fn, GetAttribute(el));
             send_face_nbr_elemdata.AddConnection(
@@ -1553,13 +1570,13 @@ void ParMesh::ExchangeFaceNbrData()
    // convert the element indices in send_face_nbr_facedata
    for (int fn = 0; fn < num_face_nbrs; fn++)
    {
-      int  num_elems  = send_face_nbr_elements.RowSize(fn);
-      int *elems      = send_face_nbr_elements.GetRow(fn);
-      int  num_verts  = send_face_nbr_vertices.RowSize(fn);
-      int *verts      = send_face_nbr_vertices.GetRow(fn);
-      int *elemdata   = send_face_nbr_elemdata.GetRow(fn);
-      int  num_sfaces = send_face_nbr_facedata.RowSize(fn)/2;
-      int *facedata   = send_face_nbr_facedata.GetRow(fn);
+      int  num_elems   = send_face_nbr_elements.RowSize(fn);
+      const int *elems = send_face_nbr_elements.GetRow(fn);
+      int  num_verts   = send_face_nbr_vertices.RowSize(fn);
+      const int *verts = send_face_nbr_vertices.GetRow(fn);
+      int *elemdata    = send_face_nbr_elemdata.GetRow(fn);
+      int  num_sfaces  = send_face_nbr_facedata.RowSize(fn)/2;
+      int *facedata    = send_face_nbr_facedata.GetRow(fn);
 
       for (int i = 0; i < num_verts; i++)
       {
@@ -1687,11 +1704,11 @@ void ParMesh::ExchangeFaceNbrData()
          break;
       }
 
-      int  elem_off   = face_nbr_elements_offset[fn];
-      int  nbr_group  = face_nbr_group[fn];
-      int  num_sfaces = gr_sface->RowSize(nbr_group-1);
-      int *sface      = gr_sface->GetRow(nbr_group-1);
-      int *facedata =
+      int elem_off     = face_nbr_elements_offset[fn];
+      int nbr_group    = face_nbr_group[fn];
+      int num_sfaces   = gr_sface->RowSize(nbr_group-1);
+      const int *sface = gr_sface->GetRow(nbr_group-1);
+      const int *facedata =
          &recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]];
 
       for (int i = 0; i < num_sfaces; i++)
@@ -2110,10 +2127,15 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
          uniform_refinement = 1;
       }
 
+      const bool have_el_to_edge = (el_to_edge != NULL);
+      const bool have_el_to_face = (el_to_face != NULL);
+
       // 1. Get table of vertex to vertex connections.
       DSTable v_to_v(NumOfVertices);
       GetVertexToVertexTable(v_to_v);
 
+      DeleteTables();
+
       // 2. Create a marker array for all edges (vertex to vertex connections).
       Array<int> middle(v_to_v.NumberOfEntries());
       middle = -1;
@@ -2312,13 +2334,15 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       while (need_refinement == 1);
 
       if (NumOfBdrElements != boundary.Size())
+      {
          mfem_error("ParMesh::LocalRefinement :"
                     " (NumOfBdrElements != boundary.Size())");
+      }
 
       // 5a. Update the groups after refinement.
-      if (el_to_face != NULL)
+      RefineGroups(v_to_v, middle);
+      if (have_el_to_face)
       {
-         RefineGroups(v_to_v, middle);
          // GetElementToFaceTable(); // Called by RefineGroups
          GenerateFaces();
       }
@@ -2340,7 +2364,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 7. Free the allocated memory.
       middle.DeleteAll();
 
-      if (el_to_edge != NULL)
+      if (have_el_to_edge)
       {
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
@@ -2356,10 +2380,14 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
          uniform_refinement = 1;
       }
 
+      const bool have_el_to_edge = (el_to_edge != NULL);
+
       // 1. Get table of vertex to vertex connections.
       DSTable v_to_v(NumOfVertices);
       GetVertexToVertexTable(v_to_v);
 
+      DeleteTables();
+
       // 2. Get edge to element connections in arrays edge1 and edge2
       int nedges  = v_to_v.NumberOfEntries();
       int *edge1  = new int[nedges];
@@ -2564,7 +2592,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       delete [] edge2;
       delete [] middle;
 
-      if (el_to_edge != NULL)
+      if (have_el_to_edge)
       {
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
          GenerateFaces();
@@ -2744,6 +2772,106 @@ void ParMesh::Rebalance()
    }
 }
 
+void ParMesh::ReorderElements(const Array<int> &ordering,
+                              bool reorder_vertices)
+{
+   if (NURBSext)
+   {
+      MFEM_WARNING("element reordering of NURBS meshes is not supported.");
+      return;
+   }
+   if (Nonconforming())
+   {
+      MFEM_WARNING("element reordering of non-conforming meshes is not"
+                   " supported.");
+      return;
+   }
+
+   const Array<int>* s2l_face = NULL;
+   Array<Pair<int,int> > sface_info;
+   if (have_face_nbr_data)
+   {
+      // Save face-neighbor data from 'faces_info'.
+      s2l_face = ((Dim == 1) ? &svert_lvert :
+                  ((Dim == 2) ? &sedge_ledge : &sface_lface));
+      sface_info.SetSize(s2l_face->Size());
+      for (int sf = 0; sf < sface_info.Size(); ++sf)
+      {
+         const FaceInfo &fi = faces_info[(*s2l_face)[sf]];
+         sface_info[sf] = Pair<int,int>(fi.Elem2No, fi.Elem2Inf);
+      }
+   }
+
+   const bool update_nodes = false;
+   Array<int> vertex_ordering;
+   ReorderElements_internal(ordering, reorder_vertices, vertex_ordering,
+                            update_nodes);
+
+   if (reorder_vertices)
+   {
+      // Replace the vertex ids in the 'shared_faces' and 'shared edges' with
+      // the reordered vertex numbers.
+      for (int sf_id = 0; sf_id < shared_faces.Size(); ++sf_id)
+      {
+         int *v = shared_faces[sf_id]->GetVertices();
+         int nv = shared_faces[sf_id]->GetNVertices();
+         for (int vi = 0; vi < nv; ++vi)
+         {
+            v[vi] = vertex_ordering[v[vi]];
+         }
+      }
+      for (int se_id = 0; se_id < shared_edges.Size(); ++se_id)
+      {
+         int *v = shared_edges[se_id]->GetVertices();
+         int nv = shared_edges[se_id]->GetNVertices();
+         for (int vi = 0; vi < nv; ++vi)
+         {
+            v[vi] = vertex_ordering[v[vi]];
+         }
+      }
+      // svert_lvert
+      for (int sv_id = 0; sv_id < svert_lvert.Size(); ++sv_id)
+      {
+         svert_lvert[sv_id] = vertex_ordering[svert_lvert[sv_id]];
+      }
+      // send_face_nbr_vertices
+      if (have_face_nbr_data)
+      {
+         int *J = send_face_nbr_vertices.GetJ();
+         const int nnz = send_face_nbr_vertices.Size_of_connections();
+         for (int j = 0; j < nnz; j++)
+         {
+            J[j] = vertex_ordering[J[j]];
+         }
+      }
+   }
+
+   // Set 'sedge_ledge' and 'sface_lface'.
+   SetSharedToLocalMaps();
+
+   if (have_face_nbr_data)
+   {
+      // Restore face-neighbor data to 'faces_info'.
+      for (int sf = 0; sf < sface_info.Size(); ++sf)
+      {
+         FaceInfo &fi = faces_info[(*s2l_face)[sf]];
+         MFEM_ASSERT(fi.Elem2No == -1 && fi.Elem2Inf == -1, "internal error");
+         fi.Elem2No = sface_info[sf].one;
+         fi.Elem2Inf = sface_info[sf].two;
+      }
+
+      // send_face_nbr_elements
+      int *J = send_face_nbr_elements.GetJ();
+      const int nnz = send_face_nbr_elements.Size_of_connections();
+      for (int j = 0; j < nnz; j++)
+      {
+         J[j] = ordering[J[j]];
+      }
+   }
+
+   if (Nodes) { Nodes->Update(); }
+}
+
 void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
 {
    int i, attr, newv[3], ind, f_ind, *v;
@@ -2782,6 +2910,7 @@ void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
    // overestimate the size of the J arrays
    if (Dim == 3)
    {
+      // Assuming only triangle shared faces
       J_group_svert = new int[group_svert.Size_of_connections()
                               + group_sedge.Size_of_connections()];
       J_group_sedge = new int[2*group_sedge.Size_of_connections()

mesh/pmesh.hpp

@@ -50,6 +50,9 @@ protected:
    /// Create from a nonconforming mesh.
    ParMesh(const ParNCMesh &pncmesh);
 
+   // Set sedge_ledge and sface_lface
+   void SetSharedToLocalMaps();
+
    // Mark all tets to ensure consistency across MPI tasks; also mark the
    // shared and boundary triangle faces using the consistently marked tets.
    virtual void MarkTetMeshForRefinement(DSTable &v_to_v);
@@ -129,6 +132,7 @@ public:
    // Local face-neighbor elements and vertices ordered by face-neighbor
    Table            send_face_nbr_elements;
    Table            send_face_nbr_vertices;
+   // Note: additional face-neighbor data is stored in 'faces_info'.
 
    ParNCMesh* pncmesh;
 
@@ -183,6 +187,9 @@ public:
    /// Load balance the mesh. NC meshes only.
    void Rebalance();
 
+   virtual void ReorderElements(const Array<int> &ordering,
+                                bool reorder_vertices = true);
+
    /** Print the part of the mesh in the calling processor adding the interface
        as boundary (for visualization purposes) using the mfem v1.0 format. */
    virtual void Print(std::ostream &out = std::cout) const;

miniapps/meshing/mesh-explorer.cpp

@@ -94,6 +94,7 @@ Mesh *read_par_mesh(int np, const char *mesh_prefix)
    return mesh;
 }
 
+
 int main (int argc, char *argv[])
 {
    int np = 0;
@@ -193,6 +194,8 @@ int main (int argc, char *argv[])
            "k) View element ratios, kappa\n"
            "x) Print sub-element stats\n"
            "p) Generate a partitioning\n"
+           "o) Reorder the mesh elements\n"
+           "n) Convert to nonconforming mesh\n"
            "S) Save\n"
            "--> " << flush;
       char mk;
@@ -399,6 +402,7 @@ int main (int argc, char *argv[])
                << "\nmax kappa    = " << max_kappa << endl;
       }
 
+      // Choices that send data to GLVis.
       if (mk == 'm' || mk == 'b' || mk == 'e' || mk == 'v' || mk == 'h' ||
           mk == 'k' || mk == 'p')
       {
@@ -420,6 +424,7 @@ int main (int argc, char *argv[])
 
          if (mk == 'e')
          {
+#if USE_CHECKERBOARD_ELEMENT_COLORING
             Array<int> coloring;
             srand(time(0));
             double a = double(rand()) / (double(RAND_MAX) + 1.);
@@ -432,10 +437,14 @@ int main (int argc, char *argv[])
                attr(i) = coloring[i];
             }
             cout << "Number of colors: " << attr.Max() + 1 << endl;
+#endif
             for (int i = 0; i < mesh->GetNE(); i++)
             {
-               // part[i] = i; // checkerboard element coloring
+#if USE_CHECKERBOARD_ELEMENT_COLORING
+               part[i] = i; // checkerboard element coloring
+#else
                attr(i) = part[i] = i; // coloring by element number
+#endif
             }
          }
 
@@ -660,6 +669,70 @@ int main (int argc, char *argv[])
          delete attr_fespace;
       }
 
+      if (mk == 'o')
+      {
+         Array<int> ordering;
+
+         cout << "Current element-to-element Table stats:\n";
+         mesh->ElementToElementTable().PrintOrderingStats();
+
+         cout << "Choose reordering algorithm:\n"
+              "c) Cuthill-McKee\n"
+              "r) reverse Cuthill-McKee\n"
+#ifdef MFEM_USE_GECKO
+              "g) Gecko\n"
+#endif
+#ifdef MFEM_USE_METIS
+              "m) Metis\n"
+#endif
+              "other char) skip reordering\n"
+              "--> " << flush;
+         char pk;
+         cin >> pk;
+         if (pk == 'c')
+         {
+            cout << "Reordering the elements using Cuthill-McKee ..." << flush;
+            mesh->GetCMElementReordering(ordering, false);
+         }
+         else if (pk == 'r')
+         {
+            cout << "Reordering the elements using reverse Cuthill-McKee ..."
+                 << flush;
+            mesh->GetCMElementReordering(ordering, true);
+         }
+#ifdef MFEM_USE_GECKO
+         else if (pk == 'g')
+         {
+            cout << "Reordering the elements using Gecko ..." << flush;
+            mesh->GetGeckoElementReordering(ordering);
+         }
+#endif
+#ifdef MFEM_USE_METIS
+         else if (pk == 'm')
+         {
+            cout << "Reordering the elements using Metis ..." << flush;
+            mesh->GetMetisElementReordering(ordering);
+         }
+#endif
+         else
+         {
+            cout << "Skipping reordering." << endl;
+            continue;
+         }
+         bool reorder_vertices = true;
+         mesh->ReorderElements(ordering, reorder_vertices);
+         cout << " done." << endl;
+
+         cout << "New element-to-element Table stats:\n";
+         mesh->ElementToElementTable().PrintOrderingStats();
+      }
+
+      if (mk == 'n')
+      {
+         bool triangles_nonconforming = true;
+         mesh->EnsureNCMesh(triangles_nonconforming);
+      }
+
       if (mk == 'S')
       {
          const char mesh_file[] = "mesh-explorer.mesh";

miniapps/performance/ex1.cpp

@@ -168,6 +168,20 @@ int main(int argc, char *argv[])
                   " the LOR preconditioner yet");
    }
 
+   // FIXME: add a command line option ...
+   Array<int> ordering;
+   cout << "Computing mesh element reordering ..." << flush;
+   // cout << " (no reordering) ..." << flush;
+   cout << " (CM) ..." << flush; mesh->GetCMElementReordering(ordering);
+   // cout << " (Metis) ..." << flush; mesh->GetMetisElementReordering(ordering);
+   // cout << " (Gecko) ..." << flush; mesh->GetGeckoElementReordering(ordering);
+   cout << " done.\nApplying the mesh element reordering ..." << flush;
+   if (ordering.Size() > 0)
+   {
+      mesh->ReorderElements(ordering);
+   }
+   cout << " done." << endl;
+
    // 5. Define a finite element space on the mesh. Here we use continuous
    //    Lagrange finite elements of the specified order. If order < 1, we
    //    instead use an isoparametric/isogeometric space.
@@ -340,6 +354,9 @@ int main(int argc, char *argv[])
    cout << " done, " << tic_toc.RealTime() << "s." << endl;
 
    // Solve with CG or PCG, depending if the matrix A_pc is available
+   cout << "Solving the linear system ..." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
    if (pc_choice != NONE)
    {
       GSSmoother M(A_pc);
@@ -349,6 +366,9 @@ int main(int argc, char *argv[])
    {
       CG(*a_oper, B, X, 1, 500, 1e-12, 0.0);
    }
+   tic_toc.Stop();
+   cout << "Solving the linear system ... done, "
+        << tic_toc.RealTime() << " sec." << endl;
 
    // 13. Recover the solution as a finite element grid function.
    if (perf && matrix_free)

miniapps/performance/ex1p.cpp

@@ -2,17 +2,18 @@
 //
 // Compile with: make ex1p
 //
-// Sample runs:  mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -mf  -pc lor
-//               mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -asm -pc ho
-//               mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -asm -pc ho -sc
-//               mpirun -np 4 ex1p -m ../../data/fichera.mesh -std  -asm -pc ho
-//               mpirun -np 4 ex1p -m ../../data/fichera.mesh -std  -asm -pc ho -sc
-//               mpirun -np 4 ex1p -m ../../data/amr-hex.mesh -perf -asm -pc ho -sc
-//               mpirun -np 4 ex1p -m ../../data/amr-hex.mesh -std  -asm -pc ho -sc
-//               mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -perf -asm -pc ho  -sc
-//               mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -std  -asm -pc ho  -sc
-//               mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -perf -mf  -pc lor
-//               mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -std  -asm -pc ho  -sc
+// Sample runs:
+//    mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -mf  -pc lor
+//    mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -asm -pc ho
+//    mpirun -np 4 ex1p -m ../../data/fichera.mesh -perf -asm -pc ho -sc
+//    mpirun -np 4 ex1p -m ../../data/fichera.mesh -std  -asm -pc ho
+//    mpirun -np 4 ex1p -m ../../data/fichera.mesh -std  -asm -pc ho -sc
+//    mpirun -np 4 ex1p -m ../../data/amr-hex.mesh -perf -asm -pc ho -sc
+//    mpirun -np 4 ex1p -m ../../data/amr-hex.mesh -std  -asm -pc ho -sc
+//    mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -perf -asm -pc ho  -sc
+//    mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -std  -asm -pc ho  -sc
+//    mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -perf -mf  -pc lor
+//    mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -std  -asm -pc ho  -sc
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem
@@ -67,20 +68,26 @@ int main(int argc, char *argv[])
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+   ParTimer timer(MPI_COMM_WORLD);
 
    // 2. Parse command-line options.
    const char *mesh_file = "../../data/fichera.mesh";
+   int el_reord_type = 0;
    int order = sol_p;
    const char *basis_type = "G"; // Gauss-Lobatto
    bool static_cond = false;
    const char *pc = "lor";
    bool perf = true;
    bool matrix_free = true;
+   int max_iter = 500;
    bool visualization = 1;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
                   "Mesh file to use.");
+   args.AddOption(&el_reord_type, "-er", "--element-reorder-type",
+                  "How to reorder the mesh elements:\n\t0 - no reordering, "
+                  "1 - Cuthill-McKee, 2 - Metis, 3 - Gecko.");
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree) or -1 for"
                   " isoparametric space.");
@@ -96,6 +103,8 @@ int main(int argc, char *argv[])
                   "ho - high-order (assembled) AMG, none.");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&max_iter, "-mi", "--max-iter",
+                  "Maximum number of iterations.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -162,8 +171,8 @@ int main(int argc, char *argv[])
       {
          if (myid == 0)
          {
-            cout << "The given mesh does not match the optimized 'geom' parameter.\n"
-                 << "Recompile with suitable 'geom' value." << endl;
+            cout << "The given mesh does not match the optimized 'geom' "
+                 << "parameter.\nRecompile with suitable 'geom' value." << endl;
          }
          delete mesh;
          MPI_Finalize();
@@ -221,6 +230,36 @@ int main(int argc, char *argv[])
                   " the LOR preconditioner yet");
    }
 
+   // Perform local reordering of the mesh elements.
+   Array<int> ordering;
+   const char *el_reord_str[] = { "no reordering", "CM", "Metis", "Gecko" };
+   if (myid == 0)
+   {
+      cout << "Computing mesh element reordering ... ("
+           << el_reord_str[el_reord_type] << ") ..." << flush;
+   }
+   switch (el_reord_type)
+   {
+      case 1: pmesh->GetCMElementReordering(ordering); break;
+#ifdef MFEM_USE_METIS
+      case 2: pmesh->GetMetisElementReordering(ordering); break;
+#endif
+#ifdef MFEM_USE_GECKO
+      case 3: pmesh->GetGeckoElementReordering(ordering); break;
+#endif
+      default: MFEM_VERIFY(el_reord_type == 0, "invalid reordering type: "
+                           << el_reord_type);
+   }
+   if (el_reord_type != 0)
+   {
+      if (myid == 0)
+      {
+         cout << " done.\nApplying the mesh element reordering ..." << flush;
+      }
+      pmesh->ReorderElements(ordering);
+   }
+   if (myid == 0) { cout << " done." << endl; }
+
    // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use continuous Lagrange finite elements of the specified order. If
    //    order < 1, we instead use an isoparametric/isogeometric space.
@@ -323,8 +362,8 @@ int main(int argc, char *argv[])
    {
       cout << "Assembling the matrix ..." << flush;
    }
-   tic_toc.Clear();
-   tic_toc.Start();
+   timer.Clear();
+   timer.Start();
    // Pre-allocate sparsity assuming dense element matrices
    a->UsePrecomputedSparsity();
 
@@ -350,10 +389,10 @@ int main(int argc, char *argv[])
          a_hpc->AssembleBilinearForm(*a); // full matrix assembly
       }
    }
-   tic_toc.Stop();
+   timer.ParStop();
    if (myid == 0)
    {
-      cout << " done, " << tic_toc.RealTime() << "s." << endl;
+      cout << " done, time: " << timer << "." << endl;
    }
 
    // 14. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
@@ -387,8 +426,8 @@ int main(int argc, char *argv[])
    {
       cout << "Assembling the preconditioning matrix ..." << flush;
    }
-   tic_toc.Clear();
-   tic_toc.Start();
+   timer.Clear();
+   timer.Start();
 
    HypreParMatrix A_pc;
    if (pc_choice == LOR)
@@ -412,18 +451,18 @@ int main(int argc, char *argv[])
          a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
       }
    }
-   tic_toc.Stop();
+   timer.ParStop();
    if (myid == 0)
    {
-      cout << " done, " << tic_toc.RealTime() << "s." << endl;
+      cout << " done, time: " << timer << "." << endl;
    }
 
    // Solve with CG or PCG, depending if the matrix A_pc is available
    CGSolver *pcg;
    pcg = new CGSolver(MPI_COMM_WORLD);
    pcg->SetRelTol(1e-6);
-   pcg->SetMaxIter(500);
-   pcg->SetPrintLevel(1);
+   pcg->SetMaxIter(max_iter);
+   pcg->SetPrintLevel(3);
 
    HypreSolver *amg = NULL;
 
@@ -434,20 +473,21 @@ int main(int argc, char *argv[])
       pcg->SetPreconditioner(*amg);
    }
 
-   tic_toc.Clear();
-   tic_toc.Start();
+   timer.Clear();
+   timer.Start();
 
    pcg->Mult(B, X);
 
-   tic_toc.Stop();
+   timer.ParStop();
    delete amg;
 
    if (myid == 0)
    {
+      cout << "Total solve time: " << timer << ".\n";
       // Note: In the pcg algorithm, the number of operator Mult() calls is
       //       N_iter and the number of preconditioner Mult() calls is N_iter+1.
-      cout << "Time per CG step: "
-           << tic_toc.RealTime() / pcg->GetNumIterations() << "s." << endl;
+      cout << "Time per CG step: max: "
+           << timer.RealTimeMax() / pcg->GetNumIterations() << "s." << endl;
    }
 
    // 15. Recover the parallel grid function corresponding to X. This is the