--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D1StaggeredUpwind.m	Thu Mar 10 16:54:26 2022 +0100
@@ -0,0 +1,99 @@
+classdef D1StaggeredUpwind < sbp.OpSet
+    % Compatible staggered and upwind operators by Ken Mattsson and Ossian O'reilly
+    properties
+        % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper.
+        % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper.
+
+        % D1_primal takes FROM dual grid TO primal grid
+        % D1_dual takes FROM primal grid TO dual grid
+
+        D1_primal % SBP operator approximating first derivative
+        D1_dual % SBP operator approximating first derivative
+
+        Dplus_primal  % Upwind operator on primal grid
+        Dminus_primal % Upwind operator on primal grid
+        Dplus_dual % Upwind operator on dual grid
+        Dminus_dual % Upwind operator on dual grid
+
+        H_primal % Norm matrix
+        H_dual % Norm matrix
+        H_primalI % H^-1
+        H_dualI % H^-1
+        e_primal_l % Left boundary operator
+        e_dual_l % Left boundary operator
+        e_primal_r % Right boundary operator
+        e_dual_r % Right boundary operator
+        m % Number of grid points.
+        m_primal % Number of grid points.
+        m_dual % Number of grid points.
+        h % Step size
+        x_primal % grid
+        x_dual % grid
+        x
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D1StaggeredUpwind(m,lim,order)
+
+          xl = lim{1};
+          xr = lim{2};
+          L = xr-xl;
+          h = L/(m-1);
+
+          m_primal = m;
+          m_dual = m+1;
+
+          switch order
+          case 2
+            [~, ~, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_2(m, L);
+          case 4
+            [~, ~, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_4(m, L);
+          case 6
+            [~, ~, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_6(m, L);
+          case 8
+            [~, ~, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_8(m, L);
+          otherwise
+           error('Invalid operator order %d.',order);
+          end
+
+          obj.m = m;
+          obj.m_primal = m_primal;
+          obj.m_dual = m_dual;
+          obj.h = h;
+
+          obj.x_primal = linspace(xl, xr, m)';
+          obj.x_dual = [xl, linspace(xl+h/2, xr-h/2, m-1), xr]';
+
+          obj.e_primal_l = sparse(m_primal,1);
+          obj.e_primal_r = sparse(m_primal,1);
+          obj.e_primal_l(1) = 1;
+          obj.e_primal_r(m_primal) = 1;
+
+          obj.e_dual_l = sparse(m_dual,1);
+          obj.e_dual_r = sparse(m_dual,1);
+          obj.e_dual_l(1) = 1;
+          obj.e_dual_r(m_dual) = 1;
+
+          obj.borrowing = [];
+          obj.x = [];
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+end