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 = (ops.Dp + ops.Dm)/2;
                DissOp = (ops.Dm - ops.Dp)/2;
            else 
                obj.D1 = ops.D1;
            end
            
            switch pde_form
                case 'skew-symmetric'
                    if (isequal(opSet, @sbp.D1Upwind))
                        D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1  + fluxSplitting(v)*DissOp)*v);
                    else
                        D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
                    end
                case 'conservative'
                    if (isequal(opSet, @sbp.D1Upwind))
                        D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
                    else
                        D = @(v) -(1/2*obj.D1*v.*v);
                    end
                otherwise
                    error('Not supported', pde_form);
            end

            obj.grid = g;

            obj.D = D;
            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','robin');
            [e, index, s] = obj.get_boundary_ops(boundary);
            switch type
                % Stable dirhclet-like boundary conditions ((u+-abs(u))*u/3)) with +- at left/right boundary
                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;
                otherwise
                    error('No such boundary condition: type = %s',type);
            end
        end

        % Returns 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, index, s] = get_boundary_ops(obj,boundary)
            switch boundary
                case {'l','L','left','Left'}
                    e = obj.e_l;
                    index = 1;
                    s = -1;
                case {'r','R','right','Right'}
                    e = obj.e_r;
                    index = 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