--- a/+scheme/Burgers1d.m	Fri Jan 18 09:02:02 2019 +0100
+++ b/+scheme/Burgers1d.m	Fri Jan 18 09:04:38 2019 +0100
@@ -1,136 +1,99 @@
 classdef Burgers1d < scheme.Scheme
     properties
-        grid % Physical grid
+        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
         
-        params
-
-        D % Non-stabalized scheme operator
-        H % Discrete norm
-        Hi % Norm inverse
+        D1
+        Hi
         e_l
         e_r
-        d_l
-        d_r
     end
 
     methods
-        function obj = Burgers1d(grid, pde_form, operator_type, order, dissipation, params)
-            assert(grid.D == 1);
-            assert(grid.size() == length(params.eps));
-            m = grid.size();
-            lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
-            switch operator_type
-                case 'narrow'
-                    ops = sbp.D4Variable(m, lim, order);
-                    D1 = ops.D1;
-                    D2 = ops.D2;
-                    if (strcmp(dissipation,'on'))
-                        DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
-                    end
-                    d_l = ops.d1_l';
-                    d_r = ops.d1_r';
-                case 'upwind-'
-                    ops = sbp.D1Upwind(m, lim, order);
-                    D1 = (ops.Dp + ops.Dm)/2;
-                    D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
-                    if (strcmp(dissipation,'on'))
-                        DissipationOp = (ops.Dp-ops.Dm)/2;
-                    end
-                    d_l = ops.e_l'*ops.Dm;
-                    d_r = ops.e_r'*ops.Dm;
-                case 'upwind+'
-                    ops = sbp.D1Upwind(m, lim, order);
-                    D1 = (ops.Dp + ops.Dm)/2;
-                    D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp;
-                    if (strcmp(dissipation,'on'))
-                        DissipationOp = (ops.Dp-ops.Dm)/2;
-                    end
-                    d_l = ops.e_l'*ops.Dp;
-                    d_r = ops.e_r'*ops.Dp;
-                case 'upwind+-'
-                    ops = sbp.D1Upwind(m, lim, order);
-                    D1 = (ops.Dp + ops.Dm)/2;
-                    D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2;
-                    if (strcmp(dissipation,'on'))
-                        DissipationOp = (ops.Dp-ops.Dm)/2;
-                    end
-                    d_l = ops.e_l'*D1;
-                    d_r = ops.e_r'*D1;
-                otherwise
-                    error('Other operator types not yet supported', operator_type);
+        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 (strcmp(dissipation,'on'))
-                        D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
+                    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, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v;
+                        D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
                     end
                 case 'conservative'
-                    if (strcmp(dissipation,'on'))
-                        D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
+                    if (isequal(opSet, @sbp.D1Upwind))
+                        D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
                     else
-                        D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
+                        D = @(v) -(1/2*obj.D1*v.*v);
                     end
                 otherwise
                     error('Not supported', pde_form);
             end
 
-            obj.grid = grid;
-            obj.order = order;
-            obj.params = params;
+            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.d_l =  d_l;
-            obj.d_r =  d_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.
-        %       data                is a function returning the data that should be applied at the boundary.
-        function [closure, penalty] = boundary_condition(obj,boundary,type,data)
+        function [closure, penalty] = boundary_condition(obj, boundary, type)
             default_arg('type','robin');
-            default_arg('data',0);
-            [e, d, i_b, s] = obj.get_boundary_ops(boundary);
+            [e, index, s] = 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*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - ((obj.params.eps(i_b) + viscosity(i_b))*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
+                % 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
 
-        % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
+        % 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, d, i_b, s] = get_boundary_ops(obj,boundary)
+        function [e, index, s] = get_boundary_ops(obj,boundary)
             switch boundary
-                case 'l'
+                case {'l','L','left','Left'}
                     e = obj.e_l;
-                    d = obj.d_l;
-                    i_b = 1;
+                    index = 1;
                     s = -1;
-                case 'r'
+                case {'r','R','right','Right'}
                     e = obj.e_r;
-                    d = obj.d_r;
-                    i_b = length(e);
+                    index = length(e);
                     s = 1;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);