diff --git a/.github/workflows/fenicsx-refs.env b/.github/workflows/fenicsx-refs.env
index 22c7e44a7e..693a736666 100644
--- a/.github/workflows/fenicsx-refs.env
+++ b/.github/workflows/fenicsx-refs.env
@@ -3,4 +3,4 @@ basix_ref=main
 ufl_repository=FEniCS/ufl
 ufl_ref=main
 ffcx_repository=FEniCS/ffcx
-ffcx_ref=main
+ffcx_ref=sclaus2/add_voidptr_cast_intrinsic
diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h
index 6bbb426e5c..0f0a0e78a6 100644
--- a/cpp/dolfinx/fem/Form.h
+++ b/cpp/dolfinx/fem/Form.h
@@ -55,6 +55,8 @@ struct integral_data
   /// @param[in] entities Indices of entities to integrate over.
   /// @param[in] coeffs Indices of the coefficients that are present
   /// (active) in `kernel`.
+  /// @param[in] custom_data Optional custom user data pointer passed to
+  /// the kernel function.
   template <typename K, typename V, typename W>
     requires std::is_convertible_v<
                  std::remove_cvref_t<K>,
@@ -64,9 +66,10 @@ struct integral_data
                                            std::vector<std::int32_t>>
                  and std::is_convertible_v<std::remove_cvref_t<W>,
                                            std::vector<int>>
-  integral_data(K&& kernel, V&& entities, W&& coeffs)
+  integral_data(K&& kernel, V&& entities, W&& coeffs,
+                std::optional<void*> custom_data = std::nullopt)
       : kernel(std::forward<K>(kernel)), entities(std::forward<V>(entities)),
-        coeffs(std::forward<W>(coeffs))
+        coeffs(std::forward<W>(coeffs)), custom_data(custom_data)
   {
   }
 
@@ -82,6 +85,11 @@ struct integral_data
   /// @brief Indices of coefficients (from the form) that are in this
   /// integral.
   std::vector<int> coeffs;
+
+  /// @brief Custom user data pointer passed to the kernel function.
+  /// This can be used to pass runtime-computed data (e.g., per-cell
+  /// quadrature rules, material properties) to the kernel.
+  std::optional<void*> custom_data = std::nullopt;
 };
 
 /// @brief A representation of finite element variational forms.
@@ -391,6 +399,30 @@ class Form
     return it->second.kernel;
   }
 
+  /// @brief Get the custom data pointer for an integral.
+  ///
+  /// The custom data pointer is passed to the kernel function during
+  /// assembly. This can be used to pass runtime-computed data to
+  /// kernels (e.g., per-cell quadrature rules, material properties).
+  ///
+  /// @warning Void pointers are inherently unsafe and cannot be type or
+  /// bounds checked. Incorrect usage of this feature can and will lead
+  /// to undefined behaviour and crashes.
+  ///
+  /// @param[in] type Integral type.
+  /// @param[in] id Integral subdomain ID.
+  /// @param[in] kernel_idx Index of the kernel (we may have multiple
+  /// kernels for a given ID in mixed-topology meshes).
+  /// @return Custom data pointer for the integral, or std::nullopt if not set.
+  std::optional<void*> custom_data(IntegralType type, int id,
+                                   int kernel_idx) const
+  {
+    auto it = _integrals.find({type, id, kernel_idx});
+    if (it == _integrals.end())
+      throw std::runtime_error("Requested integral not found.");
+    return it->second.custom_data;
+  }
+
   /// @brief Get types of integrals in the form.
   /// @return Integrals types.
   std::set<IntegralType> integral_types() const
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index 3b6148f9e9..ceb505ddb2 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/cpp/dolfinx/fem/assemble_matrix_impl.h
@@ -18,6 +18,7 @@
 #include <dolfinx/mesh/Topology.h>
 #include <functional>
 #include <iterator>
+#include <optional>
 #include <span>
 #include <tuple>
 #include <vector>
@@ -60,6 +61,7 @@ using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
 /// function mesh.
 /// @param cell_info1 Cell permutation information for the trial
 /// function mesh.
+/// @param custom_data Custom user data pointer passed to the kernel.
 template <dolfinx::scalar T>
 void assemble_cells_matrix(
     la::MatSet<T> auto mat_set, mdspan2_t x_dofmap,
@@ -74,7 +76,8 @@ void assemble_cells_matrix(
     std::span<const std::int8_t> bc1, FEkernel<T> auto kernel,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
-    std::span<const std::uint32_t> cell_info1)
+    std::span<const std::uint32_t> cell_info1,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (cells.empty())
     return;
@@ -109,7 +112,7 @@ void assemble_cells_matrix(
     // Tabulate tensor
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(c, 0), constants.data(), cdofs.data(), nullptr,
-           nullptr, nullptr);
+           nullptr, custom_data.value_or(nullptr));
 
     // Compute A = P_0 \tilde{A} P_1^T (dof transformation)
     P0(Ae, cell_info0, cell0, ndim1);  // B = P0 \tilde{A}
