def METHOD_NAME(r, element):
    """Apply an appropriate pullback to something in physical space
    :arg r: An expression wrapped in ReferenceValue.
    :arg element: The element this expression lives in.
    :returns: a pulled back expression."""
    mapping = element.mapping()
    if r.ufl_shape != element.reference_value_shape():
        raise ValueError(
            f"Expecting reference space expression with shape '{element.reference_value_shape()}', "
            f"got '{r.ufl_shape}'")
    if mapping in {"physical", "identity",
                   "contravariant Piola", "covariant Piola",
                   "double contravariant Piola", "double covariant Piola",
                   "L2 Piola", "custom"}:
        # Base case in recursion through elements. If the element
        # advertises a mapping we know how to handle, do that
        # directly.
        f = apply_known_single_pullback(r, element)
        if f.ufl_shape != element.value_shape():
            raise ValueError(f"Expecting pulled back expression with shape '{element.value_shape()}', "
                             f"got '{f.ufl_shape}'")
        return f
    elif mapping in {"symmetries", "undefined"}:
        # Need to pull back each unique piece of the reference space thing
        gsh = element.value_shape()
        rsh = r.ufl_shape
        if mapping == "symmetries":
            subelem = element.sub_elements()[0]
            fcm = element.flattened_sub_element_mapping()
            offsets = (product(subelem.reference_value_shape()) * i for i in fcm)
            elements = repeat(subelem)
        else:
            elements = sub_elements_with_mappings(element)
            # Python >= 3.8 has an initial keyword argument to
            # accumulate, but 3.7 does not.
            offsets = chain([0],
                            accumulate(product(e.reference_value_shape())
                                       for e in elements))
        rflat = as_vector([r[idx] for idx in numpy.ndindex(rsh)])
        g_components = []
        # For each unique piece in reference space, apply the appropriate pullback
        for offset, subelem in zip(offsets, elements):
            sub_rsh = subelem.reference_value_shape()
            rm = product(sub_rsh)
            rsub = [rflat[offset + i] for i in range(rm)]
            rsub = as_tensor(numpy.asarray(rsub).reshape(sub_rsh))
            rmapped = METHOD_NAME(rsub, subelem)
            # Flatten into the pulled back expression for the whole thing
            g_components.extend([rmapped[idx]
                                 for idx in numpy.ndindex(rmapped.ufl_shape)])
        # And reshape appropriately
        f = as_tensor(numpy.asarray(g_components).reshape(gsh))
        if f.ufl_shape != element.value_shape():
            raise ValueError(f"Expecting pulled back expression with shape '{element.value_shape()}', "
                             f"got '{f.ufl_shape}'")
        return f
    else:
        raise ValueError(f"Unsupported mapping type: {mapping}")