--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Burgers1D.m	Tue Sep 11 15:58:35 2018 +0200
@@ -0,0 +1,114 @@
+classdef Burgers1D < scheme.Scheme
+    properties
+        grid % Physical grid
+        order % Order accuracy for the approximation
+        
+        params
+
+        D % Non-stabalized scheme operator
+        M % Derivative norm
+        H % Discrete norm
+        Hi % Norm inverse
+        e_l
+        e_r
+        d_l
+        d_r
+    end
+
+    methods
+        function obj = Burgers1D(grid, pde_form, operator_type, order, params)
+            assert(grid.D == 1);
+            assert(grid.size() == length(params.eps));
+            m = grid.size();
+            h = grid.scaling();
+            lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
+            default_arg('pde_form','skew-symmetric');
+            default_arg('operator_type','narrow');
+
+            switch operator_type
+                case 'narrow'
+                    ops = sbp.D2Variable(m, lim, order);
+                    D1 = ops.D1;
+                    D2 = ops.D2;  
+                otherwise
+                    error('Other operator types not yet supported', operator_type);
+            end
+
+            switch pde_form
+                case 'skew-symmetric'
+                    D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(params.eps + viscosity)*v;
+                case 'conservative'
+                    D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
+            end
+
+            obj.grid = grid;
+            obj.order = order;
+            obj.params = params;
+
+            %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
+            obj.D = D;
+            obj.M = ops.M;
+            obj.H =  ops.H;
+            obj.Hi = ops.HI;
+            obj.e_l = ops.e_l;
+            obj.e_r = ops.e_r;
+            obj.d_l =  ops.d1_l;
+            obj.d_r =  ops.d1_r;
+        end
+
+        % Closure functions return the opertors applied to the own doamin to close the boundary
+        % Penalty functions return the opertors 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.
+        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
+        %       neighbour_boundary  is a string specifying which boundary to interface to.
+        function [closure, penalty] = boundary_condition(obj,boundary,type,data)
+            default_arg('type','robin');
+            default_arg('data',0);
+            [e, s, d] = obj.get_boundary_ops(boundary);
+            switch type
+                % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
+                case {'R','robin'}
+                    p = s*obj.Hi*e;
+                    closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v);
+                    switch class(data)
+                        case 'double'
+                            penalty = s*p*data;
+                        case 'function_handle'
+                            penalty = @(t) s*p*data(t);
+                        otherwise
+                            error('Wierd data argument!')
+                    end
+                % Unknown, boundary condition
+                otherwise
+                    error('No such boundary condition: type = %s',type);
+            end
+        end
+
+        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
+        % The right boundary is considered the positive boundary
+        function [e, s, d] = get_boundary_ops(obj,boundary)
+            switch boundary
+                case 'l'
+                    e = obj.e_l;
+                    s = -1;
+                    d = obj.d_l;
+                case 'r'
+                    e = obj.e_r;
+                    s = 1;
+                    d = obj.d_r;
+                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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/Rungekutta4RV.m	Tue Sep 11 15:58:35 2018 +0200
@@ -0,0 +1,54 @@
+classdef Rungekutta4RV < time.Timestepper
+    properties
+        F
+        k
+        t
+        v
+        m
+        n
+
+        % Additional members used for the RV update
+        RV
+    end
+
+
+    methods
+        function obj = Rungekutta4RV(F, k, t0, v0, RV)
+            obj.F = F;
+            obj.k = k;
+            obj.t = t0;
+            obj.v = v0;
+            obj.m = length(v0);
+            obj.n = 0;
+            obj.RV = RV;
+        end
+
+        function [v, t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function [residual, viscosity, t] = getRV(obj)
+            residual = obj.RV.getResidual();
+            viscosity = obj.RV.getViscosity();
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            v_prev = obj.v;
+            F = @(v,t) obj.F(v, t, obj.RV.getViscosity());
+            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, F);
+            obj.RV.update(obj.v, v_prev, obj.k);
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+        end
+    end
+
+
+    methods (Static)
+        function k = getTimeStep(lambda)
+            k = rk4.get_rk4_time_step(lambda);
+        end
+    end
+
+end
\ No newline at end of file