--- a/+scheme/Utux2d.m	Mon Jan 07 13:12:10 2019 +0100
+++ b/+scheme/Utux2d.m	Mon Jan 07 16:26:05 2019 +0100
@@ -4,7 +4,6 @@
         h % Grid spacing
         grid % Grid
         order % Order accuracy for the approximation
-        v0 % Initial data
 
         a % Wave speed a = [a1, a2];
           % Can either be a constant vector or a cell array of function handles.
@@ -16,9 +15,6 @@
         % Derivatives
         Dx, Dy
 
-        % Dissipation Operators
-        DissOpx, DissOpy
-
         % Boundary operators
         e_w, e_e, e_s, e_n
 
@@ -76,17 +72,21 @@
             if (isequal(opSet,@sbp.D1Upwind))
                 Dx = (ops_x.Dp + ops_x.Dm)/2;
                 Dy = (ops_y.Dp + ops_y.Dm)/2;
+                obj.Dx = kron(Dx,Iy);
+                obj.Dy = kron(Ix,Dy);
                 DissOpx = (ops_x.Dm - ops_x.Dp)/2;
                 DissOpy = (ops_y.Dm - ops_y.Dp)/2;
-                obj.Dx = kron(Dx,Iy);
-                obj.Dy = kron(Ix,Dy);
-                obj.DissOpx = kron(DissOpx,Iy);
-                obj.DissOpy = kron(Ix,DissOpy);
+                DissOpx = kron(DissOpx,Iy);
+                DissOpy = kron(Ix,DissOpy);
+
+                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*DissOpx + fluxSplitting{2}*DissOpy);
             else
                 Dx = ops_x.D1;
                 Dy = ops_y.D1;
                 obj.Dx = kron(Dx,Iy);
                 obj.Dy = kron(Ix,Dy);
+
+                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
             end
 
             % Boundary operators
@@ -99,12 +99,6 @@
             obj.h = [ops_x.h ops_y.h];
             obj.order = order;
             obj.a = a;
-
-            if (isequal(opSet,@sbp.D1Upwind))
-                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*obj.DissOpx + fluxSplitting{2}*obj.DissOpy);
-            else
-                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
-            end
         end
         % Closure functions return the opertors applied to the own domain 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 doamin.
@@ -115,36 +109,32 @@
         %       neighbour_boundary  is a string specifying which boundary to interface to.
         function [closure, penalty] = boundary_condition(obj,boundary,type)
             default_arg('type','dirichlet');
-            sigma_left = -1; % Scalar penalty parameter for left boundaries
-            sigma_right = 1; % Scalar penalty parameter for right boundaries
+            sigma_left = -1; % Scalar penalty parameter for left boundaries (West/South)
+            sigma_right = 1; % Scalar penalty parameter for right boundaries (East/North)
             switch boundary
+                % Can only specify boundary condition where there is inflow
+                % Extract the postivie resp. negative part of a, for the left
+                % resp. right boundaries, and set other values of a to zero.
+                % Then the closure will effectively only contribute to inflow boundaries
                 case {'w','W','west','West'}
-                    % Can only specify boundary condition where there is inflow
-                    % Extract positive part of a{1}
-                    a_pos = obj.a{1};
-                    a_pos(a_pos < 0) = 0;
-                    tau = sigma_left*a_pos*obj.e_w*obj.H_y;
+                    a_inflow = obj.a{1};
+                    a_inflow(a_inflow < 0) = 0;
+                    tau = sigma_left*a_inflow*obj.e_w*obj.H_y;
                     closure = obj.Hi*tau*obj.e_w';
+                case {'e','E','east','East'}
+                    a_inflow = obj.a{1};
+                    a_inflow(a_inflow > 0) = 0;
+                    tau = sigma_right*a_inflow*obj.e_e*obj.H_y;
+                    closure = obj.Hi*tau*obj.e_e';
                 case {'s','S','south','South'}
-                    % Can only specify boundary condition where there is inflow
-                    % Extract positive part of a{2}
-                    a_pos = obj.a{2};
-                    a_pos(a_pos < 0) = 0;
-                    tau = sigma_left*a_pos*obj.e_s*obj.H_x;
+                    a_inflow = obj.a{2};
+                    a_inflow(a_inflow < 0) = 0;
+                    tau = sigma_left*a_inflow*obj.e_s*obj.H_x;
                     closure = obj.Hi*tau*obj.e_s';
-                case {'e','E','east','East'}
-                    % Can only specify boundary condition where there is inflow
-                    % Extract negative part of a{1}
-                    a_neg = obj.a{1};
-                    a_neg(a_neg > 0) = 0;
-                    tau = sigma_right*a_neg*obj.e_e*obj.H_y;
-                    closure = obj.Hi*tau*obj.e_e';
                 case {'n','N','north','North'}
-                    % Can only specify boundary condition where there is inflow
-                    % Extract negative part of a{2}
-                    a_neg = obj.a{2};
-                    a_neg(a_neg > 0) = 0;
-                    tau = sigma_right*a_neg*obj.e_n*obj.H_x;
+                    a_inflow = obj.a{2};
+                    a_inflow(a_inflow > 0) = 0;
+                    tau = sigma_right*a_inflow*obj.e_n*obj.H_x;
                     closure = obj.Hi*tau*obj.e_n';
             end
             penalty = -obj.Hi*tau;