diff --git a/cpp/demo/biharmonic/CMakeLists.txt b/cpp/demo/biharmonic/CMakeLists.txt
index c373efed2ef..aff7feda83a 100644
--- a/cpp/demo/biharmonic/CMakeLists.txt
+++ b/cpp/demo/biharmonic/CMakeLists.txt
@@ -45,9 +45,7 @@ add_custom_command(
 
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-add_executable(
-  ${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/biharmonic.c
-)
+add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/biharmonic.c)
 target_link_libraries(${PROJECT_NAME} dolfinx)
 
 # Do not throw error for 'multi-line comments' (these are typical in rst which
diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
index 2f2746b29c2..a9a61d8631f 100644
--- a/cpp/demo/biharmonic/main.cpp
+++ b/cpp/demo/biharmonic/main.cpp
@@ -191,9 +191,9 @@ int main(int argc, char* argv[])
 
     //  Define variational forms
     auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}));
+        *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}));
+        fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}, {}));
 
     //  Now, the Dirichlet boundary condition ($u = 0$) can be
     //  created using the class {cpp:class}`DirichletBC`. A
diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp
new file mode 100644
index 00000000000..c1797e60e66
--- /dev/null
+++ b/cpp/demo/codim_0_assembly/main.cpp
@@ -0,0 +1,173 @@
+// # Mixed assembly with a function mesh on a subset of cells
+//
+// This demo illustrates how to:
+//
+// * Create a submesh of co-dimension 0
+// * Assemble a mixed formulation with function spaces defined on the sub mesh
+// and parent mesh
+
+#include "mixed_codim0.h"
+#include <basix/finite-element.h>
+#include <cmath>
+#include <dolfinx.h>
+#include <dolfinx/fem/Constant.h>
+#include <dolfinx/fem/petsc.h>
+#include <dolfinx/la/MatrixCSR.h>
+#include <dolfinx/la/SparsityPattern.h>
+#include <utility>
+#include <vector>
+
+using namespace dolfinx;
+using T = PetscScalar;
+using U = typename dolfinx::scalar_value_type_t<T>;
+
+int main(int argc, char* argv[])
+{
+  dolfinx::init_logging(argc, argv);
+  PetscInitialize(&argc, &argv, nullptr, nullptr);
+
+  {
+    // Create mesh and function space
+    auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
+    auto mesh = std::make_shared<mesh::Mesh<U>>(
+        mesh::create_rectangle<U>(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}},
+                                  {1, 4}, mesh::CellType::quadrilateral, part));
+
+    auto element = basix::create_element<U>(
+        basix::element::family::P, basix::cell::type::quadrilateral, 1,
+        basix::element::lagrange_variant::unset,
+        basix::element::dpc_variant::unset, false);
+
+    auto V = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace(mesh, element, {}));
+
+    // Next we find all cells of the mesh with y<0.5
+    const int tdim = mesh->topology()->dim();
+    auto marked_cells = mesh::locate_entities(
+        *mesh, tdim,
+        [](auto x)
+        {
+          using U = typename decltype(x)::value_type;
+          constexpr U eps = 1.0e-8;
+          std::vector<std::int8_t> marker(x.extent(1), false);
+          for (std::size_t p = 0; p < x.extent(1); ++p)
+          {
+            auto y = x(1, p);
+            if (std::abs(y) <= 0.5 + eps)
+              marker[p] = true;
+          }
+          return marker;
+        });
+    // We create a MeshTags object where we mark these cells with 2, and any
+    // other cell with 1
+    auto cell_map = mesh->topology()->index_map(tdim);
+    std::size_t num_cells_local
+        = mesh->topology()->index_map(tdim)->size_local()
+          + mesh->topology()->index_map(tdim)->num_ghosts();
+    std::vector<std::int32_t> cells(num_cells_local);
+    std::iota(cells.begin(), cells.end(), 0);
+    std::vector<std::int32_t> values(cell_map->size_local(), 1);
+    std::for_each(marked_cells.begin(), marked_cells.end(),
+                  [&values](auto& c) { values[c] = 2; });
+    dolfinx::mesh::MeshTags<std::int32_t> cell_marker(mesh->topology(), tdim,
+                                                      cells, values);
+
+    std::shared_ptr<mesh::Mesh<U>> submesh;
+    std::vector<std::int32_t> submesh_to_mesh;
+    {
+      auto [_submesh, _submesh_to_mesh, v_map, g_map]
+          = mesh::create_submesh(*mesh, tdim, cell_marker.find(2));
+      submesh = std::make_shared<mesh::Mesh<U>>(std::move(_submesh));
+      submesh_to_mesh = std::move(_submesh_to_mesh);
+    }
+
+    // We create the function space used for the trial space
+    auto W = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace(submesh, element, {}));
+
+    // A mixed-domain form has functions defined over different meshes. The mesh
+    // associated with the measure (dx, ds, etc.) is called the integration
+    // domain. To assemble mixed-domain forms, maps must be provided taking
+    // entities in the integration domain to entities on each mesh in the form.
+    // Since one of our forms has a measure defined over `mesh` and involves a
+    // function defined over `submesh`, we must provide a map from entities in
+    // `mesh` to entities in `submesh`. This is simply the "inverse" of
+    // `submesh_to_mesh`.
+    std::vector<std::int32_t> mesh_to_submesh(num_cells_local, -1);
+    for (std::size_t i = 0; i < submesh_to_mesh.size(); ++i)
+      mesh_to_submesh[submesh_to_mesh[i]] = i;
+
+    std::shared_ptr<const mesh::Mesh<U>> const_ptr = submesh;
+    std::map<std::shared_ptr<const mesh::Mesh<U>>,
+             std::span<const std::int32_t>>
+        entity_maps
+        = {{const_ptr, std::span<const std::int32_t>(mesh_to_submesh.data(),
+                                                     mesh_to_submesh.size())}};
+
+    // Next we compute the integration entities on the integration domain `mesh`
+    std::map<
+        fem::IntegralType,
+        std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
+        subdomain_map = {};
+    auto integration_entities = fem::compute_integration_domains(
+        fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2), tdim);
+
+    subdomain_map[fem::IntegralType::cell].push_back(
+        {3, std::span<const std::int32_t>(integration_entities.data(),
+                                          integration_entities.size())});
+
+    // We can now create the bi-linear form
+    auto a_mixed = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
+                            subdomain_map, entity_maps, V->mesh()));
+
+    la::SparsityPattern sp_mixed = fem::create_sparsity_pattern(*a_mixed);
+    sp_mixed.finalize();
+    la::MatrixCSR<double> A_mixed(sp_mixed);
+    fem::assemble_matrix(A_mixed.mat_add_values(), *a_mixed, {});
+    A_mixed.scatter_rev();
+
+    auto a = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_mixed_codim0_a, {W, W}, {}, {}, {}, {}));
+    la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
+    sp.finalize();
+
+    la::MatrixCSR<double> A(sp);
+    fem::assemble_matrix(A.mat_add_values(), *a, {});
+    A.scatter_rev();
+
+    std::vector<T> A_mixed_flattened = A_mixed.to_dense();
+    std::stringstream cc;
+    cc.precision(3);
+    cc << "A_mixed:" << std::endl;
+
+    std::size_t num_owned_rows = V->dofmap()->index_map->size_local();
+    std::size_t num_sub_cols = W->dofmap()->index_map->size_local()
+                               + W->dofmap()->index_map->num_ghosts();
+    for (std::size_t i = 0; i < num_owned_rows; i++)
+    {
+      for (std::size_t j = 0; j < num_sub_cols; j++)
+      {
+        cc << A_mixed_flattened[i * num_sub_cols + j] << " ";
+      }
+      cc << std::endl;
+    }
+
+    std::size_t num_owned_sub_rows = W->dofmap()->index_map->size_local();
+    std::vector<T> A_flattened = A.to_dense();
+    cc << "A" << std::endl;
+    for (std::size_t i = 0; i < num_owned_sub_rows; i++)
+    {
+      for (std::size_t j = 0; j < num_sub_cols; j++)
+      {
+        cc << A_flattened[i * num_sub_cols + j] << " ";
+      }
+      cc << std::endl;
+    }
+    std::cout << cc.str() << std::endl;
+  }
+
+  PetscFinalize();
+
+  return 0;
+}
diff --git a/cpp/demo/codim_0_assembly/mixed_codim0.py b/cpp/demo/codim_0_assembly/mixed_codim0.py
new file mode 100644
index 00000000000..e69660b2bc5
--- /dev/null
+++ b/cpp/demo/codim_0_assembly/mixed_codim0.py
@@ -0,0 +1,35 @@
+# This demo aims to illustrate how to assemble a matrix with a trial function
+# defined on a submesh of co-dimension 0, and a test function defined on the parent mesh
+from basix.ufl import element
+from ufl import (
+    FunctionSpace,
+    Mesh,
+    TestFunction,
+    TrialFunction,
+    dx,
+)
+
+cell = "quadrilateral"
+coord_element = element("Lagrange", cell, 1, shape=(2,))
+mesh = Mesh(coord_element)
+
+# We define the function space and test function on the full mesh
+e = element("Lagrange", cell, 1)
+V = FunctionSpace(mesh, e)
+v = TestFunction(V)
+
+# Next we define the sub-mesh
+submesh = Mesh(coord_element)
+W = FunctionSpace(submesh, e)
+p = TrialFunction(W)
+
+# And finally we define a "mass matrix" on the submesh, with the test function
+# of the parent mesh. The integration domain is the parent mesh, but we restrict integration
+# to all cells marked with subdomain_id=3, which will indicate what cells of our mesh is part
+# of the submesh
+a_mixed = p * v * dx(domain=mesh, subdomain_id=3)
+
+q = TestFunction(W)
+a = p * q * dx(domain=submesh)
+
+forms = [a_mixed, a]
diff --git a/cpp/demo/hyperelasticity/CMakeLists.txt b/cpp/demo/hyperelasticity/CMakeLists.txt
index ce398bbe127..3b2a0104817 100644
--- a/cpp/demo/hyperelasticity/CMakeLists.txt
+++ b/cpp/demo/hyperelasticity/CMakeLists.txt
@@ -40,9 +40,7 @@ else()
 
   set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-  add_executable(
-    ${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c
-  )
+  add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
   target_link_libraries(${PROJECT_NAME} dolfinx)
 
   # Do not throw error for 'multi-line comments' (these are typical in rst which
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 2ae716bce3b..9734dd640d6 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -175,10 +175,10 @@ int main(int argc, char* argv[])
     auto u = std::make_shared<fem::Function<T>>(V);
     auto a = std::make_shared<fem::Form<T>>(
         fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}));
