diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp
index 891fad97779..6ddcfacabfb 100644
--- a/cpp/demo/custom_kernel/main.cpp
+++ b/cpp/demo/custom_kernel/main.cpp
@@ -75,7 +75,7 @@ double assemble_matrix0(std::shared_ptr<const fem::FunctionSpace<T>> V,
 {
   // Kernel data (ID, kernel function, cell indices to execute over)
   std::map integrals{
-      std::pair{std::tuple{fem::IntegralType::cell, -1, 0},
+      std::pair{std::tuple{fem::IntegralType::cell, 0, 0},
                 fem::integral_data<T>(kernel, cells, std::vector<int>{})}};
 
   fem::Form<T, T> a({V, V}, integrals, V->mesh(), {}, {}, false, {});
@@ -105,7 +105,7 @@ double assemble_vector0(std::shared_ptr<const fem::FunctionSpace<T>> V,
 {
   auto mesh = V->mesh();
   std::map integrals{
-      std::pair{std::tuple{fem::IntegralType::cell, -1, 0},
+      std::pair{std::tuple{fem::IntegralType::cell, 0, 0},
                 fem::integral_data<T>(kernel, cells, std::vector<int>{})}};
   fem::Form<T> L({V}, integrals, mesh, {}, {}, false, {});
   auto dofmap = V->dofmap();
diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h
index eb0b2c77569..8a361dcf0c8 100644
--- a/cpp/dolfinx/fem/Form.h
+++ b/cpp/dolfinx/fem/Form.h
@@ -129,9 +129,12 @@ class Form
   /// @param[in] V Function spaces for the form arguments, e.g. test and
   /// trial function spaces.
   /// @param[in] integrals Integrals in the form, where
-  /// `integrals[IntegralType, domain ID, kernel index]` returns the
-  /// integral (`integral_data`) of type `IntegralType` over domain `ID`
-  /// with kernel index `kernel index`.
+  /// `integrals[IntegralType, i, kernel index]` returns the `i`th integral
+  /// (`integral_data`) of type `IntegralType` with kernel index `kernel index`.
+  /// The `i`-index refers to the position of a kernel when flattened by
+  /// sorted subdomain ids, sorted by subdomain ids. The subdomain ids can
+  /// contain duplicate entries referring to different kernels over the same
+  /// subdomain.
   /// @param[in] coefficients Coefficients in the form.
   /// @param[in] constants Constants in the form.
   /// @param[in] mesh Mesh of the domain to integrate over (the
@@ -247,7 +250,7 @@ class Form
 
     for (auto& space : _function_spaces)
     {
-      // Working map: [integral type, domain ID, kernel_idx]->entities
+      // Working map: [integral type, integral_idx, kernel_idx]->entities
       std::map<std::tuple<IntegralType, int, int>,
                std::variant<std::vector<std::int32_t>,
                             std::span<const std::int32_t>>>
@@ -269,7 +272,7 @@ class Form
         bool inverse = emap.sub_topology() == mesh0->topology();
         for (auto& [key, itg] : _integrals)
         {
-          auto [type, id, kernel_idx] = key;
+          auto [type, idx, kernel_idx] = key;
           std::vector<std::int32_t> e;
           if (type == IntegralType::cell)
             e = emap.sub_topology_to_topology(itg.entities, inverse);
@@ -299,13 +302,13 @@ class Form
 
     for (auto& [key, integral] : _integrals)
     {
-      auto [type, id, kernel_idx] = key;
+      auto [type, idx, kernel_idx] = key;
       for (int c : integral.coeffs)
       {
         if (auto mesh0 = coefficients.at(c)->function_space()->mesh();
             mesh0 == _mesh)
         {
-          _cdata.insert({{type, id, c}, std::span(integral.entities)});
+          _cdata.insert({{type, idx, c}, std::span(integral.entities)});
         }
         else
         {
@@ -331,8 +334,7 @@ class Form
           }
           else
             throw std::runtime_error("Integral type not supported.");
-
-          _cdata.insert({{type, id, c}, std::move(e)});
+          _cdata.insert({{type, idx, c}, std::move(e)});
         }
       }
     }
@@ -370,6 +372,9 @@ class Form
   }
 
   /// @brief Get the kernel function for an integral.
+  ///
+  ///
+  ///
   /// @param[in] type Integral type.
   /// @param[in] id Integral subdomain ID.
   /// @param[in] kernel_idx Index of the kernel (we may have multiple
@@ -418,31 +423,29 @@ class Form
     return it->second.coeffs;
   }
 
-  /// @brief Get the IDs for integrals (kernels) for given integral
-  /// domain type.
-  ///
-  /// The IDs correspond to the domain IDs which the integrals are
-  /// defined for in the form. `ID=-1` is the default integral over the
-  /// whole domain.
+  /// @brief Get number of integrals (kernels) for a given integral type and
+  /// kernel index.
   ///
+  /// For a form containing two integrals of `integral_a` and `integral_b`
+  /// with subdomain-ids `(1, 4)` and `(3, 4, 5)` respectively, the integrals
+  /// are stored as a flattened list, sorted by sudomain-ids
+  /// ```cpp
+  /// auto form_integrals = {integral_a, integral_b, integral_a, integral_b,
+  /// integral_b}; auto form_integral_ids = {1, 3, 4, 4, 5}.
+  /// ```
   /// @param[in] type Integral type.
-  /// @return List of IDs for given integral type.
-  std::vector<int> integral_ids(IntegralType type) const
+  /// @param[in] kernel_idx Index of the kernel (we may have multiple
+  /// kernels for a integral type in mixed-topology meshes).
+  int num_integrals(IntegralType type, int kernel_idx) const
   {
-    std::vector<int> ids;
-    for (auto& [key, integral] : _integrals)
-    {
-      auto [t, id, kernel_idx] = key;
-      if (t == type)
-        ids.push_back(id);
-    }
 
-    // IDs may be repeated in mixed-topology meshes, so remove
-    // duplicates
-    std::sort(ids.begin(), ids.end());
-    auto it = std::unique(ids.begin(), ids.end());
-    ids.erase(it, ids.end());
-    return ids;
+    int count = std::count_if(_integrals.begin(), _integrals.end(),
+                              [type, kernel_idx](auto& x)
+                              {
+                                auto [t, id, k_idx] = x.first;
+                                return t == type and k_idx == kernel_idx;
+                              });
+    return count;
   }
 
   /// @brief Mesh entity indices to integrate over for a given integral
@@ -463,16 +466,23 @@ class Form
   /// is row-major.
   ///
   /// @param[in] type Integral type.
-  /// @param[in] id Integral domain identifier.
+  /// @param[in] idx Integral index in flattened list of integral kernels.
+  /// For a form containing two integrals of `integral_a` and `integral_b`
+  /// with subdomain-ids `(1, 4)` and `(3, 4, 5)` respectively, the integrals
+  /// are stored as a flattened list, sorted by sudomain-ids
+  /// ```cpp
+  /// auto form_integrals = {integral_a, integral_b, integral_a, integral_b,
+  /// integral_b}; auto form_integral_ids = {1, 3, 4, 4, 5}.
+  /// ```
   /// @param[in] kernel_idx Index of the kernel with in the domain (we
   /// may have multiple kernels for a given ID in mixed-topology
   /// meshes).
   /// @return Entity indices in the mesh::Mesh returned by mesh() to
   /// integrate over.
-  std::span<const std::int32_t> domain(IntegralType type, int id,
+  std::span<const std::int32_t> domain(IntegralType type, int idx,
                                        int kernel_idx) const
   {
-    auto it = _integrals.find({type, id, kernel_idx});
+    auto it = _integrals.find({type, idx, kernel_idx});
     if (it == _integrals.end())
       throw std::runtime_error("Requested domain not found.");
     return it->second.entities;
@@ -501,21 +511,21 @@ class Form
   /// may exist in one domain but not another. In this case, the entity
   /// is marked with -1.
   ///
-  /// @param type Integral type.
-  /// @param rank Argument index, e.g. `0` for the test function space, `1`
+  /// @param[in] type Integral type.
+  /// @param[in] rank Argument index, e.g. `0` for the test function space, `1`
   /// for the trial function space.
-  /// @param id Integral domain identifier.
-  /// @param kernel_idx Kernel index (cell type).
+  /// @param[in] idx Integral identifier.
+  /// @param[in] kernel_idx Kernel index (cell type).
   /// @return Entity indices in the argument function space mesh that is
   /// integrated over.
   /// - For cell integrals it has shape `(num_cells,)`.
   /// - For exterior/interior facet integrals, it has shape `(num_facts, 2)`
   /// (row-major storage), where `[i, 0]` is the index of a cell and
   /// `[i, 1]` is the local index of the facet relative to the cell.
-  std::span<const std::int32_t> domain_arg(IntegralType type, int rank, int id,
+  std::span<const std::int32_t> domain_arg(IntegralType type, int rank, int idx,
                                            int kernel_idx) const
   {
-    auto it = _edata.at(rank).find({type, id, kernel_idx});
+    auto it = _edata.at(rank).find({type, idx, kernel_idx});
     if (it == _edata.at(rank).end())
       throw std::runtime_error("Requested domain for argument not found.");
     try
@@ -533,19 +543,19 @@ class Form
   /// This method is equivalent to ::domain_arg, but returns mesh entity
   /// indices for coefficient \link Function Functions. \endlink
   ///
-  /// @param type Integral type.
-  /// @param id Integral identifier index.
-  /// @param c Coefficient index.
+  /// @param[in] type Integral type.
+  /// @param[in] idx Integral identifier.
+  /// @param[in] c Coefficient index.
   /// @return Entity indices in the coefficient function space mesh that
   /// is integrated over.
   /// - For cell integrals it has shape `(num_cells,)`.
   /// - For exterior/interior facet integrals, it has shape `(num_facts, 2)`
   /// (row-major storage), where `[i, 0]` is the index of a cell and
   /// `[i, 1]` is the local index of the facet relative to the cell.
-  std::span<const std::int32_t> domain_coeff(IntegralType type, int id,
+  std::span<const std::int32_t> domain_coeff(IntegralType type, int idx,
                                              int c) const
   {
-    auto it = _cdata.find({type, id, c});
+    auto it = _cdata.find({type, idx, c});
     if (it == _cdata.end())
       throw std::runtime_error("No domain for requested integral.");
     try
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index 98f432c7f4c..90351a05b7a 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/cpp/dolfinx/fem/assemble_matrix_impl.h
@@ -589,7 +589,7 @@ void assemble_matrix(
       cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
     }
 
-    for (int i : a.integral_ids(IntegralType::cell))
+    for (int i = 0; i < a.num_integrals(IntegralType::cell, cell_type_idx); ++i)
     {
       auto fn = a.kernel(IntegralType::cell, i, cell_type_idx);
       assert(fn);
@@ -618,7 +618,8 @@ void assemble_matrix(
                          num_facets_per_cell);
     }
 
-    for (int i : a.integral_ids(IntegralType::exterior_facet))
+    for (int i = 0;
+         i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i)
     {
       if (num_cell_types > 1)
       {
@@ -649,7 +650,8 @@ void assemble_matrix(
           cell_info0, cell_info1, perms);
     }
 
-    for (int i : a.integral_ids(IntegralType::interior_facet))
+    for (int i = 0;
+         i < a.num_integrals(IntegralType::interior_facet, cell_type_idx); ++i)
     {
       if (num_cell_types > 1)
       {
diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h
index 3aea777fa3e..85c1e2cd6ee 100644
--- a/cpp/dolfinx/fem/assemble_scalar_impl.h
+++ b/cpp/dolfinx/fem/assemble_scalar_impl.h
@@ -165,7 +165,7 @@ T assemble_scalar(
   assert(mesh);
 
   T value = 0;
-  for (int i : M.integral_ids(IntegralType::cell))
+  for (int i = 0; i < M.num_integrals(IntegralType::cell, 0); ++i)
   {
     auto fn = M.kernel(IntegralType::cell, i, 0);
     assert(fn);
@@ -190,7 +190,7 @@ T assemble_scalar(
                        num_facets_per_cell);
   }
 
-  for (int i : M.integral_ids(IntegralType::exterior_facet))
+  for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i)
   {
     auto fn = M.kernel(IntegralType::exterior_facet, i, 0);
     assert(fn);
@@ -208,7 +208,7 @@ T assemble_scalar(
         perms);
   }
 
-  for (int i : M.integral_ids(IntegralType::interior_facet))
+  for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i)
   {
     auto fn = M.kernel(IntegralType::interior_facet, i, 0);
     assert(fn);
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index 01cc2e1dc08..855aedccd42 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -1012,7 +1012,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
       = element1->template dof_transformation_right_fn<T>(
           doftransform::transpose);
 
-  for (int i : a.integral_ids(IntegralType::cell))
+  for (int i = 0; i < a.num_integrals(IntegralType::cell, 0); ++i)
   {
     auto kernel = a.kernel(IntegralType::cell, i, 0);
     assert(kernel);
@@ -1057,7 +1057,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
                        num_facets_per_cell);
   }
 
-  for (int i : a.integral_ids(IntegralType::exterior_facet))
+  for (int i = 0; i < a.num_integrals(IntegralType::exterior_facet, 0); ++i)
   {
     auto kernel = a.kernel(IntegralType::exterior_facet, i, 0);
     assert(kernel);
@@ -1081,7 +1081,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
         cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
   }
 
-  for (int i : a.integral_ids(IntegralType::interior_facet))
+  for (int i = 0; i < a.num_integrals(IntegralType::interior_facet, 0); ++i)
   {
     auto kernel = a.kernel(IntegralType::interior_facet, i, 0);
     assert(kernel);
@@ -1252,7 +1252,7 @@ void assemble_vector(
       cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
     }
 
-    for (int i : L.integral_ids(IntegralType::cell))
+    for (int i = 0; i < L.num_integrals(IntegralType::cell, 0); ++i)
     {
       auto fn = L.kernel(IntegralType::cell, i, cell_type_idx);
       assert(fn);
@@ -1297,7 +1297,7 @@ void assemble_vector(
         = md::mdspan<const std::int32_t,
                      md::extents<std::size_t, md::dynamic_extent, 2>>;
 
-    for (int i : L.integral_ids(IntegralType::exterior_facet))
+    for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i)
     {
       auto fn = L.kernel(IntegralType::exterior_facet, i, 0);
       assert(fn);
@@ -1331,7 +1331,7 @@ void assemble_vector(
       }
     }
 
-    for (int i : L.integral_ids(IntegralType::interior_facet))
+    for (int i = 0; i < L.num_integrals(IntegralType::interior_facet, 0); ++i)
     {
       using mdspanx22_t
           = md::mdspan<const std::int32_t,
diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h
index 7d8a30f7ff6..01721014475 100644
--- a/cpp/dolfinx/fem/pack.h
+++ b/cpp/dolfinx/fem/pack.h
@@ -200,10 +200,10 @@ allocate_coefficient_storage(const Form<T, U>& form)
   std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> coeffs;
   for (fem::IntegralType type : form.integral_types())
   {
-    for (int id : form.integral_ids(type))
+    for (int i = 0; i < form.num_integrals(type, 0); ++i)
     {
-      coeffs.emplace_hint(coeffs.end(), std::pair{type, id},
-                          allocate_coefficient_storage(form, type, id));
+      coeffs.emplace_hint(coeffs.end(), std::pair{type, i},
+                          allocate_coefficient_storage(form, type, i));
     }
   }
 
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index d1a0e71dee0..a0e67af34a1 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -228,34 +228,33 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form<T, U>& a)
     // Create and build sparsity pattern
     for (auto type : types)
     {
-      std::vector<int> ids = a.integral_ids(type);
       switch (type)
       {
       case IntegralType::cell:
-        for (int id : ids)
+        for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
         {
           sparsitybuild::cells(pattern,
-                               {a.domain_arg(type, 0, id, cell_type_idx),
-                                a.domain_arg(type, 1, id, cell_type_idx)},
+                               {a.domain_arg(type, 0, i, cell_type_idx),
+                                a.domain_arg(type, 1, i, cell_type_idx)},
                                {{dofmaps[0], dofmaps[1]}});
         }
         break;
       case IntegralType::interior_facet:
-        for (int id : ids)
+        for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
         {
           sparsitybuild::interior_facets(
               pattern,
-              {extract_cells(a.domain_arg(type, 0, id, 0)),
-               extract_cells(a.domain_arg(type, 1, id, 0))},
+              {extract_cells(a.domain_arg(type, 0, i, 0)),
+               extract_cells(a.domain_arg(type, 1, i, 0))},
               {{dofmaps[0], dofmaps[1]}});
         }
         break;
       case IntegralType::exterior_facet:
-        for (int id : ids)
+        for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
         {
           sparsitybuild::cells(pattern,
-                               {extract_cells(a.domain_arg(type, 0, id, 0)),
-                                extract_cells(a.domain_arg(type, 1, id, 0))},
+                               {extract_cells(a.domain_arg(type, 0, i, 0)),
+                                extract_cells(a.domain_arg(type, 1, i, 0))},
                                {{dofmaps[0], dofmaps[1]}});
         }
         break;
@@ -547,7 +546,7 @@ Form<T, U> create_form_factory(
           default_cells.resize(
               topology->index_maps(tdim).at(form_idx)->size_local(), 0);
           std::iota(default_cells.begin(), default_cells.end(), 0);
-          integrals.insert({{IntegralType::cell, id, form_idx},
+          integrals.insert({{IntegralType::cell, i, form_idx},
                             {k, default_cells, active_coeffs}});
         }
         else if (sd != subdomains.end())
@@ -557,7 +556,7 @@ Form<T, U> create_form_factory(
                                              [](auto& a) { return a.first; });
           if (it != sd->second.end() and it->first == id)
           {
-            integrals.insert({{IntegralType::cell, id, form_idx},
+            integrals.insert({{IntegralType::cell, i, form_idx},
                               {k,
                                std::vector<std::int32_t>(it->second.begin(),
                                                          it->second.end()),
@@ -639,7 +638,7 @@ Form<T, U> create_form_factory(
             default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
                                       pair.end());
           }
-          integrals.insert({{IntegralType::exterior_facet, id, form_idx},
+          integrals.insert({{IntegralType::exterior_facet, i, form_idx},
                             {k, default_facets_ext, active_coeffs}});
         }
         else if (sd != subdomains.end())
@@ -649,7 +648,7 @@ Form<T, U> create_form_factory(
                                              [](auto& a) { return a.first; });
           if (it != sd->second.end() and it->first == id)
           {
-            integrals.insert({{IntegralType::exterior_facet, id, form_idx},
+            integrals.insert({{IntegralType::exterior_facet, i, form_idx},
                               {k,
                                std::vector<std::int32_t>(it->second.begin(),
                                                          it->second.end()),
@@ -757,7 +756,7 @@ Form<T, U> create_form_factory(
                   "mesh");
             }
           }
-          integrals.insert({{IntegralType::interior_facet, id, form_idx},
+          integrals.insert({{IntegralType::interior_facet, i, form_idx},
                             {k, default_facets_int, active_coeffs}});
         }
         else if (sd != subdomains.end())
@@ -766,7 +765,7 @@ Form<T, U> create_form_factory(
                                              [](auto& a) { return a.first; });
           if (it != sd->second.end() and it->first == id)
           {
-            integrals.insert({{IntegralType::interior_facet, id, form_idx},
+            integrals.insert({{IntegralType::interior_facet, i, form_idx},
                               {k,
                                std::vector<std::int32_t>(it->second.begin(),
                                                          it->second.end()),
diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py
index e6addd1b47d..5be91e0eec6 100644
--- a/python/demo/demo_static-condensation.py
+++ b/python/demo/demo_static-condensation.py
@@ -183,7 +183,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, cu
 # Prepare a Form with a condensed tabulation kernel
 formtype = form_cpp_class(PETSc.ScalarType)  # type: ignore
 cells = np.arange(msh.topology.index_map(msh.topology.dim).size_local)
-integrals = {IntegralType.cell: [(-1, tabulate_A.address, cells, np.array([], dtype=np.int8))]}
+integrals = {IntegralType.cell: [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]}
 a_cond = Form(
     formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, [], mesh=msh._cpp_object)
 )
diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index 1b5f72e11a5..2d168620c16 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -12,7 +12,6 @@
 import typing
 from collections.abc import Iterable, Sequence
 from dataclasses import dataclass
-from itertools import chain
 
 from mpi4py import MPI
 
@@ -352,29 +351,15 @@ def _form(form):
         ]
         constants = [c._cpp_object for c in form.constants()]
 
-        # Make map from integral_type to subdomain id
+        # Extract subdomain ids from ufcx_form
         subdomain_ids = {type: [] for type in sd.get(domain).keys()}
-        for integral in 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":
-                    if isinstance(integral.subdomain_id(), tuple):
-                        ids = [sid for sid in integral.subdomain_id() if sid != "everywhere"]
-                    else:
-                        ids = [integral.subdomain_id()]
-                else:
-                    ids = []
-                subdomain_ids[integral.integral_type()].append(ids)
-
-        # Chain and sort subdomain ids
-        for itg_type, marker_ids in subdomain_ids.items():
-            flattened_ids = list(chain.from_iterable(marker_ids))
-            flattened_ids.sort()
-            subdomain_ids[itg_type] = flattened_ids
-
-        # Subdomain markers (possibly empty list for some integral
-        # types)
+        integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(4)]
+        for i in range(3):
+            integral_type = IntegralType(i)
+            for j in range(integral_offsets[i], integral_offsets[i + 1]):
+                subdomain_ids[integral_type.name].append(ufcx_form.form_integral_ids[j])
+
+        # 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], subdomain_ids[key]
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 7a85513ba35..fe8a366f4ba 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -757,12 +757,7 @@ void declare_form(nb::module_& m, std::string type)
       .def_prop_ro("mesh", &dolfinx::fem::Form<T, U>::mesh)
       .def_prop_ro("function_spaces",
                    &dolfinx::fem::Form<T, U>::function_spaces)
-      .def(
-          "integral_ids",
-          [](const dolfinx::fem::Form<T, U>& self,
-             dolfinx::fem::IntegralType type)
-          { return dolfinx_wrappers::as_nbarray(self.integral_ids(type)); },
-          nb::arg("type"))
+      .def_prop_ro("num_integrals", &dolfinx::fem::Form<T, U>::num_integrals)
       .def_prop_ro("integral_types", &dolfinx::fem::Form<T, U>::integral_types)
       .def_prop_ro("needs_facet_permutations",
                    &dolfinx::fem::Form<T, U>::needs_facet_permutations)
diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py
index 894a0ebbf10..8c8caecc672 100644
--- a/python/test/unit/fem/test_assembler.py
+++ b/python/test/unit/fem/test_assembler.py
@@ -1,4 +1,4 @@
-# Copyright (C) 2018-2022 Garth N. Wells
+# Copyright (C) 2018-2025 Garth N. Wells, Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -42,7 +42,9 @@
     create_unit_cube,
     create_unit_square,
     exterior_facet_indices,
+    locate_entities,
     locate_entities_boundary,
+    meshtags,
 )
 from ufl import derivative, dS, ds, dx, inner
 from ufl.geometry import SpatialCoordinate
@@ -1507,3 +1509,54 @@ def test_vector_types():
 
     assert np.linalg.norm(x0.array - x1.array) == pytest.approx(0.0)
     assert np.linalg.norm(x0.array - x2.array) == pytest.approx(0.0, abs=1e-7)
+
+
+@dtype_parametrize
+@pytest.mark.parametrize("method", ["degree", "metadata"])
+def test_mixed_quadrature(dtype, method):
+    xtype = dtype(0).real.dtype
+    mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, dtype=xtype)
+
+    V = functionspace(mesh, ("Lagrange", 1))
+    u = Function(V, dtype=dtype)
+    u.interpolate(lambda x: x[0])
+
+    tol = 500 * np.finfo(dtype).eps
+    num_cells_local = (
+        mesh.topology.index_map(mesh.topology.dim).size_local
+        + mesh.topology.index_map(mesh.topology.dim).num_ghosts
+    )
+    values = np.full(num_cells_local, 1, dtype=np.int32)
+    left_cells = locate_entities(mesh, mesh.topology.dim, lambda x: x[0] <= 0.5 + tol)
+    values[left_cells] = 2
+    top_cells = locate_entities(mesh, mesh.topology.dim, lambda x: x[1] >= 0.5 - tol)
+    values[top_cells] = 3
+    ct = meshtags(mesh, mesh.topology.dim, np.arange(num_cells_local, dtype=np.int32), values)
+
+    dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
+
+    if method == "degree":
+        dx_1 = dx(subdomain_id=(1,), degree=1)
+        dx_2 = dx(subdomain_id=(1, 2), degree=2)
+        dx_3 = dx(subdomain_id=(2, 3), degree=3)
+    elif method == "metadata":
+        dx_1 = dx(subdomain_id=(1,), metadata={"quadrature_degree": 1})
+        dx_2 = dx(subdomain_id=(1, 2), metadata={"quadrature_degree": 2})
+        dx_3 = dx(subdomain_id=(2, 3), metadata={"quadrature_degree": 3})
+    else:
+        raise ValueError(f"Invalid method {method}")
+    form_1 = u * dx_1
+    form_2 = u * dx_2
+    form_3 = u * dx_3
+    summed_form = form_1 + form_2 + form_3
+
+    compiled_forms = form([form_1, form_2, form_3], dtype=dtype)
+    local_contributions = 0
+    for compiled_form in compiled_forms:
+        local_contributions += assemble_scalar(compiled_form)
+    global_contribution = mesh.comm.allreduce(local_contributions, op=MPI.SUM)
+
+    compiled_form = form(summed_form, dtype=dtype)
+    local_sum = assemble_scalar(compiled_form)
+    global_sum = mesh.comm.allreduce(local_sum, op=MPI.SUM)
+    assert np.isclose(global_contribution, global_sum, rtol=tol, atol=tol)
diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py
index 9fec7e6706f..0f0a0acdd2c 100644
--- a/python/test/unit/fem/test_custom_jit_kernels.py
+++ b/python/test/unit/fem/test_custom_jit_kernels.py
@@ -94,9 +94,9 @@ def test_numba_assembly(dtype):
     active_coeffs = np.array([], dtype=np.int8)
     integrals = {
         IntegralType.cell: [
-            (-1, k2.address, cells, active_coeffs),
+            (0, k2.address, cells, active_coeffs),
+            (1, k2.address, np.arange(0), active_coeffs),
             (2, k2.address, np.arange(0), active_coeffs),
-            (12, k2.address, np.arange(0), active_coeffs),
         ]
     }
     formtype = form_cpp_class(dtype)
@@ -105,7 +105,7 @@ def test_numba_assembly(dtype):
             [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object
         )
     )
-    integrals = {IntegralType.cell: [(-1, k1.address, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, k1.address, cells, active_coeffs)]}
     L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
 
     A = dolfinx.fem.assemble_matrix(a)
@@ -136,7 +136,7 @@ def test_coefficient(dtype):
     num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
     active_coeffs = np.array([0], dtype=np.int8)
     integrals = {
-        IntegralType.cell: [(1, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)]
+        IntegralType.cell: [(0, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)]
     }
     formtype = form_cpp_class(dtype)
     L = Form(
@@ -276,7 +276,7 @@ def test_cffi_assembly():
 
     ptrA = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonA"))
     active_coeffs = np.array([], dtype=np.int8)
-    integrals = {IntegralType.cell: [(-1, ptrA, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, ptrA, cells, active_coeffs)]}
     a = Form(
         _cpp.fem.Form_float64(
             [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object
@@ -284,11 +284,10 @@ def test_cffi_assembly():
     )
 
     ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL"))
-    integrals = {IntegralType.cell: [(-1, ptrL, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, ptrL, cells, active_coeffs)]}
     L = Form(
         _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object)
     )
-
     A = fem.assemble_matrix(a)
     A.scatter_reverse()
     assert np.isclose(np.sqrt(A.squared_norm()), 56.124860801609124)