diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h
index 5b3a5524a69..a9a3adc1e6c 100644
--- a/cpp/dolfinx/fem/Form.h
+++ b/cpp/dolfinx/fem/Form.h
@@ -119,11 +119,8 @@ class Form
       throw std::runtime_error("No mesh could be associated with the Form.");
 
     // Store kernels, looping over integrals by domain type (dimension)
-    for (auto& integral_type : integrals)
+    for (auto& [type, kernels] : integrals)
     {
-      const IntegralType type = integral_type.first;
-      auto& kernels = integral_type.second;
-
       // Loop over integrals kernels and set domains
       switch (type)
       {
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 0cb5cd96d5a..8891dd96138 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -172,7 +172,6 @@ class Function
     std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
     std::vector<std::int32_t> cells(num_cells, 0);
     std::iota(cells.begin(), cells.end(), 0);
-
     interpolate(v, cells, nmm_interpolation_data);
   }
 
@@ -217,6 +216,7 @@ class Function
         throw std::runtime_error(
             "Data returned by callable has wrong shape(0) size");
       }
+
       if (fshape[1] != x.size() / 3)
       {
         throw std::runtime_error(
@@ -285,6 +285,7 @@ class Function
             "Function element interpolation points has different shape to "
             "Expression interpolation points");
       }
+
       for (std::size_t i = 0; i < X0.size(); ++i)
       {
         if (std::abs(X0[i] - X1[i]) > 1.0e-10)
@@ -295,9 +296,8 @@ class Function
       }
     }
 
-    namespace stdex = std::experimental;
-
     // Array to hold evaluated Expression
+    namespace stdex = std::experimental;
     std::size_t num_cells = cells.size();
     std::size_t num_points = e.X().second[0];
     std::vector<T> fdata(num_cells * num_points * value_size);
@@ -312,7 +312,6 @@ class Function
     // value_size), i.e. xyzxyz ordering of dof values per cell per point.
     // The interpolation uses xxyyzz input, ordered for all points of each
     // cell, i.e. (value_size, num_cells*num_points)
-
     std::vector<T> fdata1(num_cells * num_points * value_size);
     stdex::mdspan<T, stdex::dextents<std::size_t, 3>> f1(
         fdata1.data(), value_size, num_cells, num_points);
@@ -371,6 +370,7 @@ class Function
       throw std::runtime_error(
           "Number of points and number of cells must be equal.");
     }
