diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt
index 44c62102e..83f928d62 100644
--- a/apf/CMakeLists.txt
+++ b/apf/CMakeLists.txt
@@ -6,6 +6,7 @@ endif()
 # Package sources
 set(SOURCES
   apf.cc
+  apfAggregateNumbering.cc
   apfCavityOp.cc
   apfElement.cc
   apfField.cc
@@ -51,6 +52,7 @@ set(SOURCES
 # Package headers
 set(HEADERS
   apf.h
+  apfAggregateNumbering.h
   apfMesh.h
   apfMesh2.h
   apfMatrix.h
diff --git a/apf/apf.cc b/apf/apf.cc
index 231d0f2fd..c9f41175f 100644
--- a/apf/apf.cc
+++ b/apf/apf.cc
@@ -453,10 +453,6 @@ bool isFrozen(Field* f)
   return f->getData()->isFrozen();
 }
 
-Function::~Function()
-{
-}
-
 Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s,
     Function* f)
 {
diff --git a/apf/apf.h b/apf/apf.h
index 44ef417db..d02c7b1bc 100644
--- a/apf/apf.h
+++ b/apf/apf.h
@@ -654,19 +654,22 @@ double* getArrayData(Field* f);
 /** \brief Initialize all nodal values with all-zero components */
 void zeroField(Field* f);
 
-/** \brief User-defined Analytic Function. */
-struct Function
+/** \brief User-defined Analytic Function for arbitrary scalar type. */
+template <class T>
+struct FunctionBase
 {
   /** \brief Possible user-defined cleanup */
-  virtual ~Function();
+  virtual ~FunctionBase() {}
   /** \brief The user's analytic function call.
-    \details For simplicity, this
-    currently only supports one node per entity.
-    \param e the entity on which the node is
-    \param result the field component values for that node
-    */
-  virtual void eval(MeshEntity* e, double* result) = 0;
+      \details For simplicity, this
+      currently only supports one node per entity.
+      \param e the entity on which the node is
+      \param result the field component values for that node
+  */
+  virtual void eval(MeshEntity * e, T * result) = 0;
 };
+typedef FunctionBase<double> Function;
+
 
 /* \brief Create a Field from a user's analytic function.
  * \details This field will use no memory, has values on all
diff --git a/apf/apfAggregateNumbering.cc b/apf/apfAggregateNumbering.cc
new file mode 100644
index 000000000..17099d815
--- /dev/null
+++ b/apf/apfAggregateNumbering.cc
@@ -0,0 +1,334 @@
+#include "apfAggregateNumbering.h"
+#include "apfAggregateNumberingClass.h"
+#include "apfField.h"
+#include "apfNumberingClass.h"
+#include "apfTagData.h"
+#include "apfUserData.h"
+#include <PCU.h>
+#include <pcu_util.h>
+#include <set>
+namespace apf
+{
+  // explicit instantiations
+  template class AggregateNumberingOf<int>;
+  template class AggregateNumberingOf<long>;
+  // this assumes that the ranks returned by a sharing
+  //  do not change when the pcu comm changes
+  //  if the sharing can handle the pcu comm changing
+  //  then this isn't needed
+  class StaticSharing : public NormalSharing
+  {
+  protected:
+    Sharing * shr;
+    // original comm
+    MPI_Comm old_cm;
+    int old_rnk;
+    // new comm
+    MPI_Comm new_cm;
+  public:
+    StaticSharing(Mesh * m,
+                  Sharing * s,
+                  MPI_Comm ncm,
+                  MPI_Comm pcm = MPI_COMM_NULL)
+      : NormalSharing(m)
+      , shr(s)
+      , old_cm(pcm)
+      , old_rnk(-1)
+      , new_cm(ncm)
+    {
+      if(old_cm == MPI_COMM_NULL)
+        old_cm = PCU_Get_Comm();
+      new_cm = PCU_Get_Comm();
+      old_rnk = PCU_Comm_Self();
+    }
+    virtual bool isOwned(MeshEntity * e)
+    {
+      int old_ownr = shr->getOwner(e);
+      if(old_ownr == old_rnk)
+        return true;
+      return false;
+    }
+    virtual int getOwner(MeshEntity * e)
+    {
+      int old_ownr = shr->getOwner(e);
+      return PCU_Foreign_To_Local(old_ownr,old_cm);
+    }
+    virtual void getCopies(MeshEntity * e, CopyArray & copies)
+    {
+      shr->getCopies(e,copies);
+      for(int ii = 0; ii < copies.getSize(); ++ii)
+        copies[ii].peer = PCU_Foreign_To_Local(copies[ii].peer,old_cm);
+    }
+    virtual bool isShared(MeshEntity * e)
+    {
+      return shr->isShared(e);
+    }
+  };
+  template <class T>
+  struct AggDataOfFunc : public FunctionBase<T>
+  {
+    AggDataOfFunc(AggregateNumberingOf<T> * n,
+                  TagDataOf<T> * i)
+      : num(n)
+      , nd_ids(i)
+    { }
+    virtual void eval(MeshEntity * e, T * dat)
+    {
+      int nds = num->getShape()->countNodesOn(num->getMesh()->getType(e));
+      int blks_per_nd = num->countBlocks();
+      int dofs_per_blk = num->blockSize();
+      int cmps = blks_per_nd * dofs_per_blk;
+      int ownr = num->getSharing()->getOwner(e);
+      int strd = num->getStride(ownr);
+      int frst = num->getFirstDof(ownr);
+      for(int nd = 0; nd < nds; ++nd)
+      {
+        T nd_id = -1;
+        nd_ids->get(e,&nd_id);
+        for(int blk = 0; blk < blks_per_nd; ++blk)
+          for(int dof = 0; dof < dofs_per_blk; ++dof)
+            dat[nd * cmps + blk * dofs_per_blk + dof] =
+              frst + (nd_id * dofs_per_blk) + (blk * strd) + dof;
+      }
+    }
+    AggregateNumberingOf<T> * num;
+    TagDataOf<T> * nd_ids;
+  };
+  template <class T>
+  T exscanT(T v);
+  template <>
+  int exscanT(int v)
+  {
+    return PCU_Exscan_Int(v);
+  }
+  template <>
+  long exscanT(long v)
+  {
+    return PCU_Exscan_Long(v);
+  }
+  template <class T>
+  void allgatherT(T i, T* o);
+  template <>
+  void allgatherT(int i, int* o)
+  {
+    PCU_Allgather_Int(i,o);
+  }
+  template <>
+  void allgatherT(long i, long* o)
+  {
+    PCU_Allgather_Long(i,o);
+  }
+  template <class T>
+  void AggregateNumberingOf<T>::init(const char * n,
+                                     Mesh * m,
+                                     Sharing * sh,
+                                     FieldShape * s,
+                                     int blks,
+                                     int bs)
+  {
+    // trick AggFieldDataOf into only allocating one int per node
+    shr = sh;
+    //FieldBase::init(n,m,s,new AggDataOf<T>(this));
+    NumberingOf<T>::components = 1;
+    FieldBase::shape = s;
+    FieldBase::mesh = m;
+    FieldBase::name = n;
+    nd_ids = new TagDataOf<T>;
+    nd_ids->init(this);
+    FieldBase::init(n,m,s,new UserDataBase<T>(new AggDataOfFunc<T>(this,static_cast<TagDataOf<T>*>(nd_ids))));
+    NumberingOf<T>::components = blks * bs;
+    blks_per_nd = blks;
+    dofs_per_blk = bs;
+    int dim = m->getDimension();
+    for(int dd = 0; dd < dim; ++dd)
+      lcl_strd += countOwnedNodes(m,dd,s,sh);
+    int sz = PCU_Comm_Peers();
+    slf = PCU_Comm_Self();
+    // do we bother setting strds and frsts since they
+    //  are reset in globalize which must be called after
+    //  init for the numbering to be valid
+    strds.resize(static_cast<unsigned>(sz));
+    // all peers in the aggregation scope necessarily
+    //  have the same stride
+    int strd = PCU_Add_Int(lcl_strd);
+    for(int ii = 0; ii < sz; ++ii)
+      strds[ii] = strd;
+    // the first dof id of the first local node is the # of
+    //  nodes before the first local node times the # of
+    //  dofs in the first block on a node
+    T dof_id_offset = lcl_strd * dofs_per_blk;
+    T frst = exscanT<T>(dof_id_offset);
+    frsts.resize(static_cast<unsigned>(sz));
+    allgatherT(frst,&frsts[0]);
+    T nd_id = 0;
+    for(int dd = 0; dd < dim; ++dd)
+    {
+      MeshIterator * it = m->begin(dd);
+      MeshEntity * e = NULL;
+      while((e = m->iterate(it)))
+      {
+        int nds = FieldBase::shape->countNodesOn(m->getType(e));
+        // only really works if nds == 1 or 0
+        //  due to the above comment ^
+        for(int nd = 0; nd < nds; ++nd)
+        {
+          if(shr->isOwned(e))
+          {
+            nd_ids->set(e,&nd_id);
+            ++nd_id;
+          }
+        }
+      }
+      m->end(it);
+    }
+  }
+  template <class T>
+  void AggregateNumberingOf<T>::init(Field * fld,
+                                     Sharing * sh,
+                                     int blks,
+                                     int bs)
+  {
+    PCU_ALWAYS_ASSERT(fld->countComponents() == blks * bs);
+    std::string nm = fld->getName();
+    nm += "_num";
+    init(nm.c_str(),
+         fld->getMesh(),
+         sh,
+         fld->getShape(),
+         blks,
+         bs);
+  }
+  template <class T>
+  void AggregateNumberingOf<T>::globalize()
+  {
+    synchronizeFieldData(nd_ids,shr,false);
+    // any rank with frst == 0 was the zero rank
+    //  in the aggregation scope under which init
+    //  was called, the offsets from these
+    //  zero-ranked processes are what we need to
+    //  sum since they include in their strides
+    //  all the nodes from the local aggregation scope
+    int strd = getScopeStride();
+    int frst = getLocalFirstDof();
+    int scp_offset = (int)(!(bool)frst) * (strd * dofs_per_blk * blks_per_nd);
+    int lcl_offset = PCU_Exscan_Int(scp_offset);
+    lcl_offset += scp_offset; // make the excan a scan
+    // remove the contribution from the local
+    //  aggregation scope, because the scope can be
+    //  any comm and only rank 0 in each scope
+    //  contributes to the above scan
+    lcl_offset -= (strd * dofs_per_blk * blks_per_nd);
+    frst += lcl_offset;
+    // strds and frsts needs to be resized/updated
+    //  since we can have adjacent parts that were
+    //  not in the aggregation scope under which
+    //  init was called
+    int sz = PCU_Comm_Peers();
+    frsts.setSize(sz);
+    allgatherT<T>(frst,&frsts[0]);
+    // each aggregation scope has its own stride
+    //  we store per-rank since that is vastly
+    //  less complicated than keeping track
+    //  of the initial scopes
+    strds.setSize(sz);
+    allgatherT<T>(strd,&strds[0]);
+  }
+  template <class T>
+  void AggregateNumberingOf<T>::getAll(MeshEntity * e, T * dat)
+  {
+    reinterpret_cast<FieldDataOf<T>*>(FieldBase::data)->get(e,dat);
+  }
+  template <class T>
+  T AggregateNumberingOf<T>::get(MeshEntity * e, int nd, int cmp)
+  {
+    int cmps = blks_per_nd * dofs_per_blk;
+    int nds = FieldBase::shape->countNodesOn(FieldBase::mesh->getType(e));
+    NewArray<T> data(nds*cmps);
+    getAll(e,&data[0]);
+    return data[nd*cmps + cmp];
+  }
+  AggNumbering * createAggNumbering(Field * f,
+                                    int blocks,
+                                    int dofs_per_block,
+                                    MPI_Comm cm,
+                                    Sharing * share)
+  {
+    bool dlt = false;
+    if(!share)
+    {
+      share = getSharing(getMesh(f));
+      dlt = true;
+    }
+    AggNumbering * n = new AggNumbering;
+    MPI_Comm pcu_cm = PCU_Get_Comm();
+    bool swtch = pcu_cm != cm;
+    Sharing * init_share = share;
+    if(swtch)
+    {
+      init_share = new StaticSharing(getMesh(f),share,cm);
+      PCU_Switch_Comm(cm);
+    }
+    n->init(f,init_share,blocks,dofs_per_block);
+    if(swtch)
+    {
+      PCU_Switch_Comm(pcu_cm);
+      delete init_share;
+    }
+    // during init the static sharing may be used which
+    //  only provides binary ownership information, not
+    //  the owner rank of unowned entities, so we replace
+    //  it with the passed-in sharing here before globalization
+    //  which requires the rank information
+    n->setSharing(share);
+    n->globalize();
+    f->getMesh()->addNumbering(n);
+    return n;
+  }
+  AggNumbering * createAggNumbering(Mesh * m,
+                                    const char * name,
+                                    FieldShape * shape,
+                                    int blocks,
+                                    int dofs_per_block,
+                                    MPI_Comm cm,
+                                    Sharing * share)
+  {
+    bool dlt = false;
+    if(!share)
+    {
+      share = getSharing(m);
+      dlt = true;
+    }
+    AggNumbering * n = new AggNumbering;
+    MPI_Comm pcu_cm = PCU_Get_Comm();
+    bool swtch = pcu_cm != cm;
+    Sharing * init_share = share;
+    if(swtch)
+    {
+      init_share = new StaticSharing(m,share,cm);
+      PCU_Switch_Comm(cm);
+    }
+    n->init(name,m,init_share,shape,blocks,dofs_per_block);
+    if(swtch)
+    {
+      PCU_Switch_Comm(pcu_cm);
+      delete init_share;
+    }
+    n->setSharing(share);
+    n->globalize();
+    m->addNumbering(n);
+    return n;
+  }
+  /* Public API */
+  Numbering * getNumbering(AggNumbering * n)
+  {
+    return static_cast<Numbering*>(n);
+  }
+  int countBlocks(AggNumbering * n)
+  {
+    return n->blockSize();
+  }
+  int countDOFsPerBlock(AggNumbering * n)
+  {
+    return n->blockSize();
+  }
+}
diff --git a/apf/apfAggregateNumbering.h b/apf/apfAggregateNumbering.h
new file mode 100644
index 000000000..ff2715fb5
--- /dev/null
+++ b/apf/apfAggregateNumbering.h
@@ -0,0 +1,42 @@
+#ifndef APF_AGGREGATE_NUMBERING_H_
+#define APF_AGGREGATE_NUMBERING_H_
+#include "apfNew.h"
+#include "apfMesh.h"
+#include "apfNumbering.h"
+#include <PCU.h>
+namespace apf
+{
+  // class declaration
+  template <class T>
+  class AggregateNumberingOf;
+  typedef AggregateNumberingOf<int> AggNumbering;
+  typedef AggregateNumberingOf<long> GlobalAggNumbering;
+  // numbering creation functions
+  AggNumbering * createAggNumbering(Field * f,
+                                    int blocks,
+                                    int dofs_per_block,
+                                    MPI_Comm cm,
+                                    Sharing * share = NULL);
+  AggNumbering * createAggNumbering(Mesh * m,
+                                    const char * name,
+                                    FieldShape * shape,
+                                    int blocks,
+                                    int dofs_per_block,
+                                    MPI_Comm cm,
+                                    Sharing * share = NULL);
+  /*
+   * @brief Since the class implementation is not public
+   *        we need an API to up-cast AggNumbering to
+   *        Numbering.
+   */
+  Numbering * getNumbering(AggNumbering * n);
+  int countBlocks(AggNumbering * n);
+  int countDOFsPerBlock(AggNumbering * n);
+}
+#endif
+/******************************************************************************
+  Copyright 2018 Scientific Computation Research Center,
+      Rensselaer Polytechnic Institute. All rights reserved.
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+*******************************************************************************/
diff --git a/apf/apfAggregateNumberingClass.h b/apf/apfAggregateNumberingClass.h
new file mode 100644
index 000000000..ae6e4c1dd
--- /dev/null
+++ b/apf/apfAggregateNumberingClass.h
@@ -0,0 +1,53 @@
+#ifndef APF_AGGREGATE_NUMBERING_CLASS_H_
+#define APF_AGGREGATE_NUMBERING_CLASS_H_
+#include "apfFieldData.h"
+#include "apfNumberingClass.h"
+namespace apf
+{
+  template <class T>
+  class AggregateNumberingOf : public NumberingOf<T>
+  {
+  public:
+    AggregateNumberingOf()
+      : NumberingOf<T>()
+      , nd_ids(NULL)
+      , shr(NULL)
+      , blks_per_nd(0)
+      , dofs_per_blk(0)
+      , slf(-1)
+      , strds()
+      , frsts()
+      , lcl_strd(0)
+    { }
+    void init(const char * n,
+              Mesh * m,
+              Sharing * sh,
+              FieldShape * s,
+              int blks,
+              int bs);
+    void init(Field * f, Sharing * sh, int blks, int bs);
+    void globalize();
+    virtual void getAll(MeshEntity * e, T * dat);
+    virtual T get(MeshEntity * e, int nd, int cmp);
+    virtual void set(MeshEntity*,int,int,T) {}
+    Sharing * getSharing() const { return shr; }
+    void setSharing(Sharing * s) { shr = s; }
+    int countBlocks() const { return blks_per_nd; }
+    int blockSize() const { return dofs_per_blk; }
+    T getLocalStride() const { return lcl_strd; }
+    T getScopeStride() const { return strds[slf]; }
+    T getLocalFirstDof() const { return frsts[slf]; }
+    T getStride(int peer) const { return strds[peer]; };
+    T getFirstDof(int peer) const { return frsts[peer]; }
+  private:
+    FieldDataOf<T> * nd_ids;
+    Sharing * shr;
+    int blks_per_nd;
+    int dofs_per_blk;
+    int slf;
+    DynamicArray<T> strds;
+    DynamicArray<T> frsts;
+    int lcl_strd;
+  };
+}
+#endif
diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc
index 1b0c3f088..dc5f5bd38 100644
--- a/apf/apfMesh.cc
+++ b/apf/apfMesh.cc
@@ -802,6 +802,31 @@ int countOwned(Mesh* m, int dim, Sharing * shr)
   return n;
 }
 
