def ic_expr(t, x, fields):
        from hedge.optemplate import CFunction
        from pymbolic.primitives import IfPositive
        from pytools.obj_array import make_obj_array

        tanh = CFunction("tanh")
        sin = CFunction("sin")

        rho = 1
        u0 = 0.05
        w = 0.05
        delta = 0.05

        from hedge.optemplate.primitives import make_common_subexpression as cse
        u = cse(make_obj_array([
            IfPositive(x[1]-1/2,
                u0*tanh(4*(3/4-x[1])/w),
                u0*tanh(4*(x[1]-1/4)/w)),
            u0*delta*sin(2*np.pi*(x[0]+1/4))]),
            "u")

        return make_obj_array([
            op.method.f_equilibrium(rho, alpha, u)
            for alpha in range(len(op.method))
            ])

def _small_mat_inverse(mat):
    m, n = mat.shape
    if m != n:
        raise ValueError("inverses only make sense for square matrices")

    if m == 1:
        return make_obj_array([1/mat[0, 0]])
    elif m == 2:
        (a, b), (c, d) = mat
        return 1/(a*d-b*c) * make_obj_array([
            [d, -b],
            [-c, a],
            ])
    else:
        raise NotImplementedError(
                "inverse formula for %dx%d matrices" % (m, n))

 def map_nodes(self, expr):
     discr = self.discr_dict[expr.where]
     from pytools.obj_array import make_obj_array
     return MultiVector(
             make_obj_array([
                 prim.NodeCoordinateComponent(i, expr.where)
                 for i in range(discr.ambient_dim)]))

def _small_mat_eigenvalues(mat):
    m, n = mat.shape
    if m != n:
        raise ValueError("eigenvalues only make sense for square matrices")

    if m == 1:
        return make_obj_array([mat[0, 0]])
    elif m == 2:
        (a, b), (c, d) = mat
        return make_obj_array([
                -(sqrt(d**2-2*a*d+4*b*c+a**2)-d-a)/2,
                 (sqrt(d**2-2*a*d+4*b*c+a**2)+d+a)/2
                ])
    else:
        raise NotImplementedError(
                "eigenvalue formula for %dx%d matrices" % (m, n))

def xyz_to_tangential(xyz_vec, where=None):
    ambient_dim = len(xyz_vec)
    tonb = tangential_onb(ambient_dim, where=where)
    return make_obj_array([
        np.dot(tonb[:, i], xyz_vec)
        for i in range(ambient_dim - 1)
        ])

    def full_output_zeros(self):
        """This includes QBX and non-QBX targets."""

        from pytools.obj_array import make_obj_array
        return make_obj_array([
                np.zeros(self.tree.ntargets, self.dtype)
                for k in self.outputs])

    def representation(self, unknown, i_domain):
        """
        :return: a symbolic expression for the representation of the PDE solution
            in domain number *i_domain*.
        """
        unk = self._structured_unknown(unknown, with_l2_weights=False)

        result = []

        for field_kind in self.field_kinds:
            if not self.is_field_present(field_kind):
                continue

            field_result = 0
            for i_interface, (i_domain_outer, i_domain_inner, interface_id) in (
                    enumerate(self.interfaces)):
                if i_domain_outer == i_domain:
                    side = self.side_out
                elif i_domain_inner == i_domain:
                    side = self.side_in
                else:
                    continue

                my_unk = unk[side, field_kind, i_interface]
                if my_unk:
                    field_result += sym.S(
                            self.kernel,
                            my_unk,
                            source=interface_id,
                            k=self.domain_K_exprs[i_domain])

            result.append(field_result)

        from pytools.obj_array import make_obj_array
        return make_obj_array(result)

    def add_div_bcs(self, tgt, bc_getter, dirichlet_tags, neumann_tags,
            stab_term, adjust_flux, flux_v, flux_arg_int,
            grad_flux_arg_count):
        from pytools.obj_array import make_obj_array, join_fields
        n_times = tgt.normal_times_flux

        def unwrap_cse(expr):
            from pymbolic.primitives import CommonSubexpression
            if isinstance(expr, CommonSubexpression):
                return expr.child
            else:
                return expr

        for tag in dirichlet_tags:
            dir_bc_w = join_fields(
                    [0]*grad_flux_arg_count,
                    [bc_getter(tag, unwrap_cse(vol_expr)) for vol_expr in
                        flux_arg_int[grad_flux_arg_count:]])
            tgt.add_boundary_flux(
                    adjust_flux(n_times(flux_v.int-stab_term)),
                    flux_arg_int, dir_bc_w, tag)

        loc_bc_vec = make_obj_array([0]*len(flux_arg_int))

        for tag in neumann_tags:
            neu_bc_w = join_fields(
                    NeumannBCGenerator(tag, bc_getter(tag, None))(tgt.operand),
                    [0]*len(flux_arg_int))

            tgt.add_boundary_flux(
                    adjust_flux(n_times(flux_v.ext)),
                    loc_bc_vec, neu_bc_w, tag)

    def __call__(self, arg):
        if isinstance(arg, int) and arg == 0:
            return 0
        from pytools.obj_array import make_obj_array
        from hedge.optemplate.primitives import OperatorBinding

        return make_obj_array([OperatorBinding(FluxOperator(f), arg) for f in self.fluxes])

    def matvec(self, x):
        if isinstance(x, np.ndarray):
            x = cl.array.to_device(self.queue, x)
            out_host = True
        else:
            out_host = False

        do_split = len(self.starts_and_ends) > 1
        from pytools.obj_array import make_obj_array

        if do_split:
            x = make_obj_array(
                    [x[start:end] for start, end in self.starts_and_ends])

        args = self.extra_args.copy()
        args[self.arg_name] = x
        result = self.bound_expr(self.queue, **args)

        if do_split:
            # re-join what was split
            joined_result = cl.array.empty(self.queue, self.total_dofs,
                    np.complex128)  # FIXME
            for res_i, (start, end) in zip(result, self.starts_and_ends):
                joined_result[start:end] = res_i
            result = joined_result

        if out_host:
            result = result.get()

        return result

    def map_int_g_ds(self, expr):
        dsource = self.rec(expr.dsource)

        ambient_dim = self.ambient_dim

        from sumpy.kernel import KernelDimensionSetter
        kernel = _insert_source_derivative_into_kernel(
                KernelDimensionSetter(ambient_dim)(expr.kernel))

        from pytools.obj_array import make_obj_array
        nabla = MultiVector(make_obj_array(
            [prim.NablaComponent(axis, None)
                for axis in range(ambient_dim)]))

        kernel_arguments = dict(
                (name, self.rec(arg_expr))
                for name, arg_expr in expr.kernel_arguments.items()
                )

        def add_dir_vec_to_kernel_args(coeff):
            result = kernel_arguments.copy()
            result[_DIR_VEC_NAME] = _get_dir_vec(coeff, ambient_dim)
            return result

        rec_operand = prim.cse(self.rec(expr.density))
        return (dsource*nabla).map(
                lambda coeff: prim.IntG(
                    kernel,
                    rec_operand, expr.qbx_forced_limit, expr.source, expr.target,
                    kernel_arguments=add_dir_vec_to_kernel_args(coeff)))

