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

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

        D2_visc % Artificial viscosity operator
    end

    methods
        function obj = AdvectionRV1D(grid, operator_type, order, waveSpeed)
            m = grid.size();
            lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
            switch operator_type
                case 'upwind+'
                    ops = sbp.D1Upwind(m, lim, order);
                    D1 = (ops.Dp + ops.Dm)/2;
                    B = ops.e_r*ops.e_r' - ops.e_l*ops.e_l';
                    obj.D2_visc = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp);
                    % max(abs()) or just abs()?
                    DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2;
                otherwise
                    error('Other operator types not yet supported', operator_type);
            end
            
            obj.D = -D1 + DissipationOp;
            
            obj.grid = grid;
            obj.order = order;

            obj.H =  ops.H;
            obj.Hi = ops.HI;
            obj.e_l = ops.e_l;
            obj.e_r = ops.e_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, s] = obj.get_boundary_ops(boundary);
            switch type
                case {'D', 'dirichlet'}
                    p = s*obj.Hi*e;
                    closure = p*e';
                otherwise
                    error('No such boundary condition: type = %s',type);
            end
            switch class(data)
                case 'double'
                    penalty = s*p*data;
                case 'function_handle'
                    penalty = @(t) s*p*data(t);
                otherwise
                    error('Wierd data argument!')
            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, s] = get_boundary_ops(obj,boundary)
            switch boundary
                case 'l'
                    e = obj.e_l;
                    s = -1;
                case 'r'
                    e = obj.e_r;
                    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