def replace_integral_domains(form, common_domain):  # TODO: Move elsewhere
    """Given a form and a domain, assign a common integration domain to
    all integrals.

    Does not modify the input form (``Form`` should always be
    immutable).  This is to support ill formed forms with no domain
    specified, sometimes occurring in pydolfin, e.g. assemble(1*dx,
    mesh=mesh).

    """
    domains = form.ufl_domains()
    if common_domain is not None:
        gdim = common_domain.geometric_dimension()
        tdim = common_domain.topological_dimension()
        if not all((gdim == domain.geometric_dimension() and
                    tdim == domain.topological_dimension())
                   for domain in domains):
            error("Common domain does not share dimensions with form domains.")

    reconstruct = False
    integrals = []
    for itg in form.integrals():
        domain = itg.ufl_domain()
        if domain != common_domain:
            itg = itg.reconstruct(domain=common_domain)
            reconstruct = True
        integrals.append(itg)
    if reconstruct:
        form = Form(integrals)
    return form

def form_info(form):
    if not isinstance(form, Form):
        error("Expecting a Form.")

    bf = form.arguments()
    cf = form.coefficients()

    s = "Form info:\n"
    s += "  rank:                          %d\n" % len(bf)
    s += "  num_coefficients:              %d\n" % len(cf)
    s += "\n"

    for f in cf:
        if f._name:
            s += "\n"
            s += "  Coefficient %d is named '%s'" % (f._count, f._name)
    s += "\n"

    integrals = form.integrals()
    integral_types = sorted(set(itg.integral_type() for itg in integrals))
    for integral_type in integral_types:
        itgs = form.integrals_by_type(integral_type)
        s += "  num_{0}_integrals:  {1}\n".format(integral_type, len(itgs))
    s += "\n"

    for integral_type in integral_types:
        itgs = form.integrals_by_type(integral_type)
        for itg in itgs:
            s += integral_info(itg)
            s += "\n"

    return s

    def extract_subelement_reference_component(self, i):
        """Extract direct subelement index and subelement relative
        reference_component index for a given reference_component index."""
        if isinstance(i, int):
            i = (i,)
        self._check_reference_component(i)

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

        # Find subelement for this index
        for sub_element_index, e in enumerate(self._sub_elements):
            sh = e.reference_value_shape()
            si = product(sh)
            if j < si:
                break
            j -= si
        if j < 0:
            error("Moved past last value reference_component!")

        # Convert index into a shape tuple
        st = shape_to_strides(sh)
        reference_component = unflatten_index(j, st)
        return (sub_element_index, reference_component)

 def math_function(self, o, df):
     # FIXME: Introduce a UserOperator type instead of this hack
     # and define user derivative() function properly
     if hasattr(o, 'derivative'):
         f, = o.ufl_operands
         return df * o.derivative()
     error("Unknown math function.")

    def __init__(self, *expressions):
        Operator.__init__(self, expressions)

        # Checks
        indexset = set(self.ufl_operands[0].ufl_free_indices)
        if not all(not (indexset ^ set(e.ufl_free_indices)) for e in self.ufl_operands):
            error("Can't combine subtensor expressions with different sets of free indices.")

