function [D, penalties] = constructFluxDiffOpWithClosures(scheme, g, order, schemeParams, opSet, BCs)
    [D, diffOp] = rv.diffops.constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
    [D, penalties] = rv.diffops.addClosuresToDiffOp(diffOp, D, BCs);
end