Very-high-level interface

geom = GeometryCollection(
    "vol": meshmode.discretization.Discretization(...),
    "bdry": pytential.QBXLayerPotentialSource(...),
    )
f = sym.Variable("f", "vol")
bc = sym.Variable("f", "bdry")
sigma = sym.Variable("sigma", "bdry")

kernel = sumpy.LaplaceKernel(3)

u_sym = (
    sym.IntG(kernel, f, source="vol", target="vol")
    + sym.D(kernel,sym.IterativeInverse(
        0.5*sigma + sym.D(kernel, sigma),
        rhs=bc - sym.IntG(kernel, f, source="vol", target="bdry")),
        variable_name="sigma"))

u = bind(geom, u_sym)(queue, f=..., bc=....)