def METHOD_NAME(grid_intervals, output_type='l2', vector_valued_output=False):
    domain = RectDomain(([-1, -1], [1, 1]))
    indicator_domain = ExpressionFunction(
        '(-2/3. <= x[0] <= -1/3.) * (-2/3. <= x[1] <= -1/3.) * 1. \
       + (-2/3. <= x[0] <= -1/3.) *  (1/3. <= x[1] <=  2/3.) * 1.',
        dim_domain=2)
    rest_of_domain = ConstantFunction(1, 2) - indicator_domain
    f = ExpressionFunction('0.5*pi*pi*cos(0.5*pi*x[0])*cos(0.5*pi*x[1])', dim_domain=2)
    parameters = {'diffusion': 2}
    thetas = [
        ExpressionParameterFunctional(
            '1.1 + sin(diffusion[0])*diffusion[1]',
            parameters,
            derivative_expressions={'diffusion': ['cos(diffusion[0])*diffusion[1]',
                                                  'sin(diffusion[0])']},
        ),
        ExpressionParameterFunctional(
            '1.1 + sin(diffusion[1])',
            parameters,
            derivative_expressions={'diffusion': ['0', 'cos(diffusion[1])']},
        ),
    ]
    diffusion = LincombFunction([rest_of_domain, indicator_domain], thetas)
    theta_J = ExpressionParameterFunctional('1 + 1/5 * diffusion[0] + 1/5 * diffusion[1]', parameters,
                                            derivative_expressions={'diffusion': ['1/5', '1/5']})
    if vector_valued_output:
        problem = StationaryProblem(
            domain, f, diffusion, outputs=[(output_type, f * theta_J), (output_type, f * 0.5 * theta_J)]
        )
    else:
        problem = StationaryProblem(domain, f, diffusion, outputs=[(output_type, f * theta_J)])
    print('Discretize ...')
    mu_bar = problem.parameters.parse([np.pi/2, np.pi/2])
    fom, _ = discretize_stationary_cg(problem, diameter=1. / grid_intervals, mu_energy_product=mu_bar)
    return fom, mu_bar