How to attatch abstract domains to expressions

Suggested design.

# Cells are mapped internally in UFL to a default domain for compatibility:
_default_domains = {
    triangle: Domain(triangle),
    tetrahedron: Domain(tetrahedron),
    #...
    }
def as_domain(D):
    if isinstance(D, Cell):
        return _default_domains[D]
    return D

# Some kind of language for defining domains:
D = Domain(triangle) # A simple wrapper around cell can be the first step
E = Domain(triangle)

# Elements are defined on domains:
U = FiniteElement("CG", D, 1)
V = FiniteElement("CG", E, 1)
W = FiniteElement("CG", D, 2)

# Coefficients are defined on elements as before, but now has a domain associated implicitly:
f = Coefficient(U)
g = Coefficient(V)
h = Coefficient(W)

# Preprocessing a form will now let us recover the set of domains in use
a = f*g*h*dx
domains = a.compute_form_data().domains()
assert domains == [D, E]

# There is a relation between the ordered list of coefficients and the ordered list of domains:
# w0 = f -> D = d0 # d0 = domain 0 just as w0 = coefficient 0
# w1 = g -> E = d1 # d1 = domain 1
# w2 = h -> D = d0
# This can be encoded as an integer mapping:
# { 0: 0, 1: 1, 2: 0 }
# These indices i:j map directly to the dofs w[i][] and cell c[j] in the upgraded ufc interface

# Coefficients in a form must be defined on the integration domain, this can be assumed and detected runtime
a = f*g*dx(D) # Assuming E subset of D such that g is defined in the entire integration domain