Added more mesh_connection methods.

data/ref-cube.mesh

@@ -11,6 +11,21 @@ MFEM mesh v1.0
 # CUBE        = 5
 # PRISM       = 6
 #
+#       7-------6
+#      /|      /|
+#     / |     / |
+#    /  |    /  |
+#   4-------5   |
+#   |   3---|---2
+#   |  /    |  /
+#   | /     | /
+#   |/      |/
+#   0-------1
+#
+#   z
+#   | y
+#   |/
+#   *--x
 
 dimension
 3

data/ref-prism.mesh

@@ -11,6 +11,25 @@ MFEM mesh v1.0
 # CUBE        = 5
 # PRISM       = 6
 #
+#
+#      5
+#     /. \
+#    / .    \
+#   3--------4
+#   |  .     |
+#   |  .     |
+#   |  .     |
+#   |  .     |
+#   |  2     |
+#   | . .    |
+#   |.     . |
+#   0--------1
+#
+#   z
+#   | y
+#   |/
+#   *--x
+
 
 dimension
 3

data/ref-pyramid.mesh

@@ -12,6 +12,22 @@ MFEM mesh v1.0
 # PRISM       = 6
 # PYRAMID     = 7
 #
+#
+#   4- _  
+#   |\ . - _
+#   | \  .   -  _
+#   |  \  3.......2
+#   |   \.       /
+#   |   .\      /
+#   |  .  \    /
+#   | .    \  /
+#   |.      \/
+#   0-------1
+#
+#   z
+#   | y
+#   |/
+#   *--x
 
 dimension
 3

data/ref-square.mesh

@@ -11,6 +11,11 @@ MFEM mesh v1.0
 # CUBE        = 5
 # PRISM       = 6
 #
+#  3----2
+#  |    |
+#  |    |
+#  0----1
+#
 
 dimension
 2

data/ref-tetrahedron.mesh

@@ -11,6 +11,24 @@ MFEM mesh v1.0
 # CUBE        = 5
 # PRISM       = 6
 #
+#   3
+#   |\
+#   |.\
+#   |  \
+#   | . \
+#   |    \
+#   |  .  \
+#   |  2   \
+#   | .  .  \
+#   |.      .\
+#   0---------1
+#
+#   z
+#   | y
+#   |/
+#   *--x
+
+
 
 dimension
 3

data/ref-triangle.mesh

@@ -11,6 +11,11 @@ MFEM mesh v1.0
 # CUBE        = 5
 # PRISM       = 6
 #
+#  2
+#  |\
+#  | \
+#  0--1
+#
 
 dimension
 2

general/array.hpp

@@ -23,6 +23,7 @@
 #include <cstring>
 #include <algorithm>
 #include <type_traits>
+#include <vector>
 
 namespace mfem
 {
@@ -92,6 +93,10 @@ public:
    template <typename CT, int N>
    explicit inline Array(const CT (&values)[N]);
 
+   /// Initilizer list contructor for syntax like Array<int> a = {1,2,3,4,5};
+   inline Array(std::initializer_list<T> list)
+   {size = list.size(); data.New(size); std::copy(list.begin(), list.end(), &(data[0]));}
+
    /// Destructor
    inline ~Array() { TypeAssert(); data.Delete(); }
 
@@ -342,7 +347,6 @@ inline bool operator!=(const Array<T> &LHS, const Array<T> &RHS)
    return !( LHS == RHS );
 }
 
-
 /// Utility function similar to std::as_const in c++17.
 template <typename T> const T &AsConst(T &a) { return a; }
 

mesh/mesh.cpp

@@ -1502,7 +1502,8 @@ void Mesh::DestroyPointers()
       FreeElement(faces[i]);
    }
 
-   DestroyTables();
+   delete connect;
+
 }
 
 void Mesh::Destroy()
@@ -2914,6 +2915,9 @@ void Mesh::FinalizeTopology(bool generate_bdr)
    // set the mesh type: 'meshgen', ...
    SetMeshGen();
 
+   if (connect) {delete connect;}
+   connect = new MeshConnections(*this);
+
    // generate the faces
    if (Dim > 2)
    {
@@ -3546,6 +3550,14 @@ void Mesh::Make1D(int n, double sx)
    bdr_attributes.Append(1); bdr_attributes.Append(2);
 }
 
+
+Mesh::Mesh()
+{
+   SetEmpty();
+   connect = new MeshConnections(*this);
+}
+
+
 Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
 {
    Dim = mesh.Dim;
@@ -3667,6 +3679,7 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
       Nodes = mesh.Nodes;
       own_nodes = 0;
    }
+   connect = new MeshConnections(*this);
 }
 
 Mesh::Mesh(Mesh &&mesh) : Mesh()
@@ -3751,6 +3764,7 @@ Mesh::Mesh(const char *filename, int generate_edges, int refine,
    {
       Load(imesh, generate_edges, refine, fix_orientation);
    }
+   connect = new MeshConnections(*this);
 }
 
 Mesh::Mesh(std::istream &input, int generate_edges, int refine,
@@ -3758,6 +3772,7 @@ Mesh::Mesh(std::istream &input, int generate_edges, int refine,
 {
    SetEmpty();
    Load(input, generate_edges, refine, fix_orientation);
+   connect = new MeshConnections(*this);
 }
 
 void Mesh::ChangeVertexDataOwnership(double *vertex_data, int len_vertex_data,
@@ -3824,6 +3839,7 @@ Mesh::Mesh(double *vertices_, int num_vertices,
    NumOfBdrElements = num_boundary_elements;
 
    FinalizeTopology();
+   connect = new MeshConnections(*this);
 }
 
 Element *Mesh::NewElement(int geom)
@@ -4305,6 +4321,8 @@ Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
       own_nodes = 1;
    }
 
+   connect = new MeshConnections(*this);
+
 #ifdef MFEM_DEBUG
    CheckElementOrientation(false);
    CheckBdrElementOrientation(false);
@@ -4316,6 +4334,7 @@ Mesh::Mesh(Mesh *orig_mesh, int ref_factor, int ref_type)
    Array<int> ref_factors(orig_mesh->GetNE());
    ref_factors = ref_factor;
    MakeRefined_(*orig_mesh, ref_factors, ref_type);
+   connect = new MeshConnections(*this);
 }
 
 void Mesh::MakeRefined_(Mesh &orig_mesh, const Array<int> ref_factors,
@@ -5917,6 +5936,17 @@ void Mesh::GetGeometries(int dim, Array<Geometry::Type> &el_geoms) const
    }
 }
 
+//TABLE_EDIT
+void Mesh::GetElementVertices(int i, Array<int> &v) const
+{ 
+   elements[i]->GetVertices(v);
+}
+
+void Mesh::GetBdrElementVertices(int i, Array<int> &v) const
+{ 
+   boundary[i]->GetVertices(v);
+}
+
 void Mesh::GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const
 {
    if (el_to_edge)
@@ -6001,6 +6031,64 @@ void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
    }
 }
 