def register_element2(family, value_rank, sobolev_space, mapping,
                      degree_range, cellnames):
    "Register new finite element family."
    if family in ufl_elements:
        error('Finite element \"%s\" has already been registered.' % family)
    ufl_elements[family] = (family, family, value_rank, sobolev_space,
                            mapping, degree_range, cellnames)

 def spatial_coordinate(self, o):
     do = self._w2v.get(o)
     # d x /d x => Argument(x.function_space())
     if do is not None:
         return do
     else:
         error("Not implemented: CoordinateDerivative found a SpatialCoordinate that is different from the one being differentiated.")

 def cofactor(self, o, A):
     # TODO: Find common subexpressions here.
     # TODO: Better implementation?
     sh = self._square_matrix_shape(A)
     if sh[0] == 2:
         return as_matrix([[A[1,1],-A[1,0]],[-A[0,1],A[0,0]]])
     elif sh[0] == 3:
         return as_matrix([
             [A[1,1]*A[2,2] - A[2,1]*A[1,2],A[2,0]*A[1,2] - A[1,0]*A[2,2], - A[2,0]*A[1,1] + A[1,0]*A[2,1]],
             [A[2,1]*A[0,2] - A[0,1]*A[2,2],A[0,0]*A[2,2] - A[2,0]*A[0,2], - A[0,0]*A[2,1] + A[2,0]*A[0,1]],
             [A[0,1]*A[1,2] - A[1,1]*A[0,2],A[1,0]*A[0,2] - A[0,0]*A[1,2], - A[1,0]*A[0,1] + A[0,0]*A[1,1]]
             ])
     elif sh[0] == 4:
         return as_matrix([
             [-A[3,1]*A[2,2]*A[1,3] - A[3,2]*A[2,3]*A[1,1] + A[1,3]*A[3,2]*A[2,1] + A[3,1]*A[2,3]*A[1,2] + A[2,2]*A[1,1]*A[3,3] - A[3,3]*A[2,1]*A[1,2],
              -A[1,0]*A[2,2]*A[3,3] + A[2,0]*A[3,3]*A[1,2] + A[2,2]*A[1,3]*A[3,0] - A[2,3]*A[3,0]*A[1,2] + A[1,0]*A[3,2]*A[2,3] - A[1,3]*A[3,2]*A[2,0],
               A[1,0]*A[3,3]*A[2,1] + A[2,3]*A[1,1]*A[3,0] - A[2,0]*A[1,1]*A[3,3] - A[1,3]*A[3,0]*A[2,1] - A[1,0]*A[3,1]*A[2,3] + A[3,1]*A[1,3]*A[2,0],
               A[3,0]*A[2,1]*A[1,2] + A[1,0]*A[3,1]*A[2,2] + A[3,2]*A[2,0]*A[1,1] - A[2,2]*A[1,1]*A[3,0] - A[3,1]*A[2,0]*A[1,2] - A[1,0]*A[3,2]*A[2,1]],
             [ A[3,1]*A[2,2]*A[0,3] + A[0,2]*A[3,3]*A[2,1] + A[0,1]*A[3,2]*A[2,3] - A[3,1]*A[0,2]*A[2,3] - A[0,1]*A[2,2]*A[3,3] - A[3,2]*A[0,3]*A[2,1],
              -A[2,2]*A[0,3]*A[3,0] - A[0,2]*A[2,0]*A[3,3] - A[3,2]*A[2,3]*A[0,0] + A[2,2]*A[3,3]*A[0,0] + A[0,2]*A[2,3]*A[3,0] + A[3,2]*A[2,0]*A[0,3],
               A[3,1]*A[2,3]*A[0,0] - A[0,1]*A[2,3]*A[3,0] - A[3,1]*A[2,0]*A[0,3] - A[3,3]*A[0,0]*A[2,1] + A[0,3]*A[3,0]*A[2,1] + A[0,1]*A[2,0]*A[3,3],
               A[3,2]*A[0,0]*A[2,1] - A[0,2]*A[3,0]*A[2,1] + A[0,1]*A[2,2]*A[3,0] + A[3,1]*A[0,2]*A[2,0] - A[0,1]*A[3,2]*A[2,0] - A[3,1]*A[2,2]*A[0,0]],
             [ A[3,1]*A[1,3]*A[0,2] - A[0,2]*A[1,1]*A[3,3] - A[3,1]*A[0,3]*A[1,2] + A[3,2]*A[1,1]*A[0,3] + A[0,1]*A[3,3]*A[1,2] - A[0,1]*A[1,3]*A[3,2],
               A[1,3]*A[3,2]*A[0,0] - A[1,0]*A[3,2]*A[0,3] - A[1,3]*A[0,2]*A[3,0] + A[0,3]*A[3,0]*A[1,2] + A[1,0]*A[0,2]*A[3,3] - A[3,3]*A[0,0]*A[1,2],
              -A[1,0]*A[0,1]*A[3,3] + A[0,1]*A[1,3]*A[3,0] - A[3,1]*A[1,3]*A[0,0] - A[1,1]*A[0,3]*A[3,0] + A[1,0]*A[3,1]*A[0,3] + A[1,1]*A[3,3]*A[0,0],
               A[0,2]*A[1,1]*A[3,0] - A[3,2]*A[1,1]*A[0,0] - A[0,1]*A[3,0]*A[1,2] - A[1,0]*A[3,1]*A[0,2] + A[3,1]*A[0,0]*A[1,2] + A[1,0]*A[0,1]*A[3,2]],
             [ A[0,3]*A[2,1]*A[1,2] + A[0,2]*A[2,3]*A[1,1] + A[0,1]*A[2,2]*A[1,3] - A[2,2]*A[1,1]*A[0,3] - A[1,3]*A[0,2]*A[2,1] - A[0,1]*A[2,3]*A[1,2],
               A[1,0]*A[2,2]*A[0,3] + A[1,3]*A[0,2]*A[2,0] - A[1,0]*A[0,2]*A[2,3] - A[2,0]*A[0,3]*A[1,2] - A[2,2]*A[1,3]*A[0,0] + A[2,3]*A[0,0]*A[1,2],
              -A[0,1]*A[1,3]*A[2,0] + A[2,0]*A[1,1]*A[0,3] + A[1,3]*A[0,0]*A[2,1] - A[1,0]*A[0,3]*A[2,1] + A[1,0]*A[0,1]*A[2,3] - A[2,3]*A[1,1]*A[0,0],
               A[1,0]*A[0,2]*A[2,1] - A[0,2]*A[2,0]*A[1,1] + A[0,1]*A[2,0]*A[1,2] + A[2,2]*A[1,1]*A[0,0] - A[1,0]*A[0,1]*A[2,2] - A[0,0]*A[2,1]*A[1,2]]
             ])
     error("Cofactor not implemented for dimension %s." % sh[0])

 def deviatoric(self, o, A):
     sh = self._square_matrix_shape(A)
     if sh[0] == 2:
         return as_matrix([[-1./2*A[1,1]+1./2*A[0,0],A[0,1]],[A[1,0],1./2*A[1,1]-1./2*A[0,0]]])
     elif sh[0] == 3:
         return as_matrix([[-1./3*A[1,1]-1./3*A[2,2]+2./3*A[0,0],A[0,1],A[0,2]],[A[1,0],2./3*A[1,1]-1./3*A[2,2]-1./3*A[0,0],A[1,2]],[A[2,0],A[2,1],-1./3*A[1,1]+2./3*A[2,2]-1./3*A[0,0]]])
     error("dev(A) not implemented for dimension %s." % sh[0])

