assemble: relax multi-domain restriction

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,

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)