+void Mesh::GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const
+{
+   MFEM_VERIFY(el_to_face != NULL, "el_to_face not generated");
+
+   el_to_face->GetRow(i, faces);
+
+   int n = faces.Size();
+   ori.SetSize(n);
+
+   for (int j = 0; j < n; j++)
+   {
+      if (faces_info[faces[j]].Elem1No == i)
+      {
+         ori[j] = faces_info[faces[j]].Elem1Inf % 64;
+      }
+      else
+      {
+         MFEM_ASSERT(faces_info[faces[j]].Elem2No == i, "internal error");
+         ori[j] = faces_info[faces[j]].Elem2Inf % 64;
+      }
+   }
+}
+
+void Mesh::GetBdrElementFace(int i, int *f, int *o) const
+{
+   const int *bv, *fv;
+
+   *f = be_to_face[i];
+   bv = boundary[i]->GetVertices();
+   fv = faces[be_to_face[i]]->GetVertices();
+
+   // find the orientation of the bdr. elem. w.r.t.
+   // the corresponding face element (that's the base)
+   switch (GetBdrElementType(i))
+   {
+      case Element::TRIANGLE:
+         *o = GetTriOrientation(fv, bv);
+         break;
+      case Element::QUADRILATERAL:
+         *o = GetQuadOrientation(fv, bv);
+         break;
+      default:
+         MFEM_ABORT("invalid geometry");
+   }
+}
+
+void Mesh::GetFaceVertices(int i, Array<int> &vert) const
+{
+   if (Dim == 1)
+   {
+      vert.SetSize(1); vert[0] = i;
+   }
+   else
+   {
+      faces[i]->GetVertices(vert);
+   }
+}
+
 void Mesh::GetEdgeVertices(int i, Array<int> &vert) const
 {
    // the two vertices are sorted: vert[0] < vert[1]
@@ -6133,52 +6221,6 @@ Table *Mesh::GetFaceToElementTable() const
    return face_elem;
 }
 
-void Mesh::GetElementFaces(int i, Array<int> &el_faces, Array<int> &ori) const
-{
-   MFEM_VERIFY(el_to_face != NULL, "el_to_face not generated");
-
-   el_to_face->GetRow(i, el_faces);
-
-   int n = el_faces.Size();
-   ori.SetSize(n);
-
-   for (int j = 0; j < n; j++)
-   {
-      if (faces_info[el_faces[j]].Elem1No == i)
-      {
-         ori[j] = faces_info[el_faces[j]].Elem1Inf % 64;
-      }
-      else
-      {
-         MFEM_ASSERT(faces_info[el_faces[j]].Elem2No == i, "internal error");
-         ori[j] = faces_info[el_faces[j]].Elem2Inf % 64;
-      }
-   }
-}
-
-void Mesh::GetBdrElementFace(int i, int *f, int *o) const
-{
-   const int *bv, *fv;
-
-   *f = be_to_face[i];
-   bv = boundary[i]->GetVertices();
-   fv = faces[be_to_face[i]]->GetVertices();
-
-   // find the orientation of the bdr. elem. w.r.t.
-   // the corresponding face element (that's the base)
-   switch (GetBdrElementType(i))
-   {
-      case Element::TRIANGLE:
-         *o = GetTriOrientation(fv, bv);
-         break;
-      case Element::QUADRILATERAL:
-         *o = GetQuadOrientation(fv, bv);
-         break;
-      default:
-         MFEM_ABORT("invalid geometry");
-   }
-}
-
 int Mesh::GetBdrElementEdgeIndex(int i) const
 {
    switch (Dim)
@@ -6191,6 +6233,37 @@ int Mesh::GetBdrElementEdgeIndex(int i) const
    return -1;
 }
 
+int GetBdrElementToFaceOrientation(int beid)
+{
+   return 0; 
+}
+
+int GetFaceOrientation(int faceid)
+{
+   return 0;
+}
+
+void GetFaceOrientation(const Array<int> &faceids, Array<int> &ori)
+{
+
+}
+
+void GetEdgeOrientationsInElement(int elemid, Array<int> &edge_ori)
+{
+
+}
+
+void GetEdgeOrientationsInFace(int faceid, Array<int> &edge_ori)
+{
+
+}
+
+void GetEdgeOrientationsInBdrElement(int beid, Array<int> &edge_ori)
+{
+
+}
+
+
 void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
 {
    int fid = GetBdrElementEdgeIndex(bdr_el);
@@ -6213,6 +6286,8 @@ void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
    info = fi.Elem1Inf + ori;
 }
 
+//TABLE_EDIT^
+
 Element::Type Mesh::GetElementType(int i) const
 {
    return elements[i]->GetType();
@@ -9220,6 +9295,7 @@ Mesh::Mesh(const NCMesh &ncmesh_)
    InitTables();
    InitFromNCMesh(ncmesh_);
    SetAttributes();
+   connect = new MeshConnections(*this);
 }
 
 void Mesh::Swap(Mesh& other, bool non_geometry)

mesh/mesh.hpp

@@ -23,6 +23,7 @@
 #include "../fem/eltrans.hpp"
 #include "../fem/coefficient.hpp"
 #include "../general/zstr.hpp"
+#include "mesh_connections.hpp"
 #ifdef MFEM_USE_ADIOS2
 #include "../general/adios2stream.hpp"
 #endif
@@ -57,11 +58,16 @@ class Mesh
 #endif
    friend class NCMesh;
    friend class NURBSExtension;
+   friend class MeshConnections;
 
 #ifdef MFEM_USE_ADIOS2
    friend class adios2stream;
 #endif
 
+public:
+
+   MeshConnections *connect = nullptr;
+
 protected:
    int Dim;
    int spaceDim;
@@ -217,6 +223,7 @@ protected:
    Array<FaceInfo> faces_info;
    Array<NCFaceInfo> nc_faces_info;
 
+   //TABLE_EDIT
    Table *el_to_edge;
    Table *el_to_face;
    Table *el_to_el;
@@ -560,7 +567,7 @@ protected:
 
 public:
 
-   Mesh() { SetEmpty(); }
+   Mesh(); 
 
    /** Copy constructor. Performs a deep copy of (almost) all data, so that the
        source mesh can be modified (e.g. deleted, refined) without affecting the
@@ -694,10 +701,10 @@ public:
 
    /** @anchor mfem_Mesh_init_ctor
        @brief _Init_ constructor: begin the construction of a Mesh object. */
-   Mesh(int Dim_, int NVert, int NElem, int NBdrElem = 0, int spaceDim_ = -1)
-   {
+   Mesh(int Dim_, int NVert, int NElem, int NBdrElem = 0, int spaceDim_ = -1) {
       if (spaceDim_ == -1) { spaceDim_ = Dim_; }
       InitMesh(Dim_, spaceDim_, NVert, NElem, NBdrElem);
+      connect = new MeshConnections(*this);
    }
 
    /** @name Methods for Mesh construction.
@@ -849,23 +856,26 @@ public:
    {
       Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
       Finalize(true); // refine = true
+      connect = new MeshConnections(*this);
    }
 
    /// Deprecated: see @a MakeCartesian2D.
    MFEM_DEPRECATED
    Mesh(int nx, int ny, Element::Type type, bool generate_edges = false,
-        double sx = 1.0, double sy = 1.0, bool sfc_ordering = true)
+        double sx = 1.0, double sy = 1.0, bool sfc_ordering = true) 
    {
       Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
       Finalize(true); // refine = true
+      connect = new MeshConnections(*this);
    }
 
    /// Deprecated: see @a MakeCartesian1D.
    MFEM_DEPRECATED
-   explicit Mesh(int n, double sx = 1.0)
+   explicit Mesh(int n, double sx = 1.0) 
    {
       Make1D(n, sx);
       // Finalize(); // reminder: not needed
+      connect = new MeshConnections(*this);
    }
 
    /** Creates mesh by reading a file in MFEM, Netgen, or VTK format. If
@@ -1098,59 +1108,147 @@ public:
       { mesh.GetGeometries(dim, *this); }
    };
 
-   /// Returns the indices of the vertices of element i.
-   void GetElementVertices(int i, Array<int> &v) const
-   { elements[i]->GetVertices(v); }
+   //TABLE_EDIT
 
-   /// Returns the indices of the vertices of boundary element i.
-   void GetBdrElementVertices(int i, Array<int> &v) const
-   { boundary[i]->GetVertices(v); }
+   /**  @brief Returns the indices of the vertices of element i.
+        Deprecated, please use MeshConnections:
+        mesh->connect.ChildrenOfEntity({MeshDim,AllIdx,i}, {0,AllIdx,{}})*/
+   MFEM_DEPRECATED void GetElementVertices(int i, Array<int> &v) const;
 
-   /// Return the indices and the orientations of all edges of element i.
-   void GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const;
+   /**  @brief Returns the indices of the vertices of bdr element i.
+        Deprecated, please use MeshConnections:
+        mesh->connect.ChildrenOfEntity({MeshDim-1,BdrIdx,i}, {0,AllIdx,{}})*/
+   MFEM_DEPRECATED void GetBdrElementVertices(int i, Array<int> &v) const;
 
-   /// Return the indices and the orientations of all edges of bdr element i.
-   void GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const;
+   /** Return the indices and the orientations of all edges of element i.
+       Works for both 2D (face=edge) and 3D faces. 
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim,AllIdx,i}, {1,AllIdx,{}})
+       mesh->GetEdgeOrientationsInElement(elemid)*/
+   MFEM_DEPRECATED void GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const;
+
+   /** Return the indices and the orientations of all edges of element i. 
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim-1,BdrIdx,i}, {1,AllIdx,{}})
+       mesh->GetEdgeOrientationsInBdrElement(i)*/
+   MFEM_DEPRECATED void GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const;
 
    /** Return the indices and the orientations of all edges of face i.
-       Works for both 2D (face=edge) and 3D faces. */
-   void GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const;
+       Works for both 2D (face=edge) and 3D faces. 
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim-1,AllIdx,i}, {1,AllIdx,{}})
+       mesh->GetEdgeOrientationsInFace(i)*/
+   MFEM_DEPRECATED void GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const;
 
