TSFC: Dual evaluation for nodal EnrichedElement (#5184)

Enable xfailed H(div) and H(curl) tests

Author: pbrubeck

tests/firedrake/regression/test_interpolate.py

@@ -175,6 +175,23 @@ def test_compound_expression():
     assert np.allclose(g.dat.data, h.dat.data)
 
 
+def test_restricted_extruded_interval():
+    mesh = ExtrudedMesh(UnitIntervalMesh(1), 1)
+    element = FiniteElement("Q", degree=4)
+    VF = FunctionSpace(mesh, element["facet"])
+    VI = FunctionSpace(mesh, element["interior"])
+    uf = Function(VF)
+    ui = Function(VI)
+
+    x = SpatialCoordinate(mesh)
+    exprs = [Constant(1), x[0], x[1], x[0]*x[1]]
+
+    for expr in exprs:
+        uf.interpolate(expr)
+        ui.interpolate(expr)
+        assert norm(expr - (uf + ui)) < 1E-12
+
+
 def test_hdiv_extruded_interval():
     mesh = ExtrudedMesh(UnitIntervalMesh(10), 10, 0.1)
     x = SpatialCoordinate(mesh)
@@ -209,8 +226,6 @@ def test_dpc_into_dq_extruded_interval():
     assert errornorm(u2, u1) < 1E-12
 
 
-# Requires the relevant FInAT or FIAT duals to be defined
-@pytest.mark.xfail(raises=NotImplementedError, reason="Requires the relevant FInAT or FIAT duals to be defined")
 def test_hdiv_2d():
     mesh = UnitCubedSphereMesh(2)
     x = SpatialCoordinate(mesh)
@@ -230,7 +245,6 @@ def test_hdiv_2d():
     assert np.allclose(g.dat.data, h.dat.data)
 
 
-@pytest.mark.xfail(raises=NotImplementedError, reason="Requires the relevant FInAT or FIAT duals to be defined")
 def test_hcurl_2d():
     mesh = UnitCubedSphereMesh(2)
     x = SpatialCoordinate(mesh)

tests/firedrake/vertexonly/test_interpolation_from_parent.py

@@ -3,7 +3,6 @@
 import numpy as np
 from functools import reduce
 from operator import add
-import subprocess
 
 
 # Utility Functions and Fixtures
@@ -70,12 +69,9 @@ def fs(request):
                         ("N2curl", 2, FunctionSpace),
                         ("N1div", 2, FunctionSpace),
                         ("N2div", 2, FunctionSpace),
-                        pytest.param(("RTCE", 2, FunctionSpace),
-                                     marks=pytest.mark.xfail(raises=(subprocess.CalledProcessError, NotImplementedError),
-                                                             reason="EnrichedElement dual basis not yet defined and FIAT duals don't have a point_dict")),
-                        pytest.param(("RTCF", 2, FunctionSpace),
-                                     marks=pytest.mark.xfail(raises=(subprocess.CalledProcessError, NotImplementedError),
-                                                             reason="EnrichedElement dual basis not yet defined and FIAT duals don't have a point_dict"))],
+                        ("RTCE", 2, FunctionSpace),
+                        ("RTCF", 2, FunctionSpace),
+                        ],
                 ids=lambda x: f"{x[2].__name__}({x[0]}{x[1]})")
 def vfs(request, parentmesh):
     family = request.param[0]

tsfc/driver.py

@@ -2,6 +2,7 @@
 import time
 import sys
 from itertools import chain
+import numpy
 from finat.physically_mapped import NeedsCoordinateMappingElement
 
 import ufl
@@ -13,6 +14,7 @@
 
 import gem
 import gem.impero_utils as impero_utils
+from gem.unconcatenate import unconcatenate
 
 import finat
 from finat.element_factory import as_fiat_cell
@@ -367,13 +369,15 @@ def predicate(index):
     # Build kernel body
     return_indices = tuple(chain.from_iterable(argument_multiindices.values()))
     return_shape = tuple(i.extent for i in return_indices)
-    return_var = gem.Variable('A', return_shape or (1,))
-    return_expr = gem.Indexed(return_var, return_indices or (0,))
+    return_var = gem.Variable('A', (numpy.prod(return_shape, dtype=int),))
+    return_expr = gem.Indexed(gem.reshape(return_var, return_shape), return_indices)
+    return_expr, = gem.optimise.remove_componenttensors([return_expr])
 
     # TODO: one should apply some GEM optimisations as in assembly,
     # but we don't for now.
     evaluation, = impero_utils.preprocess_gem([evaluation])
-    impero_c = impero_utils.compile_gem([(return_expr, evaluation)], return_indices)
+    pairs = unconcatenate([(return_expr, evaluation)])
+    impero_c = impero_utils.compile_gem(pairs, return_indices)
     index_names = {idx: f"p{i}" for (i, idx) in enumerate(basis_indices)}
     # Handle kernel interface requirements
     builder.register_requirements([evaluation])