+                            {{"B", B}, {"T", traction}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(
         fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}));
+                            {{"B", B}, {"T", traction}}, {}, {}));
 
     auto u_rotation = std::make_shared<fem::Function<T>>(V);
     u_rotation->interpolate(
diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
index 64644229b61..01f14418285 100644
--- a/cpp/demo/poisson/main.cpp
+++ b/cpp/demo/poisson/main.cpp
@@ -140,9 +140,9 @@ int main(int argc, char* argv[])
 
     // Define variational forms
     auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_poisson_a, {V, V}, {}, {{"kappa", kappa}}, {}));
+        *form_poisson_a, {V, V}, {}, {{"kappa", kappa}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_poisson_L, {V}, {{"f", f}, {"g", g}}, {}, {}));
+        *form_poisson_L, {V}, {{"f", f}, {"g", g}}, {}, {}, {}));
 
     //  Now, the Dirichlet boundary condition ($u = 0$) can be created
     //  using the class {cpp:class}`DirichletBC`. A
diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp
index d3d699388aa..31d2c6fc531 100644
--- a/cpp/demo/poisson_matrix_free/main.cpp
+++ b/cpp/demo/poisson_matrix_free/main.cpp
@@ -151,12 +151,12 @@ void solver(MPI_Comm comm)
 
   // Define variational forms
   auto L = std::make_shared<fem::Form<T, U>>(
-      fem::create_form<T>(*form_poisson_L, {V}, {}, {{"f", f}}, {}));
+      fem::create_form<T>(*form_poisson_L, {V}, {}, {{"f", f}}, {}, {}));
 
   // Action of the bilinear form "a" on a function ui
   auto ui = std::make_shared<fem::Function<T, U>>(V);
   auto M = std::make_shared<fem::Form<T, U>>(
-      fem::create_form<T>(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}));
+      fem::create_form<T>(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}, {}));
 
   // Define boundary condition
   auto u_D = std::make_shared<fem::Function<T, U>>(V);
