Redesign function/expression interpolation

[michalhabera]

Generally, there is a need for "interpolation" when
(in all cases I assume the action of target FE space dual basis on donor object is well defined)

1. object being interpolated is a FE function

   1. matching meshes
      Typical IO postprocessing, P_k -> P_1
   2. non-matching meshes
      Error norm computation, coarsening, refining

2. object being interpolated is a blackbox

   1. C++ compiled expression
      Source term preparation, method of manufactured solutions, ...

   2. UFL expression
      E.g. Scalar function computed as trace of some other tensor FE function, various internal variables updates (ODE solved explicitly), grid transfer operators in multigrid, evaluation of handcoded gradients, ... Reference to older discussion, https://bitbucket.org/fenics-project/dolfin/issues/422/allow-arbitrary-ufl-expressions-passed-to

   3. Other custom class complying to an interface.
      This would allow e.g. local newton solver be interpolated to FE function, think of return mapping in plasticity, local nonlinear system for chemical reactions, ...

---

Ideas how to handle:

1. 1. matching
      1. affine target FE
         generate cell-local interpolation matrix at compile time
      2. non-affine target FE
         generate interpolation kernel on a symbolic level - at runtime takes geometric information as pushforwards, pullbacks
   2. non-matching
      1. affine target FE
         generate interpolation kernel on a symbolic level - at runtime takes geometric information as point where to evaluate
      2. non-affine target FE
         generate interpolation kernel on a symbolic level - at runtime takes geometric information as pushforwards, pullbacks, point where to evaluate
2. 1. deprecate in favour of 2.ii.
   2. generate ufl expression kernel
      Check how this works in TSFC. Possibly employ TSFC.
   3. keep the old approach
      generate evaluate_dof which runs .eval of the blackbox.

The main difference between 1. and 2. is existence of underlying basis, linear and nodal structure.

[michalhabera]

This touches issues #7 and #2

Usecase of "pointwise" expression (expression evaluated at quadrature points) could be handled

1. (2.2 from above) Expression is expressible via UFL primitives, e.g. x[0]*x[1] with x being SpatialCoordinate, ...
   Inclusion into UFL form is trivial and is naturally evaluated at quadrature points. This is the most common case.

2. (2.3 from above) Expression is more of a blackbox, an abstract expression mentioned in issue Distinguish between abstract Expression and a FE function #7, non-expressible via UFL.
   In such case user should be forced to "interpolate" (in the sense above) this expression to FE function of Quadrature element type. It is the soft approach 1. from issue Pointwise expressions #2 . Disadvantage that "users might find it ugly to use" is valid, but only advanced users would really need to interpolate non-UFL blackbox to FE function.

[blechta]

Consider also FEniCS/ffcx#33.

[michalhabera]

To support the point 2.ii, here is simple demonstration of TSFC expression kernel generation. (sidenote: merged in one python file just for brevity)

import cffi
import importlib
import numpy as np

import ufl
import tsfc

el = ufl.FiniteElement("P", "triangle", 1)
u = ufl.Coefficient(el)

a = ufl.exp(u) + u
cexp = tsfc.compile_expression_at_points(a, [[1.0, 0]], a)

CODE_C = cexp[0].gencode()
CODE_H = """
static inline void expression_kernel (double *restrict A, const double *restrict w_0);
"""

ffi = cffi.FFI()
ffi.set_source("interpolator", CODE_C)
ffi.cdef(CODE_H)
ffi.compile()

from interpolator.lib import expression_kernel

# Target scalar value
A = np.array([0.0])

# Cell-local expansion coeffs of function
coeffs = np.array([1.0, 1.0, 1.0])

pA = ffi.cast("double *", A.ctypes.data)
pcoeffs = ffi.cast("double *", coeffs.ctypes.data)

expression_kernel(pA, pcoeffs)
print(A)

It evaluates expression exp(u)+u for u being P1 FE function at point (1.0, 0.0).

From the current implementation of interpolation in Firedrake, it is evident that only PointEvaluation target dual basis are allowed (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/interpolation.py#L146). Therefore TSFC can be currently utilized only for simple (affine) FE with point evaluated dofs and for the "matching" meshes case (i.e. point 1.i.a).

[jorgensd]

This issue should be revisited by @mscroggs once #1324 is merged.

[garth-wells]

Fixed by an number of PRs.