def METHOD_NAME(target, eval_data, order=None,
                            ls=None, nls_options=None):
    """
    Project scalar data to a scalar `target` field variable using the
    :math:`L^2` dot product.
    Parameters
    ----------
    target : FieldVariable instance
        The target variable.
    eval_data : callable or array
        Either a material-like function `eval_data()`, or an array of values in
        quadrature points that has to be reshapable to the shape required by
        `order`.
    order : int, optional
        The quadrature order. If not given, it is set to
        `2 * target.field.approx_order`.
    """
    if order is None:
       order = 2 * target.field.approx_order
    integral = Integral('i', order=order)
    un = FieldVariable('u', 'unknown', target.field)
    v = FieldVariable('v', 'test', un.field, primary_var_name=un.name)
    lhs = Term.new('dw_dot(v, %s)' % un.name, integral,
                   un.field.region, v=v, **{un.name : un})
    def _eval_data(ts, coors, mode, **kwargs):
        if mode == 'qp':
            if callable(eval_data):
                val = eval_data(ts, coors, mode, **kwargs)
            else:
                val = eval_data.reshape((coors.shape[0], 1, 1))
            return {'val' : val}
    m = Material('m', function=_eval_data)
    rhs = Term.new('dw_volume_lvf(m.val, v)', integral, un.field.region,
                   m=m, v=v)
    eq = Equation('projection', lhs - rhs)
    eqs = Equations([eq])
    if ls is None:
        ls = ScipyDirect({})
    if nls_options is None:
        nls_options = {}
    nls_status = IndexedStruct()
    nls = Newton(nls_options, lin_solver=ls, status=nls_status)
    pb = Problem('aux', equations=eqs)
    pb.set_solver(nls)
    # This sets the un variable with the projection solution.
    pb.solve(save_results=False)
    # Copy the projection solution to target.
    target.set_data(un())
    if nls_status.condition != 0:
        output('L2 projection: solver did not converge!')