-   /// Returns the indices of the vertices of face i.
-   void GetFaceVertices(int i, Array<int> &vert) const
-   {
-      if (Dim == 1)
-      {
-         vert.SetSize(1); vert[0] = i;
-      }
-      else
-      {
-         faces[i]->GetVertices(vert);
-      }
-   }
+   /** @Brief  Return the indices and the orientations of all faces of element i.
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim,AllIdx,i}, {MeshDim-1,AllIdx,{}})
+       mesh->GetFaceOrientation(faces)*/
+   MFEM_DEPRECATED void GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const;
 
-   /// Returns the indices of the vertices of edge i.
-   void GetEdgeVertices(int i, Array<int> &vert) const;
+   /** @brief Return the index and the orientation of the face of bdr element i. (3D) only
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim-1,AllIdx,i}, {0,AllIdx,{}})
+       mesh->GetBdrElementToFaceOrientation(i)*/
+   MFEM_DEPRECATED void GetBdrElementFace(int i, int *f, int *o) const;
 
-   /// Returns the face-to-edge Table (3D)
-   Table *GetFaceEdgeTable() const;
+   /** @brief Return the indices of the vertices of face i.
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({MeshDim-1,AllIdx,i}, {0,AllIdx,{}})*/
+   MFEM_DEPRECATED void GetFaceVertices(int i, Array<int> &vert) const;
 
-   /// Returns the edge-to-vertex Table (3D)
-   Table *GetEdgeVertexTable() const;
+   /** @brief Return the indices of the vertices of edge i.
+       Deprecated, please use MeshConnections:
+       mesh->connect.ChildrenOfEntity({1,AllIdx,i}, {0,AllIdx,{}})*/
+   MFEM_DEPRECATED void GetEdgeVertices(int i, Array<int> &vert) const;
 
-   /// Return the indices and the orientations of all faces of element i.
-   void GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const;
+   /** @brief Get the face to edge Table (3D)
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(2,AllIdx,1,AllIdx)*/
+   MFEM_DEPRECATED Table *GetFaceEdgeTable() const;
 
-   /// Return the index and the orientation of the face of bdr element i. (3D)
-   void GetBdrElementFace(int i, int *f, int *o) const;
+   /** @brief Get the edge to vertex Table (3D)
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(1,AllIdx,0,AllIdx)*/
+   MFEM_DEPRECATED Table *GetEdgeVertexTable() const;
 
-   /** Return the vertex index of boundary element i. (1D)
+   /** @brief  Get the element to element neightbot table
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetNeighborTable(MeshDim,AllIdx,MeshDim-1)*/
+   MFEM_DEPRECATED const Table &ElementToElementTable();
+
+   /** @brief  Get the element to face table
+       which only exists in 3D
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(0,AllIdx,MeshDim,AllIdx)*/
+   MFEM_DEPRECATED const Table &ElementToFaceTable() const;
+
+   /** @brief  Get the element to edge table
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(0,AllIdx,MeshDim,AllIdx)*/
+   MFEM_DEPRECATED const Table &ElementToEdgeTable() const;
+
+   /** @brief  The returned Table must be destroyed by the caller
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(0,AllIdx,MeshDim,AllIdx)*/
+   MFEM_DEPRECATED Table *GetVertexToElementTable();
+
+   /** @brief Return the "face"-element Table. Here "face" refers to face (3D),
+       edge (2D), or vertex (1D).
+       The returned Table must be destroyed by the caller. 
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetTable(MeshDim-1,AllIdx,MeshDim,AllIdx)*/
+   MFEM_DEPRECATED Table *GetFaceToElementTable() const;
+
+   /** @brief Return the vertex index of boundary element i. (1D)
        Return the edge index of boundary element i. (2D)
-       Return the face index of boundary element i. (3D) */
-   int GetBdrElementEdgeIndex(int i) const;
+       Return the face index of boundary element i. (3D) 
+       Deprecated, please use MeshConnections:
+       mesh->connect.GetAllIdxFromBdrIdx(i)*/
+   MFEM_DEPRECATED int GetBdrElementEdgeIndex(int i) const;
+
+   /** @brief Return the orientation # of the boundary element 
+       w.r.t. it's corresponding face.*/ 
+   int GetBdrElementToFaceOrientation(int beid);
+
+   /** @brief Return the orientation # of the face element w.r.t. to the  
+       2 elements shared by the face.  The face in the first 
+       element in the face_info has an orientation of 0
+       and the second face will have the orientation # representing the 
+       the permutation between  the face in the first element and the face
+       in the second element.*/
+   int GetFaceOrientation(int faceid);
+
+   /** @brief Return the orientation #s of the face elements w.r.t. to the  
+       2 elements shared by each face.  The face in the first 
+       element in the face_info has an orientation of 0
+       and the second face will have the orientation # representing the 
+       the permutation between  the face in the first element and the face
+       in the second element.*/
+   void GetFaceOrientation(const Array<int> &faceids, Array<int> &ori);
+
+   /** @brief Return the orientations of the edges in the given element
+       w.r.t. the default ordering (the lower vertex number is first 
+       followed by the higher vertex number).  If the edge vertices 
+       are in this order then we return 1 otherwise we return -1*/
+   void GetEdgeOrientationsInElement(int elemid, Array<int> &edge_ori);   
+
+   /** @brief Return the orientations of the edges in the given face
+       w.r.t. the default ordering (the lower vertex number is first 
+       followed by the higher vertex number).  If the edge vertices 
+       are in this order then we return 1 otherwise we return -1*/
+   void GetEdgeOrientationsInFace(int faceid, Array<int> &edge_ori);
+
+   /** @brief Return the orientations of the edges in the given boundary
+       element w.r.t. the default ordering (the lower vertex number is first 
+       followed by the higher vertex number).  If the edge vertices 
+       are in this order then we return 1 otherwise we return -1*/
+   void GetEdgeOrientationsInBdrElement(int beid, Array<int> &edge_ori);
+
+   //TABLE_EDIT
 
    /** @brief For the given boundary element, bdr_el, return its adjacent
-       element and its info, i.e. 64*local_bdr_index+bdr_orientation. */
+       element and its info, i.e. 64*local_bdr_index+bdr_orientation.*/
    void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const;
 
    /// Returns the type of element i.
@@ -1463,20 +1561,6 @@ public:
    /// Set the attribute of boundary element i.
    void SetBdrAttribute(int i, int attr) { boundary[i]->SetAttribute(attr); }
 
