Fix Real space assembly (#4331)

* Fix Real space assembly

Author: pbrubeck

tests/firedrake/regression/test_real_space.py

@@ -58,15 +58,22 @@ def test_real_nonsquare_two_form_assembly():
 
 
 @pytest.mark.skipcomplex
-def test_real_mixed_one_form_assembly():
+@pytest.mark.parametrize("coefficient", (False, True))
+def test_real_mixed_one_form_assembly(coefficient):
     mesh = UnitIntervalMesh(3)
     rfs = FunctionSpace(mesh, "Real", 0)
     cgfs = FunctionSpace(mesh, "CG", 1)
 
     mfs = cgfs*rfs
     v, q = TestFunctions(mfs)
 
-    A = assemble(conj(v) * dx + q * dx)
+    if coefficient:
+        z = Function(mfs)
+        z.subfunctions[1].assign(1)
+        u, p = split(z)
+        A = assemble(inner(u, v) * dx + inner(p, q) * dx)
+    else:
+        A = assemble(conj(v) * dx + q * dx)
 
     qq = TestFunction(rfs)
 
@@ -131,6 +138,33 @@ def test_real_mixed_empty_component_assembly():
     assemble(derivative(inner(grad(v), grad(v)) * dx, w))
 
 
+@pytest.mark.skipcomplex
+@pytest.mark.parametrize("coefficient", (False, True))
+def test_real_extruded_mixed_one_form_assembly(coefficient):
+    m = UnitIntervalMesh(3)
+    mesh = ExtrudedMesh(m, 10)
+    rfs = FunctionSpace(mesh, "Real", 0)
+    cgfs = FunctionSpace(mesh, "CG", 1)
+
+    mfs = cgfs*rfs
+    v, q = TestFunctions(mfs)
+
+    if coefficient:
+        z = Function(mfs)
+        z.subfunctions[1].assign(1)
+        u, p = split(z)
+        A = assemble(inner(u, v) * dx + inner(p, q) * dx)
+    else:
+        A = assemble(conj(v) * dx + q * dx)
+
+    qq = TestFunction(rfs)
+
+    AA = assemble(qq * dx)
+
+    np.testing.assert_almost_equal(A.dat.data[1],
+                                   AA.dat.data)
+
+
 @pytest.mark.skipcomplex
 def test_real_extruded_mixed_two_form_assembly():
     m = UnitIntervalMesh(3)

tsfc/kernel_interface/common.py

@@ -16,6 +16,7 @@
 from numpy import asarray
 from tsfc import fem, ufl_utils
 from finat.element_factory import as_fiat_cell, create_element
+from finat.ufl import MixedElement
 from tsfc.kernel_interface import KernelInterface
 from tsfc.logging import logger
 from ufl.utils.sequences import max_degree
@@ -298,10 +299,16 @@ def set_quad_rule(params, cell, integral_type, functions):
     # Check if the integral has a quad degree or quad element attached,
     # otherwise use the estimated polynomial degree attached by compute_form_data
     quad_rule = params.get("quadrature_rule", "default")
+    elements = []
+    for f in functions:
+        e = f.ufl_element()
+        if type(e) is MixedElement:
+            elements.extend(e.sub_elements)
+        else:
+            elements.append(e)
     try:
         quadrature_degree = params["quadrature_degree"]
     except KeyError:
-        elements = [f.ufl_function_space().ufl_element() for f in functions]
         quad_data = set((e.degree(), e.quadrature_scheme() or "default") for e in elements
                         if e.family() in {"Quadrature", "Boundary Quadrature"})
         if len(quad_data) == 0:
@@ -320,8 +327,7 @@ def set_quad_rule(params, cell, integral_type, functions):
     if isinstance(quad_rule, str):
         scheme = quad_rule
         fiat_cell = as_fiat_cell(cell)
-        finat_elements = set(create_element(f.ufl_element()) for f in functions
-                             if f.ufl_element().family() != "Real")
+        finat_elements = set(create_element(e) for e in elements if e.family() != "Real")
         fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements]
         fiat_cell = max_complex(fiat_cells)