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 {'mg2','multi-grid2'}
            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.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    D_flux = rv.diffops.constructFluxDiffOpWithClosures(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.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    [D_flux, ~] = rv.diffops.constructFluxDiffOp(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

function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % DiffOp for unstabilized solution vector
    D_unstable = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_t(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_unstable', D_unstable,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', penalties_scheme,...
                     'penalties_res', penalties_res);
end

function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % DiffOp for unstabilized solution vector
    D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    [D_flux, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_flux(v);
    D_t = @(v) D(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_coarse', D_coarse,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', penalties_scheme,...
                     'penalties_res', penalties_res);
end

function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(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.constructFluxDiffOpWithClosures(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

% TODO: Only works for equidistant grids
function D_coarse = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
    m = g.m();
    lim = g.lim();
    m_coarse = (m-1)/2 + 1;
    g_coarse = grid.equidistant(m_coarse, lim{1});
    D_coarse = rv.diffops.constructFluxDiffOpWithClosures(scheme, g_coarse, schemeOrder, schemeParams, opSet, BCs);
    x = g.x{1};
    x_coarse = x(1:2:end);
    % TODO: Fix interpolation
    D_coarse = @(v) interp1(x_coarse,D_coarse(v(1:2:end)),x,'spline');
end