--- a/+scheme/Hypsyst3d.m	Fri Oct 12 08:50:25 2018 +0200
+++ b/+scheme/Hypsyst3d.m	Thu Jan 17 10:44:12 2019 +0100
@@ -7,11 +7,11 @@
         X, Y, Z% Values of x and y for each grid point
         Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
         order % Order accuracy for the approximation
-        
+
         D % non-stabalized scheme operator
         A, B, C, E % Symbolic coefficient matrices
         Aevaluated,Bevaluated,Cevaluated, Eevaluated
-        
+
         H % Discrete norm
         Hx, Hy, Hz  % Norms in the x, y and z directions
         Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
@@ -19,8 +19,8 @@
         e_w, e_e, e_s, e_n, e_b, e_t
         params % Parameters for the coeficient matrice
     end
-    
-    
+
+
     methods
         % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
         function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
@@ -28,11 +28,11 @@
             xlim =  lim{1};
             ylim = lim{2};
             zlim = lim{3};
-            
+
             if length(m) == 1
                 m = [m m m];
             end
-            
+
             obj.A = A;
             obj.B = B;
             obj.C = C;
@@ -41,7 +41,7 @@
             m_y = m(2);
             m_z = m(3);
             obj.params = params;
-            
+
             switch operator
                 case 'upwind'
                     ops_x = sbp.D1Upwind(m_x,xlim,order);
@@ -52,29 +52,29 @@
                     ops_y = sbp.D2Standard(m_y,ylim,order);
                     ops_z = sbp.D2Standard(m_z,zlim,order);
             end
-            
+
             obj.x = ops_x.x;
             obj.y = ops_y.x;
             obj.z = ops_z.x;
-            
+
             obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));
             obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
             obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
-            
+
             obj.Yx = kr(obj.y,ones(m_z,1));
             obj.Zx = kr(ones(m_y,1),obj.z);
             obj.Xy = kr(obj.x,ones(m_z,1));
             obj.Zy = kr(ones(m_x,1),obj.z);
             obj.Xz = kr(obj.x,ones(m_y,1));
             obj.Yz = kr(ones(m_z,1),obj.y);
-            
+
             obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
             obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
             obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
             obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
-            
+
             obj.n = length(A(obj.params,0,0,0));
-            
+
             I_n = speye(obj.n);
             I_x = speye(m_x);
             obj.I_x = I_x;
@@ -83,31 +83,31 @@
             I_z = speye(m_z);
             obj.I_z = I_z;
             I_N = kr(I_n,I_x,I_y,I_z);
-            
+
             obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
             obj.Hx = ops_x.H;
             obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
             obj.Hy = ops_y.H;
             obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI);
             obj.Hz = ops_z.H;
-            
+
             obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
             obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
             obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
             obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
             obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
             obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
-            
+
             obj.m = m;
             obj.h = [ops_x.h ops_y.h ops_x.h];
             obj.order = order;
-            
+
             switch operator
                 case 'upwind'
                     alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end)))));
                     alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end)))));
                     alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end)))));
-                    
+
                     Ap = (obj.Aevaluated+alphaA*I_N)/2;
                     Am = (obj.Aevaluated-alphaA*I_N)/2;
                     Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
@@ -116,7 +116,7 @@
                     temp = Ap*Dmx;
                     obj.D = obj.D-temp;
                     clear Ap Am Dpx Dmx
-                    
+
                     Bp = (obj.Bevaluated+alphaB*I_N)/2;
                     Bm = (obj.Bevaluated-alphaB*I_N)/2;
                     Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
@@ -126,20 +126,20 @@
                     temp = Bp*Dmy;
                     obj.D = obj.D-temp;
                     clear Bp Bm Dpy Dmy
-                    
-                    
+
+
                     Cp = (obj.Cevaluated+alphaC*I_N)/2;
                     Cm = (obj.Cevaluated-alphaC*I_N)/2;
                     Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
                     Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
-                    
+
                     temp = Cm*Dpz;
                     obj.D = obj.D-temp;
                     temp = Cp*Dmz;
                     obj.D = obj.D-temp;
                     clear Cp Cm Dpz Dmz
                     obj.D = obj.D-obj.Eevaluated;
