classdef Burgers1D < scheme.Scheme
    properties
        grid % Physical grid
        order % Order accuracy for the approximation
        
        params

        D % Non-stabalized scheme operator
        H % Discrete norm
        Hi % Norm inverse
        e_l
        e_r
        d_l
        d_r
    end

    methods
        function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params)
            assert(grid.D == 1);
            assert(grid.size() == length(params.eps));
            m = grid.size();
            lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
            switch operator_type
                case 'narrow'
                    ops = sbp.D4Variable(m, lim, order);
                    D1 = ops.D1;
                    D2 = ops.D2;
                    if (strcmp(dissipation,'on'))
                        DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
                    end
                    d_l = ops.d1_l';
                    d_r = ops.d1_r';
                case 'upwind-'
                    ops = sbp.D1Upwind(m, lim, order);
                    D1 = (ops.Dp + ops.Dm)/2;
                    D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
                    if (strcmp(dissipation,'on'))
                        DissipationOp = (ops.Dp-ops.Dm)/2;
                    end
                    d_l = ops.e_l'*ops.Dm;
                    d_r = ops.e_r'*ops.Dm;
                case 'upwind+'
                    ops = sbp.D1Upwind(m, lim, order);
                    D1 = (ops.Dp + ops.Dm)/2;
                    D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp;
                    if (strcmp(dissipation,'on'))
                        DissipationOp = (ops.Dp-ops.Dm)/2;
                    end
                    d_l = ops.e_l'*ops.Dp;
                    d_r = ops.e_r'*ops.Dp;
                case 'upwind+-'
                    ops = sbp.D1Upwind(m, lim, order);
                    D1 = (ops.Dp + ops.Dm)/2;
                    D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2;
                    if (strcmp(dissipation,'on'))
                        DissipationOp = (ops.Dp-ops.Dm)/2;
                    end
                    d_l = ops.e_l'*D1;
                    d_r = ops.e_r'*D1;
                otherwise
                    error('Other operator types not yet supported', operator_type);
            end

            switch pde_form
                case 'skew-symmetric'
                    if (strcmp(dissipation,'on'))
                        D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
                    else
                        D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v;
                    end
                case 'conservative'
                    if (strcmp(dissipation,'on'))
                        D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
                    else
                        D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
                    end
                otherwise
                    error('Not supported', pde_form);
            end

            obj.grid = grid;
            obj.order = order;
            obj.params = params;

            obj.D = D;
            obj.H =  ops.H;
            obj.Hi = ops.HI;
            obj.e_l = ops.e_l;
            obj.e_r = ops.e_r;
            obj.d_l =  d_l;
            obj.d_r =  d_r;
        end

        % Closure functions return the operators applied to the own doamin to close the boundary
        % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
        %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
        %       type                is a string specifying the type of boundary condition if there are several.
        %       data                is a function returning the data that should be applied at the boundary.
        function [closure, penalty] = boundary_condition(obj,boundary,type,data)
            default_arg('type','robin');
            default_arg('data',0);
            [e, d, i_b, s] = obj.get_boundary_ops(boundary);
            switch type
                % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
                case {'R','robin'}
                    p = s*obj.Hi*e;
                    closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - ((obj.params.eps(i_b) + viscosity(i_b))*d*v));
                    switch class(data)
                        case 'double'
                            penalty = s*p*data;
                        case 'function_handle'
                            penalty = @(t) s*p*data(t);
                        otherwise
                            error('Wierd data argument!')
                    end
                otherwise
                    error('No such boundary condition: type = %s',type);
            end
        end

        % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
        % The right boundary is considered the positive boundary
        function [e, d, i_b, s] = get_boundary_ops(obj,boundary)
            switch boundary
                case 'l'
                    e = obj.e_l;
                    d = obj.d_l;
                    i_b = 1;
                    s = -1;
                case 'r'
                    e = obj.e_r;
                    d = obj.d_r;
                    i_b = length(e);
                    s = 1;
                otherwise
                    error('No such boundary: boundary = %s',boundary);
            end
        end

        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
            error('An interface function does not exist yet');
        end

        function N = size(obj)
            N = obj.grid.m;
        end
    end
end