diff --git a/cpp/dolfinx/common/sort.h b/cpp/dolfinx/common/sort.h
index 4ee49fd9360..6aad81ebc8d 100644
--- a/cpp/dolfinx/common/sort.h
+++ b/cpp/dolfinx/common/sort.h
@@ -67,7 +67,7 @@ struct __radix_sort
   /// @code
   /// std::array<std::int16_t, 3> a{2, 3, 1};
   /// std::array<std::int16_t, 3> i{0, 1, 2};
-  /// dolfixn::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0,
+  /// dolfinx::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0,
   /// 1} and a[i] = {1, 2, 3};
   /// @endcode
   /// @tparam R Type of range to be sorted.
diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h
index 8a361dcf0c8..6bbb426e5c2 100644
--- a/cpp/dolfinx/fem/Form.h
+++ b/cpp/dolfinx/fem/Form.h
@@ -38,9 +38,10 @@ class Function;
 enum class IntegralType : std::int8_t
 {
   cell = 0,           ///< Cell
-  exterior_facet = 1, ///< Exterior facet
+  exterior_facet = 1, ///< Facet
   interior_facet = 2, ///< Interior facet
-  vertex = 3          ///< Vertex
+  vertex = 3,         ///< Vertex
+  ridge = 4           ///< Ridge
 };
 
 /// @brief Represents integral data, containing the kernel, and a list
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index 90351a05b7a..3b6148f9e97 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/cpp/dolfinx/fem/assemble_matrix_impl.h
@@ -156,16 +156,24 @@ void assemble_cells_matrix(
   }
 }
 
-/// @brief Execute kernel over exterior facets and accumulate result in
-/// a matrix.
+/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result
+/// in a matrix.
+///
+/// Each entity is represented by (i) a cell that the entity is attached to
+/// and (ii) the local index of the entity  with respect to the cell. The
+/// kernel is executed for each entity. The kernel can access data
+/// (e.g., coefficients, basis functions) associated with the attached cell.
+/// However, entities may be attached to more than one cell. This function
+/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen
+/// from cell used to define the entity.
 ///
 /// @tparam T Matrix/form scalar type.
 /// @param[in] mat_set Function that accumulates computed entries into a
 /// matrix.
 /// @param[in] x_dofmap Dofmap for the mesh geometry.
 /// @param[in] x Mesh geometry (coordinates).
-/// @param[in] facets Facet indices (in the integration domain mesh) to
-/// execute the kernel over.
+/// @param[in] entities Integration entities (in the integration domain mesh) to
+/// execute the kernel over. These are pairs (cell, local entity index)
 /// @param[in] dofmap0 Test function (row) degree-of-freedom data
 /// holding the (0) dofmap, (1) dofmap block size and (2) dofmap cell
 /// indices.
@@ -188,17 +196,17 @@ void assemble_cells_matrix(
 /// function mesh.
 /// @param[in] cell_info1 Cell permutation information for the trial
 /// function mesh.
-/// @param[in] perms Facet permutation integer. Empty if facet
+/// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
 template <dolfinx::scalar T>
-void assemble_exterior_facets(
+void assemble_entities(
     la::MatSet<T> auto mat_set, mdspan2_t x_dofmap,
     md::mdspan<const scalar_value_t<T>,
                md::extents<std::size_t, md::dynamic_extent, 3>>
         x,
     md::mdspan<const std::int32_t,
                std::extents<std::size_t, md::dynamic_extent, 2>>
-        facets,
+        entities,
     std::tuple<mdspan2_t, int,
                md::mdspan<const std::int32_t,
                           std::extents<std::size_t, md::dynamic_extent, 2>>>
@@ -215,11 +223,11 @@ void assemble_exterior_facets(
     std::span<const std::uint32_t> cell_info1,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
 {
-  if (facets.empty())
+  if (entities.empty())
     return;
 
-  const auto [dmap0, bs0, facets0] = dofmap0;
-  const auto [dmap1, bs1, facets1] = dofmap1;
+  const auto [dmap0, bs0, entities0] = dofmap0;
+  const auto [dmap1, bs1, entities1] = dofmap1;
 
   // Data structures used in assembly
   std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
@@ -228,17 +236,17 @@ void assemble_exterior_facets(
   const int ndim0 = bs0 * num_dofs0;
   const int ndim1 = bs1 * num_dofs1;
   std::vector<T> Ae(ndim0 * ndim1);
-  assert(facets0.size() == facets.size());
-  assert(facets1.size() == facets.size());
-  for (std::size_t f = 0; f < facets.extent(0); ++f)
+  assert(entities0.size() == entities.size());
+  assert(entities1.size() == entities.size());
+  for (std::size_t f = 0; f < entities.extent(0); ++f)
   {
-    // Cell in the integration domain, local facet index relative to the
+    // Cell in the integration domain, local entity index relative to the
     // integration domain cell, and cells in the test and trial function
     // meshes
-    std::int32_t cell = facets(f, 0);
-    std::int32_t local_facet = facets(f, 1);
-    std::int32_t cell0 = facets0(f, 0);
-    std::int32_t cell1 = facets1(f, 0);
+    std::int32_t cell = entities(f, 0);
+    std::int32_t local_entity = entities(f, 1);
+    std::int32_t cell0 = entities0(f, 0);
+    std::int32_t cell1 = entities1(f, 0);
 
     // Get cell coordinates/geometry
     auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
@@ -246,12 +254,12 @@ void assemble_exterior_facets(
       std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
 
     // Permutations
-    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
+    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
 
     // Tabulate tensor
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
-           &local_facet, &perm, nullptr);
+           &local_entity, &perm, nullptr);
     P0(Ae, cell_info0, cell0, ndim1);
     P1T(Ae, cell_info1, cell1, ndim0);
 
@@ -605,7 +613,7 @@ void assemble_matrix(
           cell_info0, cell_info1);
     }
 
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
     if (a.needs_facet_permutations())
     {
       mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
@@ -614,40 +622,8 @@ void assemble_matrix(
       mesh->topology_mutable()->create_entity_permutations();
       const std::vector<std::uint8_t>& p
           = mesh->topology()->get_facet_permutations();
-      perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
-                         num_facets_per_cell);
-    }
-
-    for (int i = 0;
-         i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i)
-    {
-      if (num_cell_types > 1)
-      {
-        throw std::runtime_error("Exterior facet integrals with mixed "
-                                 "topology aren't supported yet");
-      }
-
-      using mdspanx2_t
-          = md::mdspan<const std::int32_t,
-                       md::extents<std::size_t, md::dynamic_extent, 2>>;
-
-      auto fn = a.kernel(IntegralType::exterior_facet, i, 0);
-      assert(fn);
-      auto& [coeffs, cstride]
-          = coefficients.at({IntegralType::exterior_facet, i});
-
-      std::span f = a.domain(IntegralType::exterior_facet, i, 0);
-      mdspanx2_t facets(f.data(), f.size() / 2, 2);
-      std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0);
-      mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
-      std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0);
-      mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
-      assert((facets.size() / 2) * cstride == coeffs.size());
-      impl::assemble_exterior_facets(
-          mat_set, x_dofmap, x, facets, {dofs0, bs0, facets0}, P0,
-          {dofs1, bs1, facets1}, P1T, bc0, bc1, fn,
-          md::mdspan(coeffs.data(), facets.extent(0), cstride), constants,
-          cell_info0, cell_info1, perms);
+      facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
+                               num_facets_per_cell);
     }
 
     for (int i = 0;
@@ -685,7 +661,47 @@ void assemble_matrix(
            mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
           P1T, bc0, bc1, fn,
           mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), constants,
-          cell_info0, cell_info1, perms);
+          cell_info0, cell_info1, facet_perms);
+    }
+
+    for (auto itg_type : {fem::IntegralType::exterior_facet,
+                          fem::IntegralType::vertex, fem::IntegralType::ridge})
+    {
+      md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
+          = itg_type == fem::IntegralType::exterior_facet
+                ? facet_perms
+                : md::mdspan<const std::uint8_t,
+                             md::dextents<std::size_t, 2>>{};
+
+      for (int i = 0; i < a.num_integrals(itg_type, cell_type_idx); ++i)
+      {
+        if (num_cell_types > 1)
+        {
+          throw std::runtime_error("Exterior facet integrals with mixed "
+                                   "topology aren't supported yet");
+        }
+
+        using mdspanx2_t
+            = md::mdspan<const std::int32_t,
+                         md::extents<std::size_t, md::dynamic_extent, 2>>;
+
+        auto fn = a.kernel(itg_type, i, 0);
+        assert(fn);
+        auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+
+        std::span e = a.domain(itg_type, i, 0);
+        mdspanx2_t entities(e.data(), e.size() / 2, 2);
+        std::span e0 = a.domain_arg(itg_type, 0, i, 0);
+        mdspanx2_t entities0(e0.data(), e0.size() / 2, 2);
+        std::span e1 = a.domain_arg(itg_type, 1, i, 0);
+        mdspanx2_t entities1(e1.data(), e1.size() / 2, 2);
+        assert((entities.size() / 2) * cstride == coeffs.size());
+        impl::assemble_entities(
+            mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0,
+            {dofs1, bs1, entities1}, P1T, bc0, bc1, fn,
+            md::mdspan(coeffs.data(), entities.extent(0), cstride), constants,
+            cell_info0, cell_info1, perms);
+      }
     }
   }
 }
diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h
index 1ade0a9792e..37166b46ea6 100644
--- a/cpp/dolfinx/fem/assemble_scalar_impl.h
+++ b/cpp/dolfinx/fem/assemble_scalar_impl.h
@@ -55,32 +55,41 @@ T assemble_cells(mdspan2_t x_dofmap,
   return value;
 }
 
-/// Execute kernel over exterior facets and accumulate result
+/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result
+/// in a scalar.
+///
+/// Each entity is represented by (i) a cell that the entity is attached to
+/// and (ii) the local index of the entity  with respect to the cell. The
+/// kernel is executed for each entity. The kernel can access data
+/// (e.g., coefficients, basis functions) associated with the attached cell.
+/// However, entities may be attached to more than one cell. This function
+/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen
+/// from cell used to define the entity.
 template <dolfinx::scalar T>
-T assemble_exterior_facets(
+T assemble_entities(
     mdspan2_t x_dofmap,
     md::mdspan<const scalar_value_t<T>,
                md::extents<std::size_t, md::dynamic_extent, 3>>
         x,
     md::mdspan<const std::int32_t,
                md::extents<std::size_t, md::dynamic_extent, 2>>
-        facets,
+        entities,
     FEkernel<T> auto fn, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
     std::span<scalar_value_t<T>> cdofs_b)
 {
   T value(0);
-  if (facets.empty())
+  if (entities.empty())
     return value;
 
   assert(cdofs_b.size() >= 3 * x_dofmap.extent(1));
 
   // Iterate over all facets
-  for (std::size_t f = 0; f < facets.extent(0); ++f)
+  for (std::size_t f = 0; f < entities.extent(0); ++f)
   {
-    std::int32_t cell = facets(f, 0);
-    std::int32_t local_facet = facets(f, 1);
+    std::int32_t cell = entities(f, 0);
+    std::int32_t local_entity = entities(f, 1);
 
     // Get cell coordinates/geometry
     auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
@@ -88,8 +97,8 @@ T assemble_exterior_facets(
       std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs_b.begin(), 3 * i));
 
     // Permutations
-    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
-    fn(&value, &coeffs(f, 0), constants.data(), cdofs_b.data(), &local_facet,
+    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
+    fn(&value, &coeffs(f, 0), constants.data(), cdofs_b.data(), &local_entity,
        &perm, nullptr);
   }
 
@@ -147,43 +156,6 @@ T assemble_interior_facets(
   return value;
 }
 
-/// Assemble functional over vertices
-template <dolfinx::scalar T>
-T assemble_vertices(mdspan2_t x_dofmap,
-                    md::mdspan<const scalar_value_t<T>,
-                               md::extents<std::size_t, md::dynamic_extent, 3>>
-                        x,
-                    md::mdspan<const std::int32_t,
-                               md::extents<std::size_t, md::dynamic_extent, 2>>
-                        vertices,
-                    FEkernel<T> auto fn, std::span<const T> constants,
-                    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
-                    std::span<scalar_value_t<T>> cdofs_b)
-{
-  T value(0);
-  if (vertices.empty())
-    return value;
-
-  assert(cdofs_b.size() >= 3 * x_dofmap.extent(1));
-
-  // Iterate over all cells
-  for (std::size_t index = 0; index < vertices.extent(0); ++index)
-  {
-    std::int32_t cell = vertices(index, 0);
-    std::int32_t local_vertex_index = vertices(index, 1);
-
-    // Get cell coordinates/geometry
-    auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
-    for (std::size_t i = 0; i < x_dofs.size(); ++i)
-      std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs_b.begin(), 3 * i));
-
-    fn(&value, &coeffs(index, 0), constants.data(), cdofs_b.data(),
-       &local_vertex_index, nullptr, nullptr);
-  }
-
-  return value;
-}
-
 /// Assemble functional into an scalar with provided mesh geometry.
 template <dolfinx::scalar T, std::floating_point U>
 T assemble_scalar(
@@ -216,39 +188,14 @@ T assemble_scalar(
   mesh::CellType cell_type = mesh->topology()->cell_type();
   int num_facets_per_cell
       = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
-  md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
+  md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
   if (M.needs_facet_permutations())
   {
     mesh->topology_mutable()->create_entity_permutations();
     const std::vector<std::uint8_t>& p
         = mesh->topology()->get_facet_permutations();
-    perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
-                       num_facets_per_cell);
-  }
-
-  for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i)
-  {
-    auto fn = M.kernel(IntegralType::exterior_facet, i, 0);
-    assert(fn);
-    auto& [coeffs, cstride]
-        = coefficients.at({IntegralType::exterior_facet, i});
-
-    std::span facets = M.domain(IntegralType::exterior_facet, i, 0);
-
-    // Two values per each adjacent cell (cell index and local facet
-    // index)
-    constexpr std::size_t num_adjacent_cells = 1;
-    constexpr std::size_t shape1 = 2 * num_adjacent_cells;
-
-    assert((facets.size() / 2) * cstride == coeffs.size());
-    value += impl::assemble_exterior_facets(
-        x_dofmap, x,
-        md::mdspan<const std::int32_t,
-                   md::extents<std::size_t, md::dynamic_extent, 2>>(
-            facets.data(), facets.size() / shape1, 2),
-        fn, constants,
-        md::mdspan(coeffs.data(), facets.size() / shape1, cstride), perms,
-        cdofs_b);
+    facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
+                             num_facets_per_cell);
   }
 
   for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i)