+
     if (xshape[0] != ushape[0])
     {
       throw std::runtime_error(
diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
index d6e23a8e7fe..00d43452787 100644
--- a/cpp/dolfinx/fem/assembler.h
+++ b/cpp/dolfinx/fem/assembler.h
@@ -27,7 +27,7 @@ class FunctionSpace;
 
 /// @brief Create a map of `std::span`s from a map of `std::vector`s
 template <typename T>
-std::map<std::pair<fem::IntegralType, int>, std::pair<std::span<const T>, int>>
+std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>
 make_coefficients_span(const std::map<std::pair<IntegralType, int>,
                                       std::pair<std::vector<T>, int>>& coeffs)
 {
@@ -394,7 +394,7 @@ void set_diagonal(auto set_fn, std::span<const std::int32_t> rows,
 /// @param[in] diagonal The value to add to the diagonal for rows with a
 /// boundary condition applied
 template <typename T>
-void set_diagonal(auto set_fn, const fem::FunctionSpace& V,
+void set_diagonal(auto set_fn, const FunctionSpace& V,
                   const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
                   T diagonal = 1.0)
 {
diff --git a/cpp/dolfinx/fem/interpolate.cpp b/cpp/dolfinx/fem/interpolate.cpp
index f4cf64dfd99..a4fcd970154 100644
--- a/cpp/dolfinx/fem/interpolate.cpp
+++ b/cpp/dolfinx/fem/interpolate.cpp
@@ -69,40 +69,43 @@ fem::interpolation_coords(const FiniteElement& element, const mesh::Mesh& mesh,
 
   return x;
 }
-
+//-----------------------------------------------------------------------------
 fem::nmm_interpolation_data_t fem::create_nonmatching_meshes_interpolation_data(
     const fem::FunctionSpace& Vu, const fem::FunctionSpace& Vv,
     std::span<const std::int32_t> cells)
 {
-  std::vector<double> x;
-  auto mesh = Vu.mesh();
-  auto mesh_v = Vv.mesh();
-  auto element_u = Vu.element();
 
   // Collect all the points at which values are needed to define the
   // interpolating function
+  auto element_u = Vu.element();
+  assert(element_u);
+  auto mesh = Vu.mesh();
+  assert(mesh);
   const std::vector<double> coords_b
-      = fem::interpolation_coords(*element_u, *mesh, cells);
+      = interpolation_coords(*element_u, *mesh, cells);
 
   namespace stdex = std::experimental;
   using cmdspan2_t
       = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
   using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
   cmdspan2_t coords(coords_b.data(), 3, coords_b.size() / 3);
+
   // Transpose interpolation coords
-  x.resize(coords.size());
+  std::vector<double> x(coords.size());
   mdspan2_t _x(x.data(), coords_b.size() / 3, 3);
   for (std::size_t j = 0; j < coords.extent(1); ++j)
     for (std::size_t i = 0; i < 3; ++i)
       _x(j, i) = coords(i, j);
 
   // Determine ownership of each point
-  return dolfinx::geometry::determine_point_ownership(*mesh_v, x);
+  auto mesh_v = Vv.mesh();
+  assert(mesh_v);
+  return geometry::determine_point_ownership(*mesh_v, x);
 }
-
-fem::nmm_interpolation_data_t fem::create_nonmatching_meshes_interpolation_data(
-    const dolfinx::fem::FunctionSpace& Vu,
-    const dolfinx::fem::FunctionSpace& Vv)
+//-----------------------------------------------------------------------------
+fem::nmm_interpolation_data_t
+fem::create_nonmatching_meshes_interpolation_data(const FunctionSpace& Vu,
+                                                  const FunctionSpace& Vv)
 {
   assert(Vu.mesh());
   int tdim = Vu.mesh()->topology().dim();
@@ -111,7 +114,6 @@ fem::nmm_interpolation_data_t fem::create_nonmatching_meshes_interpolation_data(
   std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
   std::vector<std::int32_t> cells(num_cells, 0);
   std::iota(cells.begin(), cells.end(), 0);
-
   return create_nonmatching_meshes_interpolation_data(Vu, Vv, cells);
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index fc663d41774..dd532b0780b 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -916,20 +916,21 @@ void interpolate(Function<T>& u, std::span<const T> f,
   }
 }
 
-/// Generate data needed to interpolate discrete functions across different
-/// meshes
+/// Generate data needed to interpolate discrete functions across
+/// different meshes.
 ///
-/// @param[out] Vu The function space of the function to interpolate into
+/// @param[out] Vu The function space of the function to interpolate
+/// into
 /// @param[in] Vv The function space of the function to interpolate from
-/// @param[in] cells Indices of the cells in the destination mesh on which to
-/// interpolate. Should be the same as the list used when calling
-/// fem::interpolation_coords.
+/// @param[in] cells Indices of the cells in the destination mesh on
+/// which to interpolate. Should be the same as the list used when
+/// calling fem::interpolation_coords.
 nmm_interpolation_data_t create_nonmatching_meshes_interpolation_data(
     const FunctionSpace& Vu, const FunctionSpace& Vv,
     std::span<const std::int32_t> cells);
 
-/// Generate data needed to interpolate discrete functions defined on different
-/// meshes. Interpolate on all cells in the mesh.
+/// Generate data needed to interpolate discrete functions defined on
+/// different meshes. Interpolate on all cells in the mesh.
 ///
 /// @param[out] Vu The function space of the function to interpolate into
 /// @param[in] Vv The function space of the function to interpolate from
@@ -937,7 +938,6 @@ nmm_interpolation_data_t
 create_nonmatching_meshes_interpolation_data(const FunctionSpace& Vu,
                                              const FunctionSpace& Vv);
 
-//----------------------------------------------------------------------------
 /// Interpolate from one finite element Function to another one
 /// @param[out] u The function to interpolate into
 /// @param[in] v The function to be interpolated
diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp
index 8de2d3a12de..d76e8181069 100644
--- a/cpp/dolfinx/fem/utils.cpp
+++ b/cpp/dolfinx/fem/utils.cpp
@@ -12,6 +12,7 @@
 #include "Function.h"
 #include "FunctionSpace.h"
 #include "dofmapbuilder.h"
+#include <algorithm>
 #include <array>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/common/Timer.h>
@@ -197,8 +198,6 @@ fem::FunctionSpace fem::create_functionspace(
     const std::function<std::vector<int>(
         const graph::AdjacencyList<std::int32_t>&)>& reorder_fn)
 {
-  assert(mesh);
-
   // Create a DOLFINx element
   auto _e = std::make_shared<FiniteElement>(e, bs);
 
@@ -220,6 +219,7 @@ fem::FunctionSpace fem::create_functionspace(
   // Create a dofmap
   ElementDofLayout layout(bs, e.entity_dofs(), e.entity_closure_dofs(), {},
                           sub_doflayout);
+  assert(mesh);
   auto dofmap = std::make_shared<const DofMap>(
       create_dofmap(mesh->comm(), layout, mesh->topology(), reorder_fn, *_e));
 
@@ -266,3 +266,128 @@ fem::FunctionSpace fem::create_functionspace(
           mesh->comm(), layout, mesh->topology(), reorder_fn, *element)));
 }
 //-----------------------------------------------------------------------------
+std::vector<std::pair<int, std::vector<std::int32_t>>>
+fem::compute_integration_domains(fem::IntegralType integral_type,
+                                 const mesh::MeshTags<int>& meshtags)
+{
+  std::shared_ptr<const mesh::Mesh> mesh = meshtags.mesh();
+  assert(mesh);
+  const mesh::Topology& topology = mesh->topology();
+  const int tdim = topology.dim();
+  const int dim = integral_type == IntegralType::cell ? tdim : tdim - 1;
+  if (dim != meshtags.dim())
+  {
+    throw std::runtime_error("Invalid MeshTags dimension: "
+                             + std::to_string(meshtags.dim()));
+  }
+
+  std::span<const std::int32_t> entities = meshtags.indices();
+  std::span<const int> values = meshtags.values();
+  {
+    assert(topology.index_map(dim));
+    auto it0 = entities.begin();
+    auto it1 = std::lower_bound(it0, entities.end(),
+                                topology.index_map(dim)->size_local());
+    entities = entities.first(std::distance(it0, it1));
+    values = values.first(std::distance(it0, it1));
+  }
+
+  std::vector<std::int32_t> entity_data;
+  std::vector<int> values1;
+  switch (integral_type)
+  {
+  case IntegralType::cell:
+    entity_data.insert(entity_data.begin(), entities.begin(), entities.end());
+    values1.insert(values1.begin(), values.begin(), values.end());
+    break;
+  default:
+    mesh->topology_mutable().create_connectivity(dim, tdim);
+    mesh->topology_mutable().create_connectivity(tdim, dim);
+    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);
+    switch (integral_type)
+    {
+    case IntegralType::exterior_facet:
+    {
+      // Create list of tagged boundary facets
+      const std::vector bfacets = mesh::exterior_facet_indices(topology);
+      std::vector<std::int32_t> facets;
+      std::set_intersection(entities.begin(), entities.end(), bfacets.begin(),
+                            bfacets.end(), std::back_inserter(facets));
+      for (auto f : facets)
+      {
+        auto index_it = std::lower_bound(entities.begin(), entities.end(), f);
+        assert(index_it != entities.end() and *index_it == f);
+        std::size_t pos = std::distance(entities.begin(), index_it);
+        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());
+        values1.push_back(values[pos]);
+      }
+    }
+    break;
+    case IntegralType::interior_facet:
+    {
+      for (std::size_t j = 0; j < entities.size(); ++j)
+      {
+        const std::int32_t f = entities[j];
+        if (f_to_c->num_links(f) == 2)
+        {
+          // Get the facet as a pair of (cell, local facet) pairs, one
+          // for each cell
+          auto facets
+              = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
+          entity_data.insert(entity_data.end(), facets.begin(), facets.end());
+          values1.push_back(values[j]);
+        }
+      }
+    }
+    break;
+    default:
+      throw std::runtime_error(
+          "Cannot compute integration domains. Integral type not supported.");
+    }
+  }
+
+  // Build permutation that sorts by meshtag value
+  std::vector<std::int32_t> perm(values1.size());
+  std::iota(perm.begin(), perm.end(), 0);
+  std::stable_sort(perm.begin(), perm.end(),
+                   [&values1](auto p0, auto p1)
+                   { return values1[p0] < values1[p1]; });
+
+  std::size_t shape = 1;
+  if (integral_type == IntegralType::exterior_facet)
+    shape = 2;
+  else if (integral_type == IntegralType::interior_facet)
+    shape = 4;
+  std::vector<std::pair<int, std::vector<std::int32_t>>> integrals;
+  {
+    // Iterator to mark the start of the group
+    auto p0 = perm.begin();
+    while (p0 != perm.end())
+    {
+      auto id0 = values1[*p0];
+      auto p1 = std::find_if_not(p0, perm.end(),
+                                 [id0, &values1](auto idx)
+                                 { return id0 == values1[idx]; });
+
+      std::vector<std::int32_t> data;
+      data.reserve(shape * std::distance(p0, p1));
+      for (auto it = p0; it != p1; ++it)
+      {
+        auto e_it0 = std::next(entity_data.begin(), (*it) * shape);
+        auto e_it1 = std::next(e_it0, shape);
+        data.insert(data.end(), e_it0, e_it1);
+      }
+
+      integrals.push_back({id0, std::move(data)});
+      p0 = p1;
+    }
+  }
+
+  return integrals;
+}
+//-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 0b05a09c8e4..0a0f82b4258 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -103,6 +103,29 @@ using scalar_value_type_t = typename scalar_value_type<T>::value_type;
 
 } // namespace impl
 
+/// @brief Given an integral type and MeshTags, compute the entities
+/// that should be integrated over.
+///
+/// This function returns as list `[(id, entities)]`, where `entities`
+/// are the entities in `meshtags` with tag value `id`. For cell
+/// integrals `entities` are the cell indices. For exterior facet
+/// integrals, `entities` is a list of `(cell_index, local_facet_index)`
+/// pairs. For interior facet integrals, `entities` is a list of
+/// `(cell_index0, local_facet_index0, cell_index1,
+/// local_facet_index1)`.
+///
+/// @note Owned mesh entities only are returned. Ghost entities are not
+/// included.
+///
+/// @param[in] integral_type Integral type
+/// @param[in] meshtags The meshtags
+/// @return A list of (integral id, entities) pairs
+/// @pre The topological dimension of the integral entity type and the
+/// topological dimension of `meshtags` must be equal.
+std::vector<std::pair<int, std::vector<std::int32_t>>>
+compute_integration_domains(IntegralType integral_type,
+                            const mesh::MeshTags<int>& meshtags);
+
 /// @brief Finite element cell kernel concept.
 ///
 /// Kernel functions that can be passed to an assembler for execution
@@ -262,6 +285,7 @@ std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
 /// @param[in] coefficients Coefficient fields in the form
 /// @param[in] constants Spatial constants in the form
 /// @param[in] subdomains Subdomain markers
+/// @pre Each value in `subdomains` must be sorted by domain id
 /// @param[in] mesh The mesh of the domain
 template <typename T>
 Form<T> create_form(
@@ -269,7 +293,10 @@ Form<T> create_form(
     const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
     const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
     const std::vector<std::shared_ptr<const Constant<T>>>& constants,
-    const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
+    const std::map<
+        IntegralType,
+        std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>&
+        subdomains,
     std::shared_ptr<const mesh::Mesh> mesh = nullptr)
 {
   if (ufcx_form.rank != (int)spaces.size())
@@ -370,29 +397,25 @@ Form<T> create_form(
       assert(k);
 
       // Build list of entities to assembler over
-      std::vector<std::int32_t> e;
       if (id == -1)
       {
         // Default kernel, operates on all (owned) cells
         assert(topology.index_map(tdim));
+        std::vector<std::int32_t> e;
         e.resize(topology.index_map(tdim)->size_local(), 0);
         std::iota(e.begin(), e.end(), 0);
+        itg.first->second.emplace_back(id, k, std::move(e));
       }
-      else if (sd != subdomains.end() and sd->second)
+      else if (sd != subdomains.end())
       {
-        assert(topology.index_map(tdim));
-        std::span<const std::int32_t> entities = sd->second->indices();
-        std::span<const int> values = sd->second->values();
-        auto it0 = entities.begin();
-        auto it1 = std::lower_bound(it0, entities.end(),
-                                    topology.index_map(tdim)->size_local());
-        entities = entities.first(std::distance(it0, it1));
-        for (std::size_t j = 0; j < entities.size(); ++j)
-          if (values[j] == id)
-            e.push_back(entities[j]);
+        // NOTE: This requires that pairs are sorted
+        auto it = std::lower_bound(sd->second.begin(), sd->second.end(), id,
+                                   [](auto& pair, auto val)
+                                   { return pair.first < val; });
+        if (it != sd->second.end() and it->first == id)
+          itg.first->second.emplace_back(id, k, it->second);
       }
 
-      itg.first->second.emplace_back(id, k, std::move(e));
       if (integral->needs_facet_permutations)
         needs_facet_permutations = true;
     }
@@ -432,7 +455,6 @@ Form<T> create_form(
       assert(k);
 
       // Build list of entities to assembler over
-      std::vector<std::int32_t> e;
       const std::vector bfacets = mesh::exterior_facet_indices(topology);
       auto f_to_c = topology.connectivity(tdim - 1, tdim);
       assert(f_to_c);
@@ -441,6 +463,7 @@ Form<T> create_form(
       if (id == -1)
       {
         // Default kernel, operates on all (owned) exterior facets
+        std::vector<std::int32_t> e;
         e.reserve(2 * bfacets.size());
         for (std::int32_t f : bfacets)
         {
@@ -449,34 +472,18 @@ Form<T> create_form(
               = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
           e.insert(e.end(), pair.begin(), pair.end());
         }
+        itg.first->second.emplace_back(id, k, std::move(e));
       }
-      else if (sd != subdomains.end() and sd->second)
+      else if (sd != subdomains.end())
       {
-        // Create list of tagged boundary facets
-        std::span<const std::int32_t> entities = sd->second->indices();
-        auto it0 = entities.begin();
-        auto it1 = std::lower_bound(it0, entities.end(),
-                                    topology.index_map(tdim - 1)->size_local());
-        entities = entities.first(std::distance(it0, it1));
-        std::span<const int> values = sd->second->values();
-        std::vector<std::int32_t> facets;
-        std::set_intersection(entities.begin(), entities.end(), bfacets.begin(),
-                              bfacets.end(), std::back_inserter(facets));
-        for (auto f : facets)
-        {
-          auto index_it = std::lower_bound(entities.begin(), entities.end(), f);
-          assert(index_it != entities.end() and *index_it == f);
-          std::size_t pos = std::distance(entities.begin(), index_it);
-          if (values[pos] == id)
-          {
-            auto facet
-                = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
-            e.insert(e.end(), facet.begin(), facet.end());
-          }
-        }
+        // NOTE: This requires that pairs are sorted
+        auto it = std::lower_bound(sd->second.begin(), sd->second.end(), id,
+                                   [](auto& pair, auto val)
+                                   { return pair.first < val; });
+        if (it != sd->second.end() and it->first == id)
+          itg.first->second.emplace_back(id, k, it->second);
       }
 
-      itg.first->second.emplace_back(id, k, std::move(e));
       if (integral->needs_facet_permutations)
         needs_facet_permutations = true;
     }
@@ -516,7 +523,6 @@ Form<T> create_form(
       assert(k);
 
       // Build list of entities to assembler over
-      std::vector<std::int32_t> e;
       auto f_to_c = topology.connectivity(tdim - 1, tdim);
       assert(f_to_c);
       auto c_to_f = topology.connectivity(tdim, tdim - 1);
@@ -524,6 +530,7 @@ Form<T> create_form(
       if (id == -1)
       {
         // Default kernel, operates on all (owned) interior facets
+        std::vector<std::int32_t> e;
         assert(topology.index_map(tdim - 1));
         std::int32_t num_facets = topology.index_map(tdim - 1)->size_local();
         e.reserve(4 * num_facets);
@@ -536,35 +543,33 @@ Form<T> create_form(
             e.insert(e.end(), pairs.begin(), pairs.end());
           }
         }
+        itg.first->second.emplace_back(id, k, std::move(e));
       }
-      else if (sd != subdomains.end() and sd->second)
+      else if (sd != subdomains.end())
       {
-        std::span<const std::int32_t> entities = sd->second->indices();
-        auto it0 = entities.begin();
-        auto it1 = std::lower_bound(it0, entities.end(),
-                                    topology.index_map(tdim - 1)->size_local());
-        entities = entities.first(std::distance(it0, it1));
-        std::span<const int> values = sd->second->values();
-        for (std::size_t j = 0; j < entities.size(); ++j)
-        {
-          const std::int32_t f = entities[j];
-          if (f_to_c->num_links(f) == 2 and values[j] == id)
-          {
-            // Get the facet as a pair of (cell, local facet) pairs, one
-            // for each cell
-            auto facets
-                = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
-            e.insert(e.end(), facets.begin(), facets.end());
-          }
-        }
+        auto it = std::lower_bound(sd->second.begin(), sd->second.end(), id,
+                                   [](auto& pair, auto val)
+                                   { return pair.first < val; });
+        if (it != sd->second.end() and it->first == id)
+          itg.first->second.emplace_back(id, k, it->second);
       }
 
-      itg.first->second.emplace_back(id, k, std::move(e));
       if (integral->needs_facet_permutations)
         needs_facet_permutations = true;
     }
   }
 
+  std::map<IntegralType,
+           std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>
+      sd;
+  for (auto& [itg, data] : subdomains)
+  {
+    std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>> x;
+    for (auto& [id, idx] : data)
+      x.emplace_back(id, std::vector(idx.data(), idx.data() + idx.size()));
+    sd.insert({itg, std::move(x)});
+  }
+
   return Form<T>(spaces, integral_data, coefficients, constants,
                  needs_facet_permutations, mesh);
 }
@@ -575,6 +580,7 @@ Form<T> create_form(
 /// @param[in] coefficients Coefficient fields in the form (by name)
 /// @param[in] constants Spatial constants in the form (by name)
 /// @param[in] subdomains Subdomain makers
+/// @pre Each value in `subdomains` must be sorted by domain id
 /// @param[in] mesh The mesh of the domain. This is required if the form
 /// has no arguments, e.g. a functional
 /// @return A Form
@@ -585,7 +591,10 @@ Form<T> create_form(
     const std::map<std::string, std::shared_ptr<const Function<T>>>&
         coefficients,
     const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
-    const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
+    const std::map<
+        IntegralType,
+        std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>&
+        subdomains,
     std::shared_ptr<const mesh::Mesh> mesh = nullptr)
 {
   // Place coefficients in appropriate order
@@ -622,6 +631,7 @@ Form<T> create_form(
 /// @param[in] coefficients Coefficient fields in the form (by name)
 /// @param[in] constants Spatial constants in the form (by name)
 /// @param[in] subdomains Subdomain markers
+/// @pre Each value in `subdomains` must be sorted by domain id
 /// @param[in] mesh The mesh of the domain. This is required if the form
 /// has no arguments, e.g. a functional.
 /// @return A Form
@@ -632,7 +642,10 @@ Form<T> create_form(
     const std::map<std::string, std::shared_ptr<const Function<T>>>&
         coefficients,
     const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
-    const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
+    const std::map<
+        IntegralType,
+        std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>&
+        subdomains,
     std::shared_ptr<const mesh::Mesh> mesh = nullptr)
 {
   ufcx_form* form = fptr();
diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index d62d7be4105..7069b72268f 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -136,17 +136,9 @@ def _form(form):
         # Extract subdomain data from UFL form
         sd = form.subdomain_data()
         domain, = list(sd.keys())  # Assuming single domain
-
-        def unwrap_mt(t):
-            """Get subdomain data for each integral type."""
-            try:
-                return t._cpp_object
-            except AttributeError:
-                return t
         # Check that subdomain data for each integral type is the same
         for data in sd.get(domain).values():
             assert all([d is data[0] for d in data])
-        subdomains = {_ufl_to_dolfinx_domain[key]: unwrap_mt(mt[0]) for (key, mt) in sd.get(domain).items()}
 
         mesh = domain.ufl_cargo()
         if mesh is None:
@@ -165,6 +157,22 @@ def unwrap_mt(t):
                                   ]._cpp_object for i in range(ufcx_form.num_coefficients)]
         constants = [c._cpp_object for c in form.constants()]
 
+        # NOTE Could remove this and let the user convert meshtags by
+        # calling compute_integration_domains themselves
+        def get_integration_domains(integral_type, subdomain):
+            """Get integration domains from subdomain data"""
+            if subdomain is None:
+                return []
+            else:
+                try:
+                    return _cpp.fem.compute_integration_domains(integral_type, subdomain._cpp_object)
+                except AttributeError:
+                    return subdomain
+
+        # Subdomain markers (possibly empty list for some integral types)
+        subdomains = {_ufl_to_dolfinx_domain[key]: get_integration_domains(
+            _ufl_to_dolfinx_domain[key], subdomain_data[0]) for (key, subdomain_data) in sd.get(domain).items()}
+
         return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)
 
     def _create_form(form):
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 6da9a7778cf..f186fdf00f0 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -414,13 +414,34 @@ void declare_form(py::module& m, const std::string& type)
                   const std::vector<std::shared_ptr<
                       const dolfinx::fem::Constant<T>>>& constants,
                   const std::map<dolfinx::fem::IntegralType,
-                                 const dolfinx::mesh::MeshTags<int>*>&
+                                 std::vector<std::pair<
+                                     std::int32_t,
+                                     py::array_t<std::int32_t,
+                                                 py::array::c_style
+                                                     | py::array::forcecast>>>>&
                       subdomains,
                   std::shared_ptr<const dolfinx::mesh::Mesh> mesh)
                {
+                 std::map<dolfinx::fem::IntegralType,
+                          std::vector<std::pair<std::int32_t,
+                                                std::vector<std::int32_t>>>>
+                     sd;
+                 for (auto& [itg, data] : subdomains)
+                 {
+                   std::vector<
+                       std::pair<std::int32_t, std::vector<std::int32_t>>>
+                       x;
+                   for (auto& [id, e] : data)
+                   {
+                     x.emplace_back(id,
+                                    std::vector(e.data(), e.data() + e.size()));
+                   }
+                   sd.insert({itg, std::move(x)});
+                 }
+
                  ufcx_form* p = reinterpret_cast<ufcx_form*>(form);
-                 return dolfinx::fem::create_form<T>(
-                     *p, spaces, coefficients, constants, subdomains, mesh);
+                 return dolfinx::fem::create_form<T>(*p, spaces, coefficients,
+                                                     constants, sd, mesh);
                }),
            py::arg("form"), py::arg("spaces"), py::arg("coefficients"),
            py::arg("constants"), py::arg("subdomains"), py::arg("mesh"),
@@ -486,13 +507,33 @@ void declare_form(py::module& m, const std::string& type)
              coefficients,
          const std::vector<std::shared_ptr<const dolfinx::fem::Constant<T>>>&
              constants,
-         const std::map<dolfinx::fem::IntegralType,
-                        const dolfinx::mesh::MeshTags<int>*>& subdomains,
+         const std::map<
+             dolfinx::fem::IntegralType,
+             std::vector<std::pair<
+                 std::int32_t,
+                 py::array_t<std::int32_t,
+                             py::array::c_style | py::array::forcecast>>>>&
+             subdomains,
          std::shared_ptr<const dolfinx::mesh::Mesh> mesh)
       {
+        std::map<
+            dolfinx::fem::IntegralType,
+            std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>
+            sd;
+        for (auto& [itg, data] : subdomains)
+        {
+          std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>> x;
+          for (auto& [id, idx] : data)
+          {
+            x.emplace_back(id,
+                           std::vector(idx.data(), idx.data() + idx.size()));
+          }
+          sd.insert({itg, std::move(x)});
+        }
+
         ufcx_form* p = reinterpret_cast<ufcx_form*>(form);
         return dolfinx::fem::create_form<T>(*p, spaces, coefficients, constants,
-                                            subdomains, mesh);
+                                            sd, mesh);
       },
       py::arg("form"), py::arg("spaces"), py::arg("coefficients"),
       py::arg("constants"), py::arg("subdomains"), py::arg("mesh"),
@@ -578,6 +619,12 @@ void fem(py::module& m)
   m.def("transpose_dofmap", &dolfinx::fem::transpose_dofmap,
         "Build the index to (cell, local index) map from a "
         "dofmap ((cell, local index ) -> index).");
+  m.def(
+      "compute_integration_domains",
+      [](dolfinx::fem::IntegralType type,
+         const dolfinx::mesh::MeshTags<int>& meshtags)
+      { return dolfinx::fem::compute_integration_domains(type, meshtags); },
+      py::arg("integral_type"), py::arg("meshtags"));
   m.def(
       "create_nonmatching_meshes_interpolation_data",
       [](const dolfinx::fem::FunctionSpace& Vu,
diff --git a/python/test/unit/fem/test_assemble_domains.py b/python/test/unit/fem/test_assemble_domains.py
index 816ee014e08..845e6b6caf4 100644
--- a/python/test/unit/fem/test_assemble_domains.py
+++ b/python/test/unit/fem/test_assemble_domains.py
@@ -7,20 +7,19 @@
 
 import numpy as np
 import pytest
-
 import ufl
-from dolfinx import cpp as _cpp
 from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
                          dirichletbc, form)
 from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                                set_bc)
-from dolfinx.mesh import (GhostMode, Mesh, create_unit_square,
+from dolfinx.mesh import (GhostMode, Mesh, create_unit_square, locate_entities,
                           locate_entities_boundary, meshtags,
                           meshtags_from_entities)
-
 from mpi4py import MPI
 from petsc4py import PETSc
 
+from dolfinx import cpp as _cpp
+
 
 @pytest.fixture
 def mesh():
@@ -229,3 +228,107 @@ def test_additivity(mode):
     assert (J1 + J3) == pytest.approx(J13)
     assert (J2 + J3) == pytest.approx(J23)
     assert (J1 + J2 + J3) == pytest.approx(J123)
+
+
+def test_manual_integration_domains():
+    """Test that specifying integration domains manually i.e.
+    by passing a list of cell indices or (cell, local facet) pairs to
+    form gives the same result as the usual approach of tagging"""
+    n = 4
+    msh = create_unit_square(MPI.COMM_WORLD, n, n)
+
+    V = FunctionSpace(msh, ("Lagrange", 1))
+    u = ufl.TrialFunction(V)
+    v = ufl.TestFunction(V)
+
+    # Create meshtags to mark some cells
+    tdim = msh.topology.dim
+    cell_map = msh.topology.index_map(tdim)
+    num_cells = cell_map.size_local + cell_map.num_ghosts
+    cell_indices = np.arange(0, num_cells)
+    cell_values = np.zeros_like(cell_indices, dtype=np.intc)
+    marked_cells = locate_entities(msh, tdim, lambda x: x[0] < 0.75)
+    cell_values[marked_cells] = 7
+    mt_cells = meshtags(msh, tdim, cell_indices, cell_values)
+
+    # Create meshtags to mark some exterior facets
+    msh.topology.create_entities(tdim - 1)
+    facet_map = msh.topology.index_map(tdim - 1)
+    num_facets = facet_map.size_local + facet_map.num_ghosts
+    facet_indices = np.arange(0, num_facets)
+    facet_values = np.zeros_like(facet_indices, dtype=np.intc)
+    marked_ext_facets = locate_entities_boundary(msh, tdim - 1, lambda x: np.isclose(x[0], 0.0))
+    marked_int_facets = locate_entities(msh, tdim - 1, lambda x: x[0] < 0.75)
+    # marked_int_facets will also contain facets on the boundary, so set
+    # these values first, followed by the values for marked_ext_facets
+    facet_values[marked_int_facets] = 3
+    facet_values[marked_ext_facets] = 6
+    mt_facets = meshtags(msh, tdim - 1, facet_indices, facet_values)
+
+    # Create measures
+    dx_mt = ufl.Measure("dx", subdomain_data=mt_cells, domain=msh)
+    ds_mt = ufl.Measure("ds", subdomain_data=mt_facets, domain=msh)
+    dS_mt = ufl.Measure("dS", subdomain_data=mt_facets, domain=msh)
+
+    g = Function(V)
+    g.interpolate(lambda x: x[1]**2)
+
+    def create_forms(dx, ds, dS):
+        a = form(ufl.inner(g * u, v) * (dx(0) + dx(7) + ds(6)) + ufl.inner(g * u("+"), v("+") + v("-")) * dS(3))
+        L = form(ufl.inner(g, v) * (dx(0) + dx(7) + ds(6)) + ufl.inner(g, v("+") + v("-")) * dS(3))
+        return (a, L)
+
+    # Create forms and assemble
+    a, L = create_forms(dx_mt, ds_mt, dS_mt)
+    A_mt = assemble_matrix(a)
+    A_mt.assemble()
+    b_mt = assemble_vector(L)
+
+    # Manually specify cells to integrate over (removing ghosts
+    # to give same result as above)
+    cell_domains = [(domain_id, cell_indices[(cell_values == domain_id)
+                                             & (cell_indices < cell_map.size_local)])
+                    for domain_id in [0, 7]]
+
+    # Manually specify exterior facets to integrate over as
+    # (cell, local facet) pairs
+    ext_facet_domain = []
+    msh.topology.create_connectivity(tdim, tdim - 1)
+    msh.topology.create_connectivity(tdim - 1, tdim)
+    c_to_f = msh.topology.connectivity(tdim, tdim - 1)
+    f_to_c = msh.topology.connectivity(tdim - 1, tdim)
+    for f in marked_ext_facets:
+        if f < facet_map.size_local:
+            c = f_to_c.links(f)[0]
+            local_f = np.where(c_to_f.links(c) == f)[0][0]
+            ext_facet_domain.append(c)
+            ext_facet_domain.append(local_f)
+    ext_facet_domains = [(6, ext_facet_domain)]
+
+    # Manually specify interior facets to integrate over
+    int_facet_domain = []
+    for f in marked_int_facets:
+        if f >= facet_map.size_local or len(f_to_c.links(f)) != 2:
+            continue
+        c_0, c_1 = f_to_c.links(f)[0], f_to_c.links(f)[1]
+        local_f_0 = np.where(c_to_f.links(c_0) == f)[0][0]
+        local_f_1 = np.where(c_to_f.links(c_1) == f)[0][0]
+        int_facet_domain.append(c_0)
+        int_facet_domain.append(local_f_0)
+        int_facet_domain.append(c_1)
+        int_facet_domain.append(local_f_1)
+    int_facet_domains = [(3, int_facet_domain)]
+
+    # Create measures
+    dx_manual = ufl.Measure("dx", subdomain_data=cell_domains, domain=msh)
+    ds_manual = ufl.Measure("ds", subdomain_data=ext_facet_domains, domain=msh)
+    dS_manual = ufl.Measure("dS", subdomain_data=int_facet_domains, domain=msh)
+
+    # Assemble forms and check
+    a, L = create_forms(dx_manual, ds_manual, dS_manual)
+    A = assemble_matrix(a)
+    A.assemble()
+    b = assemble_vector(L)
+
+    assert np.isclose((A - A_mt).norm(), 0.0)
+    assert np.isclose((b - b_mt).norm(), 0.0)