-                    
+
                 case 'standard'
                     D1_x = kr(I_n, ops_x.D1, I_y,I_z);
                     D1_y = kr(I_n, I_x, ops_y.D1,I_z);
@@ -149,7 +149,7 @@
                     error('Opperator not supported');
             end
         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 doamin.
         %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
@@ -167,15 +167,15 @@
                     error('No such boundary condition')
             end
         end
-        
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
-            error('An interface function does not exist yet');
+
+        function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
+            error('Not implemented');
         end
-        
+
         function N = size(obj)
             N = obj.m;
         end
-        
+
         function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
             params = obj.params;
             side = max(length(X),length(Y));
@@ -189,7 +189,7 @@
                 side = max(length(X),length(Y));
                 cols = cols/side;
             end
-            
+
             ret = cell(rows,cols);
             for ii = 1:rows
                 for jj = 1:cols
@@ -198,10 +198,10 @@
             end
             ret = cell2mat(ret);
         end
-        
+
         function [BM] = boundary_matrices(obj,boundary)
             params = obj.params;
-            
+
             switch boundary
                 case {'w','W','west'}
                     BM.e_ = obj.e_w;
@@ -248,7 +248,7 @@
             end
             BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
         end
-        
+
         % Characteristic bouyndary consitions
         function [closure, penalty]=boundary_condition_char(obj,BM)
             side = BM.side;
@@ -260,7 +260,7 @@
             Hi = BM.Hi;
             D = BM.D;
             e_ = BM.e_;
-            
+
             switch BM.boundpos
                 case {'l'}
                     tau = sparse(obj.n*side,pos);
@@ -276,9 +276,9 @@
                     penalty = -Hi*e_*V*tau*Vi_minus;
             end
         end
-        
+
         % General boundary condition in the form Lu=g(x)
-        function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L)           
+        function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L)
             side = BM.side;
             pos = BM.pos;
             neg = BM.neg;
@@ -288,7 +288,7 @@
             Hi = BM.Hi;
             D = BM.D;
             e_ = BM.e_;
-            
+
             switch boundary
                 case {'w','W','west'}
                     L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
@@ -303,7 +303,7 @@
                 case {'t','T','top'}
                     L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
             end
-            
+
             switch BM.boundpos
                 case {'l'}
                     tau = sparse(obj.n*side,pos);
@@ -311,7 +311,7 @@
                     Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
                     V_plus = V(:,1:pos);
                     V_minus = V(:,(pos+zeroval)+1:obj.n*side);
-                    
+
                     tau(1:pos,:) = -abs(D(1:pos,1:pos));
                     R = -inv(L*V_plus)*(L*V_minus);
                     closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
@@ -321,7 +321,7 @@
                     tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
                     Vi_plus = Vi(1:pos,:);
                     Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
-                    
+
                     V_plus = V(:,1:pos);
                     V_minus = V(:,(pos+zeroval)+1:obj.n*side);
                     R = -inv(L*V_minus)*(L*V_plus);
@@ -329,7 +329,7 @@
                     penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
             end
         end
-        
+
         % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
         % D         is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
         %                                    [d+       ]
@@ -344,13 +344,13 @@
             xs = x;
             ys = y;
             zs = z;
-            
-            
+
+
             side = max(length(x),length(y));
             Dret = zeros(obj.n,side*obj.n);
             Vret = zeros(obj.n,side*obj.n);
             Viret= zeros(obj.n,side*obj.n);
-           
+
             for ii=1:obj.n
                 for jj=1:obj.n
                     Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
@@ -358,7 +358,7 @@
                     Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
                 end
             end
-            
+
             D = sparse(Dret);
             V = sparse(Vret);
             Vi = sparse(Viret);
@@ -366,11 +366,11 @@
             Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z);
             D = obj.evaluateCoefficientMatrix(D,x,y,z);
             DD = diag(D);
-            
+
             poseig = (DD>0);
             zeroeig = (DD==0);
             negeig = (DD<0);
-            
+
             D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
             V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
             Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];