@@ -231,7 +231,7 @@ void solver(MPI_Comm comm)
   // Compute L2 error (squared) of the solution vector e = (u - u_d, u
   // - u_d)*dx
   auto E = std::make_shared<fem::Form<T>>(fem::create_form<T, U>(
-      *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, mesh));
+      *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, {}, mesh));
   T error = fem::assemble_scalar(*E);
   if (dolfinx::MPI::rank(comm) == 0)
   {
diff --git a/cpp/dolfinx/fem/DofMap.cpp b/cpp/dolfinx/fem/DofMap.cpp
index 871adce490d..5c10ac8d178 100644
--- a/cpp/dolfinx/fem/DofMap.cpp
+++ b/cpp/dolfinx/fem/DofMap.cpp
@@ -78,7 +78,8 @@ fem::DofMap build_collapsed_dofmap(const DofMap& dofmap_view,
   }
 
   // Create a map from old dofs to new dofs
-  std::vector<std::int32_t> old_to_new(dofs_view.back() + bs_view, -1);
+  std::size_t array_size = dofs_view.empty() ? 0 : dofs_view.back() + bs_view;
+  std::vector<std::int32_t> old_to_new(array_size, -1);
   for (std::size_t new_idx = 0; new_idx < sub_imap_to_imap.size(); ++new_idx)
   {
     for (int k = 0; k < bs_view; ++k)
@@ -215,7 +216,6 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
     reorder_fn = [](const graph::AdjacencyList<std::int32_t>& g)
     { return graph::reorder_gps(g); };
   }
-
   // Create new dofmap
   auto create_subdofmap = [](MPI_Comm comm, auto index_map_bs, auto& layout,
                              auto& topology, auto& reorder_fn, auto& dmap)
@@ -227,7 +227,6 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
 
       // Create new element dof layout and reset parent
       ElementDofLayout collapsed_dof_layout = layout.copy();
-
       auto [_index_map, bs, dofmaps] = build_dofmap_data(
           comm, topology, {collapsed_dof_layout}, reorder_fn);
       auto index_map
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index 9f11829526d..6687a1acff4 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -1061,8 +1061,6 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear forms, where a[j] is the form that
 /// generates A_j
-/// @param[in] x_dofmap Mesh geometry dofmap
-/// @param[in] x Mesh coordinates
 /// @param[in] constants Constants that appear in `a`
 /// @param[in] coeffs Coefficients that appear in `a`
 /// @param[in] bcs1 List of boundary conditions for each block, i.e.
@@ -1073,7 +1071,6 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
 template <dolfinx::scalar T, std::floating_point U>
 void apply_lifting(
     std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>> a,
-    mdspan2_t x_dofmap, std::span<const scalar_value_type_t<T>> x,
     const std::vector<std::span<const T>>& constants,
     const std::vector<std::map<std::pair<IntegralType, int>,
                                std::pair<std::span<const T>, int>>>& coeffs,
@@ -1081,7 +1078,6 @@ void apply_lifting(
         bcs1,
     const std::vector<std::span<const T>>& x0, T scale)
 {
-  // FIXME: make changes to reactivate this check
   if (!x0.empty() and x0.size() != a.size())
   {
     throw std::runtime_error(
@@ -1100,6 +1096,13 @@ void apply_lifting(
     std::vector<T> bc_values1;
     if (a[j] and !bcs1[j].empty())
     {
+      // Extract data from mesh
+      std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->mesh();
+      if (!mesh)
+        throw std::runtime_error("Unable to extract a mesh.");
+      mdspan2_t x_dofmap = mesh->geometry().dofmap();
+      auto x = mesh->geometry().x();
+
       assert(a[j]->function_spaces().at(0));
 
       auto V1 = a[j]->function_spaces()[1];
diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
index e50cb1b30e8..c53aef84b40 100644
--- a/cpp/dolfinx/fem/assembler.h
+++ b/cpp/dolfinx/fem/assembler.h
@@ -161,31 +161,7 @@ void apply_lifting(
   if (std::ranges::all_of(a, [](auto ptr) { return ptr == nullptr; }))
     return;
 
-  std::shared_ptr<const mesh::Mesh<U>> mesh;
-  for (auto& a_i : a)
-  {
-    if (a_i and !mesh)
-      mesh = a_i->mesh();
-    if (a_i and mesh and a_i->mesh() != mesh)
-      throw std::runtime_error("Mismatch between meshes.");
-  }
-
-  if (!mesh)
-    throw std::runtime_error("Unable to extract a mesh.");
-
-  if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
-  {
-    impl::apply_lifting<T>(b, a, mesh->geometry().dofmap(),
-                           mesh->geometry().x(), constants, coeffs, bcs1, x0,
-                           scale);
-  }
-  else
-  {
-    auto x = mesh->geometry().x();
-    std::vector<scalar_value_type_t<T>> _x(x.begin(), x.end());
-    impl::apply_lifting<T>(b, a, mesh->geometry().dofmap(), _x, constants,
-                           coeffs, bcs1, x0, scale);
-  }
+  impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, scale);
 }
 
 /// Modify b such that:
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 619a9e02d7c..3b3c717bcda 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -83,23 +83,24 @@ get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells,
 
 } // namespace impl
 
-/// @brief Given an integral type and mesh tag data, compute the
+/// @brief Given an integral type and a set of entities, 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
+/// This function returns a list `[(id, entities)]`. 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)`.
+/// `id` refers to the subdomain id used in the definition of the integration
+/// measures of the variational form.
 ///
 /// @note Owned mesh entities only are returned. Ghost entities are not
 /// included.
 ///
 /// @param[in] integral_type Integral type
 /// @param[in] topology Mesh topology
-/// @param[in] entities List of tagged mesh entities
-/// @param[in] dim Topological dimension of tagged entities
+/// @param[in] entities List of mesh entities
+/// @param[in] dim Topological dimension of entities
 /// @return List of integration entities
 /// @pre For facet integrals, the topology facet-to-cell and
 /// cell-to-facet connectivity must be computed before calling this
@@ -671,6 +672,8 @@ Form<T, U> create_form_factory(
 /// @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] entity_maps The entity maps for the form. Empty for
+/// single domain problems.
 /// @param[in] mesh Mesh of the domain. This is required if the form has
 /// no arguments, e.g. a functional.
 /// @return A Form
@@ -685,6 +688,8 @@ Form<T, U> create_form(
         IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
         subdomains,
+    const std::map<std::shared_ptr<const mesh::Mesh<U>>,
+                   std::span<const std::int32_t>>& entity_maps,
     std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
 {
   // Place coefficients in appropriate order
@@ -712,7 +717,7 @@ Form<T, U> create_form(
   }
 
   return create_form_factory(ufcx_form, spaces, coeff_map, const_map,
-                             subdomains, {}, mesh);
+                             subdomains, entity_maps, mesh);
 }
 
 /// @brief Create a Form using a factory function that returns a pointer
@@ -727,6 +732,8 @@ Form<T, U> create_form(
 /// @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] entity_maps The entity maps for the form. Empty for
+/// single domain problems.
 /// @param[in] mesh Mesh of the domain. This is required if the form has
 /// no arguments, e.g. a functional.
 /// @return A Form
@@ -741,11 +748,13 @@ Form<T, U> create_form(
         IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
         subdomains,
+    const std::map<std::shared_ptr<const mesh::Mesh<U>>,
+                   std::span<const std::int32_t>>& entity_maps,
     std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
 {
   ufcx_form* form = fptr();
   Form<T, U> L = create_form<T, U>(*form, spaces, coefficients, constants,
-                                   subdomains, mesh);
+                                   subdomains, entity_maps, mesh);
   std::free(form);
   return L;
 }
@@ -775,6 +784,10 @@ FunctionSpace<T> create_functionspace(
         "Cannot specify value shape for non-scalar base element.");
   }
 
+  if (mesh::cell_type_from_basix_type(e.cell_type())
+      != mesh->topology()->cell_type())
+    throw std::runtime_error("Cell type of element and mesh must match.");
+
   std::size_t bs = value_shape.empty()
                        ? 1
                        : std::accumulate(value_shape.begin(), value_shape.end(),
diff --git a/cpp/dolfinx/la/SparsityPattern.cpp b/cpp/dolfinx/la/SparsityPattern.cpp
index 1dab9d35d26..ac7ee1b6efc 100644
--- a/cpp/dolfinx/la/SparsityPattern.cpp
+++ b/cpp/dolfinx/la/SparsityPattern.cpp
@@ -140,6 +140,27 @@ SparsityPattern::SparsityPattern(
   }
 }
 //-----------------------------------------------------------------------------
+void SparsityPattern::insert(std::int32_t row, std::int32_t col)
+{
+  if (!_offsets.empty())
+  {
+    throw std::runtime_error(
+        "Cannot insert into sparsity pattern. It has already been finalized");
+  }
+
+  assert(_index_maps[0]);
+  const std::int32_t max_row
+      = _index_maps[0]->size_local() + _index_maps[0]->num_ghosts() - 1;
+
+  if (row > max_row or row < 0)
+  {
+    throw std::runtime_error(
+        "Cannot insert rows that do not exist in the IndexMap.");
+  }
+
+  _row_cache[row].push_back(col);
+}
+//-----------------------------------------------------------------------------
 void SparsityPattern::insert(std::span<const std::int32_t> rows,
                              std::span<const std::int32_t> cols)
 {
@@ -160,7 +181,6 @@ void SparsityPattern::insert(std::span<const std::int32_t> rows,
       throw std::runtime_error(
           "Cannot insert rows that do not exist in the IndexMap.");
     }
-
     _row_cache[row].insert(_row_cache[row].end(), cols.begin(), cols.end());
   }
 }
diff --git a/cpp/dolfinx/la/SparsityPattern.h b/cpp/dolfinx/la/SparsityPattern.h
index 33632d72813..726455a0c39 100644
--- a/cpp/dolfinx/la/SparsityPattern.h
+++ b/cpp/dolfinx/la/SparsityPattern.h
@@ -67,6 +67,19 @@ class SparsityPattern
 
   /// @brief Insert non-zero locations using local (process-wise)
   /// indices.
+  /// @param[in] row local row index
+  /// @param[in] col local column index
+  void insert(std::int32_t row, std::int32_t col);
+
+  /// @brief Insert non-zero locations using local (process-wise)
+  /// indices.
+  ///
+  /// This routine inserts non-zero locations at the outer product of rows and
+  /// cols into the sparsity pattern, i.e. adds the matrix entries at
+  ///   A[row[i], col[j]] for all i, j.
+  ///
+  /// @param[in] rows list of the local row indices 
+  /// @param[in] cols list of the local column indices
   void insert(std::span<const std::int32_t> rows,
               std::span<const std::int32_t> cols);
 
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index 53ece3f893f..e652537d80c 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -40,6 +40,7 @@ add_executable(
   common/sub_systems_manager.cpp
   common/index_map.cpp
   common/sort.cpp
+  fem/functionspace.cpp
   mesh/distributed_mesh.cpp
   common/CIFailure.cpp
   refinement/interval.cpp
diff --git a/cpp/test/fem/functionspace.cpp b/cpp/test/fem/functionspace.cpp
new file mode 100644
index 00000000000..faf3b7233dd
--- /dev/null
+++ b/cpp/test/fem/functionspace.cpp
@@ -0,0 +1,30 @@
+// Copyright (C) 2024 Paul Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <basix/finite-element.h>
+
+#include <dolfinx/fem/FunctionSpace.h>
+#include <dolfinx/fem/utils.h>
+#include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/mesh/generation.h>
+
+using namespace dolfinx;
+
+TEST_CASE("Create Function Space (mismatch of elements)", "[functionspace]")
+{
+  auto mesh = std::make_shared<mesh::Mesh<double>>(
+      dolfinx::mesh::create_rectangle<double>(
+          MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle));
+
+  auto element = basix::create_element<double>(
+      basix::element::family::P, basix::cell::type::interval, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
+  CHECK_THROWS(fem::create_functionspace<double>(mesh, element, {}));
+}
diff --git a/cpp/test/matrix.cpp b/cpp/test/matrix.cpp
index 3c7ed9dba8d..705a804768f 100644
--- a/cpp/test/matrix.cpp
+++ b/cpp/test/matrix.cpp
@@ -124,7 +124,7 @@ la::MatrixCSR<double> create_operator(MPI_Comm comm)
   auto kappa = std::make_shared<fem::Constant<double>>(2.0);
   auto a = std::make_shared<fem::Form<double, double>>(
       fem::create_form<double, double>(*form_poisson_a, {V, V}, {},
-                                       {{"kappa", kappa}}, {}));
+                                       {{"kappa", kappa}}, {}, {}));
 
   la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
   sp.finalize();
@@ -165,7 +165,7 @@ la::MatrixCSR<double> create_operator(MPI_Comm comm)
   // Define variational forms
   auto a = std::make_shared<fem::Form<double, double>>(
       fem::create_form<double, double>(*form_poisson_a, {V, V}, {},
-                                       {{"kappa", kappa}}, {}));
+                                       {{"kappa", kappa}}, {}, {}));
 
   // Create sparsity pattern
   la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
@@ -201,9 +201,9 @@ void test_matrix()
 {
   auto map0 = std::make_shared<common::IndexMap>(MPI_COMM_SELF, 8);
   la::SparsityPattern p(MPI_COMM_SELF, {map0, map0}, {1, 1});
-  p.insert(std::vector{0}, std::vector{0});
-  p.insert(std::vector{4}, std::vector{5});
-  p.insert(std::vector{5}, std::vector{4});
+  p.insert(0, 0);
+  p.insert(4, 5);
+  p.insert(5, 4);
   p.finalize();
 
   using T = float;
diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py
index 60052a2953d..a29759a9008 100644
--- a/python/dolfinx/fem/__init__.py
+++ b/python/dolfinx/fem/__init__.py
@@ -9,10 +9,12 @@
 import numpy.typing as npt
 
 from dolfinx.cpp.fem import IntegralType, transpose_dofmap
+from dolfinx.cpp.fem import compute_integration_domains as _compute_integration_domains
 from dolfinx.cpp.fem import create_interpolation_data as _create_interpolation_data
 from dolfinx.cpp.fem import create_sparsity_pattern as _create_sparsity_pattern
 from dolfinx.cpp.fem import discrete_gradient as _discrete_gradient
 from dolfinx.cpp.fem import interpolation_matrix as _interpolation_matrix
+from dolfinx.cpp.mesh import Topology
 from dolfinx.fem.assemble import (
     apply_lifting,
     assemble_matrix,
@@ -124,12 +126,45 @@ def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _Matri
     return _MatrixCSR(_interpolation_matrix(space0._cpp_object, space1._cpp_object))
 
 
+def compute_integration_domains(
+    integral_type: IntegralType, topology: Topology, entities: np.ndarray, dim: int
+):
+    """Given an integral type and a set of entities compute integration entities.
+
+    This function returns a list `[(id, entities)]`.
+    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)`.
+    `id` refers to the subdomain id used in the definition of the integration
+    measures of the variational form.
+
+    Note:
+        Owned mesh entities only are returned. Ghost entities are not included.
+
+    Note:
+        For facet integrals, the topology facet-to-cell and
+        cell-to-facet connectivity must be computed before calling this function.
+
+    Args:
+        integral_type: Integral type
+        topology: Mesh topology
+        entities: List of mesh entities
+        dim: Topological dimension of entities
+
+    Returns:
+        List of integration entities
+    """
+    return _compute_integration_domains(integral_type, topology, entities, dim)
+
+
 __all__ = [
     "Constant",
     "Expression",
     "Function",
     "ElementMetaData",
     "create_matrix",
+    "compute_integration_domains",
     "functionspace",
     "FunctionSpace",
     "create_sparsity_pattern",
diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index 6c4ed21311c..59821fd1817 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -190,7 +190,7 @@ def form(
     dtype: npt.DTypeLike = default_scalar_type,
     form_compiler_options: typing.Optional[dict] = None,
     jit_options: typing.Optional[dict] = None,
-    entity_maps: dict[Mesh, np.typing.NDArray[np.int32]] = {},
+    entity_maps: typing.Optional[dict[Mesh, np.typing.NDArray[np.int32]]] = None,
 ):
     """Create a Form or an array of Forms.
 
@@ -282,7 +282,10 @@ def _form(form):
             for (key, subdomain_data) in sd.get(domain).items()
         }
 
-        _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()}
+        if entity_maps is None:
+            _entity_maps = dict()
+        else:
+            _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()}
 
         f = ftype(
             module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)),
@@ -432,8 +435,10 @@ def create_form(
     form: CompiledForm,
     function_spaces: list[function.FunctionSpace],
     mesh: Mesh,
+    subdomains: dict[IntegralType, list[tuple[int, np.ndarray]]],
     coefficient_map: dict[ufl.Function, function.Function],
     constant_map: dict[ufl.Constant, function.Constant],
+    entity_maps: dict[Mesh, np.typing.NDArray[np.int32]] | None = None,
 ) -> Form:
     """
     Create a Form object from a data-independent compiled form
@@ -443,45 +448,24 @@ def create_form(
         function_spaces: List of function spaces associated with the form.
             Should match the number of arguments in the form.
         mesh: Mesh to associate form with
+        subdomains: A map from integral type to a list of pairs, where each pair corresponds to a
+            subdomain id and the set of of integration entities to integrate over.
+            Can be computed with {py:func}`dolfinx.fem.compute_integration_domains`.
         coefficient_map: Map from UFL coefficient to function with data
         constant_map: Map from UFL constant to constant with data
+        entity_map: A map where each key corresponds to a mesh different to the integration
+            domain `mesh`.
+        The value of the map is an array of integers, where the i-th entry is the entity in
+            the key mesh.
     """
-    sd = form.ufl_form.subdomain_data()
-    (domain,) = list(sd.keys())  # Assuming single domain
-
-    # Make map from integral_type to subdomain id
-    subdomain_ids: dict[IntegralType, list[list[int]]] = {
-        type: [] for type in sd.get(domain).keys()
-    }
-    flattened_subdomain_ids: dict[IntegralType, list[int]] = {
-        type: [] for type in sd.get(domain).keys()
-    }
-    for integral in form.ufl_form.integrals():
-        if integral.subdomain_data() is not None:
-            # Subdomain ids can be strings, its or tuples with strings and ints
-            if integral.subdomain_id() != "everywhere":
-                try:
-                    ids = [sid for sid in integral.subdomain_id() if sid != "everywhere"]
-                except TypeError:
-                    # If not tuple, but single integer id
-                    ids = [integral.subdomain_id()]
-            else:
-                ids = []
-            subdomain_ids[integral.integral_type()].append(ids)  # type: ignore
-
-        # Chain and sort subdomain ids
-        for itg_type, marker_ids in subdomain_ids.items():
-            flattened_ids: list[int] = list(chain.from_iterable(marker_ids))
-            flattened_ids.sort()
-            flattened_subdomain_ids[itg_type] = flattened_ids
+    if entity_maps is None:
+        _entity_maps = {}
+    else:
+        _entity_maps = {m._cpp_object: emap for (m, emap) in entity_maps.items()}
 
-    # 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], flattened_subdomain_ids[key]
-        )
-        for (key, subdomain_data) in sd.get(domain).items()
-    }
+    _subdomain_data = subdomains.copy()
+    for _, idomain in _subdomain_data.items():
+        idomain.sort(key=lambda x: x[0])
 
     # Extract name of ufl objects and map them to their corresponding C++ object
     ufl_coefficients = ufl.algorithms.extract_coefficients(form.ufl_form)
@@ -497,7 +481,8 @@ def create_form(
         [fs._cpp_object for fs in function_spaces],
         coefficients,
         constants,
-        subdomains,
+        _subdomain_data,
+        _entity_maps,
         mesh._cpp_object,
     )
     return Form(f, form.ufcx_form, form.code)
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 5f1a0da94bb..cf04e746fe6 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -331,7 +331,8 @@ def refine(
         mesh1 = _cpp.refinement.refine(mesh._cpp_object, redistribute)
     else:
         mesh1 = _cpp.refinement.refine(mesh._cpp_object, edges, redistribute)
-    return Mesh(mesh1, mesh._ufl_domain)
+    ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element())  # type: ignore
+    return Mesh(mesh1, ufl_domain)
 
 
 def refine_interval(
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 52495815515..4725847162a 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -797,10 +797,14 @@ void declare_form(nb::module_& m, std::string type)
          const std::map<std::string,
                         std::shared_ptr<const dolfinx::fem::Constant<T>>>&
              constants,
-         const std::map<dolfinx::fem::IntegralType,
-                        std::vector<std::pair<std::int32_t,
-                                              std::span<const std::int32_t>>>>&
+         const std::map<
+             dolfinx::fem::IntegralType,
+             std::vector<std::pair<
+                 std::int32_t, nb::ndarray<const std::int32_t, nb::c_contig>>>>&
              subdomains,
+         const std::map<std::shared_ptr<const dolfinx::mesh::Mesh<U>>,
+                        nb::ndarray<const std::int32_t, nb::c_contig>>&
+             entity_maps,
          std::shared_ptr<const dolfinx::mesh::Mesh<U>> mesh = nullptr)
       {
         std::map<
@@ -814,14 +818,18 @@ void declare_form(nb::module_& m, std::string type)
             x.emplace_back(id, std::span(idx.data(), idx.size()));
           sd.insert({itg, std::move(x)});
         }
-
+        std::map<std::shared_ptr<const dolfinx::mesh::Mesh<U>>,
+                 std::span<const int32_t>>
+            _entity_maps;
+        for (auto& [msh, map] : entity_maps)
+          _entity_maps.emplace(msh, std::span(map.data(), map.size()));
         ufcx_form* p = reinterpret_cast<ufcx_form*>(form);
-        return dolfinx::fem::create_form<T, U>(*p, spaces, coefficients,
-                                               constants, sd, mesh);
+        return dolfinx::fem::create_form<T, U>(
+            *p, spaces, coefficients, constants, sd, _entity_maps, mesh);
       },
       nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"),
-      nb::arg("constants"), nb::arg("subdomains"), nb::arg("mesh"),
-      "Create Form from a pointer to ufcx_form.");
+      nb::arg("constants"), nb::arg("subdomains"), nb::arg("entity_maps"),
+      nb::arg("mesh"), "Create Form from a pointer to ufcx_form.");
 
   m.def("create_sparsity_pattern",
         &dolfinx::fem ::create_sparsity_pattern<T, U>, nb::arg("a"),
diff --git a/python/dolfinx/wrappers/la.cpp b/python/dolfinx/wrappers/la.cpp
index e84590e31d1..5a1f7865f9a 100644
--- a/python/dolfinx/wrappers/la.cpp
+++ b/python/dolfinx/wrappers/la.cpp
@@ -8,6 +8,7 @@
 #include "caster_mpi.h"
 #include "numpy_dtype.h"
 #include <complex>
+#include <cstdint>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/la/MatrixCSR.h>
 #include <dolfinx/la/SparsityPattern.h>
@@ -265,6 +266,10 @@ void la(nb::module_& m)
                         std::span(cols.data(), cols.size()));
           },
           nb::arg("rows"), nb::arg("cols"))
+      .def("insert",
+           nb::overload_cast<int32_t, int32_t>(
+               &dolfinx::la::SparsityPattern::insert),
+           nb::arg("row"), nb::arg("col"))
       .def(
           "insert_diagonal",
           [](dolfinx::la::SparsityPattern& self,
diff --git a/python/test/unit/fem/test_assemble_mesh_independent_form.py b/python/test/unit/fem/test_assemble_mesh_independent_form.py
index 3359041f4e8..aacd934bf7e 100644
--- a/python/test/unit/fem/test_assemble_mesh_independent_form.py
+++ b/python/test/unit/fem/test_assemble_mesh_independent_form.py
@@ -47,10 +47,111 @@ def create_and_integrate(N, compiled_form):
         wh.interpolate(lambda x: x[1])
         eh = dolfinx.fem.Constant(mesh, dtype(3.0))
         ch = dolfinx.fem.Constant(mesh, dtype(2.0))
-        form = dolfinx.fem.create_form(compiled_form, [], mesh, {u: uh, w: wh}, {c: ch, e: eh})
+        form = dolfinx.fem.create_form(compiled_form, [], mesh, {}, {u: uh, w: wh}, {c: ch, e: eh})
         assert np.isclose(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(form), op=MPI.SUM), 1.5)
 
     # Create various meshes, that all uses this compiled form with a map from ufl
     # to dolfinx functions and constants
     for i in range(1, 4):
         create_and_integrate(i, compiled_form)
+
+
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(np.complex64, marks=pytest.mark.xfail_win32_complex),
+        pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex),
+    ],
+)
+def test_submesh_assembly(dtype):
+    """
+    Compile a form without an associated mesh and assemble a form over a sequence of meshes
+    """
+    real_type = dtype(0).real.dtype
+    c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,), dtype=real_type)
+    domain = ufl.Mesh(c_el)
+    el = basix.ufl.element("Lagrange", "triangle", 2, dtype=real_type)
+    V = ufl.FunctionSpace(domain, el)
+    u = ufl.TestFunction(V)
+
+    f_el = basix.ufl.element("Lagrange", "interval", 1, shape=(2,), dtype=real_type)
+    submesh = ufl.Mesh(f_el)
+    sub_el = basix.ufl.element("Lagrange", "interval", 3, dtype=real_type)
+    V_sub = ufl.FunctionSpace(submesh, sub_el)
+
+    w = ufl.Coefficient(V_sub)
+
+    subdomain_id = 3
+    J = ufl.inner(w, u) * ufl.ds(domain=domain, subdomain_id=subdomain_id)
+
+    # Compile form using dolfinx.jit.ffcx_jit
+    compiled_form = dolfinx.fem.compile_form(
+        MPI.COMM_WORLD, J, form_compiler_options={"scalar_type": dtype}
+    )
+
+    def create_and_integrate(N, compiled_form):
+        mesh = dolfinx.mesh.create_rectangle(
+            MPI.COMM_WORLD,
+            [np.array([0, 0]), np.array([2, 2])],
+            [N, N],
+            dolfinx.mesh.CellType.triangle,
+            dtype=real_type,
+        )
+        assert mesh.ufl_domain().ufl_coordinate_element() == c_el
+
+        facets = dolfinx.mesh.locate_entities_boundary(
+            mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[1], 2)
+        )
+        submesh, sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
+            mesh, mesh.topology.dim - 1, facets
+        )
+        imap = mesh.topology.index_map(mesh.topology.dim - 1)
+        num_facets = imap.size_local + imap.num_ghosts
+        parent_to_sub = np.full(num_facets, -1, dtype=np.int32)
+        parent_to_sub[sub_to_parent] = np.arange(sub_to_parent.size, dtype=np.int32)
+
+        def g(x):
+            return -3 * x[1] ** 3 + x[0]
+
+        Vh = dolfinx.fem.functionspace(mesh, u.ufl_element())
+
+        Wh = dolfinx.fem.functionspace(submesh, w.ufl_element())
+        wh = dolfinx.fem.Function(Wh, dtype=dtype)
+        wh.interpolate(g)
+
+        facet_entities = dolfinx.fem.compute_integration_domains(
+            dolfinx.fem.IntegralType.exterior_facet,
+            mesh.topology,
+            sub_to_parent,
+            mesh.topology.dim - 1,
+        )
+        subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(subdomain_id, facet_entities)]}
+        form = dolfinx.fem.create_form(
+            compiled_form, [Vh], mesh, subdomains, {w: wh}, {}, {submesh: parent_to_sub}
+        )
+
+        # Compute exact solution
+        x = ufl.SpatialCoordinate(mesh)
+        ff = dolfinx.mesh.meshtags(
+            mesh, mesh.topology.dim - 1, facets, np.full(len(facets), subdomain_id, dtype=np.int32)
+        )
+        vh = ufl.TestFunction(Vh)
+        ex_solution = dolfinx.fem.assemble_vector(
+            dolfinx.fem.form(
+                ufl.inner(g(x), vh)
+                * ufl.ds(domain=mesh, subdomain_data=ff, subdomain_id=subdomain_id),
+                dtype=dtype,
+            )
+        )
+        ex_solution.scatter_reverse(dolfinx.la.InsertMode.add)
+        bh = dolfinx.fem.assemble_vector(form)
+        bh.scatter_reverse(dolfinx.la.InsertMode.add)
+        tol = float(5e2 * np.finfo(dtype).resolution)
+        np.testing.assert_allclose(ex_solution.array, bh.array, atol=tol)
+
+    # Create various meshes, that all uses this compiled form with a map from ufl
+    # to dolfinx functions and constants
+    for i in range(1, 4):
+        create_and_integrate(i, compiled_form)
diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py
index 7f69a433f21..8764810d5cf 100644
--- a/python/test/unit/fem/test_assemble_submesh.py
+++ b/python/test/unit/fem/test_assemble_submesh.py
@@ -384,3 +384,59 @@ def coeff_expr(x):
 
 
 # TODO Test random mesh and interior facets
+
+
+@pytest.mark.petsc4py
+def test_mixed_measures():
+    """Test block assembly of forms where the integration measure in each
+    block may be different"""
+    from dolfinx.fem.petsc import assemble_vector_block
+
+    comm = MPI.COMM_WORLD
+    msh = create_unit_square(comm, 16, 21, ghost_mode=GhostMode.none)
+
+    # Create a submesh of some cells
+    tdim = msh.topology.dim
+    smsh_cells = locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
+    smsh, smsh_to_msh = create_submesh(msh, tdim, smsh_cells)[:2]
+
+    # Create function spaces over each mesh
+    V = fem.functionspace(msh, ("Lagrange", 1))
+    Q = fem.functionspace(smsh, ("Lagrange", 1))
+
+    # Define two integration measures, one over the mesh, the other over the submesh
+    dx_msh = ufl.Measure("dx", msh, subdomain_data=[(1, smsh_cells)])
+    dx_smsh = ufl.Measure("dx", smsh)
+
+    # Trial and test functions
+    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
+    p, q = ufl.TrialFunction(Q), ufl.TestFunction(Q)
+
+    # First, assemble a block vector using both dx_msh and dx_smsh
+    a = [
+        [
+            fem.form(ufl.inner(u, v) * dx_msh),
+            fem.form(ufl.inner(p, v) * dx_smsh, entity_maps={msh: smsh_to_msh}),
+        ],
+        [
+            fem.form(ufl.inner(u, q) * dx_smsh, entity_maps={msh: smsh_to_msh}),
+            fem.form(ufl.inner(p, q) * dx_smsh),
+        ],
+    ]
+    L = [fem.form(ufl.inner(2.3, v) * dx_msh), fem.form(ufl.inner(1.3, q) * dx_smsh)]
+    b0 = assemble_vector_block(L, a)
+
+    # Now, assemble the same vector using only dx_msh
+    cell_imap = msh.topology.index_map(tdim)
+    num_cells = cell_imap.size_local + cell_imap.num_ghosts
+    msh_to_smsh = np.full(num_cells, -1)
+    msh_to_smsh[smsh_to_msh] = np.arange(len(smsh_to_msh))
+    entity_maps = {smsh: msh_to_smsh}
+    L = [
+        fem.form(ufl.inner(2.3, v) * dx_msh),
+        fem.form(ufl.inner(1.3, q) * dx_msh(1), entity_maps=entity_maps),
+    ]
+    b1 = assemble_vector_block(L, a)
+
+    # Check the results are the same
+    assert np.allclose(b0.norm(), b1.norm())
diff --git a/python/test/unit/fem/test_dofmap.py b/python/test/unit/fem/test_dofmap.py
index aa1ee0b4d21..85489d82fb0 100644
--- a/python/test/unit/fem/test_dofmap.py
+++ b/python/test/unit/fem/test_dofmap.py
@@ -421,3 +421,26 @@ def test_transpose_dofmap():
     dofmap = np.array([[0, 2, 1], [3, 2, 1], [4, 3, 1]], dtype=np.int32)
     transpose = dolfinx.fem.transpose_dofmap(dofmap, 3)
     assert np.array_equal(transpose.array, [0, 2, 5, 8, 1, 4, 3, 7, 6])
+
+
+def test_empty_rank_collapse():
+    """Test that dofmap with no dofs on a rank can be collapsed"""
+    if MPI.COMM_WORLD.rank == 0:
+        nodes = np.array([[0.0], [1.0], [2.0]], dtype=np.float64)
+        cells = np.array([[0, 1], [1, 2]], dtype=np.int64)
+    else:
+        nodes = np.empty((0, 1), dtype=np.float64)
+        cells = np.empty((0, 2), dtype=np.int64)
+    c_el = element("Lagrange", "interval", 1, shape=(1,))
+
+    def self_partitioner(comm: MPI.Intracomm, n, m, topo):
+        dests = np.full(len(topo[0]) // 2, comm.rank, dtype=np.int32)
+        offsets = np.arange(len(topo[0]) // 2 + 1, dtype=np.int32)
+        return dolfinx.graph.adjacencylist(dests, offsets)
+
+    mesh = create_mesh(MPI.COMM_WORLD, cells, nodes, c_el, partitioner=self_partitioner)
+
+    el = element("Lagrange", "interval", 1, shape=(2,))
+    V = functionspace(mesh, el)
+    V_0, _ = V.sub(0).collapse()
+    assert V.dofmap.index_map.size_local == V_0.dofmap.index_map.size_local
diff --git a/python/test/unit/la/test_matrix_csr.py b/python/test/unit/la/test_matrix_csr.py
index 738b2b46b86..1f8e5c35d09 100644
--- a/python/test/unit/la/test_matrix_csr.py
+++ b/python/test/unit/la/test_matrix_csr.py
@@ -25,9 +25,9 @@ def create_test_sparsity(n, bs):
     if bs == 1:
         for i in range(2):
             for j in range(2):
-                sp.insert(np.array([2 + i]), np.array([4 + j]))
+                sp.insert(2 + i, 4 + j)
     elif bs == 2:
-        sp.insert(np.array([1]), np.array([2]))
+        sp.insert(1, 2)
     sp.finalize()
     return sp
 
@@ -126,10 +126,10 @@ def test_distributed_csr(dtype):
     sp = SparsityPattern(MPI.COMM_WORLD, [im, im], [1, 1])
     for i in range(n):
         for j in range(n + nghost):
-            sp.insert(np.array([i]), np.array([j]))
+            sp.insert(i, j)
     for i in range(n, n + nghost):
         for j in range(n, n + nghost):
-            sp.insert(np.array([i]), np.array([j]))
+            sp.insert(i, j)
     sp.finalize()
 
     mat = matrix_csr(sp, dtype=dtype)
diff --git a/python/test/unit/mesh/test_ghost_mesh.py b/python/test/unit/mesh/test_ghost_mesh.py
index 7232df8d386..489092ed40a 100644
--- a/python/test/unit/mesh/test_ghost_mesh.py
+++ b/python/test/unit/mesh/test_ghost_mesh.py
@@ -95,13 +95,11 @@ def test_ghost_connectivities(mode):
     map_f = topology.index_map(tdim - 1)
     num_facets = map_f.size_local + map_f.num_ghosts
 
-    reference = dict()
     meshR.topology.create_connectivity(tdim - 1, tdim)
     facet_mp = compute_midpoints(meshR, tdim - 1, np.arange(num_facets))
     meshR.topology.create_connectivity(tdim, tdim)
     cell_mp = compute_midpoints(meshR, tdim, np.arange(num_cells))
     reference = {tuple(row): [] for row in facet_mp}
-
     for i in range(num_facets):
         for cidx in meshR.topology.connectivity(1, 2).links(i):
             reference[tuple(facet_mp[i])].append(cell_mp[cidx].tolist())