@@ -198,6 +201,7 @@ void assemble_cells_matrix(
 /// function mesh.
 /// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
+/// @param custom_data Custom user data pointer passed to the kernel.
 template <dolfinx::scalar T>
 void assemble_entities(
     la::MatSet<T> auto mat_set, mdspan2_t x_dofmap,
@@ -221,7 +225,8 @@ void assemble_entities(
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (entities.empty())
     return;
@@ -259,7 +264,7 @@ void assemble_entities(
     // Tabulate tensor
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
-           &local_entity, &perm, nullptr);
+           &local_entity, &perm, custom_data.value_or(nullptr));
     P0(Ae, cell_info0, cell0, ndim1);
     P1T(Ae, cell_info1, cell1, ndim0);
 
@@ -338,6 +343,8 @@ void assemble_entities(
 /// function mesh.
 /// @param[in] perms Facet permutation integer. Empty if facet
 /// permutations are not required.
+/// @param[in] custom_data Optional pointer to user-supplied data passed
+/// to the kernel at runtime.
 template <dolfinx::scalar T>
 void assemble_interior_facets(
     la::MatSet<T> auto mat_set, mdspan2_t x_dofmap,
@@ -363,7 +370,8 @@ void assemble_interior_facets(
         coeffs,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (facets.empty())
     return;
@@ -440,7 +448,7 @@ void assemble_interior_facets(
                           : std::array{perms(cells[0], local_facet[0]),
                                        perms(cells[1], local_facet[1])};
     kernel(Ae.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
-           local_facet.data(), perm.data(), nullptr);
+           local_facet.data(), perm.data(), custom_data.value_or(nullptr));
 
     // Local element layout is a 2x2 block matrix with structure
     //
@@ -605,12 +613,14 @@ void assemble_matrix(
       std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
       std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx);
       auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
+      std::optional<void*> custom_data
+          = a.custom_data(IntegralType::cell, i, cell_type_idx);
       assert(cells.size() * cstride == coeffs.size());
       impl::assemble_cells_matrix(
           mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0}, P0,
           {dofs1, bs1, cells1}, P1T, bc0, bc1, fn,
           md::mdspan(coeffs.data(), cells.size(), cstride), constants,
-          cell_info0, cell_info1);
+          cell_info0, cell_info1, custom_data);
     }
 
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
@@ -646,6 +656,8 @@ void assemble_matrix(
       assert(fn);
       auto& [coeffs, cstride]
           = coefficients.at({IntegralType::interior_facet, i});
+      std::optional<void*> custom_data
+          = a.custom_data(IntegralType::interior_facet, i, 0);
 
       std::span facets = a.domain(IntegralType::interior_facet, i, 0);
       std::span facets0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0);
@@ -661,7 +673,7 @@ void assemble_matrix(
            mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
           P1T, bc0, bc1, fn,
           mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), constants,
-          cell_info0, cell_info1, facet_perms);
+          cell_info0, cell_info1, facet_perms, custom_data);
     }
 
     for (auto itg_type : {fem::IntegralType::exterior_facet,
@@ -688,6 +700,7 @@ void assemble_matrix(
         auto fn = a.kernel(itg_type, i, 0);
         assert(fn);
         auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+        std::optional<void*> custom_data = a.custom_data(itg_type, i, 0);
 
         std::span e = a.domain(itg_type, i, 0);
         mdspanx2_t entities(e.data(), e.size() / 2, 2);
@@ -700,7 +713,7 @@ void assemble_matrix(
             mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0,
             {dofs1, bs1, entities1}, P1T, bc0, bc1, fn,
             md::mdspan(coeffs.data(), entities.extent(0), cstride), constants,
-            cell_info0, cell_info1, perms);
+            cell_info0, cell_info1, perms, custom_data);
       }
     }
   }
diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h
index 37166b46ea..016c3a9bd4 100644
--- a/cpp/dolfinx/fem/assemble_scalar_impl.h
+++ b/cpp/dolfinx/fem/assemble_scalar_impl.h
@@ -17,6 +17,7 @@
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/Topology.h>
 #include <memory>
+#include <optional>
 #include <vector>
 
 namespace dolfinx::fem::impl
@@ -30,7 +31,8 @@ T assemble_cells(mdspan2_t x_dofmap,
                  std::span<const std::int32_t> cells, FEkernel<T> auto fn,
                  std::span<const T> constants,
                  md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
-                 std::span<scalar_value_t<T>> cdofs_b)
+                 std::span<scalar_value_t<T>> cdofs_b,
+                 std::optional<void*> custom_data = std::nullopt)
 {
   T value(0);
   if (cells.empty())
@@ -49,7 +51,7 @@ T assemble_cells(mdspan2_t x_dofmap,
       std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs_b.begin(), 3 * i));
 
     fn(&value, &coeffs(index, 0), constants.data(), cdofs_b.data(), nullptr,
-       nullptr, nullptr);
+       nullptr, custom_data.value_or(nullptr));
   }
 
   return value;
@@ -77,7 +79,8 @@ T assemble_entities(
     FEkernel<T> auto fn, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
-    std::span<scalar_value_t<T>> cdofs_b)
+    std::span<scalar_value_t<T>> cdofs_b,
+    std::optional<void*> custom_data = std::nullopt)
 {
   T value(0);
   if (entities.empty())
@@ -99,7 +102,7 @@ T assemble_entities(
     // Permutations
     std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
     fn(&value, &coeffs(f, 0), constants.data(), cdofs_b.data(), &local_entity,
-       &perm, nullptr);
+       &perm, custom_data.value_or(nullptr));
   }
 
   return value;
@@ -120,7 +123,8 @@ T assemble_interior_facets(
                                     md::dynamic_extent>>
         coeffs,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
-    std::span<scalar_value_t<T>> cdofs_b)
+    std::span<scalar_value_t<T>> cdofs_b,
+    std::optional<void*> custom_data = std::nullopt)
 {
   T value(0);
   if (facets.empty())
@@ -150,7 +154,7 @@ T assemble_interior_facets(
                           : std::array{perms(cells[0], local_facet[0]),
                                        perms(cells[1], local_facet[1])};
     fn(&value, &coeffs(f, 0, 0), constants.data(), cdofs_b.data(),
-       local_facet.data(), perm.data(), nullptr);
+       local_facet.data(), perm.data(), custom_data.value_or(nullptr));
   }
 
   return value;
