def METHOD_NAME(refinement=1, quadratic=False, no_exports=True):
    """
    Computes pressure gradient in a setting where bathymetry, mesh surface
    elevation, and pressure are analytical, non-trivial functions.
    """
    print_output(' ---- running refinement {:}'.format(refinement))
    # create mesh
    rho_0 = 1000.0
    physical_constants['rho0'] = rho_0
    delta_x = 120e3/refinement
    lx = 360e3
    ly = 360e3
    nx = int(lx/delta_x)
    ny = int(ly/delta_x)
    mesh2d = RectangleMesh(nx, ny, lx, ly)
    layers = 3*refinement
    # bathymetry
    P1_2d = get_functionspace(mesh2d, 'CG', 1)
    bathymetry_2d = Function(P1_2d, name='Bathymetry')
    xy = SpatialCoordinate(mesh2d)
    depth = 3600.
    bath_expr = 0.5*(depth + depth)*(1 - 0.6*tanh(4*(xy[1]-ly/2)/ly)*sin(1.5*xy[0]/ly+0.2))
    bathymetry_2d.project(bath_expr)
    mesh = extrude_mesh_sigma(mesh2d, layers, bathymetry_2d)
    bnd_len = compute_boundary_length(mesh2d)
    mesh2d.boundary_len = bnd_len
    mesh.boundary_len = bnd_len
    # make function spaces and fields
    p1 = get_functionspace(mesh, 'CG', 1)
    if quadratic:
        # NOTE for 3rd order convergence both the scalar and grad must be p2
        fs_pg = get_functionspace(mesh, 'DG', 2, 'CG', 2, vector=True, dim=2)
        fs_scalar = get_functionspace(mesh, 'DG', 2, vfamily='CG', vdegree=2)
    else:
        # the default function spaces in Thetis
        fs_pg = get_functionspace(mesh, 'DG', 1, 'CG', 2, vector=True, dim=2)
        fs_scalar = get_functionspace(mesh, 'DG', 1, vfamily='CG', vdegree=2)
    density_3d = Function(fs_scalar, name='density')
    baroc_head_3d = Function(fs_scalar, name='baroclinic head')
    int_pg_3d = Function(fs_pg, name='pressure gradient')
    elev_3d = Function(p1, name='elevation')
    bathymetry_3d = Function(p1, name='elevation')
    ExpandFunctionTo3d(bathymetry_2d, bathymetry_3d).solve()
    # analytic expressions
    xyz = SpatialCoordinate(mesh)
    elev_expr = 2000.0*sin(0.3 + 1.5*xyz[1]/ly)*cos(2*xyz[0]/lx)
    density_expr = sin((xyz[1] - 0.3)/lx)*cos(2*xyz[2]/depth)*cos(2*xyz[0]/lx)
    baroc_head_expr = -depth*sin((2*xyz[2] - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos(2*xyz[0]/lx))/depth)*sin((xyz[1] - 0.3)/lx)*cos(2*xyz[0]/lx)/2
    baroc_head_expr_dx = (depth*sin((2*xyz[2] - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos(2*xyz[0]/lx))/depth) - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos((2*xyz[2] - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos(2*xyz[0]/lx))/depth)*cos(2*xyz[0]/lx))*sin(2*xyz[0]/lx)*sin((xyz[1] - 0.3)/lx)/lx
    baroc_head_expr_dy = (-depth*ly*sin((2*xyz[2] - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos(2*xyz[0]/lx))/depth)*cos((xyz[1] - 0.3)/lx) + 6000.0*lx*sin((xyz[1] - 0.3)/lx)*cos((2*xyz[2] - 4000.0*sin((0.3*ly + 1.5*xyz[1])/ly)*cos(2*xyz[0]/lx))/depth)*cos(2*xyz[0]/lx)*cos((0.3*ly + 1.5*xyz[1])/ly))*cos(2*xyz[0]/lx)/(2*lx*ly)
    # deform mesh by elevation
    elev_3d.project(elev_expr)
    z_ref = mesh.coordinates.dat.data[:, 2]
    bath = bathymetry_3d.dat.data[:]
    eta = elev_3d.dat.data[:]
    new_z = eta*(z_ref + bath)/bath + z_ref
    mesh.coordinates.dat.data[:, 2] = new_z
    if not no_exports:
        out_density = File('density.pvd')
        out_bhead = File('baroc_head.pvd')
        out_pg = File('int_pg.pvd')
    # project initial scalar
    density_3d.project(density_expr)
    baroc_head_3d.project(baroc_head_expr)
    # compute int_pg
    fields = FieldDict()
    fields.baroc_head_3d = baroc_head_3d
    fields.int_pg_3d = int_pg_3d
    bnd_functions = {}
    int_pg_solver = InternalPressureGradientCalculator(
        fields,
        bathymetry_3d,
        bnd_functions,
        solver_parameters=None)
    int_pg_solver.solve()
    if not no_exports:
        out_density.write(density_3d)
        out_bhead.write(baroc_head_3d)
        out_pg.write(int_pg_3d)
    g_grav = physical_constants['g_grav']
    ana_sol_expr = g_grav*as_vector((
        baroc_head_expr_dx,
        baroc_head_expr_dy,))
    volume = comp_volume_3d(mesh)
    l2_err = errornorm(ana_sol_expr, int_pg_3d, degree_rise=2)/numpy.sqrt(volume)
    print_output('L2 error {:}'.format(l2_err))
    if not no_exports:
        out_density.write(density_3d)
        out_bhead.write(baroc_head_3d)
        out_pg.write(int_pg_3d.project(ana_sol_expr))
    return l2_err