classdef Burgers2d < scheme.Scheme
    properties
        grid % Physical grid
        order % Order accuracy for the approximation
        
        D % Non-stabilized scheme operator
        H % Discrete norm
        H_x, H_y % Norms in the x and y directions
        Hi % Kroneckered norm inverse
        % Boundary operators
        e_w, e_e, e_s, e_n
    end

    methods
        function obj = Burgers2d(g, order, pde_form, fluxSplitting, opSet)
            default_arg('opSet',@sbp.D2Standard);
            default_arg('fluxSplitting',{@(v)max(abs(v)),@(v)max(abs(v))});
            assertType(g, 'grid.Cartesian');

            m = g.size();
            m_x = m(1);
            m_y = m(2);
            m_tot = g.N();

            xlim = {g.x{1}(1), g.x{1}(end)};
            ylim = {g.x{2}(1), g.x{2}(end)};
            obj.grid = g;

            % Operator sets
            ops_x = opSet(m_x, xlim, order);
            ops_y = opSet(m_y, ylim, order);
            Ix = speye(m_x);
            Iy = speye(m_y);

            % Norms
            Hx = ops_x.H;
            Hy = ops_y.H;
            Hxi = ops_x.HI;
            Hyi = ops_y.HI;

            obj.H_x = Hx;
            obj.H_y = Hy;
            obj.H = kron(Hx,Hy);
            obj.Hi = kron(Hxi,Hyi);

            % Derivatives
            if (isequal(opSet,@sbp.D1Upwind))
                Dx = kron((ops_x.Dp + ops_x.Dm)/2,Iy);
                Dy = kron(Ix,(ops_y.Dp + ops_y.Dm)/2);
                DissOpx = kron((ops_x.Dp - ops_x.Dm)/2,Iy);
                DissOpy = kron(Ix,(ops_y.Dp - ops_y.Dm)/2);
                D1 = Dx + Dy;   
                switch pde_form
                    case 'skew-symmetric'
                        D = -1/3*D1;
                        switch length(fluxSplitting)
                            case 1
                                DissOp = DissOpx + DissOpy;
                                obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOp)*v;
                            case 2
                                obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v;
                        end
                    case 'conservative'
                        D = -1/2*D1;
                        switch length(fluxSplitting)
                            case 1
                                DissOp = DissOpx + DissOpy;
                                % TODO: Check if we can use fluxSplitting{1} here instead
                                obj.D = @(v) D*(v.*v) + max(abs(v))*DissOp*v;
                            case 2
                                obj.D = @(v) D*(v.*v) + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v;
                        end
                        
                end
            else
                Dx = kron(ops_x.D1,Iy);
                Dy = kron(Ix,ops_y.D1);
                D1 = Dx + Dy;
                switch pde_form
                    case 'skew-symmetric'
                        D = -1/3*D1;
                        obj.D = @(v) D*(v.*v) + spdiags(v,0,m_tot,m_tot)*D*v;
                    case 'conservative'
                        D = -1/2*D1;
                        obj.D = @(v) D*(v.*v);
                end
            end

            % Boundary operators
            obj.e_w = kr(ops_x.e_l, Iy);
            obj.e_e = kr(ops_x.e_r, Iy);
            obj.e_s = kr(Ix, ops_y.e_l);
            obj.e_n = kr(Ix, ops_y.e_r);

            obj.order = order;
        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.
        function [closure, penalty] = boundary_condition(obj,boundary,type)
            default_arg('type','dirichlet');
            s = obj.getBoundarySign(boundary);
            e = obj.getBoundaryOperator('e', boundary);
            indices = obj.getBoundaryIndices(boundary);
            H_1d = obj.getOneDirectionalNorm(boundary);
            switch type
                case {'D', 'd', 'dirichlet', 'Dirichlet'}
                    penalty_parameter = 1/3;
                    Tau = s*penalty_parameter*obj.Hi*e*H_1d/2;
                    m = obj.grid.m;
                    tau = @(v) Tau*spdiags((v(indices)-s*abs(v(indices))),0,m(1),m(2));
                    closure = @(v) Tau*((v(indices)-s*abs(v(indices))).*v(indices));
                    penalty = @(v) -tau(v);
                otherwise
                    error('No such boundary condition: type = %s',type);
            end


        end

        % Returns the boundary sign. The right boundary is considered the positive boundary
        % boundary -- string
        function s = getBoundarySign(obj, boundary)
            assertIsMember(boundary, {'w', 'e', 's', 'n'})
            switch boundary
                case {'e','n'}
                    s = 1;
                case {'w','s'}
                    s = -1;
            end
        end

        % Returns the boundary operator op for the boundary specified by the string boundary.
        % op        -- string
        % boundary  -- string
        function o = getBoundaryOperator(obj, op, boundary)
            assertIsMember(op, {'e'})
            assertIsMember(boundary, {'w', 'e', 's', 'n'})

            o = obj.([op, '_', boundary]);
        end

        % Returns square boundary quadrature matrix, of dimension
        % corresponding to the number of boundary points
        %
        % boundary -- string
        function H_b = getBoundaryQuadrature(obj, boundary)
            assertIsMember(boundary, {'w', 'e', 's', 'n'})
            H_b = obj.(['H_', boundary]);
        end

        % Returns square boundary quadrature matrix, of dimension
        % corresponding to the number of boundary points
        %
        % boundary -- string
        function H_1d = getOneDirectionalNorm(obj, boundary)
            assertIsMember(boundary, {'w', 'e', 's', 'n'})
            switch boundary
                case {'w','e'}
                    H_1d = obj.H_y;
                case {'s','n'}
                    H_1d = obj.H_x;
            end
        end

        % Returns the indices of the boundary points in the grid matrix
        % boundary -- string
        function I = getBoundaryIndices(obj, boundary)
            assertIsMember(boundary, {'w', 'e', 's', 'n'})

            ind = grid.funcToMatrix(obj.grid, 1:prod(obj.grid.m));
            switch boundary
                case 'w'
                    I = ind(1,:);
                case 'e'
                    I = ind(end,:);
                case 's'
                    I = ind(:,1)';
                case 'n'
                    I = ind(:,end)';
            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