function diffOpsStruct = constructDiffOps(method, varargin)
    switch method
        case {'standard'}
            diffOpsStruct = diffOpsStandard(varargin{:});
        case {'bdf','backwards-difference-formula'}
            diffOpsStruct = diffOpsBdf(varargin{:});
        case {'ms','multi-stage'}
            diffOpsStruct = diffOpsMultiStage(varargin{:});
        case {'mg1','multi-grid1'}
            diffOpsStruct = diffOpsMultiGrid1(varargin{:});
        case {'mg','mg2','multi-grid2'} % Default multigrid diffops
            diffOpsStruct = diffOpsMultiGrid2(varargin{:});
    end
end

function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_flux(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', {penalties_scheme});
end

function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    %% DiffOps for solution vector
    [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_flux(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                 'D_flux', D_flux,...
                 'penalties_scheme', {penalties_scheme});
end

% TODO: Remove?
function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % DiffOp for unstabilized solution vector
    [D_unstable, closures] = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    D_unstable = rv.diffops.addClosuresToDiffOp(D_unstable, closures);
    %% DiffOps for residual viscosity
    D_t = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_t(v);
    % TODO: Use residual penalties as well here?
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_unstable', D_unstable,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', {penalties_scheme});
end

function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme, D_f] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % DiffOp for unstabilized solution vector
    D_coarse = coarserGridDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_flux(v);
    D_t = @(v) D_f(v);
    % TODO: Use residual penalties as well here?
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_coarse', D_coarse,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', {penalties_scheme});
end

function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % TODO: What orders to use here?
    D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    D_flux = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder , schemeParams, opSet, BCs);
    D_flux = @(v) -D_flux(v);
    D_t = @(v) D_coarse(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', {penalties_scheme});
end

% Note: Only works for equidistant grids.
function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
    m = g.m();
    lim = g.lim();
    m_c = (m-1)/2 + 1;
    switch g.D()
        case 1
            interpOps = sbp.InterpOpsOP(m_c, m, schemeOrder, schemeOrder);
            I = interpOps.Iu2v.good;
            g_c = grid.equidistant(m_c, lim{1});
            D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
            D_c = @(v) I*D_c(v(1:2:end));
        case 2
            interpOps_x = sbp.InterpOpsOP(m_c(1), m(1), schemeOrder, schemeOrder);
            I_x = interpOps_x.Iu2v.good;
            interpOps_y = sbp.InterpOpsOP(m_c(2), m(2), schemeOrder, schemeOrder);
            I_y = interpOps_y.Iu2v.good;
            I = kron(I_x,I_y);
            g_c = grid.equidistant(m_c, lim{1}, lim{2});
            ind = grid.funcToMatrix(g, 1:g.N());
            ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
            D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
            D_c = @(v) I*D_c(v(ind_c));
    end
end

function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs)
    m = g.m();
    lim = g.lim();
    m_c = (m-1)/2 + 1;
    switch g.D()
    case 1
        g_c = grid.equidistant(m_c, lim{1});
        [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
        D_c = rv.diffops.addClosuresToDiffOp(D_c, closures);
        x = g.x{1};
        x_c = x(1:2:end);
        D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
    case 2
        g_c = grid.equidistant(m_c, lim{1}, lim{2});
        ind = grid.funcToMatrix(g, 1:g.N());
        ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
        [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
        D_c = rv.diffops.addClosuresToDiffOp(D_c, closures);
        D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(v(ind_c))),'spline')',g.N(),1);
    end
end