diff --git a/data/meshes/SplitCubeTest.g b/data/meshes/SplitCubeTest.g
new file mode 100644
index 0000000000..fb1494387b
Binary files /dev/null and b/data/meshes/SplitCubeTest.g differ
diff --git a/src/smith/numerics/functional/domain.cpp b/src/smith/numerics/functional/domain.cpp
index 9acbe708c4..fecd1fd4f0 100644
--- a/src/smith/numerics/functional/domain.cpp
+++ b/src/smith/numerics/functional/domain.cpp
@@ -622,7 +622,7 @@ Domain EntireBoundary(const mesh_t& mesh)
 }
 
 /// @brief constructs a domain from all the interior face elements in a mesh
-Domain InteriorFaces(const mesh_t& mesh)
+Domain EntireInteriorFaces(const mesh_t& mesh)
 {
   Domain output{mesh, mesh.SpaceDimension() - 1, Domain::Type::InteriorFaces};
 
@@ -660,6 +660,80 @@ Domain InteriorFaces(const mesh_t& mesh)
   return output;
 }
 
+/// @brief constructs a domain from some subset of the interior face elements in a mesh
+template <int d>
+Domain domain_of_interior_faces(const mesh_t& mesh, std::function<bool(std::vector<tensor<double, d>>, int)> predicate)
+{
+  assert(mesh.SpaceDimension() == d);
+  Domain output{mesh, d - 1, Domain::Type::InteriorFaces};
+
+  mfem::Array<int> face_id_to_bdr_id = mesh.GetFaceToBdrElMap();
+  mfem::Vector vertices;
+  mesh.GetVertices(vertices);
+
+  int edge_id = 0;
+  int tri_id = 0;
+  int quad_id = 0;
+
+  for (int f = 0; f < mesh.GetNumFaces(); f++) {
+    // discard faces with the wrong type
+    if (!mesh.GetFaceInformation(f).IsInterior()) continue;
+
+    auto geom = mesh.GetFaceGeometry(f);
+
+    mfem::Array<int> vertex_ids;
+    mesh.GetFaceVertices(f, vertex_ids);
+
+    auto x = gather<d>(vertices, vertex_ids);
+
+    int bdr_id = face_id_to_bdr_id[f];
+    int attr = (bdr_id >= 0) ? mesh.GetBdrAttribute(bdr_id) : -1;
+
+    bool add = predicate(x, attr);
+
+    switch (geom) {
+      case mfem::Geometry::SEGMENT:
+        if (add) {
+          output.edge_ids_.push_back(edge_id);
+          output.mfem_edge_ids_.push_back(f);
+        }
+        edge_id++;
+        break;
+      case mfem::Geometry::TRIANGLE:
+        if (add) {
+          output.tri_ids_.push_back(tri_id);
+          output.mfem_tri_ids_.push_back(f);
+        }
+        tri_id++;
+        break;
+      case mfem::Geometry::SQUARE:
+        if (add) {
+          output.quad_ids_.push_back(quad_id);
+          output.mfem_quad_ids_.push_back(f);
+        }
+        quad_id++;
+        break;
+      default:
+        SLIC_ERROR("unsupported element type");
+        break;
+    }
+  }
+
+  output.insert_shared_interior_face_list();
+
+  return output;
+}
+
+Domain Domain::ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func)
+{
+  return domain_of_interior_faces<2>(mesh, func);
+}
+
+Domain Domain::ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func)
+{
+  return domain_of_interior_faces<3>(mesh, func);
+}
+
 ///////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/smith/numerics/functional/domain.hpp b/src/smith/numerics/functional/domain.hpp