def _check_form_arity(preprocessed_form):
    # Check that we don't have a mixed linear/bilinear form or
    # anything like that
    # FIXME: This is slooow and should be moved to form compiler
    # and/or replaced with something faster
    if 1 != len(compute_form_arities(preprocessed_form)):
        error("All terms in form must have same rank.")

def ufl2dot(expression, formname="a", nodeoffset=0, begin=True, end=True,
            labeling="repr", object_names=None):
    if labeling == "repr":
        labeller = ReprLabeller()
    elif labeling == "compact":
        labeller = CompactLabeller(object_names or {})
        print(object_names)

    if isinstance(expression, Form):
        form = expression

        subgraphs = []
        k = 0
        for itg in form.integrals():
            prefix = "itg%d_" % k
            integralkey = "%s%s" % (itg.integral_type(), itg.subdomain_id())

            integrallabel = "%s %s" % (itg.integral_type().capitalize().replace("_", " "), "integral")
            integrallabel += " %s" % (itg.subdomain_id(),)

            integrand = itg.integrand()

            nodes = {}
            edges = []

            build_entities(integrand, nodes, edges, nodeoffset, prefix,
                           labeller)
            rootnode = nodes[id(integrand)][0]
            entitylist = format_entities(nodes, edges)
            integralnode = "%s_%s" % (formname, integralkey)
            subgraphs.append(integralgraphformat % {
                'node': integralnode,
                'label': integrallabel,
                'formname': formname,
                'root': rootnode,
                'entities': entitylist, })
            nodeoffset += len(nodes)

        s = ""
        if begin:
            s += 'digraph ufl_form\n{\n  node [shape="box"] ;\n'
        s += '  form_%s [label="Form %s"] ;' % (formname, formname)
        s += "\n".join(subgraphs)
        if end:
            s += "\n}"

    elif isinstance(expression, Expr):
        nodes = {}
        edges = []

        build_entities(expression, nodes, edges, nodeoffset, '', labeller)
        entitylist = format_entities(nodes, edges)
        s = exprgraphformat % entitylist

        nodeoffset += len(nodes)

    else:
        error("Invalid object type %s" % type(expression))

    return s, nodeoffset