def test_area_query_balls_outside_bbox(ctx_getter, dims, do_plot=False):
    """
    The input to the area query includes balls whose centers are not within
    the tree bounding box.
    """
    ctx = ctx_getter()
    queue = cl.CommandQueue(ctx)

    nparticles = 10**4
    dtype = np.float64

    particles = make_normal_particle_array(queue, nparticles, dims, dtype)

    if do_plot:
        import matplotlib.pyplot as pt
        pt.plot(particles[0].get(), particles[1].get(), "x")

    from boxtree import TreeBuilder
    tb = TreeBuilder(ctx)

    queue.finish()
    tree, _ = tb(queue, particles, max_particles_in_box=30, debug=True)

    nballs = 10**4
    from pyopencl.clrandom import PhiloxGenerator
    rng = PhiloxGenerator(ctx, seed=13)
    bbox_min = tree.bounding_box[0].min()
    bbox_max = tree.bounding_box[1].max()
    from pytools.obj_array import make_obj_array
    ball_centers = make_obj_array([
        rng.uniform(queue, nballs, dtype=dtype, a=bbox_min-1, b=bbox_max+1)
        for i in range(dims)])
    ball_radii = cl.array.empty(queue, nballs, dtype).fill(0.1)

    run_area_query_test(ctx, queue, tree, ball_centers, ball_radii)

def curl_S_volume(kernel, arg):
    from pytools import levi_civita
    from pytools.obj_array import make_obj_array

    return make_obj_array([
        sum(
            levi_civita((l, m, n)) * IntGdTarget(kernel, arg[n], m)
            for m in range(3) for n in range(3))
        for l in range(3)])

def curl(vec):
    from pytools import levi_civita
    from pytools.obj_array import make_obj_array

    return make_obj_array([
        sum(
            levi_civita((l, m, n)) * dd_axis(m, 3, vec[n])
            for m in range(3) for n in range(3))
        for l in range(3)])

def nodes(ambient_dim, where=None):
    """Return a :class:`pymbolic.geometric_algebra.MultiVector` of node
    locations.
    """

    return MultiVector(
            make_obj_array([
                NodeCoordinateComponent(i, where)
                for i in range(ambient_dim)]))

def n_cross(vec, which=None):
    nrm = normal(3, which)

    from pytools import levi_civita
    from pytools.obj_array import make_obj_array
    return make_obj_array([
        sum(
            levi_civita((i, j, k)) * nrm[j] * vec[k]
            for j in range(3) for k in range(3))
        for i in range(3)])

def cross(vec_a, vec_b):
    assert len(vec_a) == len(vec_b) == 3

    from pytools import levi_civita
    from pytools.obj_array import make_obj_array
    return make_obj_array([
        sum(
            levi_civita((i, j, k)) * vec_a[j] * vec_b[k]
            for j in range(3) for k in range(3))
        for i in range(3)])

        def make_artificial_diffusion():
            if self.artificial_viscosity_mode not in ["diffusion"]:
                return 0

            dq = self.grad_of_state()

            return make_obj_array([
                self.div(
                    to_vol_quad(self.sensor())*to_vol_quad(dq[i]),
                    to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) 
                for i in range(dq.shape[0])])

    def collision_update(self, f_bar):
        from hedge.optemplate.primitives import make_common_subexpression as cse
        rho = cse(self.rho(f_bar), "rho")
        rho_u = self.rho_u(f_bar)
        u = cse(rho_u/rho, "u")

        f_eq_func = self.method.f_equilibrium
        f_eq = make_obj_array([
            f_eq_func(rho, alpha, u) for alpha in range(len(self.method))])

        return f_bar - 1/(self.tau+1/2)*(f_bar - f_eq)

    def apply_diff(self, nabla, operand):
        from pytools.obj_array import make_obj_array, is_obj_array
        if is_obj_array(operand):
            if len(operand) != self.dimensions:
                raise ValueError("operand of apply_diff must have %d dimensions"
                        % self.dimensions)

            return sum(nabla[i](operand[i]) for i in range(self.dimensions))
        else:
            return make_obj_array(
                [nabla[i](operand) for i in range(self.dimensions)])