index a7fe5747e0..e08820d209 100644
--- a/src/smith/numerics/functional/domain.hpp
+++ b/src/smith/numerics/functional/domain.hpp
@@ -167,6 +167,18 @@ struct Domain {
   /// @overload
   static Domain ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
 
+  /**
+   * @brief create a domain from some subset of the interior faces in an a mesh
+   * @param mesh the entire mesh
+   * @param func predicate function for determining which interior faces will be
+   * included in this domain. The function's arguments are the list of vertex coordinates and
+   * an attribute index (if appropriate).
+   */
+  static Domain ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
+
+  /// @overload
+  static Domain ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
+
   /// @brief get elements by geometry type
   const std::vector<int>& get(mfem::Geometry::Type geom) const
   {
@@ -268,7 +280,7 @@ Domain EntireDomain(const mesh_t& mesh);
 Domain EntireBoundary(const mesh_t& mesh);
 
 /// @brief constructs a domain from all the interior face elements in a mesh
-Domain InteriorFaces(const mesh_t& mesh);
+Domain EntireInteriorFaces(const mesh_t& mesh);
 
 /// @brief create a new domain that is the union of `a` and `b`
 Domain operator|(const Domain& a, const Domain& b);
diff --git a/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp b/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp
index d65878dddb..c25f209129 100644
--- a/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp
+++ b/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp
@@ -61,7 +61,7 @@ void L2_index_test(std::string meshfile)
   // Construct the new functional object using the specified test and trial spaces
   Functional<test_space(trial_space)> residual(test_fespace.get(), {trial_fespace.get()});
 
-  Domain interior_faces = InteriorFaces(*mesh);
+  Domain interior_faces = EntireInteriorFaces(*mesh);
 
   // Define the integral of jumps over all interior faces
   residual.AddInteriorFaceIntegral(
diff --git a/src/smith/numerics/functional/tests/dg_restriction_operators.cpp b/src/smith/numerics/functional/tests/dg_restriction_operators.cpp
index 45a66727e2..dca0198ccd 100644
--- a/src/smith/numerics/functional/tests/dg_restriction_operators.cpp
+++ b/src/smith/numerics/functional/tests/dg_restriction_operators.cpp
@@ -298,7 +298,7 @@ void parametrized_test(int permutation)
   mfem::Mesh bmesh = generate_permuted_mesh(geom, permutation);
   auto mesh = smith::mesh::refineAndDistribute(std::move(bmesh));
 
-  Domain interior_faces = InteriorFaces(*mesh);
+  Domain interior_faces = EntireInteriorFaces(*mesh);
 
   // each one of these meshes should have two elements
   // and a single "face" that separates them
diff --git a/src/smith/numerics/functional/tests/domain_tests.cpp b/src/smith/numerics/functional/tests/domain_tests.cpp
index a73da97f0e..03582730ef 100644
--- a/src/smith/numerics/functional/tests/domain_tests.cpp
+++ b/src/smith/numerics/functional/tests/domain_tests.cpp
@@ -16,6 +16,10 @@
 #include "smith/numerics/functional/tensor.hpp"
 #include "smith/smith_config.hpp"
 #include "smith/mesh_utils/mesh_utils.hpp"
+#include "smith/numerics/functional/functional.hpp"
+#include "smith/numerics/functional/shape_aware_functional.hpp"
+#include "smith/numerics/functional/tensor.hpp"
+#include "smith/physics/state/finite_element_state.hpp"
 
 using namespace smith;
 
@@ -36,6 +40,14 @@ mfem::Mesh import_mesh(std::string meshfile)
   return mesh;
 }
 
+struct IdentityFunctor {
+  template <typename Arg1, typename Arg2>
+  SMITH_HOST_DEVICE auto operator()(Arg1, Arg2) const
+  {
+    return 1.0;
+  }
+};
+
 TEST(domain, of_edges)
 {
   {
@@ -492,6 +504,29 @@ TEST(domain, of3dBoundaryElementsFindsDofs)
   EXPECT_EQ(dof_indices.Size(), 15);
 }
 
+TEST(domain, ofInteriorBoundaryElements)
+{
+  using test_space = double;
+  using trial_space = H1<1>;
+
+  auto expectedArea = 1.0;
+  auto meshInput = buildMeshFromFile(SMITH_REPO_DIR "/data/meshes/SplitCubeTest.g");
+  auto mesh = smith::mesh::refineAndDistribute(std::move(meshInput), 0, 0);
+
+  Domain d0 = Domain::ofInteriorFaces(*mesh, by_attr<3>(2));
+
+  auto [trial_fespace, trial_fec] = smith::generateParFiniteElementSpace<trial_space>(mesh.get());
+  mfem::Vector U(trial_fespace->TrueVSize());
+
+  // Calculate the area of the internal boundary region
+  Functional<test_space(trial_space)> totalArea({trial_fespace.get()});
+
+  totalArea.AddInteriorFaceIntegral(smith::Dimension<2>{}, smith::DependsOn<>{}, IdentityFunctor{}, d0);
+  double calculatedArea = totalArea(0.0, U);
+
+  EXPECT_NEAR(calculatedArea, expectedArea, 1e-6);
+}
+
 int main(int argc, char* argv[])
 {
   ::testing::InitGoogleTest(&argc, argv);
diff --git a/src/smith/numerics/functional/tests/element_restriction_tests.cpp b/src/smith/numerics/functional/tests/element_restriction_tests.cpp
index b654093e88..0b17b296c7 100644
--- a/src/smith/numerics/functional/tests/element_restriction_tests.cpp
+++ b/src/smith/numerics/functional/tests/element_restriction_tests.cpp
@@ -82,7 +82,7 @@ TEST(patch_test_meshes, triangle_domains)
     EXPECT_EQ(L2_BER.ESize(), 4 * 3);
   }
 
-  Domain interior = InteriorFaces(*mesh);
+  Domain interior = EntireInteriorFaces(*mesh);
   EXPECT_EQ(interior.mfem_edge_ids_.size(), 4);
   EXPECT_EQ(interior.mfem_tri_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_quad_ids_.size(), 0);
@@ -152,7 +152,7 @@ TEST(patch_test_meshes, quadrilateral_domains)
     EXPECT_EQ(L2_BER.ESize(), 4 * 3);
   }
 
-  Domain interior = InteriorFaces(*mesh);
+  Domain interior = EntireInteriorFaces(*mesh);
   EXPECT_EQ(interior.mfem_edge_ids_.size(), 8);
   EXPECT_EQ(interior.mfem_tri_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_quad_ids_.size(), 0);
@@ -222,7 +222,7 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains)
     EXPECT_EQ(L2_BER.ESize(), 4 * 3);
   }
 