-   const Table &ElementToElementTable();
-
-   const Table &ElementToFaceTable() const;
-
-   const Table &ElementToEdgeTable() const;
-
-   ///  The returned Table must be destroyed by the caller
-   Table *GetVertexToElementTable();
-
-   /** Return the "face"-element Table. Here "face" refers to face (3D),
-       edge (2D), or vertex (1D).
-       The returned Table must be destroyed by the caller. */
-   Table *GetFaceToElementTable() const;
-
    /** This method modifies a tetrahedral mesh so that Nedelec spaces of order
        greater than 1 can be defined on the mesh. Specifically, we
        1) rotate all tets in the mesh so that the vertices {v0, v1, v2, v3}

mesh/mesh_connections.cpp

@@ -0,0 +1,200 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mesh_connections.hpp"
+#include "mesh.hpp"
+
+namespace mfem
+{
+
+MeshConnections::MeshConnections(const Mesh &m) :
+   mesh(m)
+{
+   T.resize(mesh.Dimension()+1);
+   NT.resize(mesh.Dimension()+1);
+   for (int i = 0; i < T.size()+1; ++i)
+   {
+      T[i].resize(T.size(), nullptr);
+      NT[i].resize(NT.size(), nullptr);
+   }
+}
+
+void MeshConnections::ConnectionsOfEntity(const EntityIndex entity, EntityIndices &connected) const
+{
+   Table *t = GetTable(entity.kind, connected.kind);
+   t->GetRow(entity.index, connected.indices);
+}
+
+void MeshConnections::ConnectionsOfEntities(const EntityIndices entities, 
+                                            EntityIndices &connected, bool covered) const
+{
+   Table *t = GetTable(entities.kind, connected.kind);
+
+   //Find all the connected entities touched by the entities
+   std::set<int> touched_set;
+   for (int i = 0; i < entities.indices.Size(); ++i)
+   {
+      int row_sz = t->RowSize(entities.indices[i]);
+      const int *row = t->GetRow(entities.indices[i]);
+      for (int j = 0; j < row_sz; ++j)
+      {
+         touched_set.insert(row[j]);
+      }
+   }
+
+   //In this case we want the entities to cover the entire connected entity to count
+   if (entities.kind.dim > connected.kind.dim && covered) 
+   {
+      Table *trev = GetTable(connected.kind, entities.kind);
+      //Put the entities into a hashing set for fast finding
+      std::unordered_set<int> entity_set(entities.indices.begin(), entities.indices.end());
+
+      connected.indices.SetSize(touched_set.size());
+      int num_covered = 0;
+      for (auto touched = touched_set.begin(); touched != touched_set.end(); ++touched)
+      {
+         int row_sz = trev->RowSize(*touched);
+         const int *row = trev->GetRow(*touched);
+         bool entity_covered = true;
+         for (int j = 0; j < row_sz; ++j)
+         {
+            if (entity_set.count(row[j]) < 1)
+            {
+               entity_covered = false;
+               break;
+            }
+         }
+
+         if (entity_covered)
+         {
+            connected.indices[num_covered] = *touched;
+            num_covered ++;
+         }      
+      }
+      connected.indices.SetSize(num_covered);      
+   } 
+   else
+   {
+      connected.indices.SetSize(touched_set.size());
+      std::copy(touched_set.begin(), touched_set.end(), connected_indices.begin());
+   }
+}
+
+bool MeshConnections::IsChild(const EntityIndex &parent, const EntityIndex &child) const
+{
+   EntityIndices children;
+   children.kind = child.kind;
+   ChildrenOfEntity(parent, children);
+   return children.Find(child.index) > -1;
+}
+
+void MeshConnections::ChildrenOfEntity(const EntityIndex &parent,
+                                       EntityIndices &children) const
+{
+   MFEM_ASSERT(parent.kind.dim > children.kind.dim, "Dimension of parent must be > child dimension.");
+   ConnectionsOfEntity(parent, children);
+}
+
+void MeshConnections::ChildrenOfEntities(const EntityIndices &parents,
+                                         EntityIndices &children) const
+{
+   MFEM_ASSERT(parents.kind.dim > children.kind.dim, "Dimension of parent must be > child dimension.");
+   ConnectionsOfEntities(parents, children, false);
+}
+
+void MeshConnections::ParentsOfEntity(const EntityIndex &child,
+                                      EntityIndices &parents) const
+{
+   MFEM_ASSERT(parents.kind.dim > child.kind.dim, "Dimension of parent must be > child dimension.");
+   ConnectionsOfEntity(child, parents);
+}
+
+void MeshConnections::ParentsOfAnyEntities(const EntityIndices &children,
+                                           EntityIndices &parents) const
+{
+   MFEM_ASSERT(parents.kind.dim > children.kind.dim, "Dimension of parent must be > child dimension.");
+   ConnectionsOfEntities(children, parents, false);
+}
+
+void MeshConnections::ParentsCoveredByEntities(const EntityIndices &children,
+                                               EntityIndices &parents) const
+{
+   MFEM_ASSERT(parents.kind.dim > children.kind.dim, "Dimension of parent must be > child dimension.");
+   ConnectionsOfEntities(children, parents, true);
+}
+
+void MeshConnections::NeighborsOfEntity(const EntityIndex &entity, EntityIndexKind shared_kind,
+                                        EntityIndices &neighbors) const
+{
+   MFEM_ASSERT(entity.kind == neighbors.kind, "Neighbors entity must have the same kind.")
+   Table *t = GetNeighborTable(entity.kind, shared_kind);
+   t->GetRow(entity.index, neighbors.indices);
+}
+
+void MeshConnections::NeighborsOfEntities(const EntityIndices &entities, EntityIndexKind shared_kind, 
+                                          EntityIndices &neighbors) const
+{
+   Table *t = GetTable(entities.kind, connected.kind);
+
+   //Find all the connected entities touched by the entities
+   //put it in a set to remove duplicates
+   std::set<int> touched_set;
+   for (int i = 0; i < entities.indices.Size(); ++i)
+   {
+      int row_sz = t->RowSize(entities.indices[i]);
+      const int *row = t->GetRow(entities.indices[i]);
+      for (int j = 0; j < row_sz; ++j)
+      {
+         touched_set.insert(row[j]);
+      }
+   }
+   neighbors.indices.SetSize(touched_set.size());
+   std::copy(touched_set.begin(), touched_set.end(), neighbors.indices.begin());   
+}
+
+int MeshConnections::GetAllIdxFromBdrIdx(int bdr_idx) const
+{
+   int all_idx;
+   return all_idx;
+}
+
+
+int MeshConnections::GetBdrIdxFromAllIdx(int all_idx) const
+{
+   int bdr_idx;
+   return bdr_idx;
+}
+
+Table* MeshConnections::GetTable(EntityIndexKind row_kind, EntityIndexKind col_kind) const
+{
+   const int mesh_dim = mesh.Dimension();
+   MFEM_ASSERT(row_kind.space != Bdr || row_kind.dim == mesh_dim - 1, "Row boundary index must be 1 less than the mesh dimension.")
+   MFEM_ASSERT(col_kind.space != Bdr || col_kind.dim == mesh_dim - 1, "Col boundary index must be 1 less than the mesh dimension.")
+   MFEM_ASSERT(row_kind != col_kind, "Row and Col kinds must be different.  For neighbor connections use GetNeighborTable.")
+
+   int row_i = row_kind.space == Bdr ? mesh_dim + 1 : row_kind.dim;
+   int col_i = col_kind.space == Bdr ? mesh_dim + 1 : col_kind.dim;
+   Table *t = T[row_i][col_i];
+   if (t != nullptr)
+   {
+      return t;
+   }
+
+   return nullptr;
+}
+
+
+Table* MeshConnections::GetNeighborTable(EntityIndexKind entity_kind, EntityIndexKind shared_kind) const
+{
+   return nullptr;
+}
+
+}

mesh/mesh_connections.hpp

@@ -0,0 +1,390 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_MESH_CONNECTIONS
+#define MFEM_MESH_CONNECTIONS
+
+#include "../general/array.hpp"
+#include <vector>
+
+namespace mfem
+{
+
+class Table;
+class Mesh;
+
+
+/** @a Boundary indices are indexing the subset of boundary entities in a mesh.  MFEM meshes 
+   currently only have one kind of Boundary entity indexed this way which is the mesh_dim-1 boundary.
+   @a All indices are indexing all of the entities of a given dimension in a mesh (including
+   those on the boundary).*/
+enum EntityIndexSpace : bool 
+{ 
+   Bdr=true,    ///Index space contains only the entities on the boundary
+   All=false    ///Index space contains all entities including those on the boundary
+};
+
+/** @brief Struct representing the kind of entity a following index will refer to.  In
+   MFEM meshes there is currently support for up to 5 entity kinds that all have their
+   own index numberings: (0D, 1D, 2D, 1D Boundary) in 2D and (0D, 1D, 2D, 3D, 2D Boundary)
+   in 3D.  These kinds are represented by the @a dim and @a boundary members.  For
+   instance a 1D boundary member on a 2D mesh will have dim=1 and space=Bdr.  It
+   should be noted that edges on the boundary of a 2D mesh have 2 numbers associated with them,
+   their 1D edge numering and their 1D boundary element numbering.  This is of course also
+   the case with 3D meshes and the 2D faces and 2D boundary elements.**/
+struct EntityIndexKind
+{
+   unsigned short dim;
+   EntityIndexSpace space;
+};
+
+/** @brief Struct representing the index number of a geometrical mesh entity.  In
+   MFEM meshes there is currently support for up to 5 entity kinds that all have their
+   own index numberings: (0D, 1D, 2D, 1D Boundary) in 2D and (0D, 1D, 2D, 3D, 2D Boundary)
+   in 3D.  These kinds are represented by the @a dim and @a boundary members.  For
+   instance a 1D boundary member on a 2D mesh will have dim=1 and space=Bdr.  It
+   should be noted that edges on the boundary of a 2D mesh have 2 numbers associated with them,
+   their 1D edge numering and their 1D boundary element numbering.  This is of course also
+   the case with 3D meshes and the 2D faces and 2D boundary elements.**/
+struct EntityIndex
+{
+  EntityIndexKind kind;
+  int index;
+};
+
+/** @brief Struct representing the plural index numbers of a geometrical mesh entities.  In
+   MFEM meshes there is currently support for up to 5 entity kinds that all have their
+   own index numberings: (0D-All, 1D=All, 2D-All, 1D-Bdr) in 2D and 
+   (0D-All, 1D-All, 2D-All, 3D-All, 2D-Bdr) in 3D.  These kinds are represented by the 
+   @a dim and @a space members.  For instance a 1D boundary member on a 2D mesh will have 
+   dim=1 and space=Bdr.  It should be noted that edges on the boundary of a 2D mesh have 
+   2 numbers associated with them, their 1D all edge numering and their 1D boundary element numbering.  This 
+   is of course also the case with 3D meshes and the 2D faces and 2D boundary elements.**/
+struct EntityIndices
+{
+   EntityIndexKind kind;
+   Array<int> indices;
+};
+
+
+/** @brief Class that exposes connections between Mesh entities such as
+    vertices, edges, faces, and volumes.  This is done without assuming
+    the elements have a particular dimension or type and will support
+    dimensions > 3 in the future.  EntityIndex numbers are described in the
+    \ref EntityIndex and \ref EntityIndices structs.  The connections are defined
+    on the Mesh object given at the time of construction, and Table objects are built
+    lazily as they are required so the first query of a paticular kind may be a lot slower
+    than the following queries.  It is worth noting that our naming convention of child
+    and parent connections actually describes ancestor/decendent relationship eg. a node can
+    be a "child" of edges, faces, and volumes.  It should also be noted that the ConnectionsOf
+    methods provide all of the functionality of the ChildrenOf/ParentsOf methods and vice-versa.
+    Both approaches are provided as a convenience to folks that are used to thiking of meshes 
+    in either context.**/
+class MeshConnections
+{
+   friend class Mesh;
+
+private:
+   /// The mesh object that the connections are defined on
+   const Mesh &mesh;
+
+   /** @brief 2D Array of tables describing connections between entity indices in the Mesh.
+       The 2D array represents tables mapping from row indices to column indices of the 4-5
+       different kinds of entities (0D/All, 1D/All, 2D/All, 3D/All, 2D/Bdr) respectively.  For example
+       the Table mapping Vertices to Boundary Elements/Faces on a 3D mesh would be at T[0][3]*/
+   mutable std::vector<std::vector<Table*>> T;
+
+   /** @brief 2D Array of tables describing neighbor connections among entity indices in the Mesh.
+       The 2D array represents tables mapping from entities to neighboring entities of different kinds
+       (0D/All, 1D/All, 2D/All, 3D/All, 2D/Bdr) across the shared dimension.  For example
+       the Table for neighboring elements shared across nodes in a 3D mesh would be NT[3][0]*/   
+   mutable std::vector<std::vector<Table*>> NT;
+
+public:
+
+
+   /** @brief Sets the Mesh object connections are defined on. */
+   /** @note When this method is called, no connection tables
+       are generated so they will all be null.  They will be generated
+       lazily as access is needed.**/
+   MeshConnections(const Mesh &m);
+
+   /** @brief General method for getting the @a connected entity indices of the given 
+       @a entity.  This method can represent both parent/child and child/parent relationships
+       by setting the dimensions of @a entity and @a connected properly and will return the 
+       same answers as the associated Children/Parents methods.  There are no self connections
+       and using this method with the @a entity.kind.dim = @a connected.kind.dim will always return an 
+       empty set of connections.  If a parent/child relationship is being used the 
+       @a connected.indices will be listed in the natural order derived from the parent's
+       element definition.  Otherwise the ordering of the connected.indices  not well 
+       defined. The @a connected struct operates as an  in/out parameter.  The kind of the 
+       connected entities is set in the @a connected.kind.dim and @a connected.kind.space members.  The 
+       connected entity index numbers are placed in @a connected.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                            entity      connected
+                          dim,bdry,id  dim,bdry,ids
+      ConnectionsOfEntity({2,AllIdx,0}, {0,AllIdx,{}}) ->  {0,AllIdx,{0,1,4,3}}, ChildrenOF relationship, nodes in element order 
+      ConnectionsOfEntity({2,AllIdx,1}, {1,Bdr,{}}) ->  {1,Bdr,{1,5}}, ChildrenOf relationship
+      ConnectionsOfEntity({1,AllIdx,7}, {2,AllIdx,{}}) ->  {2,AllIdx,{0,1}}, ParentsOf relationship
+      ConnectionsOfEntity({1,Bdr,5}, {2,AllIdx,{}}) ->  {2,AllIdx,{1}}, ParentsOf relationship
+      ConnectionsOfEntity({2,AllIdx,0}, {2,AllIdx,{}}) ->  {2,AllIdx,{}}, No connections defined of the same dimension**/
+   void ConnectionsOfEntity(const EntityIndex entity, EntityIndices &connected) const;
+
+   /** @brief General method for getting the @a connected entity indices of the given set of
+       @a entities.  This method can represent both parent/child and child/parent relationships
+       by setting the dimensions of @a entities and @a connected properly and will return the 
+       same answers as the associated Children/Parents methods.  There are no self connections
+       and using this method with the @a entity.kind.dim = @a connected.kind.dim will always return an 
+       empty set of connections.  The @a covered variable sets the behavior of which entities count
+       in ParentsOf connections.  If @a covered is false then all connected parents will be returned.  If
+       @a covered is true, the only the parents that are fully covered by the set of provided children 
+       will be returned.  The @a connected struct operates as an in/out parameter.  The kind of the 
+       connected entities is set in the @a connected.kind.dim and @a connected.kind.space members.  The 
+       connected entity index numbers are placed in @a connected.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                              entities          connected
+                            dim,bdry,id       dim,bdry,ids
+      ConnectionsOfEntities({1,AllIdx,{2,3}}, {0,AllIdx,{}}, false)      ->  {0,AllIdx,{3,4,5}}, not in natural element order
+      ConnectionsOfEntities({2,AllIdx,{0,2}}, {1,Bdr,{}}, false)      ->  {1,Bdr,{0,2,4,6}} 
+      ConnectionsOfEntities({0,AllIdx,{0,8}}, {2,AllIdx,{}}, false)      ->  {2,AllIdx,{0,3}}
+      ConnectionsOfEntities({1,AllIdx,{0,1,2,6,7}}, {2,AllIdx,{}}, true) ->  {2,AllIdx,{0}}, only elem 0 is covered                     
+    */
+   void ConnectionsOfEntities(const EntityIndices entities, EntityIndices &connected, bool covered) const;
+
+   /** @brief Returns true if the @a parent entity is indeed a
+       a parent of the @a child entity.  For the example mesh
+       with vertex, edge, face/element, and boundary element Indices
+       we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                 parent        child
+               dim,bdry,id  dim,bdry,id
+      IsChild({2,AllIdx,0}, {0,AllIdx,4})  -> true, any decendent counts as a child
+      IsChild({2,AllIdx,0}, {0,AllIdx,8})  -> false
+      IsChild({2,AllIdx,3}, {1,AllIdx,10}) -> true
+      IsChild({2,AllIdx,0}, {2,AllIdx,0})  -> false, no self children
+      IsChild({2,AllIdx,0}, {1,Bdr,4})  -> true
+      IsChild({1,Bdr,3}, {0,AllIdx,8})  -> true
+      IsChild({2,Bdr,3}, {0,AllIdx,8})  -> Error, no dim 2, boundary objects**/
+   bool IsChild(const EntityIndex &parent, const EntityIndex &child) const;
+
+   /** @brief Returns the child entity indices of the given @a parent
+       entity.  If the children have a natural order in the given element
+       this will return them in that order. The @a children struct operates as an 
+       in/out parameter.  The kind of the child entities is set in the
+       @a children.dim and @a children.space members.  The child
+       entity index numbers are placed in @a children.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                         parent       children
+                       dim,bdry,id  dim,bdry,ids
+      ChildrenOfEntity({2,AllIdx,0}, {0,AllIdx,{}}) ->  {0,AllIdx,{0,1,4,3}}, nodes in element ordering
+      ChildrenOfEntity({2,AllIdx,1}, {1,Bdr,{}}) ->  {1,Bdr,{1,5}}
+      ChildrenOfEntity({1,AllIdx,3}, {0,AllIdx,{}}) ->  {1,Bdr,{1,5}}**/
+   void ChildrenOfEntity(const EntityIndex &parent, EntityIndices &children) const;
+
+   /** @brief Returns the child entity indices of the given @a parents
+       entities.  Since this is collective call that does not duplicate entries, 
+       you cannot expect the children to be listed in their natural element
+       orders. The @a children struct operates as an 
+       in/out parameter.  The kind of the child entities is set in the
+       @a children.dim and @a children.space members.  The child
+       entity index numbers are placed in @a children.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                            parent          children
+                          dim,bdry,id     dim,bdry,ids
+      ChildrenOfEntities({1,AllIdx,{2,3}}, {0,AllIdx,{}}) ->  {0,AllIdx,{3,4,5}}, not in natural element order
+      ChildrenOfEntities({2,AllIdx,{0,2}}, {1,Bdr,{}}) ->  {1,Bdr,{0,2,4,6}}
+      ChildrenOfEntities({1,Bdr,{4,7}}, {0,AllIdx,{}}) ->  {0,AllIdx,{0,3,5,8}}**/
+   void ChildrenOfEntities(const EntityIndices &parents,EntityIndices &children) const;
+
+   /** @brief Returns the parent entity indices of the given @a child
+       entity.  The @a parents struct operates as an in/out
+       parameter.  The kind of the parent entities is set in the
+       @a parents.dim and @a parents.space members.  The parent
+       entity index numbers are placed in @a parents.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                        child         parents
+                      dim,bdry,id   dim,bdry,ids
+      ParentsOfEntity({1,AllIdx,7}, {2,AllIdx,{}}) ->  {2,AllIdx,{0,1}}
+      ParentsOfEntity({1,Bdr,5}, {2,AllIdx,{}}) ->  {2,AllIdx,{1}}**/
+   void ParentsOfEntity(const EntityIndex &child, EntityIndices &parents) const;
+
+   /** @brief Returns the parent entity indices of any of the given @a children
+       entities.  The @a parents struct operates as an in/out
+       parameter.  The kind of the parent entities is set in the
+       @a parents.dim and @a parents.space members.  The parent
+       entity index numbers are placed in @a parents.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                             children          parents
+                           dim,bdry,ids      dim,bdry,ids
+      ParentsOfAnyEntities({0,AllIdx,{0,8}}, {2,AllIdx,{}}) ->  {2,AllIdx,{0,3}}
+      ParentsOfAnyEntities({1,Bdr,{4,6}}, {2,AllIdx,{}}) ->  {2,AllIdx,{0,2}}**/
+   void ParentsOfAnyEntities(const EntityIndices &children,
+                             EntityIndices &parents) const;
+
+   /** @brief Returns the parent entity indices that have all their child entities
+       covered by the @a children entities.  The @a parents struct operates as
+       an in/out parameter.  The kind of the parent entities is set in the
+       @a parents.dim and @a parents.space members.  The parent
+       entity index numbers are placed in @a parents.indices member.  For the
+       example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                                  children                parents
+                                dim,bdry,ids            dim,bdry,ids
+      ParentsCoveredByEntities({1,AllIdx,{0,1,2,6,7}}, {2,AllIdx,{}}) ->  {2,AllIdx,{0}}, element 1 is not covered
+      ParentsCoveredByEntities({0,AllIdx,{0,1,2,3,4}}, {2,AllIdx,{}}) ->  {2,AllIdx,{0}}
+      ParentsCoveredByEntities({0,AllIdx,{0,1,2,3,6}}, {1,Bdr,{}}) ->  {1,Bdr,{0,1,4,6}}**/
+   void ParentsCoveredByEntities(const EntityIndices &children,
+                                 EntityIndices &parents) const;
+
+   /** @brief Returns the neighbor entity indices of the given @a entity.  The @a neighbors
+       struct operates as an in/out parameter.  The kind of the neighbor entities is set in the
+       @a neighbors.kind.dim and @a neighbors.kind.space members.  The neighbor
+       entity index numbers are placed in @a neighbors.indices member.  The @a shared_dim parameter
+       sets the dimension of the entities that the neighbors are sharing across.  Note:  The
+       entity.kind.dim and neighbors.kind.dim must be the same.  For example if the entity/neighbor dims are 0
+       and the shared_dim is 0 this will return the vertex that are neighbors across the edges connected
+       to the vertex entity.  For the example mesh with vertex, edge, face/element, and boundary element
+       Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                            entity                                 neighbors
+                         dim,bdry,ids                             dim,bdry,ids
+      NeighborsOfEntity({0,AllIdx,1},1,AllIdx,{0,AllIdx,{}}) ->  {0,AllIdx,{0,2,4}}
+      NeighborsOfEntity({0,AllIdx,1},2,AllIdx,{0,AllIdx,{}}) ->  {0,AllIdx,{0,2,3,4,5}}
+      NeighborsOfEntity({1,Bdr,5},0,AllIdx,{1,AllIdx,{}}) ->  {0,AllIdx,{1,11}},   edge numbering
+      NeighborsOfEntity({1,Bdr,5},0,AllIdx,{1,Bdr,{}}) ->  {0,Bdr,{1,7}},    boundary element numbering
+      NeighborsOfEntity({0,AllIdx,1},2,AllIdx,{1,AllIdx,{}}) ->  error, entity.kind.dim != neighbors.kind.dim**/
+   void NeighborsOfEntity(const EntityIndex &entity, int shared_dim, EntityIndexType shared_type,
+                          EntityIndices &neighbors) const;
+
+   /** @brief Returns the neighbor entity indices of the given @a entities.  The @a neighbors
+       struct operates as an in/out parameter.  The kind of the neighbor entities is set in the
+       @a neighbors.kind.dim and @a neighbors.kind.space members.  The neighbor
+       entity index numbers are placed in @a neighbors.ids member.  The @a shared_dim parameter
+       sets the dimension of the entities that the neighbors are sharing across.  For example if the
+       entity/neighbor dims are 0 and the shared_dim is 0 this will return the vertex that are
+       neighbors across the edges connected to the vertex entity.  Note:  The entities.dim and
+       neighbors.kind.dim must be the same.  Also, none of the index numbers from the @a entities set
+       will be returned in the @a neighbors set, even if they are neighbors of others in the
+       @a entities set.  For the example mesh with vertex, edge, face/element, and
+       boundary element Indices we have:
+
+        0,All          1,All          2,All          1,Bdr
+      6---7---8      +-4-+-5-+      +---+---+      +-2-+-3-+
+      |   |   |      9  10  11      | 2 | 3 |      6   |   7
+      3---4---5      +-2-+-3-+      +---+---+      +---+---+
+      |   |   |      6   7   8      | 0 | 1 |      4   |   5
+      0---1---2      +-0-+-1-+      +---+---+      +-0-+-1-+
+
+                               entity                                     neighbors
+                            dim,bdry,ids                                dim,bdry,ids
+      NeighborsOfEntities({0,AllIdx,{1,2}},1,AllIdx,{0,AllIdx,{}}) ->  {0,AllIdx,{0,4,5}},  2 is not included
+      NeighborsOfEntities({0,AllIdx,{0,1}},2,AllIdx,{0,AllIdx,{}}) ->  {0,AllIdx,{2,3,4,5}}
+      NeighborsOfEntities({1,Bdr,{1,5}},0,AllIdx,{1,AllIdx,{}}) ->  {1,AllIdx,{0,11}},   edge numbering
+      NeighborsOfEntities({1,Bdr,{1,5}},0,AllIdx,{1,Bdr,{}}) ->  {1,Bdr,{0,7}},    boundary element numbering
+      NeighborsOfEntities({0,AllIdx,{1,2}},2,AllIdx,{1,AllIdx,{}}) ->  error, entity.kind.dim != neighbors.kind.dim**/
+   void NeighborsOfEntities(const EntityIndices &entities, EntityIndexKind shared_kind, 
+                            EntityIndices &neighbors) const;
+
+   int GetAllIdxFromBdr(int bdr_idx) const;
+   int GetBdrFromAllIdx(int all_idx) const;
+
+   /** @brief Returns a pointer to the table with the connections from the row_kind of entity to the col_kind
+       of entity.  This method will generate the table if it currently doesn't exist.  row_kind and col_kind
+       must be different.  Use GetNeighborTable for same dimension neighbor relationships.  
+   **/
+   Table* GetTable(EntityIndexKind row_kind, EntityIndexKind col_kind) const;
+
+   /** @brief Returns a pointer to the table describing the neighbors of entities of across a shared dimension.  This
+       method will generatethe table if it currently dosen't currently exist.  row_kind and col_kind must be different.
+   **/
+   Table* GetNeighborTable(EntityIndexKind entity_kind, EntityIndexKind shared_kind);
+
+};
+
+}
+
+#endif

