Unable to interpolate into RTCF space

[tommbendall]

This might be well-known already, but @ReubenHill suggested I add this here as it might be helpful for testing dual evaluation.

Minimum example:

from firedrake import *

mesh = UnitSquareMesh(3,3,quadrilateral=True)
x, y = SpatialCoordinate(mesh)

V = FunctionSpace(mesh, "RTCF", 1)
f = Function(V)

f.interpolate(as_vector([x,y]))

gives an error like the following:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last) <ipython-input-6-b9757afeb453> in <module>
----> 1 f.interpolate(as_vector([x,y]))

~/src/firedrake/src/firedrake/firedrake/function.py in interpolate(self, expression, subset)
     353         :returns: this :class:`Function` object"""
     354         from firedrake import interpolation
--> 355         return interpolation.interpolate(expression, self,
subset=subset)
     356
     357     @FunctionMixin._ad_annotate_assign

~/src/firedrake/src/firedrake/firedrake/interpolation.py in interpolate(expr, V, subset, access)
      48        performance by using an :class:`Interpolator` instead.
      49     """
---> 50     return Interpolator(expr, V, subset=subset,
access=access).interpolate()
      51
      52

~/src/firedrake/src/firedrake/firedrake/interpolation.py in __init__(self, expr, V, subset, freeze_expr, access)
      71     def __init__(self, expr, V, subset=None, freeze_expr=False,
access=op2.WRITE):
      72         try:
---> 73             self.callable, arguments = make_interpolator(expr,
V, subset, access)
      74         except FIAT.hdiv_trace.TraceError:
      75             raise NotImplementedError("Can't interpolate onto traces sorry")

~/src/firedrake/src/firedrake/firedrake/interpolation.py in make_interpolator(expr, V, subset, access)
     187             raise NotImplementedError(
     188                 "UFL expressions for mixed functions are not yet supported.")
--> 189         loops.extend(_interpolator(V, tensor, expr, subset,
arguments, access))
     190     elif hasattr(expr, 'eval'):
     191         if len(V) > 1:

<decorator-gen-141> in _interpolator(V, tensor, expr, subset, arguments,
access)

~/src/firedrake/src/firedrake/firedrake/utils.py in wrapper(f, *args,
**kwargs)
      71         opts["type_check"] = safe
      72         try:
---> 73             return f(*args, **kwargs)
      74         finally:
      75             opts["type_check"] = check

~/src/firedrake/src/firedrake/firedrake/interpolation.py in _interpolator(V, tensor, expr, subset, arguments, access)
235 domain=V.mesh(),
236 parameters=parameters,
--> 237 coffee=False)
     238         kernel = op2.Kernel(ast, ast.name,
requires_zeroed_output_arguments=True)
     239     elif hasattr(expr, "eval"):

~/src/firedrake/src/tsfc/tsfc/driver.py in compile_expression_dual_evaluation(expression, to_element, domain, interface, parameters, coffee)
     389         expr_cache = {}         # Sharing of evaluation of the expression at points
     390         for dual in to_element.dual_basis():
--> 391             pts = tuple(sorted(dual.get_point_dict().keys()))
     392             try:
     393                 expr, point_set = expr_cache[pts]

AttributeError: 'NoneType' object has no attribute 'keys'

[ReubenHill]

I can confirm that the code snippet does run on the dual_base_fiat branch (i.e. this is not an issue in #1743 and its associated FInAT and TSFC branches). Nevertheless I was under the impression from conversations with @dham that this should work on master albeit inefficiently. Perhaps @wence- or anyone else that originally added dual evaluation can comment? (FYI I haven't looked into the exact cause of the issue here as yet.)

[wence-]

Nevertheless I was under the impression from conversations with @dham that this should work on master albeit inefficiently. Perhaps @wence- or anyone else that originally added dual evaluation can comment? (FYI I haven't looked into the exact cause of the issue here as yet.)

It looks like the RTCF elements in master don't provide a dual basis?

[dham]

They probably don't. @ReubenHill I don't think we discussed RTCF. I think we discussed RT.

[wence-]

It looks like the RTCF elements in master don't provide a dual basis?

Yeah, the implementation in FIAT of the Hdivelement does:

    # splat any PointEvaluation functionals.
    # they become a nasty mix of internal and external component DOFs
    if newelement._oldmapping == "affine":
        oldnodes = newelement.dual.nodes
        newnodes = []
        for node in oldnodes:
            if isinstance(node, functional.PointEvaluation):
                newnodes.append(functional.Functional(None, None, None, {}, "Undefined"))
            else:
                newnodes.append(node)
        newelement.dual.nodes = newnodes

So we fake nodes with no point dicts.

[ReubenHill]

See also #2077