@@ -178,11 +182,12 @@ T assemble_scalar(
     auto fn = M.kernel(IntegralType::cell, i, 0);
     assert(fn);
     auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
+    std::optional<void*> custom_data = M.custom_data(IntegralType::cell, i, 0);
     std::span<const std::int32_t> cells = M.domain(IntegralType::cell, i, 0);
     assert(cells.size() * cstride == coeffs.size());
     value += impl::assemble_cells(
         x_dofmap, x, cells, fn, constants,
-        md::mdspan(coeffs.data(), cells.size(), cstride), cdofs_b);
+        md::mdspan(coeffs.data(), cells.size(), cstride), cdofs_b, custom_data);
   }
 
   mesh::CellType cell_type = mesh->topology()->cell_type();
@@ -204,6 +209,8 @@ T assemble_scalar(
     assert(fn);
     auto& [coeffs, cstride]
         = coefficients.at({IntegralType::interior_facet, i});
+    std::optional<void*> custom_data
+        = M.custom_data(IntegralType::interior_facet, i, 0);
     std::span facets = M.domain(IntegralType::interior_facet, i, 0);
 
     constexpr std::size_t num_adjacent_cells = 2;
@@ -220,7 +227,7 @@ T assemble_scalar(
         md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                         md::dynamic_extent>>(
             coeffs.data(), facets.size() / shape1, 2, cstride),
-        facet_perms, cdofs_b);
+        facet_perms, cdofs_b, custom_data);
   }
 
   for (auto itg_type : {fem::IntegralType::exterior_facet,
@@ -236,6 +243,7 @@ T assemble_scalar(
       auto fn = M.kernel(itg_type, i, 0);
       assert(fn);
       auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+      std::optional<void*> custom_data = M.custom_data(itg_type, i, 0);
 
       std::span entities = M.domain(itg_type, i, 0);
 
@@ -248,7 +256,7 @@ T assemble_scalar(
               entities.data(), entities.size() / 2, 2),
           fn, constants,
           md::mdspan(coeffs.data(), entities.size() / 2, cstride), perms,
-          cdofs_b);
+          cdofs_b, custom_data);
     }
   }
 
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index f57e646dd9..39ab35adcf 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -76,6 +76,8 @@ using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
 /// conditions applied.
 /// @param[in] x0 Vector used in the lifting.
 /// @param[in] alpha Scaling to apply.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to the
+/// kernel at runtime.
 template <int _bs0 = -1, int _bs1 = -1, typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -92,7 +94,8 @@ void _lift_bc_cells(
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (cells.empty())
     return;
@@ -164,7 +167,7 @@ void _lift_bc_cells(
 
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
-           nullptr, nullptr, nullptr);
+           nullptr, nullptr, custom_data.value_or(nullptr));
     P0(Ae, cell_info0, c0, num_cols);
     P1T(Ae, cell_info1, c1, num_rows);
 
