def __init__(self, *cells):
        "Create a ProductCell from a given list of cells."

        self._cells = cells
        ufl_assert(len(self._cells) > 0, "Expecting at least one cell")

        self._cellname = self._cells[0].cellname()#" x ".join([c.cellname() for c in cells])
        self._topological_dimension = sum(c.topological_dimension() for c in cells)
        self._geometric_dimension = sum(c.geometric_dimension() for c in cells)

        self._repr = "ProductCell(*%r)" % list(self._cells)

        self._n = None
        self._x = SpatialCoordinate(self) # For now

        self._xi = None # ?
        self._J = None # ?
        self._Jinv = None # ?
        self._detJ = None # ?

        self._volume = None
        self._circumradius = None
        self._cellsurfacearea = None
        self._facetarea = None            # Not defined
        self._facetdiameter = None        # Not defined

 def list_tensor(self, f): # TODO: Is this right? Fix for higher order tensors too.
     "d/dx_i [x_0, ..., x_n-1] = e_i (unit vector)"
     ops = f.operands()
     n = len(ops)
     s = ops[0].shape()
     ufl_assert(s == (), "TODO: Assuming a vector, i.e. scalar operands.")
     return unit_vectors(n) # TODO: Non-scalars

 def __init__(self, arg1, arg2):
     Operator.__init__(self)
     ufl_assert(is_true_ufl_scalar(arg1), "Expecting scalar argument 1.")
     ufl_assert(is_true_ufl_scalar(arg2), "Expecting scalar argument 2.")
     self._name     = "atan_2"
     self._arg1 = arg1
     self._arg2 = arg2

