diff --git a/docker/Dockerfile b/docker/Dockerfile
index 50297c93..f3e9e5eb 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -29,8 +29,11 @@ WORKDIR /usr/lib
 RUN arch=$(uname -m) && sudo ln -s /usr/lib/${arch}-linux-gnu/libscalapack-openmpi.so.2.1.0 ./
 
 # re-install mfem with mumps, with cmake
+WORKDIR $LIB_DIR
+RUN sudo rm -rf mfem_parallel
+RUN sudo git clone https://github.com/mfem/mfem.git mfem_parallel
 WORKDIR $LIB_DIR/mfem
-RUN sudo git checkout v4.6
+RUN sudo git checkout v4.8
 WORKDIR $LIB_DIR/mfem/build
 RUN sudo -E cmake .. -DBUILD_SHARED_LIBS=YES -DMFEM_USE_MPI=YES -DMFEM_USE_GSLIB=${MG} -DMFEM_USE_METIS=YES -DMFEM_USE_METIS_5=YES -DMFEM_USE_MUMPS=YES -DGSLIB_DIR="$LIB_DIR/gslib/build" -DMUMPS_DIR="$LIB_DIR/mumps/build/local"
 RUN sudo -E make -j 16
diff --git a/include/component_topology_handler.hpp b/include/component_topology_handler.hpp
index c90222ec..0aea81d3 100644
--- a/include/component_topology_handler.hpp
+++ b/include/component_topology_handler.hpp
@@ -156,6 +156,8 @@ class ComponentTopologyHandler : public TopologyHandler
    virtual void TransferToGlobal(Array<GridFunction*> &us, Array<GridFunction*> &global_u, const int &num_var)
    { mfem_error("ComponentTopologyHandler does not yet support global grid function/mesh!\n"); }
 
+   std::map<DofTag, Array<int>> SplitDofsByBoundary(FiniteElementSpace &fes);
+
 protected:
    // Get vertex orientation of face2 (from mesh2) with respect to face1 (mesh1).
    int GetOrientation(BlockMesh *comp1, const Element::Type &be_type, const Array<int> &vtx1, const Array<int> &vtx2);
diff --git a/include/topology_handler.hpp b/include/topology_handler.hpp
index 71ee1238..62e443fe 100644
--- a/include/topology_handler.hpp
+++ b/include/topology_handler.hpp
@@ -51,6 +51,56 @@ struct TopologyData {
    Array<int> *global_bdr_attributes = NULL;
 };
 
+struct DofTag
+{
+   const Array<int> battr;
+
+   DofTag(const Array<int> &battrs)
+      : battr(battrs)
+   {
+      if (!battr.IsSorted())
+         mfem_error("DofTag: boundary attributes must be sorted beforehand!\n");
+   }
+
+   bool operator==(const DofTag &dtag) const
+   {
+      if (battr.Size() != dtag.battr.Size())
+         return false;
+
+      for (int k = 0; k < battr.Size(); k++)
+         if (battr[k] != dtag.battr[k])
+            return false;
+
+      return true;
+   }
+
+   bool operator<(const DofTag &dtag) const
+   {
+      if (battr.Size() == dtag.battr.Size())
+      {
+         for (int k = 0; k < battr.Size(); k++)
+            if (battr[k] != dtag.battr[k])
+               return (battr[k] < dtag.battr[k]);
+         
+         return false;
+      }
+      else
+         return (battr.Size() < dtag.battr.Size());
+   }
+
+   std::string print() const
+   {
+      std::string tag = "(";
+      for (int k = 0 ; k < battr.Size(); k++)
+      {
+         tag += std::to_string(battr[k]) + ",";
+      }
+      tag += ")";
+      
+      return tag;
+   }
+};
+
 const TopologyHandlerMode SetTopologyHandlerMode();
 
 class TopologyHandler
@@ -98,6 +148,8 @@ class TopologyHandler
    virtual Mesh* GetMesh(const int k) = 0;
    virtual Mesh* GetGlobalMesh() = 0;
 
+   std::map<DofTag, Array<int>> SplitBoundaryDofs(Mesh* const mesh, const FiniteElementSpace *fes) const;
+
    /*
       Methods only for ComponentTopologyHandler 
    */