tests/unit/data/quad_2x2.mesh

@@ -0,0 +1,52 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+#       dim = 0       dim = 2     Boundary Elems 
+#      6---7---8     +---+---+      +-2-+-3-+
+#      |   |   |     | 2 | 3 |      6   |   7
+#      3---4---5     +---+---+      +---+---+
+#      |   |   |     | 0 | 1 |      4   |   5
+#      0---1---2     +---+---+      +-0-+-1-+ 
+
+dimension
+2
+
+elements
+4
+1 3 0 1 4 3
+1 3 1 2 5 4
+1 3 3 4 7 6
+1 3 4 5 8 7
+
+boundary
+8
+3 1 0 1
+3 1 1 2
+3 1 7 6
+3 1 8 7
+3 1 3 0
+3 1 2 5
+3 1 6 3
+3 1 5 8
+
+vertices
+9
+2
+0 0
+1 0
+2 0
+0 1
+1 1
+2 1
+0 2
+1 2
+2 2

tests/unit/data/quadtri_2x2.mesh

@@ -0,0 +1,55 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+#       dim = 0       dim = 2     Boundary Elems 
+#      6---7---8     +---+---+      +-2-+-3-+
+#      |   | / |     | 3 |4/5|      6   | / 7
+#      3---4---5     +---+---+      +---+---+
+#      | / |   |     |0/1| 2 |      4 / |   5
+#      0---1---2     +---+---+      +-0-+-1-+ 
+
+dimension
+2
+
+elements
+6
+1 2 3 0 4
+1 2 1 4 0
+1 3 1 2 5 4
+1 3 3 4 7 6
+1 2 7 4 8
+1 2 5 8 4
+
+
+boundary
+8
+3 1 0 1
+3 1 1 2
+3 1 7 6
+3 1 8 7
+3 1 3 0
+3 1 2 5
+3 1 6 3
+3 1 5 8
+
+vertices
+9
+2
+0 0
+1 0
+2 0
+0 1
+1 1
+2 1
+0 2
+1 2
+2 2

