function [D_f, closures, penalties] = constructFluxDiffOps(scheme, g, order, schemeParams, opSet, BCs)
    diffOp = scheme(g, order, schemeParams{:}, opSet);
    if ~isa(diffOp.D, 'function_handle')
        D_f = @(v) diffOp.D*v;
    else
        D_f = diffOp.D;
	end
    penalties = cell(size(BCs));
    closures = cell(size(BCs));
    for i = 1:size(BCs,1)
        for j = 1:size(BCs,2)
            [closures{i,j}, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
        end
    end
end