From ff69532e13db3c47ff8d766bb4920a5c5f395695 Mon Sep 17 00:00:00 2001
From: ksagiyam <k.sagiyama@imperial.ac.uk>
Date: Sat, 18 Oct 2025 00:14:14 +0100
Subject: [PATCH] assemble: relax multi-domain restriction

---
 firedrake/assemble.py                         | 11 -------
 tests/firedrake/submesh/test_submesh_solve.py | 31 +++++++++++++++++++
 2 files changed, 31 insertions(+), 11 deletions(-)

diff --git a/firedrake/assemble.py b/firedrake/assemble.py
index 5db7c007ce..e151df1063 100644
--- a/firedrake/assemble.py
+++ b/firedrake/assemble.py
@@ -1090,17 +1090,6 @@ def local_kernels(self):
             each possible combination.
 
         """
-        try:
-            topology, = set(d.topology.submesh_ancesters[-1] for d in self._form.ufl_domains())
-        except ValueError:
-            raise NotImplementedError("All integration domains must share a mesh topology")
-
-        for o in itertools.chain(self._form.arguments(), self._form.coefficients()):
-            domains = extract_domains(o)
-            for domain in domains:
-                if domain is not None and domain.topology.submesh_ancesters[-1] != topology:
-                    raise NotImplementedError("Assembly with multiple meshes is not supported")
-
         if isinstance(self._form, ufl.Form):
             kernels = tsfc_interface.compile_form(
                 self._form, "form", diagonal=self.diagonal,
diff --git a/tests/firedrake/submesh/test_submesh_solve.py b/tests/firedrake/submesh/test_submesh_solve.py
index aa5e6725de..0d715d24c5 100644
--- a/tests/firedrake/submesh/test_submesh_solve.py
+++ b/tests/firedrake/submesh/test_submesh_solve.py
@@ -10,6 +10,37 @@
 cwd = abspath(dirname(__file__))
 
 
+# Solve an uncoupled system of two equations in one shot.
+def test_submesh_solve_uncoupled():
+    mesh0 = UnitTriangleMesh()
+    mesh1 = UnitSquareMesh(1, 1, quadrilateral=True)
+    V0 = FunctionSpace(mesh0, "P", 1)
+    V1 = FunctionSpace(mesh1, "Q", 1)
+    V = V0 * V1
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    u0, u1 = split(u)
+    v0, v1 = split(v)
+    dx0 = Measure(
+        "dx", domain=mesh0,
+    )
+    dx1 = Measure(
+        "dx", domain=mesh1,
+    )
+    a = (
+        inner(u0, v0) * dx0
+        + inner(u1, v1) * dx1
+    )
+    L = (
+        inner(Constant(3.), v0) * dx0
+        + inner(Constant(4.), v1) * dx1
+    )
+    solution = Function(V)
+    solve(a == L, solution)
+    assert np.allclose(solution.subfunctions[0].dat.data_ro, 3.)
+    assert np.allclose(solution.subfunctions[1].dat.data_ro, 4.)
+
+
 def _solve_helmholtz(mesh):
     V = FunctionSpace(mesh, "CG", 1)
     u = TrialFunction(V)