--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/AdvectionRV1D.m	Wed Nov 07 10:03:35 2018 -0800
@@ -0,0 +1,88 @@
+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
+    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';
+                    D2 = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp);
+                    DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2;
+                otherwise
+                    error('Other operator types not yet supported', operator_type);
+            end
+            % max(abs()) or just abs()?
+            obj.D = @(viscosity) (-D1 + D2(viscosity) + 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
\ No newline at end of file