@@ -260,6 +263,8 @@ void _lift_bc_cells(
 /// local_facet_idx)` is the permutation value for the facet attached to
 /// the cell `cell_idx` with local index `local_facet_idx` relative to
 /// the cell. Empty if facet permutations are not required.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to the
+/// kernel at runtime.
 template <typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -286,7 +291,8 @@ void _lift_bc_entities(
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
     std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (entities.empty())
     return;
@@ -345,7 +351,7 @@ void _lift_bc_entities(
     std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
-           &local_entity, &perm, nullptr);
+           &local_entity, &perm, custom_data.value_or(nullptr));
     P0(Ae, cell_info0, cell0, num_cols);
     P1T(Ae, cell_info1, cell1, num_rows);
 
@@ -415,6 +421,8 @@ void _lift_bc_entities(
 /// local_facet_idx)` is the permutation value for the facet attached to
 /// the cell `cell_idx` with local index `local_facet_idx` relative to
 /// the cell. Empty if facet permutations are not required.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to the
+/// kernel at runtime.
 template <typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -443,7 +451,8 @@ void _lift_bc_interior_facets(
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
     std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (facets.empty())
     return;
@@ -559,7 +568,7 @@ void _lift_bc_interior_facets(
                           : std::array{perms(cells[0], local_facet[0]),
                                        perms(cells[1], local_facet[1])};
     kernel(Ae.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
-           local_facet.data(), perm.data(), nullptr);
+           local_facet.data(), perm.data(), custom_data.value_or(nullptr));
 
     if (cells0[0] >= 0)
       P0(Ae, cell_info0, cells0[0], num_cols);
@@ -658,6 +667,8 @@ void _lift_bc_interior_facets(
 /// coefficient for cell `i`.
 /// @param[in] cell_info0 Cell permutation information for the test
 /// function mesh.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to
+/// the kernel at runtime.
 template <int _bs = -1, typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -671,7 +682,8 @@ void assemble_cells(
     std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
     FEkernel<T> auto kernel, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
-    std::span<const std::uint32_t> cell_info0)
+    std::span<const std::uint32_t> cell_info0,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (cells.empty())
     return;
@@ -697,8 +709,8 @@ void assemble_cells(
 
     // Tabulate vector for cell
     std::ranges::fill(be, 0);
-    kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
-           nullptr, nullptr, nullptr);
+    kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(), &c,
+           nullptr, custom_data.value_or(nullptr));
     P0(be, cell_info0, c0, 1);
 
     // Scatter cell vector to 'global' vector array
@@ -751,6 +763,8 @@ void assemble_cells(
 /// function mesh.
 /// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to
+/// the kernel at runtime.
 template <int _bs = -1, typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -769,7 +783,8 @@ void assemble_entities(
     FEkernel<T> auto kernel, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const std::uint32_t> cell_info0,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   if (entities.empty())
     return;
@@ -801,7 +816,7 @@ void assemble_entities(
     // Tabulate element vector
     std::ranges::fill(be, 0);
     kernel(be.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
-           &local_entity, &perm, nullptr);
+           &local_entity, &perm, custom_data.value_or(nullptr));
     P0(be, cell_info0, cell0, 1);
 
     // Add element vector to global vector
@@ -846,6 +861,8 @@ void assemble_entities(
 /// function mesh.
 /// @param[in] perms Facet permutation integer. Empty if facet
 /// permutations are not required.
+/// @param[in] custom_data Optional pointer to user-supplied data passed to
+/// the kernel at runtime.
 template <int _bs = -1, typename V,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
@@ -866,7 +883,8 @@ void assemble_interior_facets(
                                     md::dynamic_extent>>
         coeffs,
     std::span<const std::uint32_t> cell_info0,
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
+    std::optional<void*> custom_data = std::nullopt)
 {
   using X = scalar_value_t<T>;
 
@@ -918,7 +936,7 @@ void assemble_interior_facets(
                           : std::array{perms(cells[0], local_facet[0]),
                                        perms(cells[1], local_facet[1])};
     kernel(be.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
-           local_facet.data(), perm.data(), nullptr);
+           local_facet.data(), perm.data(), custom_data.value_or(nullptr));
 
     if (cells0[0] >= 0)
       P0(be, cell_info0, cells0[0], 1);
@@ -1026,6 +1044,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
     auto kernel = a.kernel(IntegralType::cell, i, 0);
     assert(kernel);
     auto& [_coeffs, cstride] = coefficients.at({IntegralType::cell, i});
+    std::optional<void*> custom_data = a.custom_data(IntegralType::cell, i, 0);
     std::span cells = a.domain(IntegralType::cell, i, 0);
     std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, 0);
     std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, 0);
@@ -1036,20 +1055,21 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
       _lift_bc_cells<1, 1>(b, x_dofmap, x, kernel, cells,
                            {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
                            P1T, constants, coeffs, cell_info0, cell_info1,
-                           bc_values1, bc_markers1, x0, alpha);
+                           bc_values1, bc_markers1, x0, alpha, custom_data);
     }
     else if (bs0 == 3 and bs1 == 3)
     {
       _lift_bc_cells<3, 3>(b, x_dofmap, x, kernel, cells,
                            {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
                            P1T, constants, coeffs, cell_info0, cell_info1,
-                           bc_values1, bc_markers1, x0, alpha);
+                           bc_values1, bc_markers1, x0, alpha, custom_data);
     }
     else
     {
       _lift_bc_cells(b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
                      {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
-                     cell_info1, bc_values1, bc_markers1, x0, alpha);
+                     cell_info1, bc_values1, bc_markers1, x0, alpha,
+                     custom_data);
     }
   }
 
@@ -1072,6 +1092,8 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
     assert(kernel);
     auto& [coeffs, cstride]
         = coefficients.at({IntegralType::interior_facet, i});
+    std::optional<void*> custom_data
+        = a.custom_data(IntegralType::interior_facet, i, 0);
 
     using mdspanx22_t
         = md::mdspan<const std::int32_t,
@@ -1089,7 +1111,8 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
         b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
         {dofmap1, bs1, facets1}, P1T, constants,
         mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0,
-        cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms);
+        cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms,
+        custom_data);
   }
 
   for (auto itg_type : {fem::IntegralType::exterior_facet,
@@ -1105,6 +1128,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
       auto kernel = a.kernel(itg_type, i, 0);
       assert(kernel);
       auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+      std::optional<void*> custom_data = a.custom_data(itg_type, i, 0);
 
       using mdspanx2_t
           = md::mdspan<const std::int32_t,
@@ -1120,7 +1144,7 @@ void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
           b, x_dofmap, x, kernel, entities, {dofmap0, bs0, entities0}, P0,
           {dofmap1, bs1, entities1}, P1T, constants,
           md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0,
-          cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
+          cell_info1, bc_values1, bc_markers1, x0, alpha, perms, custom_data);
     }
   }
 }
@@ -1276,24 +1300,29 @@ void assemble_vector(
       std::span cells = L.domain(IntegralType::cell, i, cell_type_idx);
       std::span cells0 = L.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
       auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
+      void* custom_data = L.custom_data(IntegralType::cell, i, cell_type_idx)
+                              .value_or(nullptr);
       assert(cells.size() * cstride == coeffs.size());
       if (bs == 1)
       {
         impl::assemble_cells<1>(
             P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
-            md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
+            md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0,
+            custom_data);
       }
       else if (bs == 3)
       {
         impl::assemble_cells<3>(
             P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
-            md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
+            md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0,
+            custom_data);
       }
       else
       {
-        impl::assemble_cells(
-            P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
-            md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
+        impl::assemble_cells(P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn,
+                             constants,
+                             md::mdspan(coeffs.data(), cells.size(), cstride),
+                             cell_info0, custom_data);
       }
     }
 
@@ -1327,6 +1356,8 @@ void assemble_vector(
       assert(fn);
       auto& [coeffs, cstride]
           = coefficients.at({IntegralType::interior_facet, i});
+      void* custom_data
+          = L.custom_data(IntegralType::interior_facet, i, 0).value_or(nullptr);
       std::span facets = L.domain(IntegralType::interior_facet, i, 0);
       std::span facets1 = L.domain_arg(IntegralType::interior_facet, 0, i, 0);
       assert((facets.size() / 4) * 2 * cstride == coeffs.size());
@@ -1339,7 +1370,7 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, facet_perms);
+            cell_info0, facet_perms, custom_data);
       }
       else if (bs == 3)
       {
@@ -1350,7 +1381,7 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, facet_perms);
+            cell_info0, facet_perms, custom_data);
       }
       else
       {
@@ -1361,7 +1392,7 @@ void assemble_vector(
              mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
             fn, constants,
             mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
-            cell_info0, facet_perms);
+            cell_info0, facet_perms, custom_data);
       }
     }
 