+int countOwnedNodes(Mesh* m, int dim, FieldShape * shp, Sharing * shr)
+{
+  if(!shp->hasNodesIn(dim))
+    return 0;
+  bool dlt = false;
+  if(shr == NULL)
+  {
+    shr = getSharing(m);
+    dlt = true;
+  }
+  MeshIterator* it = m->begin(dim);
+  MeshEntity* e;
+  int n = 0;
+  while ((e = m->iterate(it)))
+    if (shr->isOwned(e))
+    {
+      int tp = m->getType(e);
+      int nds = shp->countNodesOn(tp);
+      n += nds;
+    }
+  m->end(it);
+  if(dlt) delete shr;
+  return n;
+}
+
 void printTypes(Mesh* m)
 {
   const int dim = m->getDimension();
diff --git a/apf/apfMesh.h b/apf/apfMesh.h
index 29c3631b4..b7158cae5 100644
--- a/apf/apfMesh.h
+++ b/apf/apfMesh.h
@@ -23,6 +23,7 @@ namespace apf {
 
 class FieldShape;
 class Field;
+
 template <class T>
 class NumberingOf;
 /** \brief Numbering is meant to be a 32-bit local numbering */
@@ -568,6 +569,10 @@ int countEntitiesOn(Mesh* m, ModelEntity* me, int dim);
            the default sharing is used if none is provided */
 int countOwned(Mesh* m, int dim, Sharing * shr = NULL);
 
+/** \brief count the number of nodes on owned entities of dimension (dim) under FieldShape shape using sharing shr
+    the default sharing is used if none is provided */
+int countOwnedNodes(Mesh * m, int dim, FieldShape * shape, Sharing * shr = NULL);
+
 /** \brief print global mesh entity counts per dimension */
 void printStats(Mesh* m);
 
diff --git a/apf/apfNumbering.cc b/apf/apfNumbering.cc
index 8b05d6c3a..793c30410 100644
--- a/apf/apfNumbering.cc
+++ b/apf/apfNumbering.cc
@@ -20,6 +20,7 @@ enum { FIXED = -2, FREE_BUT_NOT_NUMBERED = -1 };
 
 template <class T>
 NumberingOf<T>::NumberingOf()
+  : FieldBase()
 {
   field = 0;
   components = 0;
diff --git a/apf/apfNumbering.h b/apf/apfNumbering.h
index f4f339e4d..f634f8dd7 100644
--- a/apf/apfNumbering.h
+++ b/apf/apfNumbering.h
@@ -68,7 +68,7 @@ int getNumber(Numbering* n, MeshEntity* e, int node, int component);
 
 /** \brief get the field being numbered
   */
-Field* getField(Numbering* n);
+Field * getField(Numbering* n);
 
 /** \brief get the FieldShape used by a Numbering */
 FieldShape* getShape(Numbering* n);
diff --git a/apf/apfNumberingClass.h b/apf/apfNumberingClass.h
index 90b150763..6fb2ce475 100644
--- a/apf/apfNumberingClass.h
+++ b/apf/apfNumberingClass.h
@@ -19,18 +19,18 @@ class NumberingOf : public FieldBase
     NumberingOf();
     virtual int countComponents() const;
     virtual int getScalarType();
-    void init(const char* n,
-              Mesh* m,
-              FieldShape* s,
+    void init(const char * n,
+              Mesh * m,
+              FieldShape * s,
               int c);
     void init(Field* f);
-    Field* getField();
-    FieldDataOf<T>* getData();
-    void getAll(MeshEntity* e, T* dat);
-    T get(MeshEntity* e, int node, int component);
-    void set(MeshEntity* e, int node, int component, T value);
-  private:
-    Field* field;
+    Field * getField();
+    FieldDataOf<T> * getData();
+    virtual void getAll(MeshEntity* e, T* dat);
+    virtual T get(MeshEntity* e, int node, int component);
+    virtual void set(MeshEntity* e, int node, int component, T value);
+  protected:
+    Field * field;
     int components;
 };
 
diff --git a/apf/apfUserData.cc b/apf/apfUserData.cc
index 15ee6a310..7f9c4872b 100644
--- a/apf/apfUserData.cc
+++ b/apf/apfUserData.cc
@@ -9,46 +9,9 @@
 
 namespace apf {
 
-UserData::UserData(Function* f)
-{
-  function = f;
-}
-
-void UserData::init(FieldBase* f)
-{
-  field = f;
-}
-
-bool UserData::hasEntity(MeshEntity* e)
-{
-  return field->getShape()->countNodesOn(field->getMesh()->getType(e)) > 0;
-}
-
-void UserData::removeEntity(MeshEntity*)
-{
-}
-
-void UserData::get(MeshEntity* e, double* data)
-{
-  function->eval(e, data);
-}
-
-void UserData::set(MeshEntity*, double const*)
-{
-}
-
-bool UserData::isFrozen()
-{
-  return false;
-}
-
-FieldData* UserData::clone()
-{
-  FieldData* newData = new UserData(function);
-  newData->init(field);
-  copyFieldData(static_cast<FieldDataOf<double>*>(newData),
-                static_cast<FieldDataOf<double>*>(field->getData()));
-  return newData;
-}
+template <> class UserDataBase<double>;
+template <> class UserDataBase<int>;
+template <> class UserDataBase<long>;
+template <> class UserDataBase<size_t>;
 
 }
diff --git a/apf/apfUserData.h b/apf/apfUserData.h
index 96980727e..041796437 100644
--- a/apf/apfUserData.h
+++ b/apf/apfUserData.h
@@ -10,26 +10,46 @@
 
 #include "apfFieldData.h"
 
-namespace apf {
-
-struct UserData : public FieldDataOf<double>
+namespace apf
+{
+template <class T>
+struct UserDataBase : public FieldDataOf<T>
 {
-  UserData(Function* f);
-  void init(FieldBase* f);
-  bool hasEntity(MeshEntity* e);
-  void removeEntity(MeshEntity* e);
-  void get(MeshEntity* e, double* data);
-  void set(MeshEntity* e, double const* data);
-  bool isFrozen();
-  virtual FieldData* clone();
-  // using const * const gives an error on gcc/7.3.0 because the return is an
-  // r-value which cannot be modified anyways
-  Function const* getFunction() const { return function; }
-  void setFunction(Function* func) { function = func; }
-  private:
-  Function* function;
+  UserDataBase(FunctionBase<T> * f)
+    : function(f)
+  { }
+  void init(FieldBase * f)
+  {
+    FieldDataOf<T>::field = f;
+  }
+  bool hasEntity(MeshEntity * e)
+  {
+    return FieldDataOf<T>::field->getShape()->countNodesOn(FieldDataOf<T>::field->getMesh()->getType(e)) > 0;
+  }
+  void removeEntity(MeshEntity * e) {}
+  void get(MeshEntity * e, T * data)
+  {
+    function->eval(e,data);
+  }
+  void set(MeshEntity * e, T const * data) {}
+  bool isFrozen() { return false; }
+  virtual FieldData* clone()
+  {
+    FieldData* newData = new UserDataBase<T>(function);
+    newData->init(FieldDataOf<T>::field);
+    copyFieldData(static_cast<FieldDataOf<T>*>(newData),
+                  static_cast<FieldDataOf<T>*>(FieldDataOf<T>::field->getData()));
+    return newData;
+  }
+  FunctionBase<T> const* getFunction() const { return function; }
+  void setFunction(FunctionBase<T>* func) { function = func; }
+private:
+  FunctionBase<T> * function;
 };
-
+typedef UserDataBase<double> UserData;
+typedef UserDataBase<int> UserDataInt;
+typedef UserDataBase<long> UserDataLong;
+typedef UserDataBase<size_t> UserDataSizeT;
 }
 
 #endif
diff --git a/pcu/PCU.h b/pcu/PCU.h
index 67b2ecccd..eaf6db268 100644
--- a/pcu/PCU.h
+++ b/pcu/PCU.h
@@ -80,6 +80,16 @@ int PCU_Max_Int(int x);
 int PCU_Or(int c);
 int PCU_And(int c);
 
+/*collective allgathers*/
+void PCU_Allgather_Doubles(double* i, size_t ni, double* o);
+void PCU_Allgather_Double(double i, double* o);
+void PCU_Allgather_Ints(int* i, size_t n, int* o);
+void PCU_Allgather_Int(int i, int* o);
+void PCU_Allgather_Longs(long* i, size_t n, long* o);
+void PCU_Allgather_Long(long i, long* o);
+void PCU_Allgather_SizeTs(size_t* i, size_t n, size_t* o);
+void PCU_Allgather_SizeT(size_t i, size_t* o);
+
 /*process-level self/peers (mpi wrappers)*/
 int PCU_Proc_Self(void);
 int PCU_Proc_Peers(void);
@@ -119,6 +129,8 @@ int PCU_Comm_Start(PCU_Method method);
 /*special MPI_Comm replacement API*/
 void PCU_Switch_Comm(MPI_Comm new_comm);
 MPI_Comm PCU_Get_Comm(void);
+int PCU_Local_To_Foreign(int local_rank, MPI_Comm foreign_comm);
+int PCU_Foreign_To_Local(int foreign_rank, MPI_Comm foreign_comm);
 
 /*stack trace helpers using GNU/Linux*/
 void PCU_Protect(void);
diff --git a/pcu/pcu.c b/pcu/pcu.c
index 5244f65e2..9755f444a 100644
--- a/pcu/pcu.c
+++ b/pcu/pcu.c
@@ -484,6 +484,82 @@ int PCU_And(int c)
   return PCU_Min_Int(c);
 }
 
+/** \brief Performs a parallel gather followed by a broadcast
+  */
+void PCU_Allgather_Doubles(double* i, size_t ni, double* o)
+{
+  if(global_state == uninit)
+    reel_fail("Allgather_Doubles called before Comm_Init");
+  int rnk = pcu_pmpi_rank();
+  int sz = pcu_pmpi_size();
+  memset(o,0,sz*ni*sizeof(double));
+  memcpy(o+(rnk*ni),i,ni*sizeof(double));
+  pcu_reduce(&(get_msg()->coll),pcu_add_doubles,o,sz*ni*sizeof(double));
+  pcu_bcast(&(get_msg()->coll),o,sz*ni*sizeof(double));
+}
+
+void PCU_Allgather_Double(double i, double* o)
+{
+  PCU_Allgather_Doubles(&i,1,o);
+}
+
+/** \brief Performs a parallel gather followed by a broadcast
+  */
+void PCU_Allgather_Ints(int* i, size_t ni, int* o)
+{
+  if(global_state == uninit)
+    reel_fail("Allgather_Ints called before Comm_Init");
+  int rnk = pcu_pmpi_rank();
+  int sz = pcu_pmpi_size();
+  memset(o,0,sz*ni*sizeof(int));
+  memcpy(o+(rnk*ni),i,ni*sizeof(int));
+  pcu_reduce(&(get_msg()->coll),pcu_add_ints,o,sz*ni*sizeof(int));
+  pcu_bcast(&(get_msg()->coll),o,sz*ni*sizeof(int));
+}
+
+void PCU_Allgather_Int(int i, int* o)
+{
+  PCU_Allgather_Ints(&i,1,o);
+}
+
+/** \brief Performs a parallel gather followed by a broadcast
+  */
+void PCU_Allgather_Longs(long* i, size_t ni, long* o)
+{
+  if(global_state == uninit)
+    reel_fail("Allgather_Longs called before Comm_Init");
+  int rnk = pcu_pmpi_rank();
+  int sz = pcu_pmpi_size();
+  memset(o,0,sz*ni*sizeof(long));
+  memcpy(o+(rnk*ni),i,ni*sizeof(long));
+  pcu_reduce(&(get_msg()->coll),pcu_add_longs,o,sz*ni*sizeof(long));
+  pcu_bcast(&(get_msg()->coll),o,sz*ni*sizeof(long));
+}
+
+void PCU_Allgather_Long(long i, long* o)
+{
+  PCU_Allgather_Longs(&i,1,o);
+}
+
+/** \brief Performs a parallel gather followed by a broadcast
+  */
+void PCU_Allgather_SizeTs(size_t* i, size_t ni, size_t *o)
+{
+  if(global_state == uninit)
+    reel_fail("Allgather_SizeTs called before Comm_Init");
+  int rnk = pcu_pmpi_rank();
+  int sz = pcu_pmpi_size();
+  memset(o,0,sz*ni*sizeof(size_t));
+  memcpy(o+(rnk*ni),i,ni*sizeof(size_t));
+  pcu_reduce(&(get_msg()->coll),pcu_add_sizets,o,sz*ni*sizeof(size_t));
+  pcu_bcast(&(get_msg()->coll),o,sz*ni*sizeof(size_t));
+}
+
+void PCU_Allgather_SizeT(size_t i, size_t* o)
+{
+  PCU_Allgather_SizeTs(&i,1,o);
+}
+
 /** \brief Returns the unique rank of the calling process.
  */
 int PCU_Proc_Self(void)
@@ -739,6 +815,20 @@ MPI_Comm PCU_Get_Comm(void)
   return pcu_pmpi_comm();
 }
 
+int PCU_Local_To_Foreign(int local_rank, MPI_Comm foreign_comm)
+{
+  if (global_state == uninit)
+    reel_fail("Local_To_Foreign called before Comm_Init");
+  return pcu_pmpi_lcl_to_frn(local_rank,foreign_comm);
+}
+
+int PCU_Foreign_To_Local(int foreign_rank, MPI_Comm foreign_comm)
+{
+  if (global_state == uninit)
+    reel_fail("Foreign_To_Local called before Comm_Init");
+  return pcu_pmpi_frn_to_lcl(foreign_rank,foreign_comm);
+}
+
 /** \brief Return the time in seconds since some time in the past
  */
 double PCU_Time(void)
diff --git a/pcu/pcu_coll.h b/pcu/pcu_coll.h
index 83ca77daa..e2ab9eae0 100644
--- a/pcu/pcu_coll.h
+++ b/pcu/pcu_coll.h
@@ -27,6 +27,7 @@
 
 typedef void pcu_merge(void* local, void* incoming, size_t size);
 void pcu_merge_assign(void* local, void* incoming, size_t size);
+void pcu_merge_gather(void* local, void* incoming, size_t size);
 void pcu_add_doubles(void* local, void* incoming, size_t size);
 void pcu_max_doubles(void* local, void* incoming, size_t size);
 void pcu_min_doubles(void* local, void* incoming, size_t size);
diff --git a/pcu/pcu_pmpi.c b/pcu/pcu_pmpi.c
index 814828a7e..fbb92f565 100644
--- a/pcu/pcu_pmpi.c
+++ b/pcu/pcu_pmpi.c
@@ -118,3 +118,30 @@ MPI_Comm pcu_pmpi_comm(void)
   return original_comm;
 }
 
+int pcu_pmpi_lcl_to_frn(int lcl_rnk, MPI_Comm frn_cm)
+{
+  MPI_Group lcl_grp;
+  MPI_Group frn_grp;
+  int frn_rnk = -1;
+  MPI_Comm_group(original_comm,&lcl_grp);
+  if(frn_cm != MPI_COMM_NULL)
+  {
+    MPI_Comm_group(frn_cm,&frn_grp);
+    MPI_Group_translate_ranks(lcl_grp, 1, &lcl_rnk, frn_grp, &frn_rnk);
+  }
+  return frn_rnk;
+}
+
+int pcu_pmpi_frn_to_lcl(int frn_rnk, MPI_Comm frn_cm)
+{
+  MPI_Group lcl_grp;
+  MPI_Group frn_grp;
+  int lcl_rnk = -1;
+  MPI_Comm_group(original_comm,&lcl_grp);
+  if(frn_cm != MPI_COMM_NULL)
+  {
+    MPI_Comm_group(frn_cm,&frn_grp);
+    MPI_Group_translate_ranks(frn_grp,1,&frn_rnk,lcl_grp,&lcl_rnk);
+  }
+  return lcl_rnk;
+}
diff --git a/pcu/pcu_pmpi.h b/pcu/pcu_pmpi.h
index c6707d23e..d23f69d84 100644
--- a/pcu/pcu_pmpi.h
+++ b/pcu/pcu_pmpi.h
@@ -26,6 +26,8 @@ bool pcu_pmpi_done(pcu_message* m);
 
 void pcu_pmpi_switch(MPI_Comm new_comm);
 MPI_Comm pcu_pmpi_comm(void);
+int  pcu_pmpi_lcl_to_frn(int lcl, MPI_Comm frn_cm);
+int  pcu_pmpi_frn_to_lcl(int frn, MPI_Comm frn_cm);
 
 extern pcu_mpi pcu_pmpi;
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c3389e565..4ec3aff08 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -124,6 +124,7 @@ if(SIMMODSUITE_SimAdvMeshing_FOUND)
 endif()
 
 # Unit tests / functionality regression tests
+test_exe_func(pcu_test pcu.cc)
 test_exe_func(qr qr.cc)
 test_exe_func(eigen_test eigen_test.cc)
 test_exe_func(integrate integrate.cc)
@@ -159,6 +160,7 @@ test_exe_func(pyramidCodeMatch ../ma/pyramidCodeMatch.cc)
 test_exe_func(newdim newdim.cc)
 test_exe_func(construct construct.cc)
 test_exe_func(test_scaling test_scaling.cc)
+test_exe_func(aggregateNumbering aggregateNumbering.cc)
 test_exe_func(mixedNumbering mixedNumbering.cc)
 test_exe_func(test_verify test_verify.cc)
 test_exe_func(hierarchic hierarchic.cc)
diff --git a/test/aggregateNumbering.cc b/test/aggregateNumbering.cc
new file mode 100644
index 000000000..d63d4c458
--- /dev/null
+++ b/test/aggregateNumbering.cc
@@ -0,0 +1,50 @@
+#include <PCU.h>
+#include <apf.h>
+#include <apfMesh2.h>
+#include <apfAggregateNumbering.h>
+#include <apfMDS.h>
+#include <gmi_mesh.h>
+#include <gmi_sim.h>
+#include <pcu_util.h>
+
+void test_numbering(apf::Mesh * m)
+{
+  apf::Field * f1 = apf::createPackedField(m, "u1", 3);
+  apf::Field * f2 = apf::createPackedField(m, "u2", 3);
+  apf::Field * f3 = apf::createPackedField(m, "u3", 3);
+  apf::Field * f4 = apf::createPackedField(m, "u4", 3);
+  apf::Sharing * def = apf::getSharing(m);
+  // global scope natural order
+  apf::AggNumbering * num_gbl_nat = apf::createAggNumbering(f1,1,3,MPI_COMM_WORLD,def);
+  // global scope block order
+  apf::AggNumbering * num_gbl_blk = apf::createAggNumbering(f2,3,1,MPI_COMM_WORLD,def);
+  // local scope natural order
+  apf::AggNumbering * num_lcl_nat = apf::createAggNumbering(f3,1,3,MPI_COMM_SELF);
+  // local scope block order
+  apf::AggNumbering * num_lcl_blk = apf::createAggNumbering(f4,3,1,MPI_COMM_SELF);
+}
+
+void write_output(apf::Mesh * m, const char * fn)
+{
+  apf::writeVtkFiles(fn,m);
+}
+
+int main(int ac, char * av[])
+{
+  PCU_ALWAYS_ASSERT(ac == 4);
+  MPI_Init(&ac,&av);
+  PCU_Comm_Init();
+  gmi_sim_start();
+  gmi_register_sim();
+  gmi_register_mesh();
+  apf::Mesh2 * m = apf::loadMdsMesh(av[1],av[2]);
+  apf::reorderMdsMesh(m);
+  test_numbering(m);
+  write_output(m,av[3]);
+  m->destroyNative();
+  apf::destroyMesh(m);
+  gmi_sim_stop();
+  PCU_Comm_Free();
+  MPI_Finalize();
+  return 0;
+}
diff --git a/test/pcu.cc b/test/pcu.cc
new file mode 100644
index 000000000..324361287
--- /dev/null
+++ b/test/pcu.cc
@@ -0,0 +1,28 @@
+#include <PCU.h>
+#include <pcu_util.h>
+
+int main(int ac, char * av[])
+{
+  MPI_Init(&ac,&av);
+  int mpi_sz = 0;
+  MPI_Comm_size(MPI_COMM_WORLD,&mpi_sz);
+  int mpi_rnk = -1;
+  MPI_Comm_rank(MPI_COMM_WORLD,&mpi_rnk);
+  PCU_Comm_Init();
+  int sz = PCU_Comm_Peers();
+  PCU_ALWAYS_ASSERT(sz == mpi_sz);
+  int rnk = PCU_Comm_Self();
+  PCU_ALWAYS_ASSERT(rnk == mpi_rnk);
+  int * gather_bfr = new int[sz];
+  int * compar_bfr = new int[sz];
+  for(int rr = 0; rr < sz; ++rr)
+    compar_bfr[rr] = rr;
+  PCU_Allgather_Int(rnk,gather_bfr);
+  for(int rr = 0; rr < sz; ++rr)
+    PCU_ALWAYS_ASSERT(gather_bfr[rr] == compar_bfr[rr]);
+  delete [] compar_bfr;
+  delete [] gather_bfr;
+  PCU_Comm_Free();
+  MPI_Finalize();
+  return 0;
+}
diff --git a/test/testing.cmake b/test/testing.cmake
index e921e9824..5e21752f0 100644
--- a/test/testing.cmake
+++ b/test/testing.cmake
@@ -18,6 +18,8 @@ function(smoke_test TESTNAME PROCS EXE)
   SET_TESTS_PROPERTIES(${tname} PROPERTIES LABELS "SMOKE_TEST" )
 endfunction(smoke_test)
 
+mpi_test(pcu_test 8 ./pcu_test)
+
 mpi_test(shapefun 1 ./shapefun)
 mpi_test(shapefun2 1 ./shapefun2)
 mpi_test(bezierElevation 1 ./bezierElevation)