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.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
%% Multigrid functions:
% TODO: Implement properly.
function v_c = coarse_v(g, g_c, v)
    V = grid.funcToMatrix(g,v);
    v_c = reshape(V(1:2:end,1:2:end)',g_c.N(),1);
end

% Note: Only works for equidistant grids.
% TODO: Change from using matlabs interp to using proper interpolation operators..
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
        g_c = grid.equidistant(m_c, lim{1});
        D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
        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});
        D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
        D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1);
    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});
        [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(coarse_v(g,g_c,v))),'spline')',g.N(),1);
    end
end