-  Domain interior = InteriorFaces(*mesh);
+  Domain interior = EntireInteriorFaces(*mesh);
   EXPECT_EQ(interior.mfem_edge_ids_.size(), 9);
   EXPECT_EQ(interior.mfem_tri_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_quad_ids_.size(), 0);
@@ -292,7 +292,7 @@ TEST(patch_test_meshes, tetrahedron_domains)
     EXPECT_EQ(L2_BER.ESize(), 12 * 6);
   }
 
-  Domain interior = InteriorFaces(*mesh);
+  Domain interior = EntireInteriorFaces(*mesh);
   EXPECT_EQ(interior.mfem_edge_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_tri_ids_.size(), 18);
   EXPECT_EQ(interior.mfem_quad_ids_.size(), 0);
@@ -362,7 +362,7 @@ TEST(patch_test_meshes, hexahedron_domains)
     EXPECT_EQ(L2_BER.ESize(), 6 * 9);
   }
 
-  Domain interior = InteriorFaces(*mesh);
+  Domain interior = EntireInteriorFaces(*mesh);
   EXPECT_EQ(interior.mfem_edge_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_tri_ids_.size(), 0);
   EXPECT_EQ(interior.mfem_quad_ids_.size(), 18);
diff --git a/src/smith/numerics/functional/tests/functional_basic_dg.cpp b/src/smith/numerics/functional/tests/functional_basic_dg.cpp
index 9adabb1569..d37ddb10de 100644
--- a/src/smith/numerics/functional/tests/functional_basic_dg.cpp
+++ b/src/smith/numerics/functional/tests/functional_basic_dg.cpp
@@ -42,7 +42,7 @@ void L2_test(std::string meshfile)
   // Construct the new functional object using the specified test and trial spaces
   Functional<test_space(trial_space)> residual(&fespace, {&fespace});
 
-  Domain interior_faces = InteriorFaces(*mesh);
+  Domain interior_faces = EntireInteriorFaces(*mesh);
 
   residual.AddInteriorFaceIntegral(
       Dimension<dim - 1>{}, DependsOn<0>{},
@@ -101,7 +101,7 @@ void L2_qoi_test(std::string meshfile)
   // Construct the new functional object using the specified test and trial spaces
   Functional<double(trial_space)> qoi({&fespace});
 
-  Domain interior_faces = InteriorFaces(*mesh);
+  Domain interior_faces = EntireInteriorFaces(*mesh);
 
   qoi.AddInteriorFaceIntegral(
       Dimension<dim - 1>{}, DependsOn<0>{},
@@ -168,7 +168,7 @@ void L2_scalar_valued_test(std::string meshfile)
   // Construct the new functional object using the specified test and trial spaces
   Functional<test_space(trial_space_0, trial_space_1)> residual(&fespace_0, {&fespace_0, &fespace_1});
 
-  Domain interior_faces = InteriorFaces(*mesh);
+  Domain interior_faces = EntireInteriorFaces(*mesh);
 
   residual.AddInteriorFaceIntegral(
       Dimension<dim - 1>{}, DependsOn<0, 1>{},
diff --git a/src/smith/physics/mesh.cpp b/src/smith/physics/mesh.cpp
index 83ed0b5257..ebdc2a6174 100644
--- a/src/smith/physics/mesh.cpp
+++ b/src/smith/physics/mesh.cpp
@@ -82,7 +82,7 @@ void Mesh::createDomains()
 {
   domains_.insert({entireBodyName(), smith::EntireDomain(*mfem_mesh_)});
   domains_.insert({entireBoundaryName(), smith::EntireBoundary(*mfem_mesh_)});
-  domains_.insert({internalBoundaryName(), smith::InteriorFaces(*mfem_mesh_)});
+  domains_.insert({internalBoundaryName(), smith::EntireInteriorFaces(*mfem_mesh_)});
 }
 
 void Mesh::errorIfDomainExists(const std::string& domain_name) const