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');
            [e, H_b, index, s] = obj.get_boundary_ops(boundary);
            switch type
                % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
                % with +- at left/right boundaries in each coordinate direction
                case {'D', 'd', 'dirichlet', 'Dirichlet'}

                    magnitude = 1/3;
                    Tau = s*magnitude*obj.Hi*e*H_b/2;
                    m = length(index);
                    tau = @(v) Tau*spdiags((v(index)-s*abs(v(index))),0,m,m);
                    closure = @(v) Tau*((v(index)-s*abs(v(index))).*v(index));
                    penalty = @(v) -tau(v);
                otherwise
                    error('No such boundary condition: type = %s',type);
            end


        end

        % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary.
        % The right boundary for each coordinate direction is considered the positive boundary
        function [e, H_b, index, s] = get_boundary_ops(obj, boundary)
            ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N());
            switch boundary
                case {'w', 'W', 'west', 'West'}
                    e = obj.e_w;
                    H_b = obj.H_y;
                    index = ind(1,:);
                    s = -1;
                case {'e', 'E', 'east', 'East'}
                    e = obj.e_e;
                    H_b = obj.H_y;
                    index = ind(end,:);
                    s = 1;
                case {'s', 'S', 'south', 'South'}
                    e = obj.e_s;
                    H_b = obj.H_x;
                    index = ind(:,1);
                    s = -1;
                case {'n', 'N', 'north', 'North'}
                    e = obj.e_n;
                    H_b = obj.H_x;
                    index = ind(:,end);
                    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