@@ -273,31 +220,36 @@ T assemble_scalar(
         md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                         md::dynamic_extent>>(
             coeffs.data(), facets.size() / shape1, 2, cstride),
-        perms, cdofs_b);
+        facet_perms, cdofs_b);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i)
+  for (auto itg_type : {fem::IntegralType::exterior_facet,
+                        fem::IntegralType::vertex, fem::IntegralType::ridge})
   {
-    auto fn = M.kernel(IntegralType::vertex, i, 0);
-    assert(fn);
-
-    auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
-
-    std::span<const std::int32_t> vertices
-        = M.domain(IntegralType::vertex, i, 0);
-    assert(vertices.size() * cstride == coeffs.size());
-
-    constexpr std::size_t num_adjacent_cells = 1;
-    // Two values per adj. cell (cell index and local vertex index).
-    constexpr std::size_t shape1 = 2 * num_adjacent_cells;
-
-    value += impl::assemble_vertices(
-        x_dofmap, x,
-        md::mdspan<const std::int32_t,
-                   md::extents<std::size_t, md::dynamic_extent, 2>>(
-            vertices.data(), vertices.size() / shape1, shape1),
-        fn, constants,
-        md::mdspan(coeffs.data(), vertices.size() / shape1, cstride), cdofs_b);
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
+        = (itg_type == fem::IntegralType::exterior_facet)
+              ? facet_perms
+              : md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>>{};
+
+    for (int i = 0; i < M.num_integrals(itg_type, 0); ++i)
+    {
+      auto fn = M.kernel(itg_type, i, 0);
+      assert(fn);
+      auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+
+      std::span entities = M.domain(itg_type, i, 0);
+
+      // Two values per each adj. cell (cell index and local entity index).
+      assert((entities.size() / 2) * cstride == coeffs.size());
+      value += impl::assemble_entities(
+          x_dofmap, x,
+          md::mdspan<const std::int32_t,
+                     md::extents<std::size_t, md::dynamic_extent, 2>>(
+              entities.data(), entities.size() / 2, 2),
+          fn, constants,
+          md::mdspan(coeffs.data(), entities.size() / 2, cstride), perms,
+          cdofs_b);
+    }
   }
 
   return value;
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index 3e8c1c2a13d..f57e646dd9b 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -229,9 +229,9 @@ void _lift_bc_cells(
 /// @param[in,out] b Vector to modify.
 /// @param[in] x_dofmap Degree-of-freedom map for the mesh geometry.
 /// @param[in] x Mesh geometry (coordinates).
-/// @param[in] kernel Kernel function to execute over each facet.
-/// @param[in] facets Facets to execute the kernel over, where for the
-/// `i`th facet `facets(i, 0)` is the attached cell and `facets(i, 1)`
+/// @param[in] kernel Kernel function to execute over each entity.
+/// @param[in] entities Entities to execute the kernel over, where for the
+/// `i`th entity `enities(i, 0)` is the attached cell and `entities(i, 1)`
 /// is the local index of the facet relative to the cell.
 /// @param[in] dofmap0 Test function (row) degree-of-freedom data
 /// holding the (0) dofmap, (1) dofmap block size and (2) dofmap
@@ -263,7 +263,7 @@ void _lift_bc_cells(
 template <typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
-void _lift_bc_exterior_facets(
+void _lift_bc_entities(
     V&& b, mdspan2_t x_dofmap,
     md::mdspan<const scalar_value_t<T>,
                md::extents<std::size_t, md::dynamic_extent, 3>>
@@ -271,7 +271,7 @@ void _lift_bc_exterior_facets(
     FEkernel<T> auto kernel,
     md::mdspan<const std::int32_t,
                md::extents<std::size_t, md::dynamic_extent, 2>>
-        facets,
+        entities,
     std::tuple<mdspan2_t, int,
                md::mdspan<const std::int32_t,
                           md::extents<std::size_t, md::dynamic_extent, 2>>>
@@ -288,11 +288,11 @@ void _lift_bc_exterior_facets(
     std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
 {
-  if (facets.empty())
+  if (entities.empty())
     return;
 
-  const auto [dmap0, bs0, facets0] = dofmap0;
-  const auto [dmap1, bs1, facets1] = dofmap1;
+  const auto [dmap0, bs0, entities0] = dofmap0;
+  const auto [dmap1, bs1, entities1] = dofmap1;
 
   const int num_rows = bs0 * dmap0.extent(1);
   const int num_cols = bs1 * dmap1.extent(1);
@@ -300,18 +300,18 @@ void _lift_bc_exterior_facets(
   // Data structures used in bc application
   std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
   std::vector<T> Ae(num_rows * num_cols), be(num_rows);
-  assert(facets0.size() == facets.size());
-  assert(facets1.size() == facets.size());
-  for (std::size_t index = 0; index < facets.extent(0); ++index)
+  assert(entities0.size() == entities.size());
+  assert(entities1.size() == entities.size());
+  for (std::size_t index = 0; index < entities.extent(0); ++index)
   {
     // Cell in integration domain, test function and trial function
     // meshes
-    std::int32_t cell = facets(index, 0);
-    std::int32_t cell0 = facets0(index, 0);
-    std::int32_t cell1 = facets1(index, 0);
+    std::int32_t cell = entities(index, 0);
+    std::int32_t cell0 = entities0(index, 0);
+    std::int32_t cell1 = entities1(index, 0);
 
-    // Local facet index
-    std::int32_t local_facet = facets(index, 1);
+    // Local entity index
+    std::int32_t local_entity = entities(index, 1);
 
     // Get dof maps for cell
     auto dofs1 = md::submdspan(dmap1, cell1, md::full_extent);
@@ -342,10 +342,10 @@ void _lift_bc_exterior_facets(
     auto dofs0 = md::submdspan(dmap0, cell0, md::full_extent);
 
     // Permutations
-    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
+    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
-           &local_facet, &perm, nullptr);
+           &local_entity, &perm, nullptr);
     P0(Ae, cell_info0, cell0, num_cols);
     P1T(Ae, cell_info1, cell1, num_rows);
 
@@ -718,7 +718,16 @@ void assemble_cells(
   }
 }
 
-/// @brief Execute kernel over cells and accumulate result in vector.
+/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result
+/// in a matrix.
+///
+/// Each entity is represented by (i) a cell that the entity is attached to
+/// and (ii) the local index of the entity  with respect to the cell. The
+/// kernel is executed for each entity. The kernel can access data
+/// (e.g., coefficients, basis functions) associated with the attached cell.
+/// However, entities may be attached to more than one cell. This function
+/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen
+/// from cell used to define the entity.
 ///
 /// @tparam T Scalar type.
 /// @tparam _bs The block size of the form test function dof map. If
@@ -730,7 +739,7 @@ void assemble_cells(
 /// @param[in,out] b The vector to accumulate into.
 /// @param[in] x_dofmap Dofmap for the mesh geometry.
 /// @param[in] x Mesh geometry (coordinates).
-/// @param[in] facets Facets (in the integration domain mesh) to execute
+/// @param[in] entities Entities (in the integration domain mesh) to execute
 /// the kernel over.
 /// @param[in] dofmap Test function (row) degree-of-freedom data holding
 /// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices.
@@ -740,19 +749,19 @@ void assemble_cells(
 /// `(cells.size(), coeffs_per_cell)`.
 /// @param[in] cell_info0 The cell permutation information for the test
 /// function mesh.
-/// @param[in] perms Facet permutation integer. Empty if facet
+/// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
 template <int _bs = -1, typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
-void assemble_exterior_facets(
+void assemble_entities(
     fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
     md::mdspan<const scalar_value_t<T>,
                md::extents<std::size_t, md::dynamic_extent, 3>>
         x,
     md::mdspan<const std::int32_t,
                std::extents<std::size_t, md::dynamic_extent, 2>>
-        facets,
+        entities,
     std::tuple<mdspan2_t, int,
                md::mdspan<const std::int32_t,
                           std::extents<std::size_t, md::dynamic_extent, 2>>>
@@ -762,24 +771,24 @@ void assemble_exterior_facets(
     std::span<const std::uint32_t> cell_info0,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
 {
-  if (facets.empty())
+  if (entities.empty())
     return;
 
-  const auto [dmap, bs, facets0] = dofmap;
+  const auto [dmap, bs, entities0] = dofmap;
   assert(_bs < 0 or _bs == bs);
 
   // Create data structures used in assembly
   const int num_dofs = dmap.extent(1);
   std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
   std::vector<T> be(bs * num_dofs);
-  assert(facets0.size() == facets.size());
-  for (std::size_t f = 0; f < facets.extent(0); ++f)
+  assert(entities0.size() == entities.size());
+  for (std::size_t f = 0; f < entities.extent(0); ++f)
   {
     // Cell in the integration domain, local facet index relative to the
     // integration domain cell, and cell in the test function mesh
-    std::int32_t cell = facets(f, 0);
-    std::int32_t local_facet = facets(f, 1);
-    std::int32_t cell0 = facets0(f, 0);
+    std::int32_t cell = entities(f, 0);
+    std::int32_t local_entity = entities(f, 1);
+    std::int32_t cell0 = entities0(f, 0);
 
     // Get cell coordinates/geometry
     auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
@@ -787,12 +796,12 @@ void assemble_exterior_facets(
       std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
 
     // Permutations
-    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
+    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
 
     // Tabulate element vector
     std::ranges::fill(be, 0);
     kernel(be.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
-           &local_facet, &perm, nullptr);
+           &local_entity, &perm, nullptr);
     P0(be, cell_info0, cell0, 1);
 
     // Add element vector to global vector
@@ -941,94 +950,6 @@ void assemble_interior_facets(
   }
 }
 
-/// @brief Execute kernel over a set of vertices and accumulate result in
-/// vector.
-///
-/// @tparam T  Scalar type
-/// @tparam _bs Block size of the form test function dof map. If less
-/// than zero the block size is determined at runtime. If `_bs` is
-/// positive the block size is used as a compile-time constant, which
-/// has performance benefits.
-/// @param[in] P0 Function that applies transformation `P0.b` in-place
-/// to `b` to transform test degrees-of-freedom.
-/// @param[in,out] b Array to accumulate into.
-/// @param[in] x_dofmap Dofmap for the mesh geometry.
-/// @param[in] x Mesh geometry (coordinates).
-/// @param[in] vertices Vertex indices `(vertices.size(), 2)` - first entry
-/// holds the index of the cell adjacent to the vertex, and the second
-/// stores the local index of the vertex within the cell.
-/// @param[in] dofmap Test function (row) degree-of-freedom data holding
-/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices.
-/// @param[in] kernel Kernel function to execute over each cell.
-/// @param[in] constants Constant coefficient data in the kernel.
-/// @param[in] coeffs Coefficient data in the kernel. It has shape
-/// `(vertices.size(), num_cell_coeffs)`. `coeffs(i, j)` is the `j`th
-/// coefficient for cell `i`.
-/// @param[in] cell_info0 Cell permutation information for the test
-/// function mesh.
-template <dolfinx::scalar T, int _bs = -1>
-void assemble_vertices(
-    fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
-    md::mdspan<const std::int32_t,
-               md::extents<std::size_t, md::dynamic_extent, 2>>
-        vertices,
-    std::tuple<mdspan2_t, int,
-               md::mdspan<const std::int32_t,
-                          md::extents<std::size_t, md::dynamic_extent, 2>>>
-        dofmap,
-    FEkernel<T> auto kernel, std::span<const T> constants,
-    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
-    std::span<const std::uint32_t> cell_info0)
-{
-  if (vertices.empty())
-    return;
-
-  const auto [dmap, bs, vertices0] = dofmap;
-  assert(_bs < 0 or _bs == bs);
-
-  // Create data structures used in assembly
-  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
-  std::vector<T> be(bs * dmap.extent(1));
-
-  // Iterate over active vertices
-  for (std::size_t index = 0; index < vertices.extent(0); ++index)
-  {
-    // Integration domain cell, local index, and test function cell
-    std::int32_t cell = vertices(index, 0);
-    std::int32_t local_index = vertices(index, 1);
-    std::int32_t c0 = vertices0(index, 0);
-
-    // Get cell coordinates/geometry
-    auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
-    for (std::size_t i = 0; i < x_dofs.size(); ++i)
-      std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
-
-    // Tabulate vector for vertex
-    std::ranges::fill(be, 0);
-    kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
-           &local_index, nullptr, nullptr);
-    P0(be, cell_info0, c0, 1);
-
-    // Scatter vertex vector to 'global' vector array
-    auto dofs = md::submdspan(dmap, c0, md::full_extent);
-    if constexpr (_bs > 0)
-    {
-      for (std::size_t i = 0; i < dofs.size(); ++i)
-        for (int k = 0; k < _bs; ++k)
-          b[_bs * dofs[i] + k] += be[_bs * i + k];
-    }
-    else
-    {
-      for (std::size_t i = 0; i < dofs.size(); ++i)
-        for (int k = 0; k < bs; ++k)
-          b[bs * dofs[i] + k] += be[bs * i + k];
-    }
-  }
-}
-
 /// Modify RHS vector to account for boundary condition such that:
 ///
 /// b <- b - alpha * A.(x_bc - x0)
@@ -1132,7 +1053,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
     }
   }
 
-  md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
+  md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
   if (a.needs_facet_permutations())
   {
     mesh::CellType cell_type = mesh->topology()->cell_type();
@@ -1141,32 +1062,8 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
     mesh->topology_mutable()->create_entity_permutations();
     const std::vector<std::uint8_t>& p
         = mesh->topology()->get_facet_permutations();
-    perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
-                       num_facets_per_cell);
-  }
-
-  for (int i = 0; i < a.num_integrals(IntegralType::exterior_facet, 0); ++i)
-  {
-    auto kernel = a.kernel(IntegralType::exterior_facet, i, 0);
-    assert(kernel);
-    auto& [coeffs, cstride]
-        = coefficients.at({IntegralType::exterior_facet, i});
-
-    using mdspanx2_t
-        = md::mdspan<const std::int32_t,
-                     md::extents<std::size_t, md::dynamic_extent, 2>>;
-    std::span f = a.domain(IntegralType::exterior_facet, i, 0);
-    mdspanx2_t facets(f.data(), f.size() / 2, 2);
-    std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0);
-    mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
-    std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0);
-    mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
-    assert(coeffs.size() == facets.extent(0) * cstride);
-    _lift_bc_exterior_facets(
-        b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
-        {dofmap1, bs1, facets1}, P1T, constants,
-        md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
-        cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
+    facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
+                             num_facets_per_cell);
   }
 
   for (int i = 0; i < a.num_integrals(IntegralType::interior_facet, 0); ++i)
@@ -1192,7 +1089,39 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
         b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
         {dofmap1, bs1, facets1}, P1T, constants,
         mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0,
-        cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
+        cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms);
+  }
+
+  for (auto itg_type : {fem::IntegralType::exterior_facet,
+                        fem::IntegralType::vertex, fem::IntegralType::ridge})
+  {
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
+        = (itg_type == fem::IntegralType::exterior_facet)
+              ? facet_perms
+              : md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>>{};
+
+    for (int i = 0; i < a.num_integrals(itg_type, 0); ++i)
+    {
+      auto kernel = a.kernel(itg_type, i, 0);
+      assert(kernel);
+      auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+
+      using mdspanx2_t
+          = md::mdspan<const std::int32_t,
+                       md::extents<std::size_t, md::dynamic_extent, 2>>;
+      std::span e = a.domain(itg_type, i, 0);
+      mdspanx2_t entities(e.data(), e.size() / 2, 2);
+      std::span e0 = a.domain_arg(itg_type, 0, i, 0);
+      mdspanx2_t entities0(e0.data(), e0.size() / 2, 2);
+      std::span e1 = a.domain_arg(itg_type, 1, i, 0);
+      mdspanx2_t entities1(e1.data(), e1.size() / 2, 2);
+      assert(coeffs.size() == entities.extent(0) * cstride);
+      _lift_bc_entities(
+          b, x_dofmap, x, kernel, entities, {dofmap0, bs0, entities0}, P0,
+          {dofmap1, bs1, entities1}, P1T, constants,
+          md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0,
+          cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
+    }
   }
 }
 
@@ -1368,7 +1297,7 @@ void assemble_vector(
       }
     }
 
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
     if (L.needs_facet_permutations())
     {
       mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
@@ -1377,48 +1306,14 @@ void assemble_vector(
       mesh->topology_mutable()->create_entity_permutations();
       const std::vector<std::uint8_t>& p
           = mesh->topology()->get_facet_permutations();
-      perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
-                         num_facets_per_cell);
+      facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
+                               num_facets_per_cell);
     }
 
     using mdspanx2_t
         = md::mdspan<const std::int32_t,
                      md::extents<std::size_t, md::dynamic_extent, 2>>;
 
-    for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i)
-    {
-      auto fn = L.kernel(IntegralType::exterior_facet, i, 0);
-      assert(fn);
-      auto& [coeffs, cstride]
-          = coefficients.at({IntegralType::exterior_facet, i});
-      std::span f = L.domain(IntegralType::exterior_facet, i, 0);
-      mdspanx2_t facets(f.data(), f.size() / 2, 2);
-      std::span f1 = L.domain_arg(IntegralType::exterior_facet, 0, i, 0);
-      mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
-      assert((facets.size() / 2) * cstride == coeffs.size());
-      if (bs == 1)
-      {
-        impl::assemble_exterior_facets<1>(
-            P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
-            md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
-            perms);
-      }
-      else if (bs == 3)
-      {
-        impl::assemble_exterior_facets<3>(
-            P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
-            md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
-            perms);
-      }
-      else
-      {
-        impl::assemble_exterior_facets(
-            P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
-            md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
-            perms);
-      }
-    }
-
     for (int i = 0; i < L.num_integrals(IntegralType::interior_facet, 0); ++i)
     {
       using mdspanx22_t
@@ -1444,7 +1339,7 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, perms);
+            cell_info0, facet_perms);
       }
       else if (bs == 3)
       {
@@ -1455,7 +1350,7 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, perms);
+            cell_info0, facet_perms);
       }
       else
       {
@@ -1466,34 +1361,52 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, perms);
+            cell_info0, facet_perms);
       }
     }
 
