extract_domains does the wrong thing with interpolate.

[dham]

Consider the following Firedrake MFE:

from firedrake import *
from firedrake.__future__ import interpolate

mesh = UnitSquareMesh(2, 2)
Q = FunctionSpace(mesh, "CG", 1)
u = Function(Q)
vom = VertexOnlyMesh(mesh, ([0.5, 0.5],))
P0DG = FunctionSpace(vom, "DG", 0)

domain = interpolate(u, P0DG).ufl_domain()

Here domain is Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 1) but it should be Mesh(VectorElement(FiniteElement('Discontinuous Lagrange', Cell(vertex, 2), 0), dim=2), 8). I.e. the domain of u is being returned but it should be the domain of P0DG. This is because extract_domains assumes that domains come from terminals, which is not true in this case. A similar issue will arise for external operator.