def _check_elements(form_data):
    for element in chain(form_data.unique_elements,
                         form_data.unique_sub_elements):
        if element.family() is None:
            error("Found element with undefined familty: %s" % repr(element))
        if element.cell() is None:
            error("Found element with undefined cell: %s" % repr(element))

def multiply_block_interior_facets(point_index, unames, ttypes, unique_tables, unique_table_num_dofs):
    rank = len(unames)
    tables = [unique_tables.get(name) for name in unames]
    num_dofs = tuple(unique_table_num_dofs[name] for name in unames)

    num_entities = max([1] + [tbl.shape[0] for tbl in tables if tbl is not None])
    ptable = numpy.zeros((num_entities,)*rank + num_dofs)
    for facets in itertools.product(*[range(num_entities)]*rank):
        vectors = []
        for i, tbl in enumerate(tables):
            if tbl is None:
                assert ttypes[i] == "ones"
                vectors.append(numpy.ones((num_dofs[i],)))
            else:
                # Some tables are compacted along entities or points
                e = 0 if tbl.shape[0] == 1 else facets[i]
                q = 0 if tbl.shape[1] == 1 else point_index
                vectors.append(tbl[e, q, :])
        if rank > 1:
            assert rank == 2
            ptable[facets[0], facets[1], ...] = numpy.outer(*vectors)
        elif rank == 1:
            ptable[facets[0], :] = vectors[0]
        else:
            error("Nothing to multiply!")

    return ptable

    def __init__(self, integrals):
        # Basic input checking (further compatibilty analysis happens
        # later)
        if not all(isinstance(itg, Integral) for itg in integrals):
            error("Expecting list of integrals.")

        # Store integrals sorted canonically to increase signature
        # stability
        self._integrals = _sorted_integrals(integrals)

        # Internal variables for caching domain data
        self._integration_domains = None
        self._domain_numbering = None

        # Internal variables for caching subdomain data
        self._subdomain_data = None

        # Internal variables for caching form argument data
        self._arguments = None
        self._coefficients = None
        self._coefficient_numbering = None

        # Internal variables for caching of hash and signature after
        # first request
        self._hash = None
        self._signature = None

        # Never use this internally in ufl!
        self._cache = {}

 def math_function(self, o, a):
     if hasattr(o, 'derivative'): # FIXME: Introduce a UserOperator type instead of this hack
         f, fp = a
         o = self.reuse_if_possible(o, f)
         op = fp * o.derivative()
         return (o, op)
     error("Unknown math function.")

def extract_arguments_and_coefficients(a):
    """Build two sorted lists of all arguments and coefficients
    in a, which can be a Form, Integral or Expr."""

    # This function is faster than extract_arguments + extract_coefficients
    # for large forms, and has more validation built in.

    # Extract lists of all form argument instances
    terminals = extract_type(a, FormArgument)
    arguments = [f for f in terminals if isinstance(f, Argument)]
    coefficients = [f for f in terminals if isinstance(f, Coefficient)]

    # Build count: instance mappings, should be one to one
    bfcounts = dict((f, f.count()) for f in arguments)
    fcounts = dict((f, f.count()) for f in coefficients)

    if len(bfcounts) != len(set(bfcounts.values())):
        msg = """\
Found different Arguments with same counts.
Did you combine test or trial functions from different spaces?
The Arguments found are:\n%s""" % "\n".join("  %s" % f for f in arguments)
        error(msg)

    if len(fcounts) != len(set(fcounts.values())):
        msg = """\
Found different coefficients with same counts.
The arguments found are:\n%s""" % "\n".join("  %s" % f for f in coefficients)
        error(msg)

    # Passed checks, so we can safely sort the instances by count
    arguments = sorted_by_count(arguments)
    coefficients = sorted_by_count(coefficients)

    return arguments, coefficients