-    for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i)
+    for (auto itg_type : {fem::IntegralType::exterior_facet,
+                          fem::IntegralType::vertex, fem::IntegralType::ridge})
     {
-      auto fn = L.kernel(IntegralType::vertex, i, 0);
-      assert(fn);
-
-      std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx);
-      std::span vertices0
-          = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx);
-
-      auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
-
-      assert(vertices.size() * cstride == coeffs.size());
-
-      impl::assemble_vertices<T>(
-          P0, b, x_dofmap, x,
-          md::mdspan<const std::int32_t,
-                     md::extents<std::size_t, md::dynamic_extent, 2>>(
-              vertices.data(), vertices.size() / 2, 2),
-          {dofs, bs,
-           md::mdspan<const std::int32_t,
-                      md::extents<std::size_t, md::dynamic_extent, 2>>(
-               vertices0.data(), vertices0.size() / 2, 2)},
-          fn, constants,
-          md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0);
+      md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
+          = (itg_type == fem::IntegralType::exterior_facet)
+                ? facet_perms
+                : md::mdspan<const std::uint8_t,
+                             md::dextents<std::size_t, 2>>{};
+      for (int i = 0; i < L.num_integrals(itg_type, 0); ++i)
+      {
+        auto fn = L.kernel(itg_type, i, 0);
+        assert(fn);
+        auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+        std::span e = L.domain(itg_type, i, 0);
+        mdspanx2_t entities(e.data(), e.size() / 2, 2);
+        std::span e1 = L.domain_arg(itg_type, 0, i, 0);
+        mdspanx2_t entities1(e1.data(), e1.size() / 2, 2);
+        assert((entities.size() / 2) * cstride == coeffs.size());
+        if (bs == 1)
+        {
+          impl::assemble_entities<1>(
+              P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
+              constants, md::mdspan(coeffs.data(), entities.extent(0), cstride),
+              cell_info0, perms);
+        }
+        else if (bs == 3)
+        {
+          impl::assemble_entities<3>(
+              P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
+              constants,
+              md::mdspan(coeffs.data(), entities.size() / 2, cstride),
+              cell_info0, perms);
+        }
+        else
+        {
+          impl::assemble_entities(
+              P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
+              constants,
+              md::mdspan(coeffs.data(), entities.size() / 2, cstride),
+              cell_info0, perms);
+        }
+      }
     }
   }
 }
diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h
index 9b08e845d81..65777820ee8 100644
--- a/cpp/dolfinx/fem/pack.h
+++ b/cpp/dolfinx/fem/pack.h
@@ -179,11 +179,8 @@ allocate_coefficient_storage(const Form<T, U>& form, IntegralType integral_type,
     const std::vector<int> offsets = form.coefficient_offsets();
     cstride = offsets.back();
     num_entities = form.domain(integral_type, id, 0).size();
-    if (integral_type == IntegralType::exterior_facet
-        or integral_type == IntegralType::interior_facet)
-    {
+    if (integral_type != IntegralType::cell)
       num_entities /= 2;
-    }
   }
 
   return {std::vector<T>(num_entities * cstride), cstride};
@@ -270,28 +267,6 @@ void pack_coefficients(const Form<T, U>& form,
         }
         break;
       }
-      case IntegralType::exterior_facet:
-      {
-        // Iterate over coefficients coefficients that are active in
-        // exterior facet integrals
-        for (int coeff : form.active_coeffs(IntegralType::exterior_facet, id))
-        {
-          auto mesh = coefficients[coeff]->function_space()->mesh();
-          std::span<const std::int32_t> facets_b
-              = form.domain_coeff(IntegralType::exterior_facet, id, coeff);
-          md::mdspan<const std::int32_t,
-                     md::extents<std::size_t, md::dynamic_extent, 2>>
-              facets(facets_b.data(), facets_b.size() / 2, 2);
-          auto cells = md::submdspan(facets, md::full_extent, 0);
-
-          std::span<const std::uint32_t> cell_info
-              = impl::get_cell_orientation_info(*coefficients[coeff]);
-          impl::pack_coefficient_entity(std::span(c), cstride,
-                                        *coefficients[coeff], cell_info, cells,
-                                        offsets[coeff]);
-        }
-        break;
-      }
       case IntegralType::interior_facet:
       {
         // Iterate over coefficients that are active in interior
@@ -322,25 +297,27 @@ void pack_coefficients(const Form<T, U>& form,
         }
         break;
       }
+      case IntegralType::exterior_facet:
       case IntegralType::vertex:
+      case IntegralType::ridge:
       {
         // Iterate over coefficients that are active in vertex integrals
-        for (int coeff : form.active_coeffs(IntegralType::vertex, id))
+        for (int coeff : form.active_coeffs(integral_type, id))
         {
           // Get coefficient mesh
           auto mesh = coefficients[coeff]->function_space()->mesh();
           assert(mesh);
 
-          std::span<const std::int32_t> vertices_b
-              = form.domain_coeff(IntegralType::vertex, id, coeff);
+          std::span<const std::int32_t> entitites_b
+              = form.domain_coeff(integral_type, id, coeff);
           md::mdspan<const std::int32_t,
                      md::extents<std::size_t, md::dynamic_extent, 2>>
-              vertices(vertices_b.data(), vertices_b.size() / 2, 2);
+              entities(entitites_b.data(), entitites_b.size() / 2, 2);
           std::span<const std::uint32_t> cell_info
               = impl::get_cell_orientation_info(*coefficients[coeff]);
           impl::pack_coefficient_entity(
               std::span(c), cstride, *coefficients[coeff], cell_info,
-              md::submdspan(vertices, md::full_extent, 0), offsets[coeff]);
+              md::submdspan(entities, md::full_extent, 0), offsets[coeff]);
         }
         break;
       }
diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp
index ec7a8501ae3..0aadd5fec66 100644
--- a/cpp/dolfinx/fem/utils.cpp
+++ b/cpp/dolfinx/fem/utils.cpp
@@ -154,6 +154,9 @@ fem::compute_integration_domains(fem::IntegralType integral_type,
   case IntegralType::vertex:
     dim = 0;
     break;
+  case IntegralType::ridge:
+    dim = tdim - 2;
+    break;
   default:
     throw std::runtime_error(
         "Cannot compute integration domains. Integral type not supported.");
@@ -199,22 +202,6 @@ fem::compute_integration_domains(fem::IntegralType integral_type,
     entity_data.insert(entity_data.begin(), entities.begin(), entities.end());
     break;
   }
-  case IntegralType::exterior_facet:
-  {
-    auto [f_to_c, c_to_f] = get_connectivities(tdim - 1);
-    // Create list of tagged boundary facets
-    const std::vector bfacets = mesh::exterior_facet_indices(topology);
-    std::vector<std::int32_t> facets;
-    std::ranges::set_intersection(entities, bfacets,
-                                  std::back_inserter(facets));
-    for (auto f : facets)
-    {
-      // Get the facet as a pair of (cell, local facet)
-      auto facet = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
-      entity_data.insert(entity_data.end(), facet.begin(), facet.end());
-    }
-    break;
-  }
   case IntegralType::interior_facet:
   {
     auto [f_to_c, c_to_f] = get_connectivities(tdim - 1);
@@ -248,17 +235,22 @@ fem::compute_integration_domains(fem::IntegralType integral_type,
     }
     break;
   }
+  case IntegralType::exterior_facet:
   case IntegralType::vertex:
+  case IntegralType::ridge:
   {
-    auto [v_to_c, c_to_v] = get_connectivities(0);
-    for (auto vertex : entities)
+    auto [e_to_c, c_to_e] = get_connectivities(dim);
+    for (auto entity : entities)
     {
-      std::array<std::int32_t, 2> pair = impl::get_cell_vertex_pairs<1>(
-          vertex, v_to_c->links(vertex), *c_to_v);
-
+      std::array<std::int32_t, 2> pair = impl::get_cell_entity_pairs<1>(
+          entity, e_to_c->links(entity), *c_to_e);
       entity_data.insert(entity_data.end(), pair.begin(), pair.end());
     }
+    break;
   }
+  default:
+    throw std::runtime_error(
+        "Cannot compute integration domains. Integral type not supported.");
   }
 
   return entity_data;
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 22436c8e50d..094cf3e55b2 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -92,17 +92,17 @@ get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells,
   return cell_local_facet_pairs;
 }
 
-/// Helper function to get an array of of (cell, local_vertex) pairs
-/// corresponding to a given vertex index.
-/// @note If the vertex is connected to multiple cells, the first one is picked.
-/// @param[in] v vertex index
-/// @param[in] cells List of cells incident to the vertex
-/// @param[in] c_to_v Cell to vertex connectivity
-/// @return Vector of (cell, local_vertex) pairs
+/// Helper function to get an array of of (cell, local_entity) pairs
+/// corresponding to a given entity index.
+/// @note If the entity is connected to multiple cells, the first one is picked.
+/// @param[in] e entity index
+/// @param[in] cells List of cells incident to the entity
+/// @param[in] c_to_e Cell to entity connectivity
+/// @return Vector of (cell, local_entity) pairs
 template <int num_cells>
 std::array<std::int32_t, 2 * num_cells>