@@ -1378,6 +1409,7 @@ void assemble_vector(
         auto fn = L.kernel(itg_type, i, 0);
         assert(fn);
         auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+        void* custom_data = L.custom_data(itg_type, i, 0).value_or(nullptr);
         std::span e = L.domain(itg_type, i, 0);
         mdspanx2_t entities(e.data(), e.size() / 2, 2);
         std::span e1 = L.domain_arg(itg_type, 0, i, 0);
@@ -1388,7 +1420,7 @@ void assemble_vector(
           impl::assemble_entities<1>(
               P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
               constants, md::mdspan(coeffs.data(), entities.extent(0), cstride),
-              cell_info0, perms);
+              cell_info0, perms, custom_data);
         }
         else if (bs == 3)
         {
@@ -1396,7 +1428,7 @@ void assemble_vector(
               P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
               constants,
               md::mdspan(coeffs.data(), entities.size() / 2, cstride),
-              cell_info0, perms);
+              cell_info0, perms, custom_data);
         }
         else
         {
@@ -1404,7 +1436,7 @@ void assemble_vector(
               P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn,
               constants,
               md::mdspan(coeffs.data(), entities.size() / 2, cstride),
-              cell_info0, perms);
+              cell_info0, perms, custom_data);
         }
       }
     }
diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py
index e6b0f6b5b2..ad0b868ca2 100644
--- a/python/demo/demo_static-condensation.py
+++ b/python/demo/demo_static-condensation.py
@@ -242,7 +242,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, cu
 
 formtype = form_cpp_class(dtype)  # type: ignore
 cells = np.arange(msh.topology.index_map(msh.topology.dim).size_local)
-integrals = {IntegralType.cell: [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]}
+integrals = {IntegralType.cell: [(0, tabulate_A.address, cells, np.array([], dtype=np.int8), None)]}
 a_cond = Form(
     formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, [], mesh=msh._cpp_object)
 )
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 9e4b365722..06cbe5ac36 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -667,8 +667,8 @@ void declare_form(nb::module_& m, std::string type)
                  std::vector<std::tuple<
                      int, std::uintptr_t,
                      nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>,
-                     nb::ndarray<const std::int8_t, nb::ndim<1>,
-                                 nb::c_contig>>>>& integrals,
+                     nb::ndarray<const std::int8_t, nb::ndim<1>, nb::c_contig>,
+                     std::optional<std::uintptr_t>>>>& integrals,
              const std::vector<std::shared_ptr<
                  const dolfinx::fem::Function<T, U>>>& coefficients,
              const std::vector<
@@ -684,17 +684,22 @@ void declare_form(nb::module_& m, std::string type)
             // Loop over kernel for each entity type
             for (auto& [type, kernels] : integrals)
             {
-              for (auto& [id, ptr, e, c] : kernels)
+              for (auto& [id, ptr, e, c, custom_data] : kernels)
               {
                 auto kn_ptr
                     = (void (*)(T*, const T*, const T*,
                                 const typename geom_type<T>::value_type*,
                                 const int*, const std::uint8_t*, void*))ptr;
+                // Convert custom_data from std::optional<uintptr_t> to
+                // std::optional<void*>
+                std::optional<void*> cd = std::nullopt;
+                if (custom_data.has_value())
+                  cd = reinterpret_cast<void*>(custom_data.value());
                 _integrals.insert(
                     {{type, id, 0},
                      {kn_ptr,
                       std::vector<std::int32_t>(e.data(), e.data() + e.size()),
-                      std::vector<int>(c.data(), c.data() + c.size())}});
+                      std::vector<int>(c.data(), c.data() + c.size()), cd}});
               }
             }
 
@@ -761,6 +766,34 @@ void declare_form(nb::module_& m, std::string type)
       .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)
