diff +sbp/D1Staggered.m @ 636:476b7e411ce3 feature/d1_staggered

Add staggered operators
author Martin Almquist <malmquist@stanford.edu>
date Tue, 07 Nov 2017 20:44:13 -0800
parents
children eead18a8964d
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D1Staggered.m	Tue Nov 07 20:44:13 2017 -0800
@@ -0,0 +1,83 @@
+classdef D1Staggered < sbp.OpSet
+    properties
+        % xp: "plus" grid with m points
+        % xm: "minus" grid with m+1 points
+
+        D1p % SBP operator approximating first derivative
+        D1m % SBP operator approximating first derivative
+        Hp % Norm matrix
+        Hm % Norm matrix
+        HpI % H^-1
+        HmI % H^-1
+        ep_l % Left boundary operator
+        em_l % Left boundary operator
+        ep_r % Right boundary operator
+        em_r % Right boundary operator
+        m % Number of grid points.
+        mp % Number of grid points.
+        mm % Number of grid points.
+        h % Step size
+        xp % grid
+        xm % grid
+        x
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D1Staggered(m,lim,order)
+
+          x_l = lim{1};
+          x_r = lim{2};
+          L = x_r-x_l;
+
+          mp = m;
+          mm = m+1;
+
+          switch order
+          case 2
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_2(m, L);
+          case 4
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_4(m, L);
+          case 6
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_6(m, L);
+          otherwise
+           error('Invalid operator order %d.',order);
+          end
+
+          obj.m = m;
+          obj.mp = mp;
+          obj.mm = mm;
+          obj.xp = x_l + xp';
+          obj.xm = x_l + xm';
+
+          D1p = Pp\Qp;
+          D1m = Pm\Qm;
+
+          obj.D1p = D1p;
+          obj.D1m = D1m;
+          pbj.Hp = Pp;
+          obj.Hm = Pm; 
+
+          obj.ep_l = sparse(mp,1);
+          obj.ep_r = sparse(mp,1);
+          obj.ep_l(1) = 1;
+          obj.ep_r(mp) = 1;
+
+          obj.em_l = sparse(mm,1);
+          obj.em_r = sparse(mm,1);
+          obj.em_l(1) = 1;
+          obj.em_r(mm) = 1;
+
+          obj.HpI = inv(obj.Hp);
+          obj.HmI = inv(obj.Hm);
+
+          obj.borrowing = [];
+          obj.x = [];
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+end