-get_cell_vertex_pairs(std::int32_t v, std::span<const std::int32_t> cells,
-                      const graph::AdjacencyList<std::int32_t>& c_to_v)
+get_cell_entity_pairs(std::int32_t e, std::span<const std::int32_t> cells,
+                      const graph::AdjacencyList<std::int32_t>& c_to_e)
 {
   static_assert(num_cells == 1); // Patch assembly not supported.
 
@@ -111,11 +111,11 @@ get_cell_vertex_pairs(std::int32_t v, std::span<const std::int32_t> cells,
   // Use first cell for assembly over by default
   std::int32_t cell = cells[0];
 
-  // Find local index of vertex within cell
-  auto cell_vertices = c_to_v.links(cell);
-  auto it = std::ranges::find(cell_vertices, v);
-  assert(it != cell_vertices.end());
-  std::int32_t local_index = std::distance(cell_vertices.begin(), it);
+  // Find local index of entity within cell
+  auto cell_entities = c_to_e.links(cell);
+  auto it = std::ranges::find(cell_entities, e);
+  assert(it != cell_entities.end());
+  std::int32_t local_index = std::distance(cell_entities.begin(), it);
 
   return {cell, local_index};
 }
@@ -146,7 +146,7 @@ get_cell_vertex_pairs(std::int32_t v, std::span<const std::int32_t> cells,
 /// @param[in] entities List of mesh entities. Depending on the `IntegralType`
 /// these are associated with different entities:
 ///     `IntegralType::cell`:             cells
-///     `IntegralType::exterior_facet`:   facets
+///     `IntegralType::exterior_facet`: facets
 ///     `IntegralType::interior_facet`:   facets
 ///     `IntegralType::vertex`:           vertices
 /// @return List of integration entity data, depending on the `IntegralType` the
@@ -291,6 +291,8 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form<T, U>& a)
         }
         break;
       case IntegralType::exterior_facet:
+      case IntegralType::ridge:
+      case IntegralType::vertex:
         for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
         {
           sparsitybuild::cells(pattern,
@@ -490,7 +492,7 @@ Form<T, U> create_form_factory(
   // integral offsets. Since the UFL forms for each type of cell should be
   // the same, I think this assumption is OK.
   const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets;
-  std::array<int, 4> num_integrals_type;
+  std::array<int, 5> num_integrals_type;
   for (std::size_t i = 0; i < num_integrals_type.size(); ++i)
     num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i];
 
@@ -502,6 +504,7 @@ Form<T, U> create_form_factory(
   }
 
   // Create facets, if required
+  // NOTE: exterior_facet and interior_facet is declared in ufcx.h
   if (num_integrals_type[exterior_facet] > 0
       or num_integrals_type[interior_facet] > 0)
   {
@@ -510,6 +513,14 @@ Form<T, U> create_form_factory(
     mesh->topology_mutable()->create_connectivity(tdim, tdim - 1);
   }
 
+  // Create ridges, if required
+  if (num_integrals_type[ridge] > 0)
+  {
+    mesh->topology_mutable()->create_entities(tdim - 2);
+    mesh->topology_mutable()->create_connectivity(tdim - 2, tdim);
+    mesh->topology_mutable()->create_connectivity(tdim, tdim - 2);
+  }
+
   // Get list of integral IDs, and load tabulate tensor into memory for
   // each
   std::map<std::tuple<IntegralType, int, int>, integral_data<T, U>> integrals;
@@ -593,78 +604,9 @@ Form<T, U> create_form_factory(
     }
   }
 
-  // Attach exterior facet kernels
-  std::vector<std::int32_t> default_facets_ext;
-  {
-    std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
-                                 + integral_offsets[exterior_facet],
-                             num_integrals_type[exterior_facet]);
-    auto sd = subdomains.find(IntegralType::exterior_facet);
-    for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
-    {
-      const ufcx_form& ufcx_form = ufcx_forms[form_idx];
-      for (int i = 0; i < num_integrals_type[exterior_facet]; ++i)
-      {
-        const int id = ids[i];
-        ufcx_integral* integral
-            = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i];
-        assert(integral);
-        check_geometry_hash(*integral, form_idx);
-
-        std::vector<int> active_coeffs;
-        for (int j = 0; j < ufcx_form.num_coefficients; ++j)
-        {
-          if (integral->enabled_coefficients[j])
-            active_coeffs.push_back(j);
-        }
-
-        impl::kernel_t<T, U> k = impl::extract_kernel<T>(integral);
-
-        // Build list of entities to assembler over
-        const std::vector bfacets = mesh::exterior_facet_indices(*topology);
-        auto f_to_c = topology->connectivity(tdim - 1, tdim);
-        assert(f_to_c);
-        auto c_to_f = topology->connectivity(tdim, tdim - 1);
-        assert(c_to_f);
-        if (id == -1)
-        {
-          // Default kernel, operates on all (owned) exterior facets
-          default_facets_ext.reserve(2 * bfacets.size());
-          for (std::int32_t f : bfacets)
-          {
-            // There will only be one pair for an exterior facet integral
-            std::array<std::int32_t, 2> pair
-                = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
-            default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
-                                      pair.end());
-          }
-          integrals.insert({{IntegralType::exterior_facet, i, form_idx},
-                            {k, default_facets_ext, active_coeffs}});
-        }
-        else if (sd != subdomains.end())
-        {
-          // NOTE: This requires that pairs are sorted
-          auto it = std::ranges::lower_bound(sd->second, id, std::less<>{},
-                                             [](auto& a) { return a.first; });
-          if (it != sd->second.end() and it->first == id)
-          {
-            integrals.insert({{IntegralType::exterior_facet, i, form_idx},
-                              {k,
-                               std::vector<std::int32_t>(it->second.begin(),
-                                                         it->second.end()),
-                               active_coeffs}});
-          }
-        }
-
-        if (integral->needs_facet_permutations)
-          needs_facet_permutations = true;
-      }
-    }
-  }
-
   // Attach interior facet kernels
-  std::vector<std::int32_t> default_facets_int;
   {
+    std::vector<std::int32_t> default_facets_int;
     std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
                                  + integral_offsets[interior_facet],
                              num_integrals_type[interior_facet]);
@@ -757,82 +699,119 @@ Form<T, U> create_form_factory(
     }
   }
 