tests/unit/data/tinyzoo-3d.mesh

@@ -0,0 +1,56 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+
+dimension
+3
+
+elements
+4
+1 5 0 6 7 1 3 9 10 4
+1 6 7 8 1 10 11 4
+1 7 4 1 8 11 5
+1 4 2 8 5 1
+
+boundary
+14
+3 3 3 9 6 0
+3 3 3 4 10 9
+3 3 0 1 4 3
+3 3 6 7 1 0
+3 3 6 9 10 7
+3 2 7 1 8
+3 2 2 8 1
+3 2 10 4 11
+3 2 5 11 4
+3 2 4 1 5
+3 2 2 5 1
+3 3 7 10 11 8
+3 2 2 8 5
+3 2 11 5 8
+
+vertices
+12
+3
+0 0 0
+1 0 0
+2 0 0
+0 1 0
+1 1 0
+2 1 0
+0 0 1
+1 0 1
+2 0 1
+0 1 1
+1 1 1
+2 1 1

tests/unit/general/test_array.cpp

@@ -18,10 +18,12 @@ TEST_CASE("Array init-list construction", "[Array]")
 {
    int ContigData[6] = {6, 5, 4, 3, 2, 1};
    Array<int> a(ContigData, 6);
-   Array<int> b({6.0, 5.0, 4.0, 3.0, 2.0, 1.0});
+   Array<int> b({6, 5, 4, 3, 2, 1});
+   Array<int> c = {6, 5, 4, 3, 2, 1};
 
    for (int i = 0; i < a.Size(); i++)
    {
       REQUIRE(a[i] == b[i]);
+      REQUIRE(a[i] == c[i]);
    }
 }