diff --git a/sketches/CMakeLists.txt b/sketches/CMakeLists.txt
index dd3d99f9..7cf6bc58 100644
--- a/sketches/CMakeLists.txt
+++ b/sketches/CMakeLists.txt
@@ -28,6 +28,7 @@ add_executable(ns_mms ns_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_dg_mms ns_dg_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_rom ns_rom.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(usns usns.cpp $<TARGET_OBJECTS:scaleupROMObj>)
+add_executable(dof_separation dof_separation.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
 file(COPY inputs/gen_interface.yml DESTINATION ${CMAKE_BINARY_DIR}/sketches/inputs)
 file(COPY meshes/2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/sketches/meshes)
diff --git a/sketches/dof_separation.cpp b/sketches/dof_separation.cpp
new file mode 100644
index 00000000..8ee1fa56
--- /dev/null
+++ b/sketches/dof_separation.cpp
@@ -0,0 +1,216 @@
+// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: MIT
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "etc.hpp"
+
+using namespace std;
+using namespace mfem;
+
+struct DofTag
+{
+   Array<int> battr;
+
+   DofTag(const Array<int> &battrs)
+      : battr(battrs)
+   { battr.Sort(); }
+
+   bool operator==(const DofTag &dtag) const
+   {
+      if (battr.Size() != dtag.battr.Size())
+         return false;
+
+      for (int k = 0; k < battr.Size(); k++)
+         if (battr[k] != dtag.battr[k])
+            return false;
+
+      return true;
+   }
+
+   bool operator<(const DofTag &dtag) const
+   {
+      if (battr.Size() == dtag.battr.Size())
+      {
+         for (int k = 0; k < battr.Size(); k++)
+            if (battr[k] != dtag.battr[k])
+               return (battr[k] < dtag.battr[k]);
+         
+         return false;
+      }
+      else
+         return (battr.Size() < dtag.battr.Size());
+   }
+
+   std::string print() const
+   {
+      std::string tag = "(";
+      for (int k = 0 ; k < battr.Size(); k++)
+      {
+         tag += std::to_string(battr[k]) + ",";
+      }
+      tag += ")";
+      
+      return tag;
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int refine = 0;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&refine, "-r", "--refine",
+                  "Number of refinements.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    largest number that gives a final mesh with no more than 10,000
+   //    elements.
+   for (int l = 0; l < refine; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 5. Define a finite element space on the mesh. Here we use the
+   //    Raviart-Thomas finite elements of the specified order.
+   // FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
+   L2_FECollection l2_coll(order, dim);
+   H1_FECollection h1_coll(order+1, dim);
+   L2_FECollection pl2_coll(order, dim);
+   H1_FECollection ph1_coll(order, dim);
+
+   FiniteElementSpace fes(&mesh, &h1_coll);
+   FiniteElementSpace ufes(&mesh, &l2_coll, dim);
+   FiniteElementSpace pfes(&mesh, &ph1_coll);
+
+   // 6. Define the BlockStructure of the problem, i.e. define the array of
+   //    offsets for each variable. The last component of the Array is the sum
+   //    of the dimensions of each block.
+   Array<int> block_offsets(dim+2); // number of variables + 1
+   block_offsets[0] = 0;
+   for (int d = 1; d <= dim; d++)
+      block_offsets[d] = ufes.GetNDofs();
+   block_offsets[dim+1] = pfes.GetVSize();
+   block_offsets.PartialSum();
+
+   std::cout << "***********************************************************\n";
+   for (int d = 1; d < block_offsets.Size(); d++)
+      printf("dim(%d) = %d\n", d, block_offsets[d] - block_offsets[d-1]);
+   printf("dim(q) = %d\n", block_offsets.Last());
+   std::cout << "***********************************************************\n";
+
+   printf("ufes vsize: %d\n", ufes.GetVSize());
+   printf("pfes vsize: %d\n", pfes.GetVSize());
+
+   Array<int> udofs(ufes.GetVSize());
+   Array<int> pdofs(pfes.GetVSize());
+   udofs = 0;
+   pdofs = 0;
+
+   int max_battr = -1;
+   for (int i = 0; i < pfes.GetNBE(); i++)
+      max_battr = max(max_battr, mesh.GetBdrAttribute(i));
+
+   Array2D<int> pbattr(pfes.GetVSize(), max_battr);
+   pbattr = 0;
+
+   Array<int> vdofs;
+   FaceElementTransformations *tr;
+   printf("pfes nbe: %d\n", pfes.GetNBE());
+   for (int i = 0; i < pfes.GetNBE(); i++)
+   {
+      const int bdr_attr = mesh.GetBdrAttribute(i);
+
+      pfes.GetBdrElementVDofs(i, vdofs);
+      for (int k = 0 ; k < vdofs.Size(); k++)
+         printf("%d\t", vdofs[k]);
+
+      tr = mesh.GetBdrFaceTransformations(i);
+      if (tr != NULL)
+      {
+         pfes.GetElementVDofs(tr->Elem1No, vdofs);
+         printf(" / ");
+         for (int k = 0 ; k < vdofs.Size(); k++)
+         {
+            printf("%d\t", vdofs[k]);
+            pdofs[vdofs[k]] += 1;
+            pbattr(vdofs[k], bdr_attr-1) = 1;
+         }
+      }
+      printf("\n");
+   }
+
+   for (int k = 0; k < pdofs.Size(); k++)
+   {
+      printf("pdofs[%d] = %d / ", k, pdofs[k]);
+      int nbattr = 0;
+      for (int d = 0; d < max_battr; d++)
+      {
+         printf("%d ", pbattr(k, d));
+         nbattr += pbattr(k, d);
+      }
+      printf(": %d\n", nbattr);
+   }
+
+   std::map<DofTag, Array<int>> pdofbattr;
+   for (int k = 0; k < pfes.GetVSize(); k++)
+   {
+      Array<int> battrs(0);
+      for (int d = 0; d < max_battr; d++)
+      {
+         if (pbattr(k, d) > 0)
+            battrs.Append(d+1);
+      }
+      battrs.Sort();
+
+      DofTag pdtag(battrs);
+      if (!pdofbattr.count(pdtag))
+         pdofbattr[pdtag] = Array<int>(0);
+
+      pdofbattr[pdtag].Append(k);
+   }
+
+   for(std::map<DofTag, Array<int>>::iterator iter = pdofbattr.begin(); iter != pdofbattr.end(); ++iter)
+   {
+      const DofTag *key = &(iter->first);
+      const Array<int> *val = &(iter->second);
+
+      printf("%s: ", key->print().c_str());
+      for (int k = 0 ; k < val->Size(); k++)
+      {
+         printf("%d ", (*val)[k]);
+      }
+      printf("\n");
+   }
+
+   return 0;
+}
\ No newline at end of file
diff --git a/sketches/domain_decomposition.cpp b/sketches/domain_decomposition.cpp
index 4def42cd..bc1462f6 100644
--- a/sketches/domain_decomposition.cpp
+++ b/sketches/domain_decomposition.cpp
@@ -132,7 +132,7 @@ int main(int argc, char* argv[])
     for (int m = 0; m < 1; m++){
       printf("Mesh %d\n", m);
       for (int i = 0; i < meshes[m]->GetNBE(); i++) {
-        int fn = meshes[m]->GetBdrFace(i);
+        int fn = meshes[m]->GetBdrElementFaceIndex(i);
         printf("Boundary Element %d : Face %d\n", i, fn);
       }
       // for (int i = 0; i < meshes[m]->GetNBE(); i++) {
@@ -146,7 +146,7 @@ int main(int argc, char* argv[])
       //   printf("Boundary element %d - Adjacent element %d\n", i, elem_id);
       //   printf("Face index: %d, face orientation : %d\n", face_info / 64, face_info % 64);
       //
-      //   int fn = meshes[m]->GetBdrFace(i);
+      //   int fn = meshes[m]->GetBdrElementFaceIndex(i);
       //   int face_inf[2];
       //   meshes[m]->GetFaceInfos(fn, &face_inf[0], &face_inf[1]);
       //   printf("From faces_info - Face index: %d, face orientation : %d\n", face_inf[0] / 64, face_inf[0] % 64);
diff --git a/sketches/generate_interface.cpp b/sketches/generate_interface.cpp
index 314d5cfb..e06e7341 100644
--- a/sketches/generate_interface.cpp
+++ b/sketches/generate_interface.cpp
@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
          for (int be = 0; be < meshes[m]->GetNBE(); be++)
          {
             int battr = meshes[m]->GetBdrAttribute(be);
-            int face = meshes[m]->GetBdrFace(be);
+            int face = meshes[m]->GetBdrElementFaceIndex(be);
             Array<int> vtx;
             // comp.GetFaceVertices(face, vtx);
             meshes[m]->GetBdrElementVertices(be, vtx);
@@ -115,7 +115,7 @@ int main(int argc, char *argv[])
       for (int be = 0; be < comp.GetNBE(); be++)
       {
          int battr = comp.GetBdrAttribute(be);
-         int face = comp.GetBdrFace(be);
+         int face = comp.GetBdrElementFaceIndex(be);
          Array<int> vtx;
          // comp.GetFaceVertices(face, vtx);
          comp.GetBdrElementVertices(be, vtx);
@@ -168,7 +168,7 @@ int main(int argc, char *argv[])
       // These will be stored in a 'interface pair' file, with the comp/battr pair above.
       Array<std::map<int,int>*> vtx_2to1(n_interface_pairs);
       // Array<std::map<int,int>*> be_2to1(n_interface_pairs);
-      Array<Array<std::pair<int,int>>*> be_pairs(n_interface_pairs);
+      // Array<Array<std::pair<int,int>>*> be_pairs(n_interface_pairs);
       for (int k = 0; k < n_interface_pairs; k++)
       {
          int bidx1 = comp.bdr_attributes.Find(if_battr1[k]);
@@ -177,18 +177,18 @@ int main(int argc, char *argv[])
          assert(bdr_vtx[bidx1]->Size() == bdr_vtx[bidx2]->Size());
 
          vtx_2to1[k] = new std::map<int,int>;
-         be_pairs[k] = new Array<std::pair<int,int>>(0);
+         // be_pairs[k] = new Array<std::pair<int,int>>(0);
       }
 
       (*vtx_2to1[0])[1] = 3;
       (*vtx_2to1[0])[0] = 2;
       // (*be_pairs[0])[0] = 1;
-      be_pairs[0]->Append({1, 0});
+      // be_pairs[0]->Append({1, 0});
 
       (*vtx_2to1[1])[0] = 1;
       (*vtx_2to1[1])[2] = 3;
       // (*be_pairs[0])[2] = 3;
-      be_pairs[1]->Append({3, 2});
+      // be_pairs[1]->Append({3, 2});
 
       // These are the informations that will be stored in the 'global topology' file.
       Array<int> global_comps(topol_data.numSub);
@@ -252,8 +252,8 @@ int main(int argc, char *argv[])
       //          tmp.BE2 = (*be_pair)[be].second;
 
       //          // use the face index from each component mesh.
-      //          int f1 = meshes[tmp.Mesh1]->GetBdrFace(tmp.BE1);
-      //          int f2 = meshes[tmp.Mesh2]->GetBdrFace(tmp.BE2);
+      //          int f1 = meshes[tmp.Mesh1]->GetBdrElementFaceIndex(tmp.BE1);
+      //          int f2 = meshes[tmp.Mesh2]->GetBdrElementFaceIndex(tmp.BE2);
       //          int Inf1, Inf2, dump;
       //          meshes[tmp.Mesh1]->GetFaceInfos(f1, &Inf1, &dump);
       //          meshes[tmp.Mesh2]->GetFaceInfos(f2, &Inf2, &dump);
diff --git a/sketches/multidomain.cpp b/sketches/multidomain.cpp
index 71aa43b8..171e2779 100644
--- a/sketches/multidomain.cpp
+++ b/sketches/multidomain.cpp
@@ -93,18 +93,21 @@ int main(int argc, char *argv[])
    // interface attribute starts after the parent mesh boundary attributes.
    int if_attr = mesh.bdr_attributes.Max() + 1;
    Array<InterfaceInfo> interface_infos(0);
+
+   // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes.
+   const int generated_battr = mesh.bdr_attributes.Max() + 1;
    
    for (int i = 0; i < numSub; i++) {
       // printf("Submesh %d\n", i);
       for (int ib = 0; ib < submeshes[i]->GetNBE(); ib++) {
-         if (submeshes[i]->GetBdrAttribute(ib) != SubMesh::GENERATED_ATTRIBUTE) continue;
+         if (submeshes[i]->GetBdrAttribute(ib) != generated_battr) continue;
 
-         int parent_face_i = (*parent_face_map_2d[i])[submeshes[i]->GetBdrFace(ib)];
+         int parent_face_i = (*parent_face_map_2d[i])[submeshes[i]->GetBdrElementFaceIndex(ib)];
          for (int j = i+1; j < numSub; j++) {
             for (int jb = 0; jb < submeshes[j]->GetNBE(); jb++) {
-               int parent_face_j = (*parent_face_map_2d[j])[submeshes[j]->GetBdrFace(jb)];
+               int parent_face_j = (*parent_face_map_2d[j])[submeshes[j]->GetBdrElementFaceIndex(jb)];
                if (parent_face_i == parent_face_j) {
-                  MFEM_ASSERT(submeshes[j]->GetBdrAttribute(jb) == SubMesh::GENERATED_ATTRIBUTE,
+                  MFEM_ASSERT(submeshes[j]->GetBdrAttribute(jb) == generated_battr,
                               "This interface element has been already set!");
                   if (interface_attributes[i][j] <= 0) {
                      interface_attributes[i][j] = if_attr;
@@ -169,15 +172,15 @@ int main(int argc, char *argv[])
          int interface_attr = submeshes[i]->GetBdrAttribute(ib);
          if (interface_attr <= mesh.bdr_attributes.Max()) continue;
 
-         int parent_face_i = (*parent_face_map_2d[i])[submeshes[i]->GetBdrFace(ib)];
+         int parent_face_i = (*parent_face_map_2d[i])[submeshes[i]->GetBdrElementFaceIndex(ib)];
          
          for (int j = 0; j < numSub; j++) {
             if (i == j) continue;
             for (int jb = 0; jb < submeshes[j]->GetNBE(); jb++) {
-               int parent_face_j = (*parent_face_map_2d[j])[submeshes[j]->GetBdrFace(jb)];
+               int parent_face_j = (*parent_face_map_2d[j])[submeshes[j]->GetBdrElementFaceIndex(jb)];
                if (parent_face_i == parent_face_j) {
                   printf("(BE %d, face %d) - parent face %d, attr %d - Submesh %d (BE %d, face %d)\n",
-                         ib, submeshes[i]->GetBdrFace(ib), parent_face_i, interface_attr, j, jb, submeshes[j]->GetBdrFace(jb));
+                         ib, submeshes[i]->GetBdrElementFaceIndex(ib), parent_face_i, interface_attr, j, jb, submeshes[j]->GetBdrElementFaceIndex(jb));
                }
             }
          }
@@ -247,6 +250,9 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
    if (parent_face_map == NULL)
       parent_face_map = new Array<int>(BuildFaceMap2D(pm, sm));
 
+   // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes.
+   const int generated_battr = pm.bdr_attributes.Max() + 1;
+
    // Setting boundary element attribute of submesh for 2D.
    // This does not support 2D.
    // Array<int> parent_face_to_be = mesh.GetFaceToBdrElMap();
@@ -254,10 +260,10 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
    parent_face_to_be = -1;
    for (int i = 0; i < pm.GetNBE(); i++)
    {
-      parent_face_to_be[pm.GetBdrElementEdgeIndex(i)] = i;
+      parent_face_to_be[pm.GetBdrElementFaceIndex(i)] = i;
    }
    for (int k = 0; k < sm.GetNBE(); k++) {
-      int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrFace(k)]];
+      int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrElementFaceIndex(k)]];
       if (pbeid != -1)
       {
          int attr = pm.GetBdrElement(pbeid)->GetAttribute();
@@ -268,7 +274,7 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
          // This case happens when a domain is extracted, but the root parent
          // mesh didn't have a boundary element on the surface that defined
          // it's boundary. It still creates a valid mesh, so we allow it.
-         sm.GetBdrElement(k)->SetAttribute(SubMesh::GENERATED_ATTRIBUTE);
+         sm.GetBdrElement(k)->SetAttribute(generated_battr);
       }
    }
 
diff --git a/sketches/poisson_dg.cpp b/sketches/poisson_dg.cpp
index 947d0a3a..3f0b7169 100644
--- a/sketches/poisson_dg.cpp
+++ b/sketches/poisson_dg.cpp
@@ -257,20 +257,23 @@ int main(int argc, char *argv[])
     Array2D<int> interface_attributes(numSub, numSub);
     interface_attributes = -1;
 
+    // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes.
+    const int generated_battr = mesh.bdr_attributes.Max() + 1;
+
     // interface attribute starts after the parent mesh boundary attributes.
     Array<InterfaceInfo> interface_infos(0);
     int if_attr = mesh.bdr_attributes.Max() + 1;
     for (int i = 0; i < numSub; i++) {
         // printf("Submesh %d\n", i);
         for (int ib = 0; ib < meshes[i]->GetNBE(); ib++) {
-          if (meshes[i]->GetBdrAttribute(ib) != SubMesh::GENERATED_ATTRIBUTE) continue;
+          if (meshes[i]->GetBdrAttribute(ib) != generated_battr) continue;
 
-          int parent_face_i = (*parent_face_map_2d[i])[meshes[i]->GetBdrFace(ib)];
+          int parent_face_i = (*parent_face_map_2d[i])[meshes[i]->GetBdrElementFaceIndex(ib)];
           for (int j = i+1; j < numSub; j++) {
               for (int jb = 0; jb < meshes[j]->GetNBE(); jb++) {
-                int parent_face_j = (*parent_face_map_2d[j])[meshes[j]->GetBdrFace(jb)];
+                int parent_face_j = (*parent_face_map_2d[j])[meshes[j]->GetBdrElementFaceIndex(jb)];
                 if (parent_face_i == parent_face_j) {
-                    MFEM_ASSERT(meshes[j]->GetBdrAttribute(jb) == SubMesh::GENERATED_ATTRIBUTE,
+                    MFEM_ASSERT(meshes[j]->GetBdrAttribute(jb) == generated_battr,
                                 "This interface element has been already set!");
                     if (interface_attributes[i][j] <= 0) {
                       interface_attributes[i][j] = if_attr;
@@ -323,16 +326,16 @@ int main(int argc, char *argv[])
           int interface_attr = meshes[i]->GetBdrAttribute(ib);
           if (interface_attr <= mesh.bdr_attributes.Max()) continue;
 
-          int parent_face_i = (*parent_face_map_2d[i])[meshes[i]->GetBdrFace(ib)];
+          int parent_face_i = (*parent_face_map_2d[i])[meshes[i]->GetBdrElementFaceIndex(ib)];
           
           for (int j = 0; j < numSub; j++) {
               if (i == j) continue;
               for (int jb = 0; jb < meshes[j]->GetNBE(); jb++) {
-                int parent_face_j = (*parent_face_map_2d[j])[meshes[j]->GetBdrFace(jb)];
+                int parent_face_j = (*parent_face_map_2d[j])[meshes[j]->GetBdrElementFaceIndex(jb)];
                 if (parent_face_i == parent_face_j) {
                   interface_parent.Append(parent_face_i);
                   printf("(BE %d, face %d) - parent face %d, attr %d - Submesh %d (BE %d, face %d)\n",
-                        ib, meshes[i]->GetBdrFace(ib), parent_face_i, interface_attr, j, jb, meshes[j]->GetBdrFace(jb));
+                        ib, meshes[i]->GetBdrElementFaceIndex(ib), parent_face_i, interface_attr, j, jb, meshes[j]->GetBdrElementFaceIndex(jb));
                 }
               }
           }
@@ -353,7 +356,7 @@ int main(int argc, char *argv[])
     cout << "Number of unknowns in each subdomain: " << total_num_dofs << endl;
 
     for (int i = 0; i < meshes[0]->GetNBE(); i++) {
-      int fn = meshes[0]->GetBdrFace(i);
+      int fn = meshes[0]->GetBdrElementFaceIndex(i);
       int bdrAttr = meshes[0]->GetBdrAttribute(i);
 
       int elem_id, face_info;
@@ -578,7 +581,7 @@ int main(int argc, char *argv[])
       // Correcting the local face transformation if orientation needs correction.
       {
         int faceInf1, faceInf2;
-        int face1 = mesh1->GetBdrFace(interface_infos[bn].BE1);
+        int face1 = mesh1->GetBdrElementFaceIndex(interface_infos[bn].BE1);
         mesh1->GetFaceInfos(face1, &faceInf1, &faceInf2);
         if (faceInf1 != interface_infos[bn].Inf1)
         {
@@ -594,7 +597,7 @@ int main(int argc, char *argv[])
                                             (*tr1).Loc1.Transf, interface_infos[bn].Inf1);
         }
 
-        int face2 = mesh2->GetBdrFace(interface_infos[bn].BE2);
+        int face2 = mesh2->GetBdrElementFaceIndex(interface_infos[bn].BE2);
         mesh2->GetFaceInfos(face2, &faceInf2, &faceInf1);
         if (faceInf2 != interface_infos[bn].Inf2)
         {
@@ -618,8 +621,8 @@ int main(int argc, char *argv[])
         printf("tr2\t%d\t%d\n", tr2->Elem1No, tr2->Elem2No);
         printf("\n");
 
-        int face1 = mesh1->GetBdrFace(interface_infos[bn].BE1);
-        int face2 = mesh2->GetBdrFace(interface_infos[bn].BE2);
+        int face1 = mesh1->GetBdrElementFaceIndex(interface_infos[bn].BE1);
+        int face2 = mesh2->GetBdrElementFaceIndex(interface_infos[bn].BE2);
 
         printf("mesh1 face info\n");
 
@@ -875,6 +878,9 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
   if (parent_face_map == NULL)
     parent_face_map = new Array<int>(BuildFaceMap2D(pm, sm));
 
+  // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes.
+  const int generated_battr = pm.bdr_attributes.Max() + 1;
+
   // Setting boundary element attribute of submesh for 2D.
   // This does not support 2D.
   // Array<int> parent_face_to_be = mesh.GetFaceToBdrElMap();
@@ -882,10 +888,10 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
   parent_face_to_be = -1;
   for (int i = 0; i < pm.GetNBE(); i++)
   {
-    parent_face_to_be[pm.GetBdrElementEdgeIndex(i)] = i;
+    parent_face_to_be[pm.GetBdrElementFaceIndex(i)] = i;
   }
   for (int k = 0; k < sm.GetNBE(); k++) {
-    int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrFace(k)]];
+    int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrElementFaceIndex(k)]];
     if (pbeid != -1)
     {
       int attr = pm.GetBdrElement(pbeid)->GetAttribute();
@@ -896,7 +902,7 @@ void BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, Array<int> *parent_face
       // This case happens when a domain is extracted, but the root parent
       // mesh didn't have a boundary element on the surface that defined
       // it's boundary. It still creates a valid mesh, so we allow it.
-      sm.GetBdrElement(k)->SetAttribute(SubMesh::GENERATED_ATTRIBUTE);
+      sm.GetBdrElement(k)->SetAttribute(generated_battr);
     }
   }
 
diff --git a/sketches/poisson_serial.cpp b/sketches/poisson_serial.cpp
index 3a625479..d46d0bbc 100644
--- a/sketches/poisson_serial.cpp
+++ b/sketches/poisson_serial.cpp
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
    //     printf("Boundary element %d - Adjacent element %d\n", i, elem_id);
    //     printf("Face index: %d, face orientation : %d\n", face_info / 64, face_info % 64);
    //
-   //     int fn = mesh.GetBdrFace(i);
+   //     int fn = mesh.GetBdrElementFaceIndex(i);
    //     int face_inf[2];
    //     mesh.GetFaceInfos(fn, &face_inf[0], &face_inf[1]);
    //     printf("From faces_info - Face index: %d, face orientation : %d\n", face_inf[0] / 64, face_inf[0] % 64);
@@ -250,7 +250,7 @@ int main(int argc, char *argv[])
      cout << "Number of unknowns in each subdomain: " << total_num_dofs << endl;
 
      for (int i = 0; i < meshes[0]->GetNBE(); i++) {
-       int fn = meshes[0]->GetBdrFace(i);
+       int fn = meshes[0]->GetBdrElementFaceIndex(i);
        int bdrAttr = meshes[0]->GetBdrAttribute(i);
 
        int elem_id, face_info;
diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp
index 25a2a863..693e5604 100644
--- a/src/component_topology_handler.cpp
+++ b/src/component_topology_handler.cpp
@@ -110,6 +110,71 @@ void ComponentTopologyHandler::ExportInfo(Array<Mesh*> &mesh_ptrs, TopologyData
    topol_data.global_bdr_attributes = &bdr_attributes;
 }
 
+std::map<DofTag, Array<int>> ComponentTopologyHandler::SplitDofsByBoundary(FiniteElementSpace &fes)
+{
+   Mesh *mesh = fes.GetMesh();
+
+   /* Get max boundary attributes. */
+   int max_battr = -1;
+   for (int i = 0; i < fes.GetNBE(); i++)
+      max_battr = max(max_battr, mesh->GetBdrAttribute(i));
+
+   /*
+      Table that stores boundary attributes each dof belongs to.
+      Size of (number of vdofs x max_battr).
+   */
+   Array2D<bool> dof_battr(fes.GetVSize(), max_battr);
+   dof_battr = false;
+
+   /*
+      Loop over boundary elements,
+      Log all dofs that belong to the boundary element.
+   */
+   Array<int> vdofs;
+   FaceElementTransformations *tr;
+   for (int i = 0; i < fes.GetNBE(); i++)
+   {
+      const int bdr_attr = mesh->GetBdrAttribute(i);
+
+      fes.GetBdrElementVDofs(i, vdofs);
+
+      tr = mesh->GetBdrFaceTransformations(i);
+      if (tr != NULL)
+      {
+         fes.GetElementVDofs(tr->Elem1No, vdofs);
+         for (int k = 0 ; k < vdofs.Size(); k++)
+            dof_battr(vdofs[k], bdr_attr-1) = true;
+      }
+   }
+
+   /*
+      Sort out dofs by same boundary elements.
+      Dofs can belong to multiple boundary elements,
+      representing edge, face, or vertex.
+   */
+   std::map<DofTag, Array<int>> dof_dict;
+   for (int k = 0; k < fes.GetVSize(); k++)
+   {
+      Array<int> battrs(0);
+      for (int d = 0; d < max_battr; d++)
+      {
+         if (dof_battr(k, d))
+            battrs.Append(d+1);
+      }
+      battrs.Sort();
+
+      /* Creat DofTag if it does not exist */
+      DofTag dtag(battrs);
+      if (!dof_dict.count(dtag))
+         dof_dict[dtag] = Array<int>(0);
+
+      /* Append dof to the corresponding DofTag */
+      dof_dict[dtag].Append(k);
+   }
+
+   return dof_dict;
+}
+
 void ComponentTopologyHandler::SetupComponents()
 {
    assert(num_comp > 0);
@@ -660,8 +725,8 @@ void ComponentTopologyHandler::SetupReferenceInterfaces()
          tmp.BE2 = pair[1];
 
          // use the face index from each component mesh.
-         int f1 = comp1->GetBdrFace(tmp.BE1);
-         int f2 = comp2->GetBdrFace(tmp.BE2);
+         int f1 = comp1->GetBdrElementFaceIndex(tmp.BE1);
+         int f2 = comp2->GetBdrElementFaceIndex(tmp.BE2);
          int Inf1, Inf2, dump;
          comp1->GetFaceInfos(f1, &Inf1, &dump);
          comp2->GetFaceInfos(f2, &Inf2, &dump);
diff --git a/src/topology_handler.cpp b/src/topology_handler.cpp
index 3f5670a4..e0178d52 100644
--- a/src/topology_handler.cpp
+++ b/src/topology_handler.cpp
@@ -50,6 +50,56 @@ TopologyHandler::TopologyHandler(const TopologyHandlerMode &input_type)
    }
 }
 
+std::map<DofTag, Array<int>> TopologyHandler::SplitBoundaryDofs(
+   Mesh* const mesh, const FiniteElementSpace *fes) const
+{
+   std::map<DofTag, Array<int>> battr2dofs;
+
+   // get max boundary attribute.
+   int max_battr = -1;
+   for (int i = 0; i < fes->GetNBE(); i++)
+      max_battr = max(max_battr, mesh->GetBdrAttribute(i));
+
+   // 2d array to store boundary attributes that each dof belongs to.
+   Array2D<int> dof_battr(fes->GetVSize(), max_battr);
+   dof_battr = 0;
+
+   Array<int> vdofs;
+   FaceElementTransformations *tr;
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      const int bdr_attr = mesh->GetBdrAttribute(i);
+
+      tr = mesh->GetBdrFaceTransformations(i);
+      if (tr != NULL)
+      {
+         fes->GetElementVDofs(tr->Elem1No, vdofs);
+         for (int k = 0 ; k < vdofs.Size(); k++)
+            dof_battr(vdofs[k], bdr_attr-1) = 1;
+      }
+   }
+
+   for (int k = 0; k < fes->GetVSize(); k++)
+   {
+      Array<int> battrs(0);
+      for (int d = 0; d < max_battr; d++)
+      {
+         if (dof_battr(k, d) > 0)
+         // boundary attribute starts from 1.
+            battrs.Append(d+1);
+      }
+      battrs.Sort();
+
+      DofTag pdtag(battrs);
+      if (!battr2dofs.count(pdtag))
+         battr2dofs[pdtag] = Array<int>(0);
+
+      battr2dofs[pdtag].Append(k);
+   }
+
+   return battr2dofs;
+}
+
 void TopologyHandler::GetInterfaceTransformations(Mesh *m1, Mesh *m2, const InterfaceInfo *if_info,
                                                    FaceElementTransformations* &tr1,
                                                    FaceElementTransformations* &tr2)
@@ -63,7 +113,7 @@ void TopologyHandler::GetInterfaceTransformations(Mesh *m1, Mesh *m2, const Inte
 
    // Correcting the local face1 transformation if orientation needs correction.
    int faceInf1, faceInf2;
-   int face1 = m1->GetBdrFace(if_info->BE1);
+   int face1 = m1->GetBdrElementFaceIndex(if_info->BE1);
    m1->GetFaceInfos(face1, &faceInf1, &faceInf2);
    if (faceInf1 != if_info->Inf1)
    {
@@ -80,7 +130,7 @@ void TopologyHandler::GetInterfaceTransformations(Mesh *m1, Mesh *m2, const Inte
    }
 
    // Correcting the local face1 transformation if orientation needs correction.
-   int face2 = m2->GetBdrFace(if_info->BE2);
+   int face2 = m2->GetBdrElementFaceIndex(if_info->BE2);
    m2->GetFaceInfos(face2, &faceInf2, &faceInf1);
    if (faceInf2 != if_info->Inf2)
    {
@@ -332,11 +382,11 @@ void SubMeshTopologyHandler::BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm,
    parent_face_to_be = -1;
    for (int i = 0; i < pm.GetNBE(); i++)
    {
-      parent_face_to_be[pm.GetBdrElementEdgeIndex(i)] = i;
+      parent_face_to_be[pm.GetBdrElementFaceIndex(i)] = i;
    }
    for (int k = 0; k < sm.GetNBE(); k++)
    {
-      int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrFace(k)]];
+      int pbeid = parent_face_to_be[(*parent_face_map)[sm.GetBdrElementFaceIndex(k)]];
       if (pbeid != -1)
       {
          int attr = pm.GetBdrElement(pbeid)->GetAttribute();
@@ -374,14 +424,14 @@ void SubMeshTopologyHandler::BuildInterfaceInfos()
       for (int ib = 0; ib < meshes[i]->GetNBE(); ib++)
       {
          if (meshes[i]->GetBdrAttribute(ib) != generated_battr) continue;
-         int parent_face_i = (*parent_face_map[i])[meshes[i]->GetBdrFace(ib)];
+         int parent_face_i = (*parent_face_map[i])[meshes[i]->GetBdrElementFaceIndex(ib)];
 
          // Loop over each subdomain, each boundary element, to find the match.
          for (int j = i+1; j < numSub; j++)
          {
             for (int jb = 0; jb < meshes[j]->GetNBE(); jb++)
             {
-               int parent_face_j = (*parent_face_map[j])[meshes[j]->GetBdrFace(jb)];
+               int parent_face_j = (*parent_face_map[j])[meshes[j]->GetBdrElementFaceIndex(jb)];
                if (parent_face_i != parent_face_j) continue;
 
                assert(meshes[j]->GetBdrAttribute(jb) == generated_battr);