-  // Attach vertex kernels
+  // Attach exterior entity integrals
   {
-    for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
+    for (IntegralType itg_type : {IntegralType::exterior_facet,
+                                  IntegralType::vertex, IntegralType::ridge})
     {
-      const ufcx_form& form = ufcx_forms[form_idx];
-
-      std::span<const int> ids(form.form_integral_ids
-                                   + integral_offsets[vertex],
-                               num_integrals_type[vertex]);
-      auto sd = subdomains.find(IntegralType::vertex);
-      for (int i = 0; i < num_integrals_type[vertex]; ++i)
+      std::size_t dim;
+      switch (itg_type)
       {
-        const int id = ids[i];
-        ufcx_integral* integral
-            = form.form_integrals[integral_offsets[vertex] + i];
-        assert(integral);
-        check_geometry_hash(*integral, form_idx);
+      case IntegralType::exterior_facet:
+      {
+        dim = tdim - 1;
+        break;
+      }
+      case IntegralType::ridge:
+      {
+        dim = tdim - 2;
+        break;
+      }
+      case IntegralType::vertex:
+      {
+        dim = 0;
+        break;
+      }
+      default:
+        throw std::runtime_error("Unsupported integral type");
+      }
 
-        std::vector<int> active_coeffs;
-        for (int j = 0; j < form.num_coefficients; ++j)
+      const std::function<std::vector<std::int32_t>(const mesh::Topology&,
+                                                    IntegralType)>
+          get_default_integration_entities
+          = [dim](const mesh::Topology& topology, IntegralType itg_type)
+      {
+        if (itg_type == IntegralType::exterior_facet)
         {
-          if (integral->enabled_coefficients[j])
-            active_coeffs.push_back(j);
+          // Integrate over all owned exterior facets
+          return mesh::exterior_facet_indices(topology);
+        }
+        else
+        {
+          // Integrate over all owned entities
+          std::int32_t num_entities = topology.index_map(dim)->size_local();
+          std::vector<std::int32_t> entities(num_entities);
+          std::iota(entities.begin(), entities.end(), 0);
+          return entities;
         }
+      };
 
-        impl::kernel_t<T, U> k = impl::extract_kernel<T>(integral);
-        assert(k);
+      std::vector<std::int32_t> default_entities_ext;
 
-        // Build list of entities to assembler over
-        auto v_to_c = topology->connectivity(0, tdim);
-        assert(v_to_c);
-        auto c_to_v = topology->connectivity(tdim, 0);
-        assert(c_to_v);
-
-        // pack for a range of vertices a flattened list of cell index c_i and
-        // local vertex index l_i:
-        //  [c_0, l_0, ..., c_n, l_n]
-        auto get_cells_and_vertices = [v_to_c, c_to_v](auto vertices_range)
+      std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
+                                   + integral_offsets[(std::int8_t)itg_type],
+                               num_integrals_type[(std::int8_t)itg_type]);
+      auto sd = subdomains.find(itg_type);
+      for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
+      {
+        const ufcx_form& ufcx_form = ufcx_forms[form_idx];
+        for (int i = 0; i < num_integrals_type[(std::int8_t)itg_type]; ++i)
         {
-          std::vector<std::int32_t> cell_and_vertex;
-          cell_and_vertex.reserve(2 * vertices_range.size());
-          for (std::int32_t vertex : vertices_range)
+          const int id = ids[i];
+          ufcx_integral* integral
+              = ufcx_form.form_integrals[integral_offsets[(std::int8_t)itg_type]
+                                         + i];
+          assert(integral);
+          check_geometry_hash(*integral, form_idx);
+
+          std::vector<int> active_coeffs;
+          for (int j = 0; j < ufcx_form.num_coefficients; ++j)
           {
-            std::array<std::int32_t, 2> pair = impl::get_cell_vertex_pairs<1>(
-                vertex, v_to_c->links(vertex), *c_to_v);
-
-            cell_and_vertex.insert(cell_and_vertex.end(), pair.begin(),
-                                   pair.end());
+            if (integral->enabled_coefficients[j])
+              active_coeffs.push_back(j);
           }
-          assert(cell_and_vertex.size() == 2 * vertices_range.size());
-          return cell_and_vertex;
-        };
 
-        if (id == -1)
-        {
-          // Default vertex kernel operates on all (owned) vertices
-          std::int32_t num_vertices = topology->index_map(0)->size_local();
-          std::vector<std::int32_t> cells_and_vertices = get_cells_and_vertices(
-              std::ranges::views::iota(0, num_vertices));
+          impl::kernel_t<T, U> k = impl::extract_kernel<T>(integral);
 
-          integrals.insert({{IntegralType::vertex, i, form_idx},
-                            {k, cells_and_vertices, active_coeffs}});
-        }
-        else if (sd != subdomains.end())
-        {
-          // NOTE: This requires that pairs are sorted
-          auto it = std::ranges::lower_bound(sd->second, id, std::less<>{},
-                                             [](auto& a) { return a.first; });
-          if (it != sd->second.end() and it->first == id)
+          // Build list of entities to assembler over
+          auto e_to_c = topology->connectivity(dim, tdim);
+          assert(e_to_c);
+          auto c_to_e = topology->connectivity(tdim, dim);
+          assert(c_to_e);
+          if (id == -1)
           {
-            integrals.insert({{IntegralType::vertex, i, form_idx},
-                              {k,
-                               std::vector<std::int32_t>(it->second.begin(),
-                                                         it->second.end()),
-                               active_coeffs}});
+            std::vector default_entities
+                = get_default_integration_entities(*topology, itg_type);
+            // Default kernel
+            default_entities_ext.reserve(2 * default_entities.size());
+            for (std::int32_t e : default_entities)
+            {
+              // There will only be one pair for an exterior facet integral
+              std::array<std::int32_t, 2> pair = impl::get_cell_entity_pairs<1>(
+                  e, e_to_c->links(e), *c_to_e);
+              default_entities_ext.insert(default_entities_ext.end(),
+                                          pair.begin(), pair.end());
+            }
+            integrals.insert({{itg_type, i, form_idx},
+                              {k, default_entities_ext, active_coeffs}});
+          }
+          else if (sd != subdomains.end())
+          {
+            // NOTE: This requires that pairs are sorted
+            auto it = std::ranges::lower_bound(sd->second, id, std::less<>{},
+                                               [](auto& a) { return a.first; });
+            if (it != sd->second.end() and it->first == id)
+            {
+              integrals.insert({{itg_type, i, form_idx},
+                                {k,
+                                 std::vector<std::int32_t>(it->second.begin(),
+                                                           it->second.end()),
+                                 active_coeffs}});
+            }
           }
+
+          if (integral->needs_facet_permutations)
+            needs_facet_permutations = true;
         }
       }
     }
diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index c5499a9b3b4..ec39950d340 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -167,6 +167,11 @@ def get_integration_domains(
                 subdomain._cpp_object.topology.create_connectivity(0, tdim)
                 subdomain._cpp_object.topology.create_connectivity(tdim, 0)
 
+            if integral_type is IntegralType.ridge:
+                tdim = subdomain.topology.dim
+                subdomain._cpp_object.topology.create_connectivity(tdim - 2, tdim)
+                subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 2)
+
             # Compute integration domains only for each subdomain id in
             # the integrals. If a process has no integral entities,
             # insert an empty array.
@@ -220,6 +225,7 @@ def form_cpp_class(
     "exterior_facet": IntegralType.exterior_facet,
     "interior_facet": IntegralType.interior_facet,
     "vertex": IntegralType.vertex,
+    "ridge": IntegralType.ridge,
 }
 
 
@@ -373,7 +379,7 @@ def _form(form):
 
         # Extract subdomain ids from ufcx_form
         subdomain_ids = {type: [] for type in sd.get(domain).keys()}
-        integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)]
+        integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(6)]
         for i in range(len(integral_offsets) - 1):
             integral_type = IntegralType(i)
             for j in range(integral_offsets[i], integral_offsets[i + 1]):
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index eda2ce5a6d0..428d91d0026 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -1273,7 +1273,8 @@ void fem(nb::module_& m)
              "exterior facet integral")
       .value("interior_facet", dolfinx::fem::IntegralType::interior_facet,
              "exterior facet integral")
-      .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral");
+      .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral")
+      .value("ridge", dolfinx::fem::IntegralType::ridge, "ridge integral");
 
   declare_function_space<float>(m, "float32");
   declare_function_space<double>(m, "float64");
diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py
index f13cf68d37c..ce0932a406a 100644
--- a/python/test/unit/fem/test_assembler.py
+++ b/python/test/unit/fem/test_assembler.py
@@ -48,7 +48,7 @@
     locate_entities_boundary,
     meshtags,
 )
-from ufl import derivative, dS, ds, dx, inner
+from ufl import derivative, dP, dr, dS, ds, dx, inner
 from ufl.geometry import SpatialCoordinate
 
 dtype_parametrize = pytest.mark.parametrize(
@@ -235,7 +235,32 @@ def test_assembly_bcs(self, mode):
         mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode)
         V = functionspace(mesh, ("Lagrange", 1))
         u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
-        a = form(inner(u, v) * dx + inner(u, v) * ds)
+
+        exterior_ridges = locate_entities_boundary(
+            mesh, mesh.topology.dim - 2, lambda x: np.full(x.shape[1], True, dtype=bool)
+        )
+        ridge_tags = dolfinx.mesh.meshtags(
+            mesh,
+            mesh.topology.dim - 2,
+            exterior_ridges,
+            np.ones_like(exterior_ridges, dtype=np.int32),
+        )
+        dr_exterior = dr(subdomain_data=ridge_tags, subdomain_id=1)
+        exterior_vertices = locate_entities_boundary(
+            mesh, 0, lambda x: np.isclose(x[0], 0.0) | np.isclose(x[1], 1.0)
+        )
+        vertex_tags = dolfinx.mesh.meshtags(
+            mesh, 0, exterior_vertices, np.ones_like(exterior_vertices, dtype=np.int32)
+        )
+        dP_exterior = dP(subdomain_data=vertex_tags, subdomain_id=1)
+        x = ufl.SpatialCoordinate(mesh)
+        f = ufl.exp(x[0] + x[1])
+        a = form(
+            inner(u, v) * dx
+            + inner(u, v) * ds
+            - inner(u, v) * dr_exterior
+            + f * inner(u, v) * dP_exterior
+        )
         L = form(inner(1.0, v) * dx)
 
         bdofsV = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0))