def diag(A):
    """UFL operator: Take the diagonal part of rank 2 tensor A _or_
    make a diagonal rank 2 tensor from a rank 1 tensor.

    Always returns a rank 2 tensor. See also diag_vector."""

    # TODO: Make a compound type or two for this operator

    # Get and check dimensions
    r = A.rank()
    if r == 1:
        n, = A.shape()
    elif r == 2:
        m, n = A.shape()
        ufl_assert(m == n, "Can only take diagonal of square tensors.")
    else:
        error("Expecting rank 1 or 2 tensor.")

    # Build matrix row by row
    rows = []
    for i in range(n):
        row = [0]*n
        row[i] = A[i] if r == 1 else A[i,i]
        rows.append(row)
    return as_matrix(rows)

 def _square_matrix_shape(self, A):
     sh = A.shape()
     if self._dim is not None:
         sh = complete_shape(sh, self._dim) # FIXME: Does this ever happen? Pretty sure other places assume shapes are always "complete".
     ufl_assert(sh[0] == sh[1], "Expecting square matrix.")
     ufl_assert(sh[0] is not None, "Unknown dimension.")
     return sh

 def _variable_derivative(self, o, f, v):
     f, fp = f
     v, vp = v
     ufl_assert(isinstance(vp, Zero), "TODO: What happens if vp != 0, i.e. v depends the differentiation variable?")
     # Are there any issues with indices here? Not sure, think through it...
     oprime = o.reconstruct(fp, v)
     return (o, oprime)

 def restricted(self, o):
     ufl_assert(self.current_restriction is None,
         "Not expecting twice restricted expression.")
     self.current_restriction = o._side
     e, = o.operands()
     self.visit(e)
     self.current_restriction = None

 def component_tensor(self, f):
     x, i = f.operands()
     s = f.shape()
     ufl_assert(len(s) == 1, "TODO: Assuming a vector, i.e. scalar operands.")
     n, = s
     d = ListTensor([1]*n) # TODO: Non-scalars
     return (d, None)

    def __init__(self, *elements):
        "Create TensorProductElement from a given list of elements."

        self._sub_elements = list(elements)
        ufl_assert(len(self._sub_elements) > 0,
                   "Cannot create TensorProductElement from empty list.")
        self._repr = "TensorProductElement(*%r)" % self._sub_elements
        family = "TensorProductElement"

        # Define cell as the product of each subcell
        cell = ProductCell(*[e.cell() for e in self._sub_elements])
        domain = as_domain(cell) # FIXME: figure out what this is supposed to mean :)

        # Define polynomial degree as the maximal of each subelement
        degree = max(e.degree() for e in self._sub_elements)

        # No quadrature scheme defined
        quad_scheme = None

        # For now, check that all subelements have the same value
        # shape, and use this.
        # TODO: Not sure if this makes sense, what kind of product is used to build the basis?
        value_shape = self._sub_elements[0].value_shape()
        ufl_assert(all(e.value_shape() == value_shape
                       for e in self._sub_elements),
                   "All subelements in must have same value shape")

        super(TensorProductElement, self).__init__(family, domain, degree,
                                                   quad_scheme, value_shape)

    def extract_subelement_component(self, i):
        """Extract direct subelement index and subelement relative
        component index for a given component index"""
        if isinstance(i, int):
            i = (i,)
        self._check_component(i)

        # Select between indexing modes
        if len(self.value_shape()) == 1:
            # Indexing into a long vector of flattened subelement shapes
            j, = i

            # Find subelement for this index
            for k, e in enumerate(self._sub_elements):
                sh = e.value_shape()
                si = product(sh)
                if j < si:
                    break
                j -= si
            ufl_assert(j >= 0, "Moved past last value component!")

            # Convert index into a shape tuple
            j = index_to_component(j, sh)
        else:
            # Indexing into a multidimensional tensor
            # where subelement index is first axis
            k = i[0]
            ufl_assert(k < len(self._sub_elements),
                       "Illegal component index (dimension %d)." % k)
            j = i[1:]
        return (k, j)

 def __add__(self, other):
     "Add two elements, creating an enriched element"
     ufl_assert(isinstance(other, FiniteElementBase),
                "Can't add element and %s." % other.__class__)
     warning_blue("WARNING: Creating an EnrichedElement,\n         " +\
                  "if you intended to create a MixedElement use '*' instead of '+'.")
     from ufl.finiteelement import EnrichedElement
     return EnrichedElement(self, other)

 def __init__(self, spatial_dim, coefficients, arguments, coefficient_derivatives, cache=None):
     ForwardAD.__init__(self, spatial_dim, var_shape=(), var_free_indices=(),
                        var_index_dimensions={}, cache=cache)
     self._v = arguments
     self._w = coefficients
     self._cd = coefficient_derivatives
     ufl_assert(isinstance(self._w, Tuple), "Expecting a Tuple.")
     ufl_assert(isinstance(self._v, Tuple), "Expecting a Tuple.")

    def d(self):
        """The dimension of the cell.

        Only valid if the geometric and topological dimensions are the same."""
        ufl_assert(self._topological_dimension == self._geometric_dimension,
                   "Cell.d is undefined when geometric and"+\
                   "topological dimensions are not the same.")
        return self._geometric_dimension

 def __init__(self, integrals):
     self._dintegrals = integral_sequence_to_dict(integrals)
     ufl_assert(all(isinstance(itg, Integral) for itg in integrals),
                "Expecting list of integrals.")
     self._signature = None
     self._hash = None
     self._form_data = None
     self._is_preprocessed = False

 def reconstruct(self, **kwargs):
     kwargs["family"] = kwargs.get("family", self.family())
     kwargs["domain"] = kwargs.get("domain", self.domain())
     kwargs["degree"] = kwargs.get("degree", self.degree())
     ufl_assert("dim" not in kwargs, "Cannot change dim in reconstruct.")
     kwargs["dim"] = len(self._sub_elements)
     kwargs["quad_scheme"] = kwargs.get("quad_scheme", self.quadrature_scheme())
     return VectorElement(**kwargs)

 def reconstruct_from_elements(self, *elements):
     "Reconstruct a mixed element from new subelements."
     if all(a == b for (a,b) in izip(elements, self._sub_elements)):
         return self
     ufl_assert(all(a.value_shape() == b.value_shape()
                    for (a,b) in izip(elements, self._sub_elements)),
         "Expecting new elements to have same value shape as old ones.")
     return MixedElement(*elements, value_shape=self.value_shape())

 def conditional(self, f): # TODO: Is this right? What about non-scalars?
     c, a, b = f.operands()
     s = f.shape()
     ufl_assert(s == (), "TODO: Assuming scalar valued expressions.")
     _0 = Zero()
     _1 = IntValue(1)
     da = conditional(c, _1, _0)
     db = conditional(c, _0, _1)
     return (None, da, db)

 def __new__(cls, A):
     sh = A.shape()
     ufl_assert(len(sh) == 2, "Symmetric part of tensor with rank != 2 is undefined.")
     if sh[0] != sh[1]:
         error("Cannot take symmetric part of rectangular matrix with dimensions %s." % repr(sh))
     ufl_assert(not A.free_indices(), "Not expecting free indices in Sym.")
     if isinstance(A, Zero):
         return Zero(A.shape(), A.free_indices(), A.index_dimensions())
     return CompoundTensorOperator.__new__(cls)

    def scalar_value(self, x):
        if len(x.shape()) != len(self.component()):
            self.print_visit_stack()
        ufl_assert(len(x.shape()) == len(self.component()), "Component size mismatch.")

        s = set(x.free_indices()) - set(self._index2value.keys())
        if s: error("Free index set mismatch, these indices have no value assigned: %s." % str(s))

        return x._uflclass(x.value())

def extract_element_map(elements):
    "Build map from elements to element index in ordered tuple."
    element_map = {}
    unique_elements = unique_tuple(elements)
    for element in elements:
        indices = [i for (i, e) in enumerate(unique_elements) if e == element]
        ufl_assert(len(indices) == 1, "Unable to find unique index for element.")
        element_map[element] = i
    return element_map