conflict

Author: pbrubeck

firedrake/assemble.py

@@ -1107,17 +1107,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/regression/test_multiple_domains.py

@@ -58,11 +58,6 @@ def test_mismatching_meshes_real_space(mesh1, mesh3):
         project(donor, target)
 
 
-def test_mismatching_topologies(mesh1, mesh3):
-    with pytest.raises(NotImplementedError):
-        assemble(1*dx(domain=mesh1) + 2*dx(domain=mesh3))
-
-
 def test_functional(mesh1, mesh2):
     c = Constant(1)
 
@@ -123,3 +118,90 @@ def test_two_form(mesh1, mesh2, form, expect):
     val = assemble(form).M.values
 
     assert np.allclose(val, expect)
+
+
+def test_multi_domain_solve():
+    mesh1 = UnitSquareMesh(7, 7, quadrilateral=True)
+    x1, y1 = SpatialCoordinate(mesh1)
+    mesh2 = UnitSquareMesh(8, 8)
+    x2, y2 = SpatialCoordinate(mesh2)
+    V1 = FunctionSpace(mesh1, "Q", 3)
+    V2 = FunctionSpace(mesh2, "CG", 2)
+    V = V1 * V2
+
+    u1, u2 = TrialFunctions(V)
+    v1, v2 = TestFunctions(V)
+
+    a = (
+        inner(grad(u1), grad(v1))*dx(domain=mesh1)
+        + inner(grad(u2), grad(v2))*dx(domain=mesh2)
+    )
+
+    u_exact_expr1 = sin(pi * x1) * sin(pi * y1)
+    u_exact_expr2 = x2 * y2 * (1 - x2) * (1 - y2)
+    f1 = -div(grad(u_exact_expr1))
+    f2 = -div(grad(u_exact_expr2))
+
+    L = (
+        inner(f1, v1)*dx(domain=mesh1)
+        + inner(f2, v2)*dx(domain=mesh2)
+    )
+
+    bc1 = DirichletBC(V.sub(0), 0, "on_boundary")
+    bc2 = DirichletBC(V.sub(1), 0, "on_boundary")
+    u_sol = Function(V)
+    solve(a == L, u_sol, bcs=[bc1, bc2])
+    u1_sol, u2_sol = u_sol.subfunctions
+
+    u_exact = Function(V)
+    u1_exact, u2_exact = u_exact.subfunctions
+    u1_exact.interpolate(u_exact_expr1)
+    u2_exact.interpolate(u_exact_expr2)
+
+    err1 = errornorm(u1_exact, u1_sol)
+    assert err1 < 1e-5
+    err2 = errornorm(u2_exact, u2_sol)
+    assert err2 < 1e-5
+
+
+def test_multi_domain_assemble():
+    mesh1 = UnitSquareMesh(1, 1, quadrilateral=True)
+    mesh2 = UnitSquareMesh(2, 2)
+    V1 = FunctionSpace(mesh1, "Q", 1)
+    V2 = FunctionSpace(mesh2, "CG", 1)
+    V = V1 * V2
+
+    u = TrialFunctions(V)
+    v = TestFunctions(V)
+    f = split(Function(V))
+
+    for i, j in [(0, 1), (1, 0)]:
+        a1 = inner(u[i], v[j])*dx(domain=mesh1)
+        with pytest.raises(NotImplementedError):
+            assemble(a1)
+        a2 = inner(u[i], v[j])*dx(domain=mesh2)
+        with pytest.raises(NotImplementedError):
+            assemble(a2)
+        l1 = inner(f[i], v[j])*dx(domain=mesh1)
+        with pytest.raises(NotImplementedError):
+            assemble(l1)
+        l2 = inner(f[i], v[j])*dx(domain=mesh2)
+        with pytest.raises(NotImplementedError):
+            assemble(l2)
+
+    for i, j in [(0, 0), (1, 1)]:
+        a = inner(u[i], v[j])*dx(domain=mesh1)
+        if i == 1:
+            with pytest.raises(NotImplementedError):
+                assemble(a)
+            continue
+        A = assemble(a)
+        assert A.M.values.shape == (V.dim(), V.dim())
+
+    a = inner(u[0], v[0])*dx(domain=mesh1) + inner(u[0], v[1])*dx(domain=mesh2)
+    with pytest.raises(NotImplementedError):
+        assemble(a)
+
+    a = inner(u[0], v[0])*dx(domain=mesh1) + inner(u[1], v[1])*dx(domain=mesh2)
+    A = assemble(a)
+    assert A.M.values.shape == (V.dim(), V.dim())

tsfc/driver.py

@@ -145,6 +145,11 @@ def compile_integral(integral_data, form_data, prefix, parameters, *, diagonal=F
     domain_integral_type_map.update(dict.fromkeys(coefficient_meshes, "cell"))
     domain_integral_type_map.update(integral_data.domain_integral_type_map)
 
+    for arg in arguments:
+        if domain_integral_type_map[extract_unique_domain(arg)] is None:
+            raise NotImplementedError("Assembly of forms over unrelated meshes is not supported. "
+                                      "Try using Submeshes or cross-mesh interpolation.")
+
     integral_data_info = TSFCIntegralDataInfo(
         domain=integral_data.domain,
         integral_type=integral_data.integral_type,
@@ -156,6 +161,7 @@ def compile_integral(integral_data, form_data, prefix, parameters, *, diagonal=F
         coefficient_split=coefficient_split,
         coefficient_numbers=coefficient_numbers,
     )
+
     builder = firedrake_interface_loopy.KernelBuilder(
         integral_data_info,
         scalar_type,