@@ -1877,3 +1902,254 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False):
     assert expected_value_l == pytest.approx(
         value_l.array[: expected_value_l.size], abs=5e4 * np.finfo(rdtype).eps
     )
+
+
+@pytest.mark.parametrize(
+    "cell_type",
+    [
+        mesh.CellType.triangle,
+        mesh.CellType.quadrilateral,
+    ],
+)
+@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet])
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(
+            np.complex64,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+        pytest.param(
+            np.complex128,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+    ],
+)
+def test_ridge_integrals_rank1_2D(cell_type, ghost_mode, dtype):
+    comm = MPI.COMM_WORLD
+    rdtype = np.real(dtype(0)).dtype
+
+    msh = None
+    msh = mesh.create_unit_square(
+        comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype
+    )
+    gdim = msh.geometry.dim
+    V = dolfinx.fem.functionspace(msh, ("Lagrange", 3))
+
+    x = ufl.SpatialCoordinate(msh)
+    u = ufl.TestFunction(V)
+
+    integrand = ufl.conj(u) * x[gdim - 1]
+    dr = ufl.Measure("dr", domain=msh)
+    F = dolfinx.fem.form(integrand * dr, dtype=dtype)
+    b = dolfinx.fem.assemble_vector(F)
+
+    dP = ufl.Measure("dP", domain=msh)
+    Fp = dolfinx.fem.form(integrand * dP, dtype=dtype)
+    b_p = dolfinx.fem.assemble_vector(Fp)
+    tol = np.finfo(rdtype).eps
+    np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol)
+
+
+@pytest.mark.parametrize(
+    "cell_type",
+    [
+        mesh.CellType.triangle,
+        mesh.CellType.quadrilateral,
+        mesh.CellType.tetrahedron,
+        # mesh.CellType.pyramid,
+        # mesh.CellType.prism,
+        mesh.CellType.hexahedron,
+    ],
+)
+@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet])
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(
+            np.complex64,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+        pytest.param(
+            np.complex128,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+    ],
+)
+def test_ridge_integrals_rank0(cell_type, ghost_mode, dtype):
+    comm = MPI.COMM_WORLD
+    rdtype = np.real(dtype(0)).dtype
+
+    msh = None
+    cell_dim = mesh.cell_dim(cell_type)
+    if cell_dim == 2:
+        msh = mesh.create_unit_square(
+            comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype
+        )
+    elif cell_dim == 3:
+        msh = mesh.create_unit_cube(
+            comm, 4, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype
+        )
+    else:
+        raise RuntimeError("Bad dimension")
+    gdim = msh.geometry.dim
+
+    x = ufl.SpatialCoordinate(msh)
+
+    def marked_ridges(x):
+        return np.isclose(x[0], 1) & np.isclose(x[1], 0.5)
+
+    exterior_ridges = dolfinx.mesh.locate_entities_boundary(
+        msh, msh.topology.dim - 2, marked_ridges
+    )
+    et = dolfinx.mesh.meshtags(
+        msh, msh.topology.dim - 2, exterior_ridges, np.full_like(exterior_ridges, 33)
+    )
+    dr = ufl.Measure("dr", domain=msh, subdomain_data=et, subdomain_id=33)
+    integrand = x[0] + x[1] + x[gdim - 1] ** 2
+    J_compiled = dolfinx.fem.form(integrand * dr, dtype=dtype)
+    J = dolfinx.fem.assemble_scalar(J_compiled)
+    J_sum = msh.comm.allreduce(J, op=MPI.SUM)
+    if cell_dim == 2:
+        ref_sol = 1 + 1 / 2 + (1 / 2) ** 2
+    elif cell_dim == 3:
+        ref_sol = 1 + 1 / 2 + 1 / 3
+    else:
+        raise RuntimeError("Bad dimension")
+    tol = 50 * np.finfo(rdtype).eps
+    assert np.isclose(J_sum, ref_sol, atol=tol, rtol=tol)
+
+
+@pytest.mark.parametrize("coefficient", [True, False])
+@pytest.mark.parametrize(
+    "cell_type",
+    [
+        mesh.CellType.tetrahedron,
+        # mesh.CellType.pyramid,
+        # mesh.CellType.prism,
+        mesh.CellType.hexahedron,
+    ],
+)
+@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet])
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(
+            np.complex64,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+        pytest.param(
+            np.complex128,
+            marks=pytest.mark.skipif(
+                os.name == "nt", reason="win32 platform does not support C99 _Complex numbers"
+            ),
+        ),
+    ],
+)
+def test_ridge_integrals_rank1_3D(cell_type, ghost_mode, dtype, coefficient):
+    comm = MPI.COMM_WORLD
+    rdtype = np.real(dtype(0)).dtype
+
+    N = 4
+    msh = mesh.create_unit_cube(
+        comm, N, N, N, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype
+    )
+    gdim = msh.geometry.dim
+
+    el = ("Lagrange", 3, (gdim,))
+
+    volume_element = basix.ufl.element(el[0], msh.basix_cell(), el[1], shape=el[2], dtype=rdtype)
+    V = dolfinx.fem.functionspace(msh, volume_element)
+
+    x = ufl.SpatialCoordinate(msh)
+    v = ufl.TestFunction(V)
+
+    def marked_ridges(x):
+        return np.isclose(x[0], 1) & np.isclose(x[1], 0.5)
+
+    exterior_ridges = dolfinx.mesh.locate_entities_boundary(
+        msh, msh.topology.dim - 2, marked_ridges
+    )
+    et = dolfinx.mesh.meshtags(
+        msh, msh.topology.dim - 2, exterior_ridges, np.full_like(exterior_ridges, 33)
+    )
+
+    def f(mod, x):
+        return x[0] + mod.sin(x[1]), x[2] + x[1], x[0] ** 2 - 3 * x[1]
+
+    def integrand(x, v):
+        if coefficient:
+            Z = dolfinx.fem.functionspace(v.ufl_function_space().mesh, ("Lagrange", 1, (gdim,)))
+            z = dolfinx.fem.Function(Z, dtype=dtype)
+            z.interpolate(lambda x: f(np, x))
+        else:
+            z = ufl.as_vector(f(ufl, x))
+        return ufl.inner(z, v)
+
+    metadata = {"quadrature_degree": 10}
+    dr = ufl.Measure("dr", domain=msh, subdomain_data=et, subdomain_id=33, metadata=metadata)
+    F = dolfinx.fem.form(integrand(x, v) * dr, dtype=dtype)
+    b = dolfinx.fem.assemble_vector(F)
+    b.scatter_reverse(la.InsertMode.add)
+    b.scatter_forward()
+
+    # Create reference solution on unit interval
+    assert gdim == 3
+    if comm.rank == 0:
+        nodes = np.zeros((N + 1, gdim), dtype=rdtype)
+        nodes[:, 0] = 1
+        nodes[:, 1] = 0.5
+        nodes[:, 2] = np.linspace(0, 1, N + 1)
+        connectivity = (
+            np.repeat(np.arange(nodes.shape[0]), 2)[1:-1]
+            .reshape(nodes.shape[0] - 1, 2)
+            .astype(np.int64)
+        )
+    else:
+        nodes = np.zeros((0, gdim), dtype=rdtype)
+        connectivity = np.zeros((0, 2), dtype=np.int64)
+
+    c_el = ufl.Mesh(
+        basix.ufl.element(
+            "Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],), dtype=rdtype
+        ),
+    )
+    line_mesh = dolfinx.mesh.create_mesh(
+        MPI.COMM_WORLD,
+        x=nodes,
+        cells=connectivity,
+        e=c_el,
+        partitioner=dolfinx.mesh.create_cell_partitioner(ghost_mode),
+    )
+
+    line_element = basix.ufl.element(
+        el[0], line_mesh.basix_cell(), el[1], shape=el[2], dtype=rdtype
+    )
+    Q = dolfinx.fem.functionspace(line_mesh, line_element)
+    q = ufl.TestFunction(Q)
+
+    x_e = ufl.SpatialCoordinate(line_mesh)
+    F_ref = dolfinx.fem.form(
+        integrand(x_e, q) * ufl.dx(domain=line_mesh, metadata=metadata), dtype=dtype
+    )
+    b_ref = dolfinx.fem.assemble_vector(F_ref)
+    b_ref.scatter_reverse(la.InsertMode.add)
+    b_ref.scatter_forward()
+    tol = 10 * np.finfo(rdtype).eps
+    assert np.isclose(la.norm(b), la.norm(b_ref), rtol=tol, atol=tol)