+      .def(
+          "custom_data",
+          [](const dolfinx::fem::Form<T, U>& self,
+             dolfinx::fem::IntegralType type, int id,
+             int kernel_idx) -> std::optional<std::uintptr_t>
+          {
+            /* We use std::uintptr_t to map void* into Python to allow
+            seamless interop with numpy.ctypes.data.
+            Example:
+            # Create data to pass to kernel:
+            data = np.array([2.0], dtype=np.float64)
+            # Pass pointer at Form creation via integrals tuple:
+            integrals = {IntegralType.cell: [(id, kernel_ptr, cells,
+            coeffs, data.ctypes.data)]}
+            formtype = form_cpp_class(dtype)
+            L = Form(formtype([V._cpp_object], integrals, [], [], False, [],
+            mesh=mesh._cpp_object)) Important: You must keep the NumPy array
+            alive in Python while the form is being used; if the array is
+            garbage collected, the pointer becomes invalid! */
+            auto cd = self.custom_data(type, id, kernel_idx);
+            return cd ? std::optional<std::uintptr_t>((std::uintptr_t)*cd)
+                      : std::nullopt;
+          },
+          nb::arg("type"), nb::arg("id"), nb::arg("kernel_idx"),
+          "Get custom data pointer for an integral, or None if not set. "
+          "Warning: void-pointers are inherently unsafe and cannot be type "
+          "or bounds checked. Incorrect usage of this feature can and will "
+          "lead to undefined behaviour and crashes.")
       .def(
           "domains",
           [](const dolfinx::fem::Form<T, U>& self,
diff --git a/python/test/unit/fem/test_custom_data.py b/python/test/unit/fem/test_custom_data.py
new file mode 100644
index 0000000000..4df6675592
--- /dev/null
+++ b/python/test/unit/fem/test_custom_data.py
@@ -0,0 +1,415 @@
+"""Unit tests for custom_data functionality in assembly."""
+
+# Copyright (C) 2025 Susanne Claus
+#
+# This file is part of DOLFINx (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+from mpi4py import MPI
+
+import numpy as np
+import pytest
+
+import dolfinx
+from dolfinx import la
+from dolfinx.fem import Form, IntegralType, form_cpp_class, functionspace
+from dolfinx.mesh import create_unit_square
+from ffcx.codegeneration.utils import (
+    numba_ufcx_kernel_signature,
+    voidptr_to_float64_ptr,
+)
+
+numba = pytest.importorskip("numba")
+ufcx_signature = numba_ufcx_kernel_signature
+
+
+def tabulate_rank1_with_custom_data(dtype, xdtype):
+    """Kernel that reads a scaling factor from custom_data.
+
+    Note: custom_data must be set to a valid pointer before assembly.
+    """
+
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate(b_, w_, c_, coords_, local_index, orientation, custom_data):
+        b = numba.carray(b_, (3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        # Cast void* to float64* and read the scale value
+        typed_ptr = voidptr_to_float64_ptr(custom_data)
+        scale = typed_ptr[0]
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        # 2x Element area Ae
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        b[:] = scale * Ae / 6.0
+
+    return tabulate
+
+
+def tabulate_rank2_with_custom_data(dtype, xdtype):
+    """Kernel that reads a scaling factor from custom_data for matrix assembly.
+
+    Note: custom_data must be set to a valid pointer before assembly.
+    """
+
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate(A_, w_, c_, coords_, entity_local_index, cell_orientation, custom_data):
+        A = numba.carray(A_, (3, 3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        # Cast void* to float64* and read the scale value
+        typed_ptr = voidptr_to_float64_ptr(custom_data)
+        scale = typed_ptr[0]
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        # 2x Element area Ae
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        B = np.array([y1 - y2, y2 - y0, y0 - y1, x2 - x1, x0 - x2, x1 - x0], dtype=dtype).reshape(
+            2, 3
+        )
+        A[:, :] = scale * np.dot(B.T, B) / (2 * Ae)
+
+    return tabulate
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_vector_assembly(dtype):
+    """Test that custom_data is correctly passed to kernels during vector assembly."""
+    xdtype = np.real(dtype(0)).dtype
+    k1 = tabulate_rank1_with_custom_data(dtype, xdtype)
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 13, 13, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    tdim = mesh.topology.dim
+    num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
+    cells = np.arange(num_cells, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Create custom_data with scale=1.0
+    scale_value = np.array([1.0], dtype=dtype)
+    scale_ptr = scale_value.ctypes.data
+
+    # Pass custom_data at form creation time via the 5th element of the integrals tuple
+    integrals = {IntegralType.cell: [(0, k1.address, cells, active_coeffs, scale_ptr)]}
+    formtype = form_cpp_class(dtype)
+    L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
+
+    # Assemble with scale=1.0
+    b1 = dolfinx.fem.assemble_vector(L)
+    b1.scatter_reverse(la.InsertMode.add)
+    norm1 = la.norm(b1)
+
+    # Verify we can read back the custom_data pointer
+    assert L._cpp_object.custom_data(IntegralType.cell, 0, 0) == scale_ptr
+
+    # Update custom_data to scale=2.0 (by modifying the underlying array)
+    scale_value[0] = 2.0
+    b2 = dolfinx.fem.assemble_vector(L)
+    b2.scatter_reverse(la.InsertMode.add)
+    norm2 = la.norm(b2)
+
+    # The norm with scale=2 should be 2x the norm with scale=1
+    assert np.isclose(norm2, 2.0 * norm1)
+
+    # Test with scale=3.0
+    scale_value[0] = 3.0
+    b3 = dolfinx.fem.assemble_vector(L)
+    b3.scatter_reverse(la.InsertMode.add)
+    norm3 = la.norm(b3)
+
+    assert np.isclose(norm3, 3.0 * norm1)
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_matrix_assembly(dtype):
+    """Test that custom_data is correctly passed to kernels during matrix assembly."""
+    xdtype = np.real(dtype(0)).dtype
+    k2 = tabulate_rank2_with_custom_data(dtype, xdtype)
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 13, 13, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Create custom_data with scale=1.0
+    scale_value = np.array([1.0], dtype=dtype)
+
+    # Pass custom_data at form creation time via the 5th element
+    integrals = {
+        IntegralType.cell: [(0, k2.address, cells, active_coeffs, scale_value.ctypes.data)]
+    }
+    formtype = form_cpp_class(dtype)
+    a = Form(
+        formtype(
+            [V._cpp_object, V._cpp_object],
+            integrals,
+            [],
+            [],
+            False,
+            [],
+            mesh=mesh._cpp_object,
+        )
+    )
+
+    # Assemble with scale=1.0
+    A1 = dolfinx.fem.assemble_matrix(a)
+    A1.scatter_reverse()
+    norm1 = np.sqrt(A1.squared_norm())
+
+    # Update custom_data to scale=2.0 (by modifying the underlying array)
+    scale_value[0] = 2.0
+    A2 = dolfinx.fem.assemble_matrix(a)
+    A2.scatter_reverse()
+    norm2 = np.sqrt(A2.squared_norm())
+
+    # The norm with scale=2 should be 2x the norm with scale=1
+    assert np.isclose(norm2, 2.0 * norm1)
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_default_nullptr(dtype):
+    """Test that custom_data defaults to nullptr (None) when not provided."""
+    xdtype = np.real(dtype(0)).dtype
+
+    # Define a simple kernel that doesn't use custom_data
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate_simple(b_, w_, c_, coords_, local_index, orientation, custom_data):
+        b = numba.carray(b_, (3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        b[:] = Ae / 6.0
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 5, 5, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    tdim = mesh.topology.dim
+    num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
+    cells = np.arange(num_cells, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Pass None as custom_data (5th element) to get std::nullopt
+    integrals = {IntegralType.cell: [(0, tabulate_simple.address, cells, active_coeffs, None)]}
+    formtype = form_cpp_class(dtype)
+    L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
+
+    # custom_data should be None (std::nullopt) when passed as None
+    assert L._cpp_object.custom_data(IntegralType.cell, 0, 0) is None
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_struct(dtype):
+    """Test passing a struct with multiple values through custom_data."""
+    xdtype = np.real(dtype(0)).dtype
+
+    # Define a kernel that reads two values from custom_data
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate_with_struct(b_, w_, c_, coords_, local_index, orientation, custom_data):
+        b = numba.carray(b_, (3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        # Cast void* to float64* and read two values: [scale, offset]
+        typed_ptr = voidptr_to_float64_ptr(custom_data)
+        scale = typed_ptr[0]
+        offset = typed_ptr[1]
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        b[:] = scale * Ae / 6.0 + offset
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 5, 5, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    tdim = mesh.topology.dim
+    # Only iterate over local cells (not ghosts) to avoid double-counting
+    # contributions after scatter_reverse
+    num_local_cells = mesh.topology.index_map(tdim).size_local
+    cells = np.arange(num_local_cells, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Test 1: scale=1.0, offset=0.0 (baseline)
+    struct_data = np.array([1.0, 0.0], dtype=dtype)
+
+    # Pass custom_data at form creation time
+    integrals = {
+        IntegralType.cell: [
+            (0, tabulate_with_struct.address, cells, active_coeffs, struct_data.ctypes.data)
+        ]
+    }
+    formtype = form_cpp_class(dtype)
+    L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
+
+    b_baseline = dolfinx.fem.assemble_vector(L)
+    b_baseline.scatter_reverse(la.InsertMode.add)
+    norm_baseline = la.norm(b_baseline)
+
+    # Test 2: scale=2.0, offset=0.0 - should double the norm
+    struct_data[0] = 2.0
+    struct_data[1] = 0.0
+    b_scaled = dolfinx.fem.assemble_vector(L)
+    b_scaled.scatter_reverse(la.InsertMode.add)
+    norm_scaled = la.norm(b_scaled)
+    assert np.isclose(norm_scaled, 2.0 * norm_baseline)
+
+    # Test 3: scale=0.0, offset=1.0 - pure offset contribution
+    struct_data[0] = 0.0
+    struct_data[1] = 1.0
+    b_offset = dolfinx.fem.assemble_vector(L)
+    b_offset.scatter_reverse(la.InsertMode.add)
+    # With offset=1.0, each DOF gets contribution from each cell it touches
+    # Interior nodes touch 6 cells, edge nodes touch 3-4, corner nodes touch 1-2
+    # The sum of all contributions equals 3 * num_local_cells (3 DOFs per cell, offset=1.0 each)
+    # Sum only local DOFs and gather across processes
+    local_sum = np.sum(b_offset.array[: V.dofmap.index_map.size_local * V.dofmap.index_map_bs])
+    total_contribution = mesh.comm.allreduce(local_sum, op=MPI.SUM)
+    total_cells = mesh.comm.allreduce(num_local_cells, op=MPI.SUM)
+    assert np.isclose(total_contribution, 3.0 * total_cells)
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_multiple_parameters(dtype):
+    """Test custom_data with multiple parameters demonstrating complex data passing.
+
+    This test shows how to pass multiple values through custom_data, simulating
+    the use case of passing runtime-computed parameters like material properties
+    or integration parameters.
+    """
+    xdtype = np.real(dtype(0)).dtype
+
+    # Kernel that uses three parameters: coefficient, exponent base, and additive term
+    # Computes: coeff * base^area + additive
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate_with_params(b_, w_, c_, coords_, local_index, orientation, custom_data):
+        b = numba.carray(b_, (3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        # Read three parameters from custom_data
+        typed_ptr = voidptr_to_float64_ptr(custom_data)
+        coeff = typed_ptr[0]
+        power = typed_ptr[1]
+        additive = typed_ptr[2]
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        # Use power as a simple multiplier (avoiding actual power function complexity)
+        val = coeff * power * Ae / 6.0 + additive
+        b[:] = val
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 4, 4, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    tdim = mesh.topology.dim
+    num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
+    cells = np.arange(num_cells, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Test with specific parameters: coeff=2, power=3, additive=0
+    params = np.array([2.0, 3.0, 0.0], dtype=dtype)
+
+    # Pass custom_data at form creation time
+    integrals = {
+        IntegralType.cell: [
+            (0, tabulate_with_params.address, cells, active_coeffs, params.ctypes.data)
+        ]
+    }
+    formtype = form_cpp_class(dtype)
+    L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
+
+    b1 = dolfinx.fem.assemble_vector(L)
+    b1.scatter_reverse(la.InsertMode.add)
+    norm1 = la.norm(b1)
+
+    # Change parameters: coeff=1, power=6, additive=0 (should give same result: 1*6 = 2*3)
+    params[0] = 1.0
+    params[1] = 6.0
+    b2 = dolfinx.fem.assemble_vector(L)
+    b2.scatter_reverse(la.InsertMode.add)
+    norm2 = la.norm(b2)
+
+    assert np.isclose(norm1, norm2)
+
+
+@pytest.mark.parametrize("dtype", [np.float64])
+def test_custom_data_global_parameter_update(dtype):
+    """Test updating custom_data between assemblies for parameter studies.
+
+    This demonstrates a key use case: running multiple assemblies with
+    different parameter values without recreating the form. The custom_data
+    points to a parameter that can be modified between assembly calls.
+    """
+    xdtype = np.real(dtype(0)).dtype
+
+    # Kernel that reads a diffusion coefficient from custom_data
+    @numba.cfunc(ufcx_signature(dtype, xdtype), nopython=True)
+    def tabulate_diffusion(b_, w_, c_, coords_, local_index, orientation, custom_data):
+        b = numba.carray(b_, (3), dtype=dtype)
+        coordinate_dofs = numba.carray(coords_, (3, 3), dtype=xdtype)
+
+        # Read diffusion coefficient from custom_data
+        typed_ptr = voidptr_to_float64_ptr(custom_data)
+        kappa = typed_ptr[0]  # Diffusion coefficient
+
+        x0, y0 = coordinate_dofs[0, :2]
+        x1, y1 = coordinate_dofs[1, :2]
+        x2, y2 = coordinate_dofs[2, :2]
+
+        # 2x Element area Ae
+        Ae = abs((x0 - x1) * (y2 - y1) - (y0 - y1) * (x2 - x1))
+        # Simple load vector scaled by diffusion coefficient
+        b[:] = kappa * Ae / 6.0
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 8, 8, dtype=xdtype)
+    V = functionspace(mesh, ("Lagrange", 1))
+
+    tdim = mesh.topology.dim
+    num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
+    cells = np.arange(num_cells, dtype=np.int32)
+    active_coeffs = np.array([], dtype=np.int8)
+
+    # Parameter array - this can be updated between assemblies
+    kappa = np.array([1.0], dtype=dtype)
+
+    integrals = {
+        IntegralType.cell: [
+            (0, tabulate_diffusion.address, cells, active_coeffs, kappa.ctypes.data)
+        ]
+    }
+    formtype = form_cpp_class(dtype)
+    L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
+
+    # Store results for different kappa values
+    results = []
+    kappa_values = [0.1, 1.0, 10.0, 100.0]
+
+    for k in kappa_values:
+        kappa[0] = k
+        b = dolfinx.fem.assemble_vector(L)
+        b.scatter_reverse(la.InsertMode.add)
+        results.append(la.norm(b))
+
+    # Verify linear scaling: norm should scale linearly with kappa
+    # norm(kappa=k) / norm(kappa=1) should equal k
+    base_norm = results[1]  # kappa=1.0
+    for i, k in enumerate(kappa_values):
+        expected_ratio = k / 1.0
+        actual_ratio = results[i] / base_norm
+        assert np.isclose(actual_ratio, expected_ratio, rtol=1e-10)
diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py
index 242e867469..3f798f66c1 100644
--- a/python/test/unit/fem/test_custom_jit_kernels.py
+++ b/python/test/unit/fem/test_custom_jit_kernels.py
@@ -93,9 +93,9 @@ def test_numba_assembly(dtype):
     active_coeffs = np.array([], dtype=np.int8)
     integrals = {
         IntegralType.cell: [
-            (0, k2.address, cells, active_coeffs),
-            (1, k2.address, np.arange(0), active_coeffs),
-            (2, k2.address, np.arange(0), active_coeffs),
+            (0, k2.address, cells, active_coeffs, None),
+            (1, k2.address, np.arange(0), active_coeffs, None),
+            (2, k2.address, np.arange(0), active_coeffs, None),
         ]
     }
     formtype = form_cpp_class(dtype)
@@ -104,7 +104,7 @@ def test_numba_assembly(dtype):
             [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object
         )
     )
-    integrals = {IntegralType.cell: [(0, k1.address, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, k1.address, cells, active_coeffs, None)]}
     L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object))
 
     A = dolfinx.fem.assemble_matrix(a)
@@ -135,7 +135,9 @@ 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: [(0, 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, None)
+        ]
     }
     formtype = form_cpp_class(dtype)
     L = Form(
@@ -275,7 +277,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: [(0, ptrA, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, ptrA, cells, active_coeffs, None)]}
     a = Form(
         _cpp.fem.Form_float64(
             [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object
@@ -283,7 +285,7 @@ def test_cffi_assembly():
     )
 
     ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL"))
-    integrals = {IntegralType.cell: [(0, ptrL, cells, active_coeffs)]}
+    integrals = {IntegralType.cell: [(0, ptrL, cells, active_coeffs, None)]}
     L = Form(
         _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object)
     )