def apply_nested_forward_ad(expr, dim):
    if isinstance(expr, Terminal):
        # A terminal needs no differentiation applied
        return expr
    elif not isinstance(expr, Derivative):
        # Apply AD recursively to children
        preops = expr.operands()
        postops = tuple(apply_nested_forward_ad(o, dim) for o in preops)
        # Reconstruct if necessary
        need_reconstruct = not (preops == postops) # FIXME: Is this efficient? O(n)?
        if need_reconstruct:
            expr = expr.reconstruct(*postops)
        return expr
    elif isinstance(expr, Grad):
        # Apply AD recursively to children
        f, = expr.operands()
        f = apply_nested_forward_ad(f, dim)
        # Apply Grad-specialized AD to expanded child
        return compute_grad_forward_ad(f, dim)
    elif isinstance(expr, VariableDerivative):
        # Apply AD recursively to children
        f, v = expr.operands()
        f = apply_nested_forward_ad(f, dim)
        # Apply Variable-specialized AD to expanded child
        return compute_variable_forward_ad(f, v, dim)
    elif isinstance(expr, CoefficientDerivative):
        # Apply AD recursively to children
        f, w, v, cd = expr.operands()
        f = apply_nested_forward_ad(f, dim)
        # Apply Coefficient-specialized AD to expanded child
        return compute_coefficient_forward_ad(f, w, v, cd, dim)
    else:
        error("Unknown type.")

def compute_form_with_arity(form, arity, arguments=None):
    """Compute parts of form of given arity."""

    # Extract all arguments in form
    if arguments is None:
        arguments = form.arguments()

    parts = [arg.part() for arg in arguments]
    if set(parts) - {None}:
        error("compute_form_with_arity cannot handle parts.")

    if len(arguments) < arity:
        warning("Form has no parts with arity %d." % arity)
        return 0*form

    # Assuming that the form is not a sum of terms
    # that depend on different arguments, e.g. (u+v)*dx
    # would result in just v*dx. But that doesn't make
    # any sense anyway.
    sub_arguments = set(arguments[:arity])
    pe = PartExtracter(sub_arguments)

    def _transform(e):
        e, provides = pe.visit(e)
        if provides == sub_arguments:
            return e
        return Zero()
    return map_integrands(_transform, form)

    def equals_func_with_collision_measuring(self, other):
        # Call equals
        equal = equals_func(self, other)

        # Get properties
        st = type(self)
        ot = type(other)
        sn = st.__name__
        on = ot.__name__
        sh = hash(self)
        oh = hash(other)
        key = (sn, on)

        # If hashes are the same but objects are not equal, we have a
        # collision
        hash_total[key] += 1
        if sh == oh and not equal:
            hash_collisions[key] += 1
        elif sh != oh and equal:
            error("Equal objects must always have the same hash! Objects are:\n{0}\n{1}".format(self, other))
        elif sh == oh and equal:
            hash_equals[key] += 1
        elif sh != oh and not equal:
            hash_notequals[key] += 1

        return equal

def compute_form_action(form, coefficient):
    """Compute the action of a form on a Coefficient.

    This works simply by replacing the last Argument
    with a Coefficient on the same function space (element).
    The form returned will thus have one Argument less
    and one additional Coefficient at the end if no
    Coefficient has been provided.
    """
    # TODO: Check whatever makes sense for coefficient

    # Extract all arguments
    arguments = form.arguments()

    parts = [arg.part() for arg in arguments]
    if set(parts) - {None}:
        error("compute_form_action cannot handle parts.")

    # Pick last argument (will be replaced)
    u = arguments[-1]

    fs = u.ufl_function_space()
    if coefficient is None:
        coefficient = Coefficient(fs)
    elif coefficient.ufl_function_space() != fs:
        debug("Computing action of form on a coefficient in a different function space.")
    return replace(form, {u: coefficient})