classdef Burgers1d < scheme.Scheme
    properties
        m % Number of points in each direction, possibly a vector
        h % Grid spacing
        grid % Grid
        order % Order accuracy for the approximation

        H % Discrete norm
        D
        
        D1
        Hi
        e_l
        e_r
    end

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

            m = g.size();
            xl = g.getBoundary('l');
            xr = g.getBoundary('r');
            xlim = {xl, xr};

            ops = opSet(m, xlim, order);

            if (isequal(opSet, @sbp.D1Upwind))
                obj.D1 = sparse((ops.Dp + ops.Dm)/2);
                DissOp = sparse((ops.Dm - ops.Dp)/2);
                switch pde_form
                    case 'quasi-linear'
                        obj.D = @(v) -((spdiag(v)*obj.D1  + fluxSplitting(v)*DissOp)*v);
                    case 'skew-symmetric'
                        obj.D = @(v) -(1/3*obj.D1*(v.*v) + (1/3*spdiag(v)*obj.D1  + fluxSplitting(v)*DissOp)*v);
                    case 'conservative'
                        obj.D = @(v) -(1/2*obj.D1*(v.*v) + fluxSplitting(v)*DissOp*v);
                end
            else 
                obj.D1 = sparse(ops.D1);
                switch pde_form
                    case 'quasi-linear'
                        obj.D = @(v) -(spdiag(v)*obj.D1*v);
                    case 'skew-symmetric'
                        obj.D = @(v) -(1/3*obj.D1*(v.*v) + 1/3*spdiag(v)*obj.D1*v);
                    case 'conservative'
                        obj.D = @(v) -1/2*obj.D1*(v.*v);
                end
            end
            obj.grid = g;

            obj.H =  ops.H;
            obj.Hi = ops.HI;

            obj.e_l = ops.e_l;
            obj.e_r = ops.e_r;

            obj.m = m;
            obj.h = ops.h;
            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);
            index = obj.getBoundaryIndex(boundary);
            switch type
                % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
                % with +- at left/right boundaries
                case {'D', 'd', 'dirichlet', 'Dirichlet'}
                    % tau = s*e;
                    % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
                    % penalty = -obj.Hi*tau;

                    penalty_parameter = 1/3;
                    tau = @(v) s*penalty_parameter*obj.Hi*e*(v(index)-s*abs(v(index)))/2;
                    closure = @(v) tau(v)*v(index);
                    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, {'l', 'r'})

            switch boundary
                case {'r'}
                    s = 1;
                case {'l'}
                    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, {'l', 'r'})

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

        % Returns square boundary quadrature matrix, of dimension
        % corresponding to the number of boundary points
        %
        % boundary -- string
        % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
        function H_b = getBoundaryQuadrature(obj, boundary)
            assertIsMember(boundary, {'l', 'r'})
            H_b = 1;
        end

        % Returns the boundary index. The right boundary has the last index
        % boundary -- string
        function index = getBoundaryIndex(obj, boundary)
            assertIsMember(boundary, {'l', 'r'})
            switch boundary
                case {'r'}
                    index = length(obj.e_r);
                case {'l'}
                    index = 1;
            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