tests/unit/mesh/test_mesh_connections.cpp

@@ -0,0 +1,283 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+using namespace mfem;
+
+#include "unit_tests.hpp"
+
+bool cmp_set(const Array<int> &a, const Array<int> &b);
+
+
+TEST_CASE("Expression Construction", "[EntityIndices]")
+{
+   Array<int> a = {1,2,3};
+   EntityIndices idx = {1, AllIdx, a};
+   EntityIndices idx2 = {1, AllIdx, Array<int>{1,2,3}};
+   REQUIRE(idx.indices[2] == 3);
+   REQUIRE(idx2.indices[2] == 3);
+}
+
+TEST_CASE("Matching Hand Picked Indices", "[MeshConnections]")
+{
+   SECTION("2x2 Quad Mesh")
+   {
+      Mesh mesh("./data/quad_2x2.mesh");
+      EntityIndices result;
+
+      SECTION("Cells to Vertices")
+      {
+         result = {0, AllIdx, Array<int>{}};
+         mesh.connect->ChildrenOfEntity({2,AllIdx,0}, result);
+         REQUIRE(result.indices == Array<int>{0,1,4,3});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,1}, result);
+         REQUIRE(result.indices == Array<int>{1,2,5,4});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,2}, result);
+         REQUIRE(result.indices == Array<int>{3,4,7,6});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,3}, result);
+         REQUIRE(result.indices == Array<int>{4,5,8,7});
+      }
+
+      SECTION("Cells to Boundary Edges")
+      {
+         result = {1, BdrIdx, Array<int>{}};
+         mesh.connect->ChildrenOfEntity({2,AllIdx,0}, result);
+         REQUIRE(result.indices == Array<int>{4,0});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,1}, result);
+         REQUIRE(result.indices == Array<int>{1,5});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,2}, result);
+         REQUIRE(result.indices == Array<int>{2,6});
+         mesh.connect->ChildrenOfEntity({2,AllIdx,3}, result);
+         REQUIRE(result.indices == Array<int>{7,3});
+      }
+
+      SECTION("Boundary Edges to Verts")
+      {
+         result = {0, AllIdx, Array<int>{}};
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,0}, result);
+         REQUIRE(result.indices == Array<int>{0,1});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,1}, result);
+         REQUIRE(result.indices == Array<int>{1,2});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,2}, result);
+         REQUIRE(result.indices == Array<int>{7,6});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,3}, result);
+         REQUIRE(result.indices == Array<int>{8,7});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,4}, result);
+         REQUIRE(result.indices == Array<int>{3,0});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,5}, result);
+         REQUIRE(result.indices == Array<int>{2,5});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,6}, result);
+         REQUIRE(result.indices == Array<int>{6,3});
+         mesh.connect->ChildrenOfEntity({1,BdrIdx,7}, result);
+         REQUIRE(result.indices == Array<int>{5,8});         
+      }
+
+      SECTION("Verts to Cells")
+      {
+         result = {2, AllIdx, Array<int>{}};
+         mesh.connect->ParentsOfEntity({0,AllIdx,0}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,1}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,2}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,3}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,2}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,4}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,2,3}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,5}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,3}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,6}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,7}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,3}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,8}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3}));
+      }
+
+      SECTION("Verts to Boundary Edges")
+      {
+         result = {1, BdrIdx, Array<int>{}};
+         mesh.connect->ParentsOfEntity({0,AllIdx,0}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,4}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,1}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,2}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,5}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,3}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{4,6}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,4}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,5}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{5,7}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,6}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,6}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,7}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,3}));
+         mesh.connect->ParentsOfEntity({0,AllIdx,8}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,7}));
+      }
+
+      SECTION("Boundary Edges to Cells")
+      {
+         result = {2, AllIdx, Array<int>{}};
+         mesh.connect->ParentsOfEntity({1,BdrIdx,0}, result);
+         REQUIRE(result.indices == Array<int>{0});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,1}, result);
+         REQUIRE(result.indices == Array<int>{1});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,2}, result);
+         REQUIRE(result.indices == Array<int>{2});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,3}, result);
+         REQUIRE(result.indices == Array<int>{3});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,4}, result);
+         REQUIRE(result.indices == Array<int>{0});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,5}, result);
+         REQUIRE(result.indices == Array<int>{1});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,6}, result);
+         REQUIRE(result.indices == Array<int>{2});
+         mesh.connect->ParentsOfEntity({1,BdrIdx,7}, result);
+         REQUIRE(result.indices == Array<int>{3});         
+      }      
+
+      SECTION("Cell Neighbors Across Vertices")
+      {
+         result = {2, AllIdx, Array<int>{}};
+         mesh.connect->NeighborsOfEntity({2,AllIdx,0}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,2,3}));
+         mesh.connect->NeighborsOfEntity({2,AllIdx,1}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,2,3}));
+         mesh.connect->NeighborsOfEntity({2,AllIdx,2}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,3}));
+         mesh.connect->NeighborsOfEntity({2,AllIdx,3}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,2}));
+      }
+
+      SECTION("Boundary Edge Neighbors Across Vertices")
+      {
+         result = {1, BdrIdx, Array<int>{}};
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,0}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,4}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,1}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,5}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,2}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,6}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,3}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,7}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,4}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,6}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,5}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,7}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,6}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,4}));
+         mesh.connect->NeighborsOfEntity({1,BdrIdx,7}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,5}));
+      }
+
+      SECTION("Vertex Neighbors Across Cells")
+      {
+         result = {0, AllIdx, Array<int>{}};
+         mesh.connect->NeighborsOfEntity({0,AllIdx,0}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,3,4}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,1}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,2,3,4,5}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,2}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,4,5}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,3}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,4,6,7}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,4}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,2,3,5,6,7,8}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,5}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,2,4,7,8}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,6}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,4,7}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,7}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,4,5,6,8}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,8}, 0, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{4,5,7}));  
+      }
+
+      SECTION("Vertex Neighbors Across Edges")
+      {
+         result = {0, AllIdx, Array<int>{}};
+         mesh.connect->NeighborsOfEntity({0,AllIdx,0}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,3}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,1}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,2,4}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,2}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,5}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,3}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,4,6}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,4}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,3,5,7}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,5}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,4,8}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,6}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{3,7}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,7}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{4,6,8}));
+         mesh.connect->NeighborsOfEntity({0,AllIdx,8}, 1, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{5,7}));  
+      }
+
+      SECTION("Children of Multiple Entities")
+      {
+         result = {0, AllIdx, Array<int>{}};
+         mesh.connect->ChildrenOfEntities({1,BdrIdx,{2,3}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{6,7,8}));
+         mesh.connect->ChildrenOfEntities({2,AllIdx,{0,1}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{0,1,2,3,4,5}));
+
+         result = {1, BdrIdx, Array<int>{}};
+         mesh.connect->ChildrenOfEntities({2,AllIdx,{2,3}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,3,6,7}));
+      }
+
+      SECTION("Parents of Any Entities")
+      {
+         result = {2, AllIdx, Array<int>{}};
+         mesh.connect->ParentsOfAnyEntities({1,BdrIdx,{1,2}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,2}));
+         mesh.connect->ParentsOfAnyEntities({0,AllIdx,{5,6,7,8}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{1,2,3}));
+
+         result = {1, BdrIdx, Array<int>{}};
+         mesh.connect->ParentsOfAnyEntities({0,AllIdx,{3,4,6,7}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2,3,4,6}));
+      }
+
+      SECTION("Parents Covered By Entities")
+      {
+         result = {2, AllIdx, Array<int>{}};
+         mesh.connect->ParentsCoveredByEntities({0,AllIdx,{3,4,5,6,7}}, result);
+         REQUIRE(cmp_set(result.indices, Array<int>{2}));
+      }
+   }
+}
+
+
+bool cmp_set(const Array<int> &a, const Array<int> &b)
+{
+   if (a.Size() != b.Size())
+   {
+      return false;
+   }
+
+   bool same_set = true;
+   for (int i = 0; i < a.Size(); ++i)
+   {
+      if (std::find(b.begin(), b.end(), a[i]) == b.end())
+      {
+         same_set = false;
+         break;
+      }
+   }
+   return same_set;
+}