changeset 756:f891758ad7a4 feature/d1_staggered

Merge with feature/utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Jun 2018 14:30:45 -0700
parents 14f0058356f2 (current diff) 0be9b4d6737b (diff)
children 179f234f6cbf
files +multiblock/Def.m
diffstat 40 files changed, 2930 insertions(+), 142 deletions(-) [+]
line wrap: on
line diff
--- a/+grid/Cartesian.m	Fri Jun 15 18:10:26 2018 -0700
+++ b/+grid/Cartesian.m	Sat Jun 16 14:30:45 2018 -0700
@@ -5,6 +5,7 @@
         m % Number of points in each direction
         x % Cell array of vectors with node placement for each dimension.
         h % Spacing/Scaling
+        lim % Cell array of left and right boundaries for each dimension.
     end
 
     % General d dimensional grid with n points
@@ -27,6 +28,7 @@
             end
 
             obj.h = [];
+            obj.lim = [];
         end
         % n returns the number of points in the grid
         function o = N(obj)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Circle.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,98 @@
+classdef Circle < multiblock.DefCurvilinear
+    properties
+        r, c
+
+        hs
+        r_arc
+        omega
+    end
+
+    methods
+        function obj = Circle(r, c, hs)
+            default_arg('r', 1);
+            default_arg('c', [0; 0]);
+            default_arg('hs', 0.435);
+
+
+            % alpha = 0.75;
+            % hs = alpha*r/sqrt(2);
+
+            % Square should not be a square, it should be an arc. The arc radius
+            % is chosen so that the three angles of the meshes are all equal.
+            % This gives that the (half)arc opening angle of should be omega = pi/12
+            omega = pi/12;
+            r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); %  = hs* 1/sin(omega)
+            c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0];
+
+            cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]);
+
+            c2 = cir(0);
+            c3 = cir(1);
+
+            s1 = [-hs; -hs];
+            s2 = [ hs; -hs];
+            s3 = [ hs;  hs];
+            s4 = [-hs;  hs];
+
+            sp2 = parametrization.Curve.line(s2,c2);
+            sp3 = parametrization.Curve.line(s3,c3);
+
+            Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]);
+            Se2 = Se1.rotate(c,pi/2);
+            Se3 = Se2.rotate(c,pi/2);
+            Se4 = Se3.rotate(c,pi/2);
+
+
+            S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1);
+
+            A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse);
+            B = A.rotate(c,1*pi/2).rotate_edges(-1);
+            C = A.rotate(c,2*pi/2).rotate_edges(-1);
+            D = A.rotate(c,3*pi/2).rotate_edges(0);
+
+            blocks = {S,A,B,C,D};
+            blocksNames = {'S','A','B','C','D'};
+
+            conn = cell(5,5);
+            conn{1,2} = {'e','w'};
+            conn{1,3} = {'n','s'};
+            conn{1,4} = {'w','s'};
+            conn{1,5} = {'s','w'};
+
+            conn{2,3} = {'n','e'};
+            conn{3,4} = {'w','e'};
+            conn{4,5} = {'w','s'};
+            conn{5,2} = {'n','s'};
+
+            boundaryGroups = struct();
+            boundaryGroups.E = multiblock.BoundaryGroup({2,'e'});
+            boundaryGroups.N = multiblock.BoundaryGroup({3,'n'});
+            boundaryGroups.W = multiblock.BoundaryGroup({4,'n'});
+            boundaryGroups.S = multiblock.BoundaryGroup({5,'e'});
+            boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}});
+
+            obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames);
+
+            obj.r     = r;
+            obj.c     = c;
+            obj.hs    = hs;
+            obj.r_arc = r_arc;
+            obj.omega = omega;
+        end
+
+        function ms = getGridSizes(obj, m)
+            m_S = m;
+
+            % m_Radial
+            s = 2*obj.hs;
+            innerArc = obj.r_arc*obj.omega;
+            outerArc = obj.r*pi/2;
+            shortSpoke = obj.r-s/sqrt(2);
+            x = (1/(2-sqrt(3))-1)*obj.hs;
+            longSpoke =  (obj.r+x)-obj.r_arc;
+            m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2);
+
+            ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]};
+        end
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Rectangle.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,188 @@
+classdef Rectangle < multiblock.Definition
+    properties
+
+    blockTi % Transfinite interpolation objects used for plotting
+    xlims
+    ylims
+    blockNames % Cell array of block labels
+    nBlocks
+    connections % Cell array specifying connections between blocks
+    boundaryGroups % Structure of boundaryGroups
+
+    end
+
+
+    methods
+        % Creates a divided rectangle
+        % x and y are vectors of boundary and interface positions.
+        % blockNames: cell array of labels. The id is default.
+        function obj = Rectangle(x,y,blockNames)
+            default_arg('blockNames',[]);
+
+            n = length(y)-1; % number of blocks in the y direction.
+            m = length(x)-1; % number of blocks in the x direction.
+            N = n*m; % number of blocks
+
+            if ~issorted(x)
+                error('The elements of x seem to be in the wrong order');
+            end
+            if ~issorted(flip(y))
+                error('The elements of y seem to be in the wrong order');
+            end
+
+            % Dimensions of blocks and number of points
+            blockTi = cell(N,1);
+            xlims = cell(N,1);
+            ylims = cell(N,1);
+            for i = 1:n
+                for j = 1:m
+                    p1 = [x(j), y(i+1)];
+                    p2 = [x(j+1), y(i)];
+                    I = flat_index(m,j,i);
+                    blockTi{I} = parametrization.Ti.rectangle(p1,p2);
+                    xlims{I} = {x(j), x(j+1)};
+                    ylims{I} = {y(i+1), y(i)};
+                end
+            end
+
+            % Interface couplings
+            conn = cell(N,N);
+            for i = 1:n
+                for j = 1:m
+                    I = flat_index(m,j,i);
+                    if i < n
+                        J = flat_index(m,j,i+1);
+                        conn{I,J} = {'s','n'};
+                    end
+
+                    if j < m
+                        J = flat_index(m,j+1,i);
+                        conn{I,J} = {'e','w'};
+                    end
+                end
+            end
+
+            % Block names (id number as default)
+            if isempty(blockNames)
+                obj.blockNames = cell(1, N);
+                for i = 1:N
+                    obj.blockNames{i} = sprintf('%d', i);
+                end
+            else
+                assert(length(blockNames) == N);
+                obj.blockNames = blockNames;
+            end
+            nBlocks = N;
+
+            % Boundary groups
+            boundaryGroups = struct();
+            nx = m;
+            ny = n;
+            E = cell(1,ny);
+            W = cell(1,ny);
+            S = cell(1,nx);
+            N = cell(1,nx);
+            for i = 1:ny
+                E_id = flat_index(m,nx,i);
+                W_id = flat_index(m,1,i);
+                E{i} = {E_id,'e'};
+                W{i} = {W_id,'w'};
+            end
+            for j = 1:nx
+                S_id = flat_index(m,j,ny);
+                N_id = flat_index(m,j,1);
+                S{j} = {S_id,'s'};
+                N{j} = {N_id,'n'};
+            end  
+            boundaryGroups.E = multiblock.BoundaryGroup(E);
+            boundaryGroups.W = multiblock.BoundaryGroup(W);
+            boundaryGroups.S = multiblock.BoundaryGroup(S);
+            boundaryGroups.N = multiblock.BoundaryGroup(N);
+            boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]);
+            boundaryGroups.WS = multiblock.BoundaryGroup([W,S]);
+            boundaryGroups.WN = multiblock.BoundaryGroup([W,N]);
+            boundaryGroups.ES = multiblock.BoundaryGroup([E,S]);
+            boundaryGroups.EN = multiblock.BoundaryGroup([E,N]);
+
+            obj.connections = conn;
+            obj.nBlocks = nBlocks;
+            obj.boundaryGroups = boundaryGroups;
+            obj.blockTi = blockTi;
+            obj.xlims = xlims;
+            obj.ylims = ylims;
+
+        end
+
+
+        % Returns a multiblock.Grid given some parameters
+        % ms: cell array of [mx, my] vectors
+        % For same [mx, my] in every block, just input one vector.
+        function g = getGrid(obj, ms, varargin)
+
+            default_arg('ms',[21,21])
+
+            % Extend ms if input is a single vector
+            if (numel(ms) == 2) && ~iscell(ms)
+                m = ms;
+                ms = cell(1,obj.nBlocks);
+                for i = 1:obj.nBlocks
+                    ms{i} = m;
+                end
+            end
+
+            grids = cell(1, obj.nBlocks);
+            for i = 1:obj.nBlocks
+                grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i});
+            end
+
+            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
+        end
+
+        % label is the type of label used for plotting,
+        % default is block name, 'id' show the index for each block.
+        function show(obj, label, gridLines, varargin)
+            default_arg('label', 'name')
+            default_arg('gridLines', false);
+
+            if isempty('label') && ~gridLines
+                for i = 1:obj.nBlocks
+                    obj.blockTi{i}.show(2,2);
+                end
+                axis equal
+                return
+            end
+
+            if gridLines
+                m = 10;
+                for i = 1:obj.nBlocks
+                    obj.blockTi{i}.show(m,m);
+                end
+            end
+
+
+            switch label
+                case 'name'
+                    labels = obj.blockNames;
+                case 'id'
+                    labels = {};
+                    for i = 1:obj.nBlocks
+                        labels{i} = num2str(i);
+                    end
+                otherwise
+                    axis equal
+                    return
+            end
+
+            for i = 1:obj.nBlocks
+                parametrization.Ti.label(obj.blockTi{i}, labels{i});
+            end
+
+            axis equal
+        end
+
+        % Returns the grid size of each block in a cell array
+        % The input parameters are determined by the subclass
+        function ms = getGridSizes(obj, varargin)
+        end
+    end
+end
--- a/+multiblock/Def.m	Fri Jun 15 18:10:26 2018 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,101 +0,0 @@
-classdef Def
-    properties
-        nBlocks
-        blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation
-        blockNames
-        connections % Cell array specifying connections between blocks
-        boundaryGroups % Structure of boundaryGroups
-    end
-
-    methods
-        % Defines a multiblock setup for transfinite interpolation blocks
-        % TODO: How to bring in plotting of points?
-        function obj = Def(blockMaps, connections, boundaryGroups, blockNames)
-            default_arg('boundaryGroups', struct());
-            default_arg('blockNames',{});
-
-            nBlocks = length(blockMaps);
-
-            obj.nBlocks = nBlocks;
-
-            obj.blockMaps = blockMaps;
-
-            assert(all(size(connections) == [nBlocks, nBlocks]));
-            obj.connections = connections;
-
-
-            if isempty(blockNames)
-                obj.blockNames = cell(1, nBlocks);
-                for i = 1:length(blockMaps)
-                    obj.blockNames{i} = sprintf('%d', i);
-                end
-            else
-                assert(length(blockNames) == nBlocks);
-                obj.blockNames = blockNames;
-            end
-
-            obj.boundaryGroups = boundaryGroups;
-        end
-
-        function g = getGrid(obj, varargin)
-            ms = obj.getGridSizes(varargin{:});
-
-            grids = cell(1, obj.nBlocks);
-            for i = 1:obj.nBlocks
-                grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i});
-            end
-
-            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
-        end
-
-        function show(obj, label, gridLines, varargin)
-            default_arg('label', 'name')
-            default_arg('gridLines', false);
-
-            if isempty('label') && ~gridLines
-                for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(2,2);
-                end
-                axis equal
-                return
-            end
-
-            if gridLines
-                ms = obj.getGridSizes(varargin{:});
-                for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
-                end
-            end
-
-
-            switch label
-                case 'name'
-                    labels = obj.blockNames;
-                case 'id'
-                    labels = {};
-                    for i = 1:obj.nBlocks
-                        labels{i} = num2str(i);
-                    end
-                otherwise
-                    axis equal
-                    return
-            end
-
-            for i = 1:obj.nBlocks
-                parametrization.Ti.label(obj.blockMaps{i}, labels{i});
-            end
-
-            axis equal
-        end
-    end
-
-    methods (Abstract)
-        % Returns the grid size of each block in a cell array
-        % The input parameters are determined by the subclass
-        ms = getGridSizes(obj, varargin)
-        % end
-    end
-
-end
-
-
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/DefCurvilinear.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,101 @@
+classdef DefCurvilinear < multiblock.Definition
+    properties
+        nBlocks
+        blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation
+        blockNames
+        connections % Cell array specifying connections between blocks
+        boundaryGroups % Structure of boundaryGroups
+    end
+
+    methods
+        % Defines a multiblock setup for transfinite interpolation blocks
+        % TODO: How to bring in plotting of points?
+        function obj = DefCurvilinear(blockMaps, connections, boundaryGroups, blockNames)
+            default_arg('boundaryGroups', struct());
+            default_arg('blockNames',{});
+
+            nBlocks = length(blockMaps);
+
+            obj.nBlocks = nBlocks;
+
+            obj.blockMaps = blockMaps;
+
+            assert(all(size(connections) == [nBlocks, nBlocks]));
+            obj.connections = connections;
+
+
+            if isempty(blockNames)
+                obj.blockNames = cell(1, nBlocks);
+                for i = 1:length(blockMaps)
+                    obj.blockNames{i} = sprintf('%d', i);
+                end
+            else
+                assert(length(blockNames) == nBlocks);
+                obj.blockNames = blockNames;
+            end
+
+            obj.boundaryGroups = boundaryGroups;
+        end
+
+        function g = getGrid(obj, varargin)
+            ms = obj.getGridSizes(varargin{:});
+
+            grids = cell(1, obj.nBlocks);
+            for i = 1:obj.nBlocks
+                grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i});
+            end
+
+            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
+        end
+
+        function show(obj, label, gridLines, varargin)
+            default_arg('label', 'name')
+            default_arg('gridLines', false);
+
+            if isempty('label') && ~gridLines
+                for i = 1:obj.nBlocks
+                    obj.blockMaps{i}.show(2,2);
+                end
+                axis equal
+                return
+            end
+
+            if gridLines
+                ms = obj.getGridSizes(varargin{:});
+                for i = 1:obj.nBlocks
+                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
+                end
+            end
+
+
+            switch label
+                case 'name'
+                    labels = obj.blockNames;
+                case 'id'
+                    labels = {};
+                    for i = 1:obj.nBlocks
+                        labels{i} = num2str(i);
+                    end
+                otherwise
+                    axis equal
+                    return
+            end
+
+            for i = 1:obj.nBlocks
+                parametrization.Ti.label(obj.blockMaps{i}, labels{i});
+            end
+
+            axis equal
+        end
+    end
+
+    methods (Abstract)
+        % Returns the grid size of each block in a cell array
+        % The input parameters are determined by the subclass
+        ms = getGridSizes(obj, varargin)
+        % end
+    end
+
+end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/Definition.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,11 @@
+classdef Definition
+    methods (Abstract)
+
+        % Returns a multiblock.Grid given some parameters
+        g = getGrid(obj, varargin)
+
+        % label is the type of label used for plotting,
+        % default is block name, 'id' show the index for each block.
+        show(obj, label, gridLines, varargin)
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_2.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,50 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<3)
+        error(['Operator requires at least ' num2str(3) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+    % D1 operator
+    diags   = -1:1;
+    stencil = [-1/2 0 1/2];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    scheme_width = 3;
+    scheme_radius = (scheme_width-1)/2;
+    
+    r = 1:m;
+    offset = scheme_width;
+    r = r + offset;
+
+    function D2 = D2_fun(c)
+        c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
+
+        Mm1 = -c(r-1)/2 - c(r)/2;
+        M0  =  c(r-1)/2 + c(r)   + c(r+1)/2;
+        Mp1 =            -c(r)/2 - c(r+1)/2;
+
+        vals = [Mm1,M0,Mp1];
+        diags = -scheme_radius : scheme_radius;
+        M = spdiagsVariablePeriodic(vals,diags); 
+
+        M=M/h;
+        D2=HI*(-M );
+    end
+    D2 = @D2_fun;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_4.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,57 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<5)
+        error(['Operator requires at least ' num2str(5) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+    S = d1_l*d1_l' + d1_r*d1_r';
+
+    % D1 operator
+    stencil = [1/12 -2/3 0 2/3 -1/12];
+    diags = -2:2;
+    Q = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
+
+
+    scheme_width = 5;
+    scheme_radius = (scheme_width-1)/2;
+    
+    r = 1:m;
+    offset = scheme_width;
+    r = r + offset;
+
+    function D2 = D2_fun(c)
+        c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
+
+        % Note: these coefficients are for -M.
+        Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r);
+        Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1);
+        M0  = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2);
+        Mp1  = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2);
+        Mp2  = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2);
+
+        vals = -[Mm2,Mm1,M0,Mp1,Mp2];
+        diags = -scheme_radius : scheme_radius;
+        M = spdiagsVariablePeriodic(vals,diags); 
+
+        M=M/h;
+        D2=HI*(-M );
+
+    end
+    D2 = @D2_fun;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_6.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,58 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<7)
+        error(['Operator requires at least ' num2str(7) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+
+    % D1 operator
+    diags   = -3:3;
+    stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    % D2 operator
+    scheme_width = 7;
+    scheme_radius = (scheme_width-1)/2;
+
+    r = 1:m;
+    offset = scheme_width;
+    r = r + offset;
+
+    function D2 = D2_fun(c)
+        c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
+
+        Mm3 =  c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r);
+        Mm2 =  c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2);
+        Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1);
+        M0 =  c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r);
+        Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1);
+        Mp2 =  c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2);
+        Mp3 =  c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3);
+
+        vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3];
+        diags = -scheme_radius : scheme_radius;
+        M = spdiagsVariablePeriodic(vals,diags); 
+
+        M=M/h;
+        D2=HI*(-M );
+    end
+    D2 = @D2_fun;
+
+    
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,47 @@
+function [IC2F,IF2C,Hc,Hf] = IntOp_orders_2to2_ratio2to1(mc,hc,ACC)
+
+% ACC is a string.
+% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C.
+% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F.
+ratio = 2;
+mf = ratio*mc-1; 
+hf = hc/ratio;
+
+switch ACC
+    case 'F2C'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2;
+    case 'C2F'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1;
+end
+
+stencil_width = length(stencil_F2C);
+stencil_hw = (stencil_width-1)/2;
+[BC_rows,BC_cols] = size(BC_F2C);
+
+%%% Norm matrices %%%
+Hc = speye(mc,mc);
+HcUm = length(HcU); 
+Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm);
+Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm);
+Hc = Hc*hc;
+
+Hf = speye(mf,mf);
+HfUm = length(HfU);
+Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm);
+Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm);
+Hf = Hf*hf;
+%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IF2C from stencil and BC
+IF2C = sparse(mc,mf);
+for i = BC_rows+1 : mc-BC_rows
+    IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX>
+end
+IF2C(1:BC_rows,1:BC_cols) = BC_F2C;
+IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IC2F using symmetry condition %%%%
+IC2F = Hf\IF2C.'*Hc;
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,17 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2
+%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:36:07
+
+stencil_F2C = [-1.0./8.0,1.0./4.0,3.0./4.0,1.0./4.0,-1.0./8.0];
+if nargout > 1
+    BC_F2C = reshape([6.288191560529559e-1,-6.440957802647795e-2,6.855471317807184e-1,1.572264341096408e-1,-2.438028498638558e-1,7.469014249319279e-1,-8.431231982626708e-2,2.921561599131335e-1,1.374888185644862e-2,-1.318744409282243e-1],[2,5]);
+end
+if nargout > 2
+    HcU = 1.0./2.0;
+end
+if nargout > 3
+    HfU = 1.0./2.0;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,17 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1
+%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:36:08
+
+stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0];
+if nargout > 1
+    BC_F2C = [1.0./2.0,1.0./2.0];
+end
+if nargout > 2
+    HcU = 1.0./2.0;
+end
+if nargout > 3
+    HfU = 1.0./2.0;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,47 @@
+function [IC2F,IF2C,Hc,Hf] = IntOp_orders_4to4_ratio2to1(mc,hc,ACC)
+
+% ACC is a string.
+% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C.
+% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F.
+ratio = 2;
+mf = ratio*mc-1; 
+hf = hc/ratio;
+
+switch ACC
+    case 'F2C'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3;
+    case 'C2F'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2;
+end
+
+stencil_width = length(stencil_F2C);
+stencil_hw = (stencil_width-1)/2;
+[BC_rows,BC_cols] = size(BC_F2C);
+
+%%% Norm matrices %%%
+Hc = speye(mc,mc);
+HcUm = length(HcU); 
+Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm);
+Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm);
+Hc = Hc*hc;
+
+Hf = speye(mf,mf);
+HfUm = length(HfU);
+Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm);
+Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm);
+Hf = Hf*hf;
+%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IF2C from stencil and BC
+IF2C = sparse(mc,mf);
+for i = BC_rows+1 : mc-BC_rows
+    IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX>
+end
+IF2C(1:BC_rows,1:BC_cols) = BC_F2C;
+IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IC2F using symmetry condition %%%%
+IC2F = Hf\IF2C.'*Hc;
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3
+%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:36:01
+
+stencil_F2C = [7.0./2.56e2,-1.0./3.2e1,-7.0./6.4e1,9.0./3.2e1,8.5e1./1.28e2,9.0./3.2e1,-7.0./6.4e1,-1.0./3.2e1,7.0./2.56e2];
+if nargout > 1
+    BC_F2C = reshape([7.523257802630956e-2,2.447812262221267e-1,-1.679313063616916e-1,1.290510950666589,6.315723344677289e-3,1.67178747937954e-1,1.982667903557025,-7.554383893379468e-1,7.215271362899867e-1,-2.820478831383137,1.807034411305536,-7.589683751979843e-1,-7.685973268458095e-1,3.965751544535173e-1,4.119789638051451e-1,1.574000556785898e-1,-1.113618964927466e-1,3.631000220124657e-1,1.639694219982991,-9.114074742274862e-1,4.952715520862987e-1,-4.524162151353456e-1,2.65481378213589e-1,-2.268273408674622e-1,3.455365956773523e-1,-2.009765975089827e-1,1.722179564305811e-1,-5.52933488235368e-1,3.18109976271229e-1,-2.171481232558432e-1,1.033835580108027e-1,-5.911351224351343e-2,3.960076712054991e-2],[3,11]);
+end
+if nargout > 2
+    t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2
+%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:36:05
+
+stencil_F2C = [-1.0./2.56e2,-1.0./3.2e1,1.0./6.4e1,9.0./3.2e1,6.1e1./1.28e2,9.0./3.2e1,1.0./6.4e1,-1.0./3.2e1,-1.0./2.56e2];
+if nargout > 1
+    BC_F2C = reshape([1.0./2.0,0.0,0.0,1.77e2./2.72e2,3.0./8.0,-5.9e1./6.88e2,1.125919117647059e-2,3.546742584745763e-1,1.335392441860465e-2,-4.9e1./5.44e2,2.335805084745763e-1,3.204941860465116e-1,-1.194852941176471e-2,1.350635593220339e-2,5.308866279069767e-1,-9.0./5.44e2,-1.11228813559322e-2,2.943313953488372e-1,-2.803308823529412e-2,2.105402542372881e-2,-1.580668604651163e-2,-9.0./5.44e2,1.430084745762712e-2,-5.450581395348837e-2,-9.191176470588235e-4,7.944915254237288e-4,-5.450581395348837e-3,1.0./5.44e2,-1.588983050847458e-3,2.180232558139535e-3,2.297794117647059e-4,-1.986228813559322e-4,2.725290697674419e-4],[3,11]);
+end
+if nargout > 2
+    t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,47 @@
+function [IC2F,IF2C,Hc,Hf] = IntOp_orders_6to6_ratio2to1(mc,hc,ACC)
+
+% ACC is a string.
+% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C.
+% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F.
+ratio = 2;
+mf = ratio*mc-1; 
+hf = hc/ratio;
+
+switch ACC
+    case 'F2C'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4;
+    case 'C2F'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3;
+end
+
+stencil_width = length(stencil_F2C);
+stencil_hw = (stencil_width-1)/2;
+[BC_rows,BC_cols] = size(BC_F2C);
+
+%%% Norm matrices %%%
+Hc = speye(mc,mc);
+HcUm = length(HcU); 
+Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm);
+Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm);
+Hc = Hc*hc;
+
+Hf = speye(mf,mf);
+HfUm = length(HfU);
+Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm);
+Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm);
+Hf = Hf*hf;
+%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IF2C from stencil and BC
+IF2C = sparse(mc,mf);
+for i = BC_rows+1 : mc-BC_rows
+    IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX>
+end
+IF2C(1:BC_rows,1:BC_cols) = BC_F2C;
+IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IC2F using symmetry condition %%%%
+IC2F = Hf\IF2C.'*Hc;
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4
+%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:35:47
+
+stencil_F2C = [-5.95703125e-3,3.0./5.12e2,3.57421875e-2,-2.5e1./5.12e2,-8.935546875e-2,7.5e1./2.56e2,3.17e2./5.12e2,7.5e1./2.56e2,-8.935546875e-2,-2.5e1./5.12e2,3.57421875e-2,3.0./5.12e2,-5.95703125e-3];
+if nargout > 1
+    BC_F2C = reshape([5.233890618131365e-1,-1.594459192645467e-2,3.532688727637403e-2,8.021234957689208e-1,3.90683205173562e-1,-1.73222951632239e-1,-8.87662686483442e-2,2.90091235796637e-1,-1.600356115709148e-1,-1.044025375027475e-1,2.346179009198368e-1,6.091329306528956e-1,1.561275522703128e-1,-1.168382445709856e-1,1.040987887887311,-2.387061036980731e-1,1.363504965974361e-1,1.082611928255256e-1,-5.745310654054326e-1,3.977694249198785e-1,-9.376217911402619e-1,3.554518646054656e-2,5.609157787396987e-5,-1.564625311232018e-1,6.107907974027401e-1,-3.786608696698368e-1,7.02265869951125e-1,2.054294270642538e-1,-1.302300112378257e-1,2.478941407690889e-1,-2.657085326191479e-1,1.568761445933572e-1,-2.632906518005349e-1,-1.228047556139644e-1,7.182193248980271e-2,-1.1291238242346e-1,5.811258780405158e-2,-3.466364400805378e-2,5.683680203338252e-2,1.781337575097077e-2,-1.030572999042704e-2,1.577197067502767e-2,-1.443802115768554e-2,8.391316308374261e-3,-1.29534117369052e-2,-1.548018627738296e-3,8.794183904936319e-4,-1.298961407229805e-3,1.573818938200601e-3,-8.940753636685258e-4,1.320610764016968e-3],[3,17]);
+end
+if nargout > 2
+    t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3
+%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:35:56
+
+stencil_F2C = [1.07421875e-3,3.0./5.12e2,-6.4453125e-3,-2.5e1./5.12e2,1.611328125e-2,7.5e1./2.56e2,2.45e2./5.12e2,7.5e1./2.56e2,1.611328125e-2,-2.5e1./5.12e2,-6.4453125e-3,3.0./5.12e2,1.07421875e-3];
+if nargout > 1
+    BC_F2C = reshape([1.0./2.0,0.0,0.0,0.0,6.876076086160158e-1,1.5e1./3.2e1,-3.461879841386942e-1,3.502577439820862e-2,3.099721991995751e-3,2.228547003766753e-1,9.363652999354482e-3,-3.15791046603844e-3,-1.05789143206462e-1,2.35563165945226e-1,6.070385985337514e-1,-4.847496617839149e-2,-4.809232246867902e-3,5.154694068405061e-3,7.049223649022501e-1,1.016749349808733e-2,3.460117842882262e-2,-4.996621498584866e-2,5.210650763094799e-1,2.233592875303228e-1,-2.239024779745769e-3,4.254668119953384e-4,9.213872590833641e-3,3.910524351558127e-1,-9.048711215107334e-2,8.597375629526346e-2,-3.468219752858724e-1,3.250682992395969e-1,-1.309107466664224e-1,1.199249722306043e-1,-4.071552384959425e-1,1.474417123938701e-1,-3.02863877390285e-2,3.046016554982103e-2,-1.081315128181483e-1,1.120705238850532e-2,7.977722870036266e-2,-6.772545887788229e-2,2.00270418837606e-1,-5.427183680490763e-2,6.524845570188292e-2,-5.637610037043203e-2,1.711235764016968e-1,-4.203646902407166e-2,4.639548341453586e-3,-4.055884445496545e-3,1.258630924935448e-2,-2.776451559292779e-3,-1.023784366070774e-2,8.817108288936985e-3,-2.659664906860937e-2,7.14445034054861e-3,-1.435466025441424e-3,1.240274208149505e-3,-3.764718680837329e-3,1.028716295017727e-3,1.032012418492197e-3,-8.794183904936319e-4,2.597922814459609e-3,-6.571159498040679e-4,1.892022767235695e-4,-1.612267049238325e-4,4.76285849317595e-4,-1.204712574640791e-4],[4,17]);
+end
+if nargout > 2
+    t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,47 @@
+function [IC2F,IF2C,Hc,Hf] = IntOp_orders_8to8_ratio2to1(mc,hc,ACC)
+
+% ACC is a string.
+% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C.
+% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F.
+ratio = 2;
+mf = ratio*mc-1; 
+hf = hc/ratio;
+
+switch ACC
+    case 'F2C'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5;
+    case 'C2F'
+        [stencil_F2C,BC_F2C,HcU,HfU] = ...
+        sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4;
+end
+
+stencil_width = length(stencil_F2C);
+stencil_hw = (stencil_width-1)/2;
+[BC_rows,BC_cols] = size(BC_F2C);
+
+%%% Norm matrices %%%
+Hc = speye(mc,mc);
+HcUm = length(HcU); 
+Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm);
+Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm);
+Hc = Hc*hc;
+
+Hf = speye(mf,mf);
+HfUm = length(HfU);
+Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm);
+Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm);
+Hf = Hf*hf;
+%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IF2C from stencil and BC
+IF2C = sparse(mc,mf);
+for i = BC_rows+1 : mc-BC_rows
+    IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX>
+end
+IF2C(1:BC_rows,1:BC_cols) = BC_F2C;
+IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Create IC2F using symmetry condition %%%%
+IC2F = Hf\IF2C.'*Hc;
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5
+%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:35:02
+
+stencil_F2C = [1.324608212425595e-3,-1.220703125e-3,-1.059686569940476e-2,1.1962890625e-2,3.708902994791667e-2,-5.9814453125e-2,-7.417805989583333e-2,2.99072265625e-1,5.927225748697917e-1,2.99072265625e-1,-7.417805989583333e-2,-5.9814453125e-2,3.708902994791667e-2,1.1962890625e-2,-1.059686569940476e-2,-1.220703125e-3,1.324608212425595e-3];
+if nargout > 1
+    BC_F2C = reshape([2.087365925359584,-1.227221822236823,1.090916714284468e1,-1.041312149906314,1.134214373258162,-1.269686811479259,2.075368250436277,-1.520771409696474e1,1.389750473974189,-1.48485748783569,-1.040579640776707,8.899721508624742e-1,-7.177691661032869,6.882712043721976e-1,-7.599347599660409e-1,-1.26440526850463,1.160658043364166,-5.391509178445742,6.679923135342851e-1,-7.521725463922262e-1,3.752349810676949,-2.906684049793639,2.672876913566556e1,-2.493145660873,2.782825829667436,1.443588084644871e-1,-1.307230271348211e-1,2.216804748506809,1.437719804483573e-1,-1.145830952112369e-1,3.2149320801083e-1,-2.34361439272989e-1,1.855837403850803,1.286992417665136e-1,-5.758238484341108e-2,-1.301103772185027,9.973408313121463e-1,-8.837520107531879,9.534930382604288e-1,-1.429199518808459e-2,-2.384453600787036,1.818345997558567,-1.577342441239303e1,1.422846302346762,3.154159322410927e-3,3.589377915408905e-1,-2.109524979864723e-1,9.363227187483732e-1,6.301238593470628e-2,5.732390257964316e-2,2.793582225299658,-2.033528056927806,1.621324028555783e1,-1.189222718426265,6.189671387692164e-2,5.342131492864642e-1,-4.222231296997605e-1,3.787253511209165,-3.350064226195165e-1,1.245325350855226e-2,-2.184146687275183,1.551331905163865,-1.19387312118926e1,8.349589004472998e-1,-1.35383591594438e-1,-8.014367493869704e-1,5.668415377001197e-1,-4.302842643863972,2.84183319528176e-1,3.83514596570461e-2,9.45794880458861e-1,-6.732638898386457e-1,5.234390273661413,-3.83498639699941e-1,1.440453495443018e-1,5.645065377027613e-1,-4.039220855868618e-1,3.170156487564899,-2.383074337879712e-1,1.216617375214008e-1,-1.937863224145611e-1,1.358891032584724e-1,-1.023378722711107,6.830621267493578e-2,-4.169135139842575e-3,2.270306426219975e-2,-2.38399209865252e-2,2.928388757260058e-1,-4.02060181933195e-2,6.880595419697505e-2,9.638828970166005e-2,-6.98778025332229e-2,5.609317796857676e-1,-4.429526993510467e-2,2.882554468024363e-2,3.582431391759518e-2,-2.671910644265286e-2,2.251335440088556e-1,-1.966486325944329e-2,1.791756577091828e-2,-3.350790201878844e-1,2.581411034605136e-1,-2.283079071439122,2.162064132948073e-1,-2.321676317817421e-1,6.356644993195555e-2,-4.921907141164279e-2,4.384534411523264e-1,-4.198474091936761e-2,4.597203884996652e-2,2.261662512481835e-2,-1.740429407604355e-2,1.537038474929879e-1,-1.452709057733876e-2,1.556101844926456e-2,3.097679325854209e-2,-2.394872918869654e-2,2.128879105995857e-1,-2.032077838507769e-2,2.213372706946953e-2],[5,24]);
+end
+if nargout > 2
+    t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,18 @@
+function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4
+%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24
+%    [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24
+
+%    This function was generated by the Symbolic Math Toolbox version 8.0.
+%    21-May-2018 15:35:39
+
+stencil_F2C = [-2.564929780505952e-4,-1.220703125e-3,2.051943824404762e-3,1.1962890625e-2,-7.181803385416667e-3,-5.9814453125e-2,1.436360677083333e-2,2.99072265625e-1,4.820454915364583e-1,2.99072265625e-1,1.436360677083333e-2,-5.9814453125e-2,-7.181803385416667e-3,1.1962890625e-2,2.051943824404762e-3,-1.220703125e-3,-2.564929780505952e-4];
+if nargout > 1
+    BC_F2C = reshape([-1.673327822994457e1,1.665420585653232e1,-1.973927473454954e2,2.826257908914756e1,-6.156813484052936e1,3.974966126696054,2.880639373222355e1,-2.660797294399895e1,3.202303847857831e2,-4.598960974033412e1,1.003152507870138e2,-6.481221721577338,8.061874893005659,-7.706608975219419,9.234201483040626e1,-1.322147613066874e1,2.880209985918895e1,-1.859523138300886,-2.19528557898034e1,2.13763239391502e1,-2.476320412310993e2,3.572924955195535e1,-7.795290920161567e1,5.036097395109802,-1.287197329845376,1.244100160910233,-1.394684346537857e1,2.112320869857597,-4.60441856527933,2.978264879234002e-1,1.935199269380359,-1.88580031845933,2.330807783420723e1,-2.91749306966336,6.644427567385864,-4.302395355488564e-1,-1.148984264037686,1.109982778464842,-1.314799057941425e1,2.137017065505111,-4.084082382904703,2.606486561845403e-1,1.565462612018858,-1.508309105065004,1.772442871757061e1,-2.40060893717649,6.31422909250156,-4.052444871368164e-1,3.179380137695096e-1,-3.071741522636364e-1,3.636264382483117,-5.183996466001354e-1,2.322192602480252,-6.467931586900172e-2,6.277081831671831e-1,-6.107890751034906e-1,7.335850436798604,-1.091016872160133,3.089218691662174,6.988892111841499e-2,1.639756762532723,-1.583883992641009,1.876236237708381e1,-2.685390559486947,5.859855615166555,3.919312696253764e-3,1.611803934956753,-1.540778869185931,1.795955382224485e1,-2.495154500270907,5.003492494640851,-4.157919873785453e-2,-6.213413520816026e-1,6.242581523412814e-1,-7.823535324240331,1.222418137438074,-3.120142966701915,2.975831327541895e-1,2.12230710100161,-2.0451058741188,2.412395225004063e1,-3.423544496956294,7.325914828631551,-4.793445914188063e-1,-4.31809447483846,4.158975059021769,-4.906211424150023e1,6.975058620124194,-1.501259943091044e1,9.418983630154896e-1,-6.138797742037704e-1,5.810969368378602e-1,-6.6821474226667,9.113864351863538e-1,-1.814400062103318,1.002574935832871e-1,-9.860044240994036e-2,9.448685284232805e-2,-1.106188847377502,1.552858946803229e-1,-3.265426230113952e-1,2.062042480203908e-2,-2.745034502159738,2.655226942873441,-3.151193386341271e1,4.520932327069998,-9.883338555855938,6.423439846307804e-1,2.139968493382333,-2.067738462547628,2.450223844197562e1,-3.50699574523559,7.635005983192545,-4.923960464980415e-1,1.036406872394752,-1.002008838408551,1.188339891704641e1,-1.70302259679446,3.715789696036691,-2.407301455417931e-1,1.606428513214277,-1.552524985439685,1.840245870449254e1,-2.635131753653705,5.741508438443495,-3.708346905830083e-1,5.659562678739624e-2,-5.465656201303744e-2,6.471932864664631e-1,-9.253169145609999e-2,2.010909570941804e-1,-1.292046400706276e-2,5.435391241988039e-1,-5.252672844681129e-1,6.225561585258754,-8.91344337215917e-1,1.941632399616401,-1.253426955105431e-1,-1.552116972709218,1.499962759958308,-1.777819805127299e1,2.545472086708348,-5.545140384142844,3.580057322157566e-1],[6,24]);
+end
+if nargout > 2
+    t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158];
+    HcU = t2;
+end
+if nargout > 3
+    HfU = t2;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpMC_orders_2to2_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,35 @@
+function [IC2F,IF2C] = intOpMC_orders_2to2_ratio2to1(mc)
+
+mf = 2*(mc-1) + 1;
+
+stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0];
+stencil_width = length(stencil_F2C);
+stencil_halfwidth = (stencil_width-1)/2;
+
+BC_F2C = [1.0./2.0,1.0./2.0];
+
+Hc = speye(mc,mc);
+Hc(1,1) = 1/2;
+Hc(end,end) = 1/2;
+
+Hf = 1/2*speye(mf,mf);
+Hf(1,1) = 1/4;
+Hf(end,end) = 1/4; 
+
+IF2C = sparse(mc,mf);
+[BCrows, BCcols] = size(BC_F2C);
+IF2C(1:BCrows, 1:BCcols) = BC_F2C;
+IF2C(mc-BCrows+1:mc, mf-BCcols+1:mf) = rot90(BC_F2C,2);
+
+for i = BCrows+1 : mc-BCrows
+	IF2C(i,(2*i-stencil_halfwidth-1) :(2*i+stencil_halfwidth-1))...
+		 = stencil_F2C;
+end
+
+IC2F = Hf\(IF2C'*Hc);
+
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpMC_orders_4to4_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,66 @@
+% Marks New interpolation operators
+% 4th order to 2nd order accurate (diagonal norm) 
+% M=9 is the minimum amount of points on the coarse mesh
+
+function [IC2F,IF2C] = intOpMC_orders_4to4_ratio2to1(M_C)
+
+M_F=M_C*2-1;
+
+% Coarse to fine
+I1=zeros(M_F,M_C);
+t1=[  2047/2176 , 129/1088 , -129/2176 , 0 , 0 , 0 , 0 ,...
+      0  ;...
+    429/944 , 279/472 , -43/944 , 0 , 0 , 0 , 0 , 0 ;...
+    111/2752 , 4913/5504 , 3/32 , -147/5504 ,...
+      0 , 0 , 0 , 0  ;...
+    -103/784 , 549/784 , 387/784 , -1/16 , 0 ,...
+      0 , 0 , 0  ;...
+    -335/3072 , 205/768 , 2365/3072 , 49/512 ,...
+      -3/128 , 0 , 0 , 0  ;...
+    -9/256 , 5/256 , 129/256 , 147/256 ,...
+      -1/16 , 0 , 0 , 0  ;...
+    5/192 , -59/1024 , 43/512 , 2695/3072 ,...
+      3/32 , -3/128 , 0 , 0 ;...
+    23/768 , -37/768 , -43/768 , 147/256 ,...
+      9/16 , -1/16 , 0 , 0  ;...
+    13/2048 , -11/1024 , -43/2048 , 49/512 ,...
+      55/64 , 3/32 , -3/128 , 0 ;...
+    -1/384 , 1/256 , 0 , -49/768 , 9/16 ,...
+      9/16 , -1/16 , 0  ;...
+    -1/1024 , 3/2048 , 0 , -49/2048 , 3/32 ,...
+      55/64 , 3/32 , -3/128];
+
+
+t2=[-3/128 3/32 55/64 3/32 -3/128];
+t3=[-1/16 9/16 9/16 -1/16];
+I1(1:11,1:8)=t1;
+I1(M_F-10:M_F,M_C-7:M_C)=fliplr(flipud(t1));
+I1(12,5:8)=t3;
+for i=13:2:M_F-12
+    j=(i-1)/2;
+    I1(i,j-1:j+3)=t2;
+    I1(i+1,j:j+3)=t3;
+end
+
+% Fine to coarse
+I2=zeros(M_C,M_F);
+
+t1=[ 2047/4352 , 429/544 , 111/2176 , -103/544 ...
+      , -335/2176 , -27/544 , 5/136 , 23/544 ,...
+      39/4352 , -1/272 , -3/2176 ;...
+    129/7552 , 279/944 , 4913/15104 , 549/1888 ...
+      , 205/1888 , 15/1888 , -3/128 , -37/1888 ...
+      , -33/7552 , 3/1888 , 9/15104 ];
+t2=[-3/256 -1/32 3/64 9/32 55/128 9/32 3/64 -1/32 -3/256];
+
+I2(1:2,1:11)=t1;      
+
+I2(M_C-1:M_C,M_F-10:M_F)=fliplr(flipud(t1));
+
+for i=3:M_C-2
+    j=2*(i-3)+1;
+    I2(i,j:j+8)=t2;
+end
+
+IC2F = sparse(I1);
+IF2C = sparse(I2);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpMC_orders_6to6_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,115 @@
+% Marks New interpolation operators
+% 6th order accurate (diagonal norm) 
+% M=19 is the minimum amount of points on the coarse mesh
+
+function [IC2F,IF2C] = intOpMC_orders_6to6_ratio2to1(M_C)
+
+M_F=M_C*2-1;
+
+% Coarse to fine
+I1=zeros(M_F,M_C);
+
+t1=    [6854313/6988288 , 401925/6988288 , -401925/6988288 ...
+      , 133975/6988288 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0  ; ...
+    560547/1537664 , 1201479/1537664 , -240439/1537664 ...
+      , 16077/1537664 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ...
+    203385/5552128 , 1225647/1388032 , 364155/2776064 ...
+      , -80385/1388032 , 39385/5552128 , 0 , 0 , 0 , 0 , 0 , 0 , 0  ; ...
+    -145919/2743808 , 721527/1371904 , 105687/171488 , ...
+      -25/256 , 23631/2743808 , 0 , 0 , 0 , 0 , 0 , 0 , 0  ; ...
+    -178863/4033024 , 1178085/8066048 , ...
+      1658587/2016512 , 401925/4033024 , -15/512 , ...
+      43801/8066048 , 0 , 0 , 0 , 0 , 0 , 0 ; ...
+    -1668147/11213056 , 4193225/11213056 , ...
+      375675/2803264 , 2009625/2803264 , ...
+      -984625/11213056 , 3/256 , 0 , 0 , 0 , 0 , 0 , 0 ; ...
+    -561187/2949120 , 831521/1474560 , -788801/1474560 ...
+      , 412643/368640 , 39385/589824 , -43801/1474560 ...
+      , 5/1024 , 0 , 0 , 0 , 0 , 0 ; ...
+    23/1024 , 23/147456 , -43435/221184 , ...
+      26795/36864 , 39385/73728 , -43801/442368 , ...
+      3/256 , 0 , 0 , 0 , 0 , 0  ; ...
+    79379/368640 , -1664707/2949120 , 284431/737280 , ...
+      26795/294912 , 606529/737280 , 43801/589824 , ...
+      -15/512 , 5/1024 , 0 , 0 , 0 , 0 ; ...
+    3589/27648 , -2225/6144 , 22939/73728 , ...
+      -26795/221184 , 39385/73728 , 43801/73728 , ...
+      -25/256 , 3/256 , 0 , 0 , 0 , 0 ; ...
+    -720623/14745600 , 10637/92160 , -89513/1474560 , ...
+      -5359/147456 , 39385/589824 , 3372677/3686400 , ...
+      75/1024 , -15/512 , 5/1024 , 0 , 0 , 0 ; ...
+    -6357/81920 , 55219/276480 , -8707/61440 , ...
+      5359/368640 , -39385/442368 , 43801/73728 , ...
+      75/128 , -25/256 , 3/256 , 0 , 0 , 0 ; ...
+    -13315/884736 , 2589/65536 , -479/16384 , ...
+      5359/884736 , -7877/294912 , 43801/589824 , ...
+      231/256 , 75/1024 , -15/512 , 5/1024 , 0 , 0  ; ...
+    8299/737280 , -7043/245760 , 5473/276480 , 0 , ...
+      7877/737280 , -43801/442368 , 75/128 , ...
+      75/128 , -25/256 , 3/256 , 0 , 0 ; ...
+    11027/2949120 , -8461/884736 , 655/98304 , 0 , ...
+      7877/1769472 , -43801/1474560 , 75/1024 , ...
+      231/256 , 75/1024 , -15/512 , 5/1024 , 0 ; ...
+    -601/614400 , 601/245760 , -601/368640 , 0 , 0 , ...
+      43801/3686400 , -25/256 , 75/128 , ...
+      75/128 , -25/256 , 3/256 , 0 ; ...
+    -601/1474560 , 601/589824 , -601/884736 , 0 , 0 , ...
+      43801/8847360 , -15/512 , 75/1024 , ...
+      231/256 , 75/1024 , -15/512 , 5/1024 ] ;
+  
+  t2=  [5/1024 , -15/512 , 75/1024, 231/256 , 75/1024 , -15/512 , 5/1024];
+      
+  t3=  [3/256 , -25/256 , 75/128 , 75/128 , -25/256 , 3/256];     
+
+I1(1:17,1:12)=t1;
+I1(M_F-16:M_F,M_C-11:M_C)=fliplr(flipud(t1));
+I1(18,7:12)=t3;
+for i=19:2:M_F-18
+    j=(i-3)/2;
+    I1(i,j-1:j+5)=t2;
+    I1(i+1,j:j+5)=t3;
+end
+
+% Fine to coarse
+I2=zeros(M_C,M_F);
+
+t1=[6854313/13976576 , 2802735/3494144 , ...
+      1016925/27953152 , -729595/6988288 , ...
+      -894315/13976576 , -1668147/6988288 , ...
+      -8417805/27953152 , 15525/436768 , ...
+      1190685/3494144 , 89725/436768 , ...
+      -2161869/27953152 , -858195/6988288 , ...
+      -332875/13976576 , 124485/6988288 , ...
+      165405/27953152 , -5409/3494144 , -9015/13976576;  ...
+    80385/12301312 , 1201479/3075328 , 1225647/6150656 ...
+      , 721527/3075328 , 1178085/24602624 , ...
+      838645/6150656 , 60843/300032 , 345/6150656 , ...
+      -4994121/24602624 , -100125/768832 , ...
+      31911/768832 , 55219/768832 , 349515/24602624 , ...
+      -63387/6150656 , -42305/12301312 , 5409/6150656 ...
+      , 9015/24602624  ; ...
+    -80385/5552128 , -240439/1388032 , 364155/5552128 ...
+      , 105687/173504 , 1658587/2776064 , 75135/694016 ...
+      , -2366403/5552128 , -217175/1388032 , ...
+      853293/2776064 , 344085/1388032 , ...
+      -268539/5552128 , -78363/694016 , -64665/2776064 ...
+      , 5473/347008 , 29475/5552128 , -1803/1388032 , ...
+      -3005/5552128];
+
+
+  
+ t2=[  5/2048 , 3/512 , -15/1024 , -25/512 , ...
+      75/2048 , 75/256 , 231/512 , 75/256 , ...
+      75/2048 , -25/512 , -15/1024 , 3/512 , 5/2048];
+  
+I2(1:3,1:17)=t1;      
+
+I2(M_C-2:M_C,M_F-16:M_F)=fliplr(flipud(t1));
+
+for i=4:M_C-3
+    j=2*(i-4)+1;
+    I2(i,j:j+12)=t2;
+end
+IC2F = sparse(I1);
+IF2C = sparse(I2);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/intOpMC_orders_8to8_ratio2to1.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,53 @@
+% Marks New interpolation operators
+% 8th order accurate (diagonal norm) 
+% M=19 is the minimum amount of points on the coarse mesh
+
+function [IC2F,IF2C] = intOpMC_orders_8to8_ratio2to1(M_C)
+
+M_F=M_C*2-1;
+
+% Coarse to fine
+I1=zeros(M_F,M_C);
+
+t1=    [0.99340647972821530406e0 0.26374081087138783775e-1 -0.39561121630708175663e-1 0.26374081087138783775e-1 -0.65935202717846959438e-2 0 0 0 0 0 0 0 0 0 0 0; 0.31183959860287729600e0 0.94014160558849081601e0 -0.31646240838273622401e0 0.65141605588490816007e-1 -0.66040139712270400169e-3 0 0 0 0 0 0 0 0 0 0 0; -0.33163591467432411328e-1 0.11092588191512759415e1 -0.10539936193077965253e0 -0.77189144409925778387e-2 0.60418595406382404105e-1 -0.23395546718453703858e-1 0 0 0 0 0 0 0 0 0 0; -0.63951987538041216890e-1 0.56657207530367360435e0 0.56073157416571775151e0 -0.67107298938782711711e-1 0.54915118559238359570e-2 -0.17358748484912632117e-2 0 0 0 0 0 0 0 0 0 0; 0.22970968995770386894e0 -0.84424237330911973043e0 0.20693327723706358111e1 -0.42910115554542331920e0 -0.13191478386190563918e0 0.11675567167691342983e0 -0.10539821288804421121e-1 0 0 0 0 0 0 0 0 0; 0.10195433776615307267e0 -0.45344410547614678949e0 0.11049699103276728065e1 0.26297465812771521534e0 -0.38617448080010680341e-1 0.23925781250000000000e-1 -0.17631339153836246986e-2 0 0 0 0 0 0 0 0 0; -0.33882356049961665258e0 0.12718893449497745294e1 -0.16822328941798626751e1 0.17813590728082409984e1 0.11793036905675500906e0 -0.18266200597575250398e0 0.37689938246258754051e-1 -0.51502644057974593120e-2 0 0 0 0 0 0 0 0; -0.50400650482311136546e0 0.19901277318227816880e1 -0.29708252195230746436e1 0.23722122500767090037e1 0.24457622727717445252e0 -0.15152936311741758892e0 0.21886284536938453771e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0 0; 0.37882732714731866331e-2 0.12115606818363056270e0 -0.53924884156527999826e0 0.88886598075786086512e0 0.27660232216640139169e0 0.33730204543181760125e0 -0.12179633685076087365e0 0.38041730885639608773e-1 -0.47112422807823442564e-2 0 0 0 0 0 0 0; 0.41123781086486062886e0 -0.14418811327145423418e1 0.16793759240044487847e1 -0.57156511000645013503e0 0.24685906775203751929e0 0.76471858554416232639e0 -0.11045284035765094522e0 0.24149101163134162809e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0; 0.36463669929508579110e0 -0.13698743320477484358e1 0.18268133133896437122e1 -0.93074264690533336964e0 0.10888458847499176143e0 0.85685706622610101431e0 0.24359267370152174731e0 -0.13314605809973863070e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0 0; 0.25541802394537278164e0 -0.10497853549910180363e1 0.15921730465359157986e1 -0.96615555845660927855e0 -0.49371813550407503858e-1 0.76471858554416232639e0 0.55226420178825472608e0 -0.12074550581567081404e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0 0; -0.93469656691661424206e-1 0.26634337856674539612e0 -0.17694629839462736800e0 -0.64947940656920645005e-1 -0.54442294237495880713e-1 0.33730204543181760125e0 0.61880473767909428853e0 0.26629211619947726141e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0; -0.49445400002561969522e0 0.18168098496802059221e1 -0.23462477366129557292e1 0.11091140417405116705e1 0.98743627100815007716e-2 -0.15294371710883246528e0 0.55226420178825472608e0 0.60372752907835407022e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0; -0.24633538738293361369e0 0.93021448723433316305e0 -0.12527182680719820488e1 0.63698038058706249354e0 0.15554941210713108775e-1 -0.16865102271590880063e0 0.24359267370152174731e0 0.67646871560981190145e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0; 0.23150628239599695492e0 -0.82682792555928440531e0 0.10237569354829334077e1 -0.45129113643047460593e0 -0.10075880316409694665e-2 0.30588743421766493056e-1 -0.11045284035765094522e0 0.60372752907835407022e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0; 0.19518658723021413499e0 -0.70515100787561674409e0 0.88870680011413432419e0 -0.40458631782898657259e0 -0.19443676513391385969e-2 0.48186006490259657322e-1 -0.12179633685076087365e0 0.26629211619947726141e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0; -0.43403699494753775353e-1 0.15442805550926787092e0 -0.18997683853068679729e0 0.82584189359473172950e-1 0 -0.31213003491598462302e-2 0.22090568071530189043e-1 -0.12074550581567081404e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0; -0.58783367481931040778e-1 0.20994477957758583150e0 -0.25976152530269196017e0 0.11403438083569734975e0 0 -0.60232508112824571652e-2 0.34798953385931678187e-1 -0.13314605809973863070e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0; 0.63390647713259204145e-2 -0.22373993350504987875e-1 0.27185871992161665013e-1 -0.11561529976981026786e-1 0 0 -0.22541395991357335758e-2 0.24149101163134162809e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0; 0.10649583863025939259e-1 -0.37634916628131671889e-1 0.45812371547331598336e-1 -0.19540204529147604911e-1 0 0 -0.43498691732414597734e-2 0.38041730885639608773e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0; -0.45575492476359756866e-3 0.15951422366725914903e-2 -0.19141706840071097884e-2 0.79757111833629574515e-3 0 0 0 -0.24641939962381798784e-2 0.23925781250000000000e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0; -0.87948159845213680356e-3 0.30781855945824788125e-2 -0.36938227134989745750e-2 0.15390927972912394062e-2 0 0 0 -0.47552163607049510966e-2 0.37689938246258754051e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2;];
+
+  
+  t2=  [-0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.1856585148354920384923865861096125662293072684152233190798249652677391616531107981e82 / 0.2770141785660523684734143913470637620565185365865679672509451726827594519874538713e82 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83;];
+
+  t3=  [-0.5e1 / 0.2048e4 0.49e2 / 0.2048e4 -0.245e3 / 0.2048e4 0.1225e4 / 0.2048e4 0.1225e4 / 0.2048e4 -0.245e3 / 0.2048e4 0.49e2 / 0.2048e4 -0.5e1 / 0.2048e4;];
+I1(1:23,1:16)=t1;
+I1(M_F-22:M_F,M_C-15:M_C)=fliplr(flipud(t1));
+I1(24,9:16)=t3;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DEBUGGING
+
+for i=25:2:M_F-24
+    j=(i-3)/2;
+    I1(i,j-2:j+6)=t2;
+    I1(i+1,j-1:j+6)=t3;
+end
+
+% Fine to coarse: A
+I2=zeros(M_C,M_F);
+
+t1=[0.49670323986410765203e0 0.80670591743192512511e0 -0.14476656476698073533e-1 -0.19497555250083395132e0 0.16074268813765884450e0 0.22100911221277072755e0 -0.53042418939472250452e0 -0.86254139670456354517e0 0.64231825172866668299e-2 0.69727164011248914487e0 0.61825742343094006838e0 0.43307239695721032895e0 -0.15848187861199173328e0 -0.83836831742920925562e0 -0.41767239062238727391e0 0.39252899650233765024e0 0.33094736964907844809e0 -0.73592865087013788440e-1 -0.99669762780958210515e-1 0.10748160731101219580e-1 0.18056833808814782786e-1 -0.77275234787125894193e-3 -0.14911993994710636483e-2; 0.25487859584304862664e-2 0.47007080279424540800e0 0.93589176759289320118e-1 0.33386224041716468423e0 -0.11418395501435496457e0 -0.18998278818813064037e0 0.38484431284491986603e0 0.65828018436035862232e0 0.39704539050575728856e-1 -0.47252462545568042557e0 -0.44892698918501097924e0 -0.34402935194949605213e0 0.87284452472801643396e-1 0.59539401290875351190e0 0.30484430526276345964e0 -0.27096308216867871783e0 -0.23108785344796321535e0 0.50608234918774219796e-1 0.68801842319351676214e-1 -0.73322707316320135247e-2 -0.12333488857215226757e-1 0.52275043402033040521e-3 0.10087644967132781708e-2; -0.22656973277393401539e-1 -0.93771184228725468159e0 -0.52699680965389826267e-1 0.19581430555012001670e1 0.16586148101136731602e1 0.27435837109255836264e1 -0.30164708462284474624e1 -0.58235016129647971089e1 -0.10472767830023645070e1 0.32615209891555982371e1 0.35478595826728208891e1 0.30921640208240511054e1 -0.34364793368678654581e0 -0.45566547247355317663e1 -0.24329078834671892590e1 0.19882413967858906122e1 0.17259601262271516763e1 -0.36895458453625989435e0 -0.50448363278283993228e0 0.52797763052066775846e-1 0.88972343374038343284e-1 -0.37175165926095403240e-2 -0.71737841052106668688e-2; 0.21626748627943672418e-2 0.27636709246281031633e-1 -0.55259484658035070394e-3 -0.33553649469391355855e-1 -0.49244245327794891233e-1 0.93489376222103426899e-1 0.45734620580442859300e0 0.66579609152388478842e0 0.24716623315221892947e0 -0.15893464065423852815e0 -0.25881084331023040874e0 -0.26865808253702445366e0 -0.18060020510041282529e-1 0.30841043055726240020e0 0.17712461121229459766e0 -0.12549015561536110372e0 -0.11250298507032009025e0 0.22964117700293735857e-1 0.31709446610808019222e-1 -0.32149051440245356510e-2 -0.54335286230388551025e-2 0.22177994574852164638e-3 0.42797426992744435493e-3;];
+
+
+  
+ t2=[-0.23556211403911721282e-2 -0.12207031250000000000e-2 0.18844969123129377026e-1 0.11962890625000000000e-1 -0.65957391930952819589e-1 -0.59814453125000000000e-1 0.13191478386190563918e0 0.29907226562500000000e0 0.33510652017261795103e0 0.29907226562500000000e0 0.13191478386190563918e0 -0.59814453125000000000e-1 -0.65957391930952819589e-1 0.11962890625000000000e-1 0.18844969123129377026e-1 -0.12207031250000000000e-2 -0.23556211403911721282e-2;];
+
+  
+I2(1:4,1:23)=t1;      
+
+I2(M_C-3:M_C,M_F-22:M_F)=fliplr(flipud(t1));
+
+for i=5:M_C-4
+    j=2*(i-5)+1;
+    I2(i,j:j+16)=t2;
+end
+
+IC2F = sparse(I1);
+IF2C = sparse(I2);
+
+
--- a/+sbp/D2Variable.m	Fri Jun 15 18:10:26 2018 -0700
+++ b/+sbp/D2Variable.m	Sat Jun 16 14:30:45 2018 -0700
@@ -26,22 +26,39 @@
             obj.x = linspace(x_l,x_r,m)';
 
             switch order
+
+                case 6
+
+                    [obj.H, obj.HI, obj.D1, obj.D2, ...
+                    ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,...
+                     obj.d1_l, obj.d1_r] = ...
+                        sbp.implementations.d4_variable_6(m, obj.h);
+                    obj.borrowing.M.d1 = 0.1878;
+                    obj.borrowing.R.delta_D = 0.3696;
+                    % Borrowing e^T*D1 - d1 from R
+
                 case 4
                     [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
                         obj.e_r, obj.d1_l, obj.d1_r] = ...
                         sbp.implementations.d2_variable_4(m,obj.h);
                     obj.borrowing.M.d1 = 0.2505765857;
+
+                    obj.borrowing.R.delta_D = 0.577587500088313;
+                    % Borrowing e^T*D1 - d1 from R
                 case 2
                     [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
                         obj.e_r, obj.d1_l, obj.d1_r] = ...
                         sbp.implementations.d2_variable_2(m,obj.h);
                     obj.borrowing.M.d1 = 0.3636363636; 
                     % Borrowing const taken from Virta 2014
+
+                    obj.borrowing.R.delta_D = 1.000000538455350;
+                    % Borrowing e^T*D1 - d1 from R
                     
                 otherwise
                     error('Invalid operator order %d.',order);
             end
-
+            obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h,
             obj.m = m;
             obj.M = [];
         end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D2VariablePeriodic.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,71 @@
+classdef D2VariablePeriodic < sbp.OpSet
+    properties
+        D1 % SBP operator approximating first derivative
+        H % Norm matrix
+        HI % H^-1
+        Q % Skew-symmetric matrix
+        e_l % Left boundary operator
+        e_r % Right boundary operator
+        D2 % SBP operator for second derivative
+        M % Norm matrix, second derivative
+        d1_l % Left boundary first derivative
+        d1_r % Right boundary first derivative
+        m % Number of grid points.
+        h % Step size
+        x % grid
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D2VariablePeriodic(m,lim,order)
+
+            x_l = lim{1};
+            x_r = lim{2};
+            L = x_r-x_l;
+            obj.h = L/m;
+            x = linspace(x_l,x_r,m+1)';
+            obj.x = x(1:end-1);
+
+            switch order
+
+                case 6
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                        obj.e_r, obj.d1_l, obj.d1_r] = ...
+                        sbp.implementations.d2_variable_periodic_6(m,obj.h);
+                    obj.borrowing.M.d1 = 0.1878;
+                    obj.borrowing.R.delta_D = 0.3696;
+                    % Borrowing e^T*D1 - d1 from R
+
+                case 4
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                        obj.e_r, obj.d1_l, obj.d1_r] = ...
+                        sbp.implementations.d2_variable_periodic_4(m,obj.h);
+                    obj.borrowing.M.d1 = 0.2505765857;
+
+                    obj.borrowing.R.delta_D = 0.577587500088313;
+                    % Borrowing e^T*D1 - d1 from R
+                case 2
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                        obj.e_r, obj.d1_l, obj.d1_r] = ...
+                        sbp.implementations.d2_variable_periodic_2(m,obj.h);
+                    obj.borrowing.M.d1 = 0.3636363636; 
+                    % Borrowing const taken from Virta 2014
+
+                    obj.borrowing.R.delta_D = 1.000000538455350;
+                    % Borrowing e^T*D1 - d1 from R
+                    
+                otherwise
+                    error('Invalid operator order %d.',order);
+            end
+            obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h,
+
+            obj.m = m;
+            obj.M = [];
+        end
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+
+
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/InterpAWW.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,67 @@
+classdef InterpAWW < sbp.InterpOps
+    properties
+
+        % Interpolation operators
+        IC2F
+        IF2C
+
+        % Orders used on coarse and fine sides
+        order_C
+        order_F
+
+        % Grid points, refinement ratio.
+        ratio
+        m_C
+        m_F
+
+        % Boundary accuracy of IC2F and IF2C.
+        acc_C2F
+        acc_F2C
+    end
+
+    methods
+        % accOp : String, 'C2F' or 'F2C'. Specifies which of the operators
+        % should have higher accuracy.
+        function obj = InterpAWW(m_C,m_F,order_C,order_F,accOp)
+
+            ratio = (m_F-1)/(m_C-1);
+            h_C = 1;
+
+            assert(order_C == order_F,...
+                    'Error: Different orders of accuracy not available');
+
+            switch ratio
+            case 2
+                switch order_C
+                case 2
+                    [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_2to2_ratio2to1(m_C, h_C, accOp);
+                case 4
+                    [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_4to4_ratio2to1(m_C, h_C, accOp);
+                case 6
+                    [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_6to6_ratio2to1(m_C, h_C, accOp);
+                case 8
+                    [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_8to8_ratio2to1(m_C, h_C, accOp);
+                otherwise
+                    error(['Order ' num2str(order_C) ' not available.']);
+                end
+            otherwise
+                error(['Grid ratio ' num2str(ratio) ' not available']);
+            end
+
+            obj.IC2F = IC2F;
+            obj.IF2C = IF2C;
+            obj.order_C = order_C;
+            obj.order_F = order_F;
+            obj.ratio = ratio;
+            obj.m_C = m_C;
+            obj.m_F = m_F;
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_orders' num2str(obj.order_F) 'to'...
+            num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1'];
+        end
+
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/InterpMC.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,61 @@
+classdef InterpMC < sbp.InterpOps
+    properties
+
+        % Interpolation operators
+        IC2F
+        IF2C
+
+        % Orders used on coarse and fine sides
+        order_C
+        order_F
+
+        % Grid points, refinement ratio.
+        ratio
+        m_C
+        m_F
+    end
+
+    methods
+        function obj = InterpMC(m_C,m_F,order_C,order_F)
+
+            ratio = (m_F-1)/(m_C-1);
+
+            assert(order_C == order_F,...
+                    'Error: Different orders of accuracy not available');
+
+            switch ratio
+            case 2
+                switch order_C
+                case 2
+                    [IC2F,IF2C] = sbp.implementations.intOpMC_orders_2to2_ratio2to1(m_C);
+                case 4
+                    [IC2F,IF2C] = sbp.implementations.intOpMC_orders_4to4_ratio2to1(m_C);
+                case 6
+                    [IC2F,IF2C] = sbp.implementations.intOpMC_orders_6to6_ratio2to1(m_C);
+                case 8
+                    [IC2F,IF2C] = sbp.implementations.intOpMC_orders_8to8_ratio2to1(m_C);
+                otherwise
+                    error(['Order ' num2str(order_C) ' not available.']);
+                end
+            otherwise
+                error(['Grid ratio ' num2str(ratio) ' not available']);
+            end
+
+            obj.IC2F = IC2F;
+            obj.IF2C = IF2C;
+            obj.order_C = order_C;
+            obj.order_F = order_F;
+            obj.ratio = ratio;
+            obj.m_C = m_C;
+            obj.m_F = m_F;
+
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_orders' num2str(obj.order_F) 'to'...
+            num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1'];
+        end
+
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/InterpOps.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,14 @@
+classdef (Abstract) InterpOps
+    properties (Abstract)
+        % C and F may refer to coarse and fine, but it's not a must.
+        IC2F % Interpolation operator from "C" to "F"
+        IF2C % Interpolation operator from "F" to "C"
+
+    end
+
+    methods (Abstract)
+        % Returns a string representation of the type of operator.
+        str = string(obj)
+    end
+
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Elastic2dVariable.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,420 @@
+classdef Elastic2dVariable < scheme.Scheme
+
+% Discretizes the elastic wave equation:
+% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i 
+% opSet should be cell array of opSets, one per dimension. This
+% is useful if we have periodic BC in one direction.
+
+    properties
+        m % Number of points in each direction, possibly a vector
+        h % Grid spacing
+
+        grid
+        dim
+
+        order % Order of accuracy for the approximation
+
+        % Diagonal matrices for varible coefficients
+        LAMBDA % Variable coefficient, related to dilation
+        MU     % Shear modulus, variable coefficient
+        RHO, RHOi % Density, variable
+
+        D % Total operator
+        D1 % First derivatives
+
+        % Second derivatives
+        D2_lambda
+        D2_mu
+
+        % Traction operators used for BC
+        T_l, T_r
+        tau_l, tau_r
+
+        H, Hi % Inner products
+        phi % Borrowing constant for (d1 - e^T*D1) from R
+        gamma % Borrowing constant for d1 from M
+        H11 % First element of H
+        e_l, e_r
+        d1_l, d1_r % Normal derivatives at the boundary
+        E % E{i}^T picks out component i
+        
+        H_boundary % Boundary inner products
+
+        % Kroneckered norms and coefficients
+        RHOi_kron
+        Hi_kron
+    end
+
+    methods
+
+        function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
+            default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
+            default_arg('lambda_fun', @(x,y) 0*x+1);
+            default_arg('mu_fun', @(x,y) 0*x+1);
+            default_arg('rho_fun', @(x,y) 0*x+1);
+            dim = 2;
+
+            assert(isa(g, 'grid.Cartesian'))
+
+            lambda = grid.evalOn(g, lambda_fun);
+            mu = grid.evalOn(g, mu_fun);
+            rho = grid.evalOn(g, rho_fun);
+            m = g.size();
+            m_tot = g.N();
+
+            h = g.scaling();
+            lim = g.lim;
+
+            % 1D operators
+            ops = cell(dim,1);
+            for i = 1:dim
+                ops{i} = opSet{i}(m(i), lim{i}, order);
+            end
+
+            % Borrowing constants
+            for i = 1:dim
+                beta = ops{i}.borrowing.R.delta_D;
+                obj.H11{i} = ops{i}.borrowing.H11;
+                obj.phi{i} = beta/obj.H11{i};
+                obj.gamma{i} = ops{i}.borrowing.M.d1;
+            end
+
+            I = cell(dim,1);
+            D1 = cell(dim,1);
+            D2 = cell(dim,1);
+            H = cell(dim,1);
+            Hi = cell(dim,1);
+            e_l = cell(dim,1);
+            e_r = cell(dim,1);
+            d1_l = cell(dim,1);
+            d1_r = cell(dim,1);
+
+            for i = 1:dim
+                I{i} = speye(m(i));
+                D1{i} = ops{i}.D1;
+                D2{i} = ops{i}.D2;
+                H{i} =  ops{i}.H;
+                Hi{i} = ops{i}.HI;
+                e_l{i} = ops{i}.e_l;
+                e_r{i} = ops{i}.e_r;
+                d1_l{i} = ops{i}.d1_l;
+                d1_r{i} = ops{i}.d1_r;
+            end
+
+            %====== Assemble full operators ========
+            LAMBDA = spdiag(lambda);
+            obj.LAMBDA = LAMBDA;
+            MU = spdiag(mu);
+            obj.MU = MU;
+            RHO = spdiag(rho);
+            obj.RHO = RHO;
+            obj.RHOi = inv(RHO);
+
+            obj.D1 = cell(dim,1);
+            obj.D2_lambda = cell(dim,1);
+            obj.D2_mu = cell(dim,1);
+            obj.e_l = cell(dim,1);
+            obj.e_r = cell(dim,1);
+            obj.d1_l = cell(dim,1);
+            obj.d1_r = cell(dim,1);
+
+            % D1
+            obj.D1{1} = kron(D1{1},I{2});
+            obj.D1{2} = kron(I{1},D1{2});
+
+            % Boundary operators
+            obj.e_l{1} = kron(e_l{1},I{2});
+            obj.e_l{2} = kron(I{1},e_l{2});
+            obj.e_r{1} = kron(e_r{1},I{2});
+            obj.e_r{2} = kron(I{1},e_r{2});
+
+            obj.d1_l{1} = kron(d1_l{1},I{2});
+            obj.d1_l{2} = kron(I{1},d1_l{2});
+            obj.d1_r{1} = kron(d1_r{1},I{2});
+            obj.d1_r{2} = kron(I{1},d1_r{2});
+
+            % D2
+            for i = 1:dim
+                obj.D2_lambda{i} = sparse(m_tot);
+                obj.D2_mu{i} = sparse(m_tot);
+            end
+            ind = grid.funcToMatrix(g, 1:m_tot);
+
+            for i = 1:m(2)
+                D_lambda = D2{1}(lambda(ind(:,i)));
+                D_mu = D2{1}(mu(ind(:,i)));
+
+                p = ind(:,i);
+                obj.D2_lambda{1}(p,p) = D_lambda;
+                obj.D2_mu{1}(p,p) = D_mu;
+            end
+
+            for i = 1:m(1)
+                D_lambda = D2{2}(lambda(ind(i,:)));
+                D_mu = D2{2}(mu(ind(i,:)));
+
+                p = ind(i,:);
+                obj.D2_lambda{2}(p,p) = D_lambda;
+                obj.D2_mu{2}(p,p) = D_mu;
+            end
+
+            % Quadratures
+            obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
+            obj.H_boundary = cell(dim,1);
+            obj.H_boundary{1} = H{2};
+            obj.H_boundary{2} = H{1};
+
+            % E{i}^T picks out component i.
+            E = cell(dim,1);
+            I = speye(m_tot,m_tot);
+            for i = 1:dim
+                e = sparse(dim,1);
+                e(i) = 1;
+                E{i} = kron(I,e);
+            end
+            obj.E = E;
+
+            % Differentiation matrix D (without SAT)
+            D2_lambda = obj.D2_lambda;
+            D2_mu = obj.D2_mu;
+            D1 = obj.D1;
+            D = sparse(dim*m_tot,dim*m_tot);
+            d = @kroneckerDelta;    % Kronecker delta
+            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+            for i = 1:dim
+                for j = 1:dim
+                    D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +...
+                                            db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
+                                          );
+                    D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +...
+                                            db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
+                                            D2_mu{j}*E{i}' ...
+                                          );
+                end
+            end
+            obj.D = D;
+            %=========================================%
+
+            % Numerical traction operators for BC.
+            % Because d1 =/= e0^T*D1, the numerical tractions are different
+            % at every boundary.
+            T_l = cell(dim,1);
+            T_r = cell(dim,1);
+            tau_l = cell(dim,1);
+            tau_r = cell(dim,1);
+            % tau^{j}_i = sum_k T^{j}_{ik} u_k
+
+            d1_l = obj.d1_l;
+            d1_r = obj.d1_r;
+            e_l = obj.e_l;
+            e_r = obj.e_r;
+            D1 = obj.D1;
+
+            % Loop over boundaries
+            for j = 1:dim
+                T_l{j} = cell(dim,dim);
+                T_r{j} = cell(dim,dim);
+                tau_l{j} = cell(dim,1);
+                tau_r{j} = cell(dim,1);
+
+                % Loop over components
+                for i = 1:dim
+                    tau_l{j}{i} = sparse(m_tot,dim*m_tot);
+                    tau_r{j}{i} = sparse(m_tot,dim*m_tot);
+                    for k = 1:dim
+                        T_l{j}{i,k} = ... 
+                        -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
+                        -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... 
+                        -d(i,k)*MU*e_l{j}*d1_l{j}';
+
+                        T_r{j}{i,k} = ... 
+                        d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
+                        +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... 
+                        +d(i,k)*MU*e_r{j}*d1_r{j}';
+
+                        tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
+                        tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
+                    end
+
+                end
+            end
+            obj.T_l = T_l;
+            obj.T_r = T_r;
+            obj.tau_l = tau_l;
+            obj.tau_r = tau_r;
+
+            % Kroneckered norms and coefficients
+            I_dim = speye(dim);
+            obj.RHOi_kron = kron(obj.RHOi, I_dim);
+            obj.Hi_kron = kron(obj.Hi, I_dim);
+
+            % Misc.
+            obj.m = m;
+            obj.h = h;
+            obj.order = order;
+            obj.grid = g;
+            obj.dim = dim;
+
+        end
+
+
+        % Closure functions return the operators applied to the own domain 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 doamin.
+        %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
+        %       type                is a cell array of strings specifying the type of boundary condition for each component.
+        %       data                is a function returning the data that should be applied at the boundary.
+        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
+        %       neighbour_boundary  is a string specifying which boundary to interface to.
+        function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
+            default_arg('type',{'free','free'});
+            default_arg('parameter', []);
+
+            % j is the coordinate direction of the boundary
+            % nj: outward unit normal component. 
+            % nj = -1 for west, south, bottom boundaries
+            % nj = 1  for east, north, top boundaries
+            [j, nj] = obj.get_boundary_number(boundary);
+            switch nj
+            case 1
+                e = obj.e_r;
+                d = obj.d1_r;
+                tau = obj.tau_r{j};
+                T = obj.T_r{j};
+            case -1
+                e = obj.e_l;
+                d = obj.d1_l;
+                tau = obj.tau_l{j};
+                T = obj.T_l{j};
+            end
+
+            E = obj.E;
+            Hi = obj.Hi;
+            H_gamma = obj.H_boundary{j};
+            LAMBDA = obj.LAMBDA;
+            MU = obj.MU;
+            RHOi = obj.RHOi;
+
+            dim = obj.dim;
+            m_tot = obj.grid.N();
+
+            RHOi_kron = obj.RHOi_kron;
+            Hi_kron = obj.Hi_kron;
+
+            % Preallocate
+            closure = sparse(dim*m_tot, dim*m_tot);
+            penalty = cell(dim,1);
+            for k = 1:dim
+                penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
+            end
+
+            % Loop over components that we (potentially) have different BC on
+            for k = 1:dim
+                switch type{k}
+
+                % Dirichlet boundary condition
+                case {'D','d','dirichlet','Dirichlet'}
+
+                    tuning = 1.2;
+                    phi = obj.phi{j};
+                    h = obj.h(j);
+                    h11 = obj.H11{j}*h;
+                    gamma = obj.gamma{j};
+
+                    a_lambda = dim/h11 + 1/(h11*phi);
+                    a_mu_i = 2/(gamma*h);
+                    a_mu_ij = 2/h11 + 1/(h11*phi);
+
+                    d = @kroneckerDelta;  % Kronecker delta
+                    db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+                    alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
+                                          + d(i,j)* a_mu_i*MU ...
+                                          + db(i,j)*a_mu_ij*MU ); 
+
+                    % Loop over components that Dirichlet penalties end up on
+                    for i = 1:dim
+                        C = T{k,i};
+                        A = -d(i,k)*alpha(i,j);
+                        B = A + C;
+                        closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); 
+                        penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
+                    end 
+
+                % Free boundary condition
+                case {'F','f','Free','free','traction','Traction','t','T'}
+                        closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); 
+                        penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
+
+                % Unknown boundary condition
+                otherwise
+                    error('No such boundary condition: type = %s',type);
+                end
+            end
+        end
+
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+            tuning = 1.2;
+            % tuning = 20.2;
+            error('Interface not implemented');
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [j, nj] = get_boundary_number(obj, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch boundary
+                case {'w','W','west','West','s','S','south','South'}
+                    nj = -1;
+                case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                    nj = 1;
+            end
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [return_op] = get_boundary_operator(obj, op, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch op
+                case 'e'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.e_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.e_r{j};
+                    end
+                case 'd'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.d1_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.d1_r{j};
+                    end
+                otherwise
+                    error(['No such operator: operatr = ' op]);
+            end
+
+        end
+
+        function N = size(obj)
+            N = prod(obj.m);
+        end
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Heat2dVariable.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,262 @@
+classdef Heat2dVariable < scheme.Scheme
+
+% Discretizes the Laplacian with variable coefficent,
+% In the Heat equation way (i.e., the discretization matrix is not necessarily 
+% symmetric)
+% u_t = div * (kappa * grad u ) 
+% opSet should be cell array of opSets, one per dimension. This
+% is useful if we have periodic BC in one direction.
+
+    properties
+        m % Number of points in each direction, possibly a vector
+        h % Grid spacing
+
+        grid
+        dim
+
+        order % Order of accuracy for the approximation
+
+        % Diagonal matrix for variable coefficients
+        KAPPA % Variable coefficient
+
+        D % Total operator
+        D1 % First derivatives
+
+        % Second derivatives
+        D2_kappa
+
+        H, Hi % Inner products
+        e_l, e_r
+        d1_l, d1_r % Normal derivatives at the boundary
+        
+        H_boundary % Boundary inner products
+
+    end
+
+    methods
+
+        function obj = Heat2dVariable(g ,order, kappa_fun, opSet)
+            default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
+            default_arg('kappa_fun', @(x,y) 0*x+1);
+            dim = 2;
+
+            assert(isa(g, 'grid.Cartesian'))
+
+            kappa = grid.evalOn(g, kappa_fun);
+            m = g.size();
+            m_tot = g.N();
+
+            h = g.scaling();
+            lim = g.lim;
+
+            % 1D operators
+            ops = cell(dim,1);
+            for i = 1:dim
+                ops{i} = opSet{i}(m(i), lim{i}, order);
+            end
+
+            I = cell(dim,1);
+            D1 = cell(dim,1);
+            D2 = cell(dim,1);
+            H = cell(dim,1);
+            Hi = cell(dim,1);
+            e_l = cell(dim,1);
+            e_r = cell(dim,1);
+            d1_l = cell(dim,1);
+            d1_r = cell(dim,1);
+
+            for i = 1:dim
+                I{i} = speye(m(i));
+                D1{i} = ops{i}.D1;
+                D2{i} = ops{i}.D2;
+                H{i} =  ops{i}.H;
+                Hi{i} = ops{i}.HI;
+                e_l{i} = ops{i}.e_l;
+                e_r{i} = ops{i}.e_r;
+                d1_l{i} = ops{i}.d1_l;
+                d1_r{i} = ops{i}.d1_r;
+            end
+
+            %====== Assemble full operators ========
+            KAPPA = spdiag(kappa);
+            obj.KAPPA = KAPPA;
+
+            obj.D1 = cell(dim,1);
+            obj.D2_kappa = cell(dim,1);
+            obj.e_l = cell(dim,1);
+            obj.e_r = cell(dim,1);
+            obj.d1_l = cell(dim,1);
+            obj.d1_r = cell(dim,1);
+
+            % D1
+            obj.D1{1} = kron(D1{1},I{2});
+            obj.D1{2} = kron(I{1},D1{2});
+
+            % Boundary operators
+            obj.e_l{1} = kron(e_l{1},I{2});
+            obj.e_l{2} = kron(I{1},e_l{2});
+            obj.e_r{1} = kron(e_r{1},I{2});
+            obj.e_r{2} = kron(I{1},e_r{2});
+
+            obj.d1_l{1} = kron(d1_l{1},I{2});
+            obj.d1_l{2} = kron(I{1},d1_l{2});
+            obj.d1_r{1} = kron(d1_r{1},I{2});
+            obj.d1_r{2} = kron(I{1},d1_r{2});
+
+            % D2
+            for i = 1:dim
+                obj.D2_kappa{i} = sparse(m_tot);
+            end
+            ind = grid.funcToMatrix(g, 1:m_tot);
+
+            for i = 1:m(2)
+                D_kappa = D2{1}(kappa(ind(:,i)));
+                p = ind(:,i);
+                obj.D2_kappa{1}(p,p) = D_kappa;
+            end
+
+            for i = 1:m(1)
+                D_kappa = D2{2}(kappa(ind(i,:)));
+                p = ind(i,:);
+                obj.D2_kappa{2}(p,p) = D_kappa;
+            end
+
+            % Quadratures
+            obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
+            obj.H_boundary = cell(dim,1);
+            obj.H_boundary{1} = H{2};
+            obj.H_boundary{2} = H{1};
+
+            % Differentiation matrix D (without SAT)
+            D2_kappa = obj.D2_kappa;
+            D1 = obj.D1;
+            D = sparse(m_tot,m_tot);
+            for i = 1:dim
+                D = D + D2_kappa{i};
+            end
+            obj.D = D;
+            %=========================================%
+
+            % Misc.
+            obj.m = m;
+            obj.h = h;
+            obj.order = order;
+            obj.grid = g;
+            obj.dim = dim;
+
+        end
+
+
+        % Closure functions return the operators applied to the own domain 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 doamin.
+        %       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.
+        %       data                is a function returning the data that should be applied at the boundary.
+        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
+        %       neighbour_boundary  is a string specifying which boundary to interface to.
+        function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
+            default_arg('type','Neumann');
+            default_arg('parameter', []);
+
+            % j is the coordinate direction of the boundary
+            % nj: outward unit normal component. 
+            % nj = -1 for west, south, bottom boundaries
+            % nj = 1  for east, north, top boundaries
+            [j, nj] = obj.get_boundary_number(boundary);
+            switch nj
+            case 1
+                e = obj.e_r;
+                d = obj.d1_r;
+            case -1
+                e = obj.e_l;
+                d = obj.d1_l;
+            end
+
+            Hi = obj.Hi;
+            H_gamma = obj.H_boundary{j};
+            KAPPA = obj.KAPPA;
+            kappa_gamma = e{j}'*KAPPA*e{j}; 
+
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
+                    closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 
+                    penalty =  nj*Hi*d{j}*kappa_gamma*H_gamma;
+
+            % Free boundary condition
+            case {'N','n','neumann','Neumann'}
+                    closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 
+                    penalty =  nj*Hi*e{j}*kappa_gamma*H_gamma; 
+
+            % Unknown boundary condition
+            otherwise
+                error('No such boundary condition: type = %s',type);
+            end
+        end
+
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+            error('Interface not implemented');
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [j, nj] = get_boundary_number(obj, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch boundary
+                case {'w','W','west','West','s','S','south','South'}
+                    nj = -1;
+                case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                    nj = 1;
+            end
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [return_op] = get_boundary_operator(obj, op, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch op
+                case 'e'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.e_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.e_r{j};
+                    end
+                case 'd'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.d1_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.d1_r{j};
+                    end
+                otherwise
+                    error(['No such operator: operatr = ' op]);
+            end
+
+        end
+
+        function N = size(obj)
+            N = prod(obj.m);
+        end
+    end
+end
--- a/+scheme/LaplaceCurvilinear.m	Fri Jun 15 18:10:26 2018 -0700
+++ b/+scheme/LaplaceCurvilinear.m	Sat Jun 16 14:30:45 2018 -0700
@@ -38,22 +38,29 @@
         du_n, dv_n
         gamm_u, gamm_v
         lambda
+
+        interpolation_type
     end
 
     methods
         % Implements  a*div(b*grad(u)) as a SBP scheme
         % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
 
-        function obj = LaplaceCurvilinear(g ,order, a, b, opSet)
+        function obj = LaplaceCurvilinear(g ,order, a, b, opSet, interpolation_type)
             default_arg('opSet',@sbp.D2Variable);
             default_arg('a', 1);
             default_arg('b', 1);
+            default_arg('interpolation_type','AWW');
 
             if b ~=1
                 error('Not implemented yet')
             end
 
-            assert(isa(g, 'grid.Curvilinear'))
+            % assert(isa(g, 'grid.Curvilinear'))
+            if isa(a, 'function_handle')
+                a = grid.evalOn(g, a);
+                a = spdiag(a);
+            end
 
             m = g.size();
             m_u = m(1);
@@ -209,6 +216,7 @@
             obj.h = [h_u h_v];
             obj.order = order;
             obj.grid = g;
+            obj.interpolation_type = interpolation_type;
 
             obj.a = a;
             obj.b = b;
@@ -269,13 +277,70 @@
         end
 
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+
+            [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
+            [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            Hi = obj.Hi;
+            a = obj.a;
+
+            m_u = length(H_b_u);
+            m_v = length(H_b_v);
+
+            grid_ratio = m_u/m_v;
+             if grid_ratio ~= 1
+
+                [ms, index] = sort([m_u, m_v]);
+                orders = [obj.order, neighbour_scheme.order];
+                orders = orders(index);
+
+                switch obj.interpolation_type
+                case 'MC'
+                    interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
+                    if grid_ratio < 1
+                        I_v2u_good = interpOpSet.IF2C;
+                        I_v2u_bad = interpOpSet.IF2C;
+                        I_u2v_good = interpOpSet.IC2F;
+                        I_u2v_bad = interpOpSet.IC2F;
+                    elseif grid_ratio > 1
+                        I_v2u_good = interpOpSet.IC2F;
+                        I_v2u_bad = interpOpSet.IC2F;
+                        I_u2v_good = interpOpSet.IF2C;
+                        I_u2v_bad = interpOpSet.IF2C;
+                    end
+                case 'AWW'
+                    %String 'C2F' indicates that ICF2 is more accurate.
+                    interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
+                    interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 
+                    if grid_ratio < 1 
+                        % Local is coarser than neighbour
+                        I_v2u_good = interpOpSetF2C.IF2C;
+                        I_v2u_bad = interpOpSetC2F.IF2C;
+                        I_u2v_good = interpOpSetC2F.IC2F;
+                        I_u2v_bad = interpOpSetF2C.IC2F;
+                    elseif grid_ratio > 1
+                        % Local is finer than neighbour 
+                        I_v2u_good = interpOpSetC2F.IC2F;
+                        I_v2u_bad = interpOpSetF2C.IC2F;
+                        I_u2v_good = interpOpSetF2C.IF2C;
+                        I_u2v_bad = interpOpSetC2F.IF2C;
+                    end
+                otherwise
+                    error(['Interpolation type ' obj.interpolation_type ...
+                         ' is not available.' ]);
+                end
+
+             else 
+                % No interpolation required
+                I_v2u_good = speye(m_u,m_u);
+                I_v2u_bad = speye(m_u,m_u);
+                I_u2v_good = speye(m_u,m_u);
+                I_u2v_bad = speye(m_u,m_u);
+            end
+
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
             tuning = 1.2;
             % tuning = 20.2;
-            [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
-            [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-
             u = obj;
             v = neighbour_scheme;
 
@@ -284,18 +349,24 @@
             b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
             b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
 
-            tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
-            tau1 = tuning * spdiag(tau1);
-            tau2 = 1/2;
+            tau_u = -1./(4*b1_u) -1./(4*b2_u);
+            tau_v = -1./(4*b1_v) -1./(4*b2_v);
+
+            tau_u = tuning * spdiag(tau_u);
+            tau_v = tuning * spdiag(tau_v);
+            beta_u = tau_v;
 
-            sig1 = -1/2;
-            sig2 = 0;
+            closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
+                      a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ...
+                      a*1/2*Hi*d_u*H_b_u*e_u' + ...
+                      -a*1/2*Hi*e_u*H_b_u*d_u';
 
-            tau = (e_u*tau1 + tau2*d_u)*H_b_u;
-            sig = (sig1*e_u + sig2*d_u)*H_b_u;
+            penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ...
+                      -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ...
+                      -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ...
+                      -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v';
+            
 
-            closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
-            penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
         end
 
         % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Schrodinger2d.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,374 @@
+classdef Schrodinger2d < scheme.Scheme
+
+% Discretizes the Laplacian with constant coefficent,
+% in the Schrödinger equation way (i.e., the discretization matrix is not necessarily 
+% definite)
+% u_t = a*i*Laplace u 
+% opSet should be cell array of opSets, one per dimension. This
+% is useful if we have periodic BC in one direction.
+
+    properties
+        m % Number of points in each direction, possibly a vector
+        h % Grid spacing
+
+        grid
+        dim
+
+        order % Order of accuracy for the approximation
+
+        % Diagonal matrix for variable coefficients
+        a % Constant coefficient
+
+        D % Total operator
+        D1 % First derivatives
+
+        % Second derivatives
+        D2
+
+        H, Hi % Inner products
+        e_l, e_r
+        d1_l, d1_r % Normal derivatives at the boundary
+        e_w, e_e, e_s, e_n
+        d_w, d_e, d_s, d_n
+        
+        H_boundary % Boundary inner products
+
+        interpolation_type % MC or AWW
+
+    end
+
+    methods
+
+        function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type)
+            default_arg('interpolation_type','AWW');
+            default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
+            default_arg('a',1);
+            dim = 2;
+
+            assert(isa(g, 'grid.Cartesian'))
+            if isa(a, 'function_handle')
+                a = grid.evalOn(g, a);
+                a = spdiag(a);
+            end
+
+            m = g.size();
+            m_tot = g.N();
+
+            h = g.scaling();
+            xlim = {g.x{1}(1), g.x{1}(end)};
+            ylim = {g.x{2}(1), g.x{2}(end)};
+            lim = {xlim, ylim};
+
+            % 1D operators
+            ops = cell(dim,1);
+            for i = 1:dim
+                ops{i} = opSet{i}(m(i), lim{i}, order);
+            end
+
+            I = cell(dim,1);
+            D1 = cell(dim,1);
+            D2 = cell(dim,1);
+            H = cell(dim,1);
+            Hi = cell(dim,1);
+            e_l = cell(dim,1);
+            e_r = cell(dim,1);
+            d1_l = cell(dim,1);
+            d1_r = cell(dim,1);
+
+            for i = 1:dim
+                I{i} = speye(m(i));
+                D1{i} = ops{i}.D1;
+                D2{i} = ops{i}.D2;
+                H{i} =  ops{i}.H;
+                Hi{i} = ops{i}.HI;
+                e_l{i} = ops{i}.e_l;
+                e_r{i} = ops{i}.e_r;
+                d1_l{i} = ops{i}.d1_l;
+                d1_r{i} = ops{i}.d1_r;
+            end
+
+            % Constant coeff D2
+            for i = 1:dim
+                D2{i} = D2{i}(ones(m(i),1));
+            end
+
+            %====== Assemble full operators ========
+            obj.D1 = cell(dim,1);
+            obj.D2 = cell(dim,1);
+            obj.e_l = cell(dim,1);
+            obj.e_r = cell(dim,1);
+            obj.d1_l = cell(dim,1);
+            obj.d1_r = cell(dim,1);
+
+            % D1
+            obj.D1{1} = kron(D1{1},I{2});
+            obj.D1{2} = kron(I{1},D1{2});
+
+            % Boundary operators
+            obj.e_l{1} = kron(e_l{1},I{2});
+            obj.e_l{2} = kron(I{1},e_l{2});
+            obj.e_r{1} = kron(e_r{1},I{2});
+            obj.e_r{2} = kron(I{1},e_r{2});
+
+            obj.d1_l{1} = kron(d1_l{1},I{2});
+            obj.d1_l{2} = kron(I{1},d1_l{2});
+            obj.d1_r{1} = kron(d1_r{1},I{2});
+            obj.d1_r{2} = kron(I{1},d1_r{2});
+
+            % D2
+            obj.D2{1} = kron(D2{1},I{2});
+            obj.D2{2} = kron(I{1},D2{2});
+
+            % Quadratures
+            obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
+            obj.H_boundary = cell(dim,1);
+            obj.H_boundary{1} = H{2};
+            obj.H_boundary{2} = H{1};
+
+            % Differentiation matrix D (without SAT)
+            D2 = obj.D2;
+            D = sparse(m_tot,m_tot);
+            for j = 1:dim
+                D = D + a*1i*D2{j};
+            end
+            obj.D = D;
+            %=========================================%
+
+            % Misc.
+            obj.m = m;
+            obj.h = h;
+            obj.order = order;
+            obj.grid = g;
+            obj.dim = dim;
+            obj.a = a;
+            obj.e_w = obj.e_l{1};
+            obj.e_e = obj.e_r{1};
+            obj.e_s = obj.e_l{2};
+            obj.e_n = obj.e_r{2};
+            obj.d_w = obj.d1_l{1};
+            obj.d_e = obj.d1_r{1};
+            obj.d_s = obj.d1_l{2};
+            obj.d_n = obj.d1_r{2};
+            obj.interpolation_type = interpolation_type;
+
+        end
+
+
+        % Closure functions return the operators applied to the own domain 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 doamin.
+        %       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.
+        %       data                is a function returning the data that should be applied at the boundary.
+        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
+        %       neighbour_boundary  is a string specifying which boundary to interface to.
+        function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
+            default_arg('type','Neumann');
+            default_arg('parameter', []);
+
+            % j is the coordinate direction of the boundary
+            % nj: outward unit normal component. 
+            % nj = -1 for west, south, bottom boundaries
+            % nj = 1  for east, north, top boundaries
+            [j, nj] = obj.get_boundary_number(boundary);
+            switch nj
+            case 1
+                e = obj.e_r;
+                d = obj.d1_r;
+            case -1
+                e = obj.e_l;
+                d = obj.d1_l;
+            end
+
+            Hi = obj.Hi;
+            H_gamma = obj.H_boundary{j};
+            a = e{j}'*obj.a*e{j};
+
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
+                    closure =  nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); 
+                    penalty = -nj*Hi*d{j}*a*1i*H_gamma;
+
+            % Free boundary condition
+            case {'N','n','neumann','Neumann'}
+                    closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); 
+                    penalty =  nj*Hi*e{j}*a*1i*H_gamma; 
+
+            % Unknown boundary condition
+            otherwise
+                error('No such boundary condition: type = %s',type);
+            end
+        end
+
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+            % Get neighbour boundary operator
+
+            [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary);
+            [coord, n] = get_boundary_number(obj, boundary);
+            switch n_nei
+            case 1
+                % North or east boundary
+                e_neighbour = neighbour_scheme.e_r;
+                d_neighbour = neighbour_scheme.d1_r;
+            case -1
+                % South or west boundary
+                e_neighbour = neighbour_scheme.e_l;
+                d_neighbour = neighbour_scheme.d1_l;
+            end
+
+            e_neighbour = e_neighbour{coord_nei};
+            d_neighbour = d_neighbour{coord_nei};
+            H_gamma = obj.H_boundary{coord};
+            Hi = obj.Hi;
+            a = obj.a;
+
+            switch coord_nei
+            case 1
+                m_neighbour = neighbour_scheme.m(2);
+            case 2
+                m_neighbour = neighbour_scheme.m(1);
+            end
+
+            switch coord
+            case 1
+                m = obj.m(2);
+            case 2
+                m = obj.m(1);
+            end
+
+           switch n
+            case 1
+                % North or east boundary
+                e = obj.e_r;
+                d = obj.d1_r;
+            case -1
+                % South or west boundary
+                e = obj.e_l;
+                d = obj.d1_l;
+            end
+            e = e{coord};
+            d = d{coord}; 
+
+            Hi = obj.Hi;
+            sigma = -n*1i*a/2;
+            tau = -n*(1i*a)'/2;
+
+            grid_ratio = m/m_neighbour;
+             if grid_ratio ~= 1
+
+                [ms, index] = sort([m, m_neighbour]);
+                orders = [obj.order, neighbour_scheme.order];
+                orders = orders(index);
+
+                switch obj.interpolation_type
+                case 'MC'
+                    interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
+                    if grid_ratio < 1
+                        I_neighbour2local_e = interpOpSet.IF2C;
+                        I_neighbour2local_d = interpOpSet.IF2C;
+                        I_local2neighbour_e = interpOpSet.IC2F;
+                        I_local2neighbour_d = interpOpSet.IC2F;
+                    elseif grid_ratio > 1
+                        I_neighbour2local_e = interpOpSet.IC2F;
+                        I_neighbour2local_d = interpOpSet.IC2F;
+                        I_local2neighbour_e = interpOpSet.IF2C;
+                        I_local2neighbour_d = interpOpSet.IF2C;
+                    end
+                case 'AWW'
+                    %String 'C2F' indicates that ICF2 is more accurate.
+                    interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
+                    interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 
+                    if grid_ratio < 1 
+                        % Local is coarser than neighbour
+                        I_neighbour2local_e = interpOpSetF2C.IF2C;
+                        I_neighbour2local_d = interpOpSetC2F.IF2C;
+                        I_local2neighbour_e = interpOpSetC2F.IC2F;
+                        I_local2neighbour_d = interpOpSetF2C.IC2F;
+                    elseif grid_ratio > 1
+                        % Local is finer than neighbour 
+                        I_neighbour2local_e = interpOpSetC2F.IC2F;
+                        I_neighbour2local_d = interpOpSetF2C.IC2F;
+                        I_local2neighbour_e = interpOpSetF2C.IF2C;
+                        I_local2neighbour_d = interpOpSetC2F.IF2C;
+                    end
+                otherwise
+                    error(['Interpolation type ' obj.interpolation_type ...
+                         ' is not available.' ]);
+                end
+
+             else 
+                % No interpolation required
+                I_neighbour2local_e = speye(m,m);
+                I_neighbour2local_d = speye(m,m);
+                I_local2neighbour_e = speye(m,m);
+                I_local2neighbour_d = speye(m,m);
+            end
+
+            closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
+            penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ...
+                      -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; 
+             
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [j, nj] = get_boundary_number(obj, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch boundary
+                case {'w','W','west','West','s','S','south','South'}
+                    nj = -1;
+                case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                    nj = 1;
+            end
+        end
+
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [return_op] = get_boundary_operator(obj, op, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch op
+                case 'e'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.e_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.e_r{j};
+                    end
+                case 'd'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.d1_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.d1_r{j};
+                    end
+                otherwise
+                    error(['No such operator: operator = ' op]);
+            end
+
+        end
+
+        function N = size(obj)
+            N = prod(obj.m);
+        end
+    end
+end
--- a/+scheme/Utux.m	Fri Jun 15 18:10:26 2018 -0700
+++ b/+scheme/Utux.m	Sat Jun 16 14:30:45 2018 -0700
@@ -2,7 +2,7 @@
    properties
         m % Number of points in each direction, possibly a vector
         h % Grid spacing
-        x % Grid
+        grid % Grid
         order % Order accuracy for the approximation
 
         H % Discrete norm
@@ -17,41 +17,40 @@
 
 
     methods 
-         function obj = Utux(m,xlim,order,operator)
-             default_arg('a',1);
-           
-           %Old operators  
-           % [x, h] = util.get_grid(xlim{:},m);
-           %ops = sbp.Ordinary(m,h,order);
-           
-           
+         function obj = Utux(g ,order, operator)
+             default_arg('operator','Standard');
+             
+             m = g.size();
+             xl = g.getBoundary('l');
+             xr = g.getBoundary('r');
+             xlim = {xl, xr};
+             
            switch operator
-               case 'NonEquidistant'
-              ops = sbp.D1Nonequidistant(m,xlim,order);
-              obj.D1 = ops.D1;
+%                case 'NonEquidistant'
+%               ops = sbp.D1Nonequidistant(m,xlim,order);
+%               obj.D1 = ops.D1;
                case 'Standard'
               ops = sbp.D2Standard(m,xlim,order);
               obj.D1 = ops.D1;
-               case 'Upwind'
-              ops = sbp.D1Upwind(m,xlim,order);
-              obj.D1 = ops.Dm;
+%                case 'Upwind'
+%               ops = sbp.D1Upwind(m,xlim,order);
+%               obj.D1 = ops.Dm;
                otherwise
                    error('Unvalid operator')
            end
-              obj.x=ops.x;
+           
+            obj.grid = g;
 
-            
             obj.H =  ops.H;
             obj.Hi = ops.HI;
         
             obj.e_l = ops.e_l;
             obj.e_r = ops.e_r;
-            obj.D=obj.D1;
+            obj.D = -obj.D1;
 
             obj.m = m;
             obj.h = ops.h;
             obj.order = order;
-            obj.x = ops.x;
 
         end
         % Closure functions return the opertors applied to the own doamin to close the boundary
@@ -61,17 +60,27 @@
         %       data                is a function returning the data that should be applied at the boundary.
         %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
         %       neighbour_boundary  is a string specifying which boundary to interface to.
-        function [closure, penalty] = boundary_condition(obj,boundary,type,data)
-            default_arg('type','neumann');
-            default_arg('data',0);
+        function [closure, penalty] = boundary_condition(obj,boundary,type)
+            default_arg('type','dirichlet');
             tau =-1*obj.e_l;  
             closure = obj.Hi*tau*obj.e_l';       
-            penalty = 0*obj.e_l;
+            penalty = -obj.Hi*tau;
                 
          end
           
          function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
-          error('An interface function does not exist yet');
+             switch boundary
+                 % Upwind coupling
+                 case {'l','left'}
+                     tau = -1*obj.e_l;
+                     closure = obj.Hi*tau*obj.e_l';       
+                     penalty = -obj.Hi*tau*neighbour_scheme.e_r';
+                 case {'r','right'}
+                     tau = 0*obj.e_r;
+                     closure = obj.Hi*tau*obj.e_r';       
+                     penalty = -obj.Hi*tau*neighbour_scheme.e_l';
+             end
+                 
          end
       
         function N = size(obj)
@@ -81,9 +90,9 @@
     end
 
     methods(Static)
-        % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
+        % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
         % and bound_v of scheme schm_v.
-        %   [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
+        %   [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
         function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
             [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
             [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Utux2D.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,286 @@
+classdef Utux2D < scheme.Scheme
+   properties
+        m % Number of points in each direction, possibly a vector
+        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.
+
+        H % Discrete norm
+        H_x, H_y % Norms in the x and y directions
+        Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
+
+        % Derivatives
+        Dx, Dy
+        
+        % Boundary operators
+        e_w, e_e, e_s, e_n
+        
+        D % Total discrete operator
+
+        % String, type of interface coupling
+        % Default: 'upwind'
+        % Other:   'centered'
+        coupling_type 
+
+        % String, type of interpolation operators
+        % Default: 'AWW' (Almquist Wang Werpers)
+        % Other:   'MC' (Mattsson Carpenter)
+        interpolation_type
+
+        
+        % Cell array, damping on upwstream and downstream sides.
+        interpolation_damping
+
+    end
+
+
+    methods 
+         function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping)
+            
+            default_arg('interpolation_damping',{0,0});
+            default_arg('interpolation_type','AWW'); 
+            default_arg('coupling_type','upwind'); 
+            default_arg('a',1/sqrt(2)*[1, 1]); 
+            default_arg('opSet',@sbp.D2Standard);
+
+            assert(isa(g, 'grid.Cartesian'))
+            if iscell(a)
+                a1 = grid.evalOn(g, a{1});
+                a2 = grid.evalOn(g, a{2});
+                a = {spdiag(a1), spdiag(a2)};
+            else
+                a = {a(1), a(2)};
+            end
+             
+            m = g.size();
+            m_x = m(1);
+            m_y = m(2);
+            m_tot = g.N();
+
+            xlim = {g.x{1}(1), g.x{1}(end)};
+            ylim = {g.x{2}(1), g.x{2}(end)};
+            obj.grid = g;
+
+            % Operator sets
+            ops_x = opSet(m_x, xlim, order);
+            ops_y = opSet(m_y, ylim, order);
+            Ix = speye(m_x);
+            Iy = speye(m_y);
+            
+            % Norms
+            Hx = ops_x.H;
+            Hy = ops_y.H;
+            Hxi = ops_x.HI;
+            Hyi = ops_y.HI;
+            
+            obj.H_x = Hx;
+            obj.H_y = Hy;
+            obj.H = kron(Hx,Hy);
+            obj.Hi = kron(Hxi,Hyi);
+            obj.Hx = kron(Hx,Iy);
+            obj.Hy = kron(Ix,Hy);
+            obj.Hxi = kron(Hxi,Iy);
+            obj.Hyi = kron(Ix,Hyi);
+            
+            % Derivatives
+            Dx = ops_x.D1;
+            Dy = ops_y.D1;
+            obj.Dx = kron(Dx,Iy);
+            obj.Dy = kron(Ix,Dy);
+           
+            % Boundary operators
+            obj.e_w = kr(ops_x.e_l, Iy);
+            obj.e_e = kr(ops_x.e_r, Iy);
+            obj.e_s = kr(Ix, ops_y.e_l);
+            obj.e_n = kr(Ix, ops_y.e_r);
+
+            obj.m = m;
+            obj.h = [ops_x.h ops_y.h];
+            obj.order = order;
+            obj.a = a;
+            obj.coupling_type = coupling_type;
+            obj.interpolation_type = interpolation_type;
+            obj.interpolation_damping = interpolation_damping;
+            obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
+
+        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.
+        %       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.
+        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
+        %       neighbour_boundary  is a string specifying which boundary to interface to.
+        function [closure, penalty] = boundary_condition(obj,boundary,type)
+            default_arg('type','dirichlet');
+            
+            sigma = -1; % Scalar penalty parameter
+            switch boundary
+                case {'w','W','west','West'}
+                    tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
+                    closure = obj.Hi*tau*obj.e_w';
+                    
+                case {'s','S','south','South'}
+                    tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
+                    closure = obj.Hi*tau*obj.e_s';
+            end  
+            penalty = -obj.Hi*tau;
+                
+         end
+          
+         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+             
+             % Get neighbour boundary operator
+             switch neighbour_boundary
+                 case {'e','E','east','East'}
+                     e_neighbour = neighbour_scheme.e_e;
+                     m_neighbour = neighbour_scheme.m(2);
+                 case {'w','W','west','West'}
+                     e_neighbour = neighbour_scheme.e_w;
+                     m_neighbour = neighbour_scheme.m(2);
+                 case {'n','N','north','North'}
+                     e_neighbour = neighbour_scheme.e_n;
+                     m_neighbour = neighbour_scheme.m(1);
+                 case {'s','S','south','South'}
+                     e_neighbour = neighbour_scheme.e_s;
+                     m_neighbour = neighbour_scheme.m(1);
+             end
+             
+             switch obj.coupling_type
+             
+             % Upwind coupling (energy dissipation)
+             case 'upwind'
+                 sigma_ds = -1; %"Downstream" penalty
+                 sigma_us = 0; %"Upstream" penalty
+
+             % Energy-preserving coupling (no energy dissipation)
+             case 'centered'
+                 sigma_ds = -1/2; %"Downstream" penalty
+                 sigma_us = 1/2; %"Upstream" penalty
+
+             otherwise
+                error(['Interface coupling type ' coupling_type ' is not available.'])
+             end
+
+             % Check grid ratio for interpolation
+             switch boundary
+                 case {'w','W','west','West','e','E','east','East'}
+                     m = obj.m(2);       
+                 case {'s','S','south','South','n','N','north','North'}
+                     m = obj.m(1);
+             end
+             grid_ratio = m/m_neighbour;
+             if grid_ratio ~= 1
+
+                [ms, index] = sort([m, m_neighbour]);
+                orders = [obj.order, neighbour_scheme.order];
+                orders = orders(index);
+
+                switch obj.interpolation_type
+                case 'MC'
+                    interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
+                    if grid_ratio < 1
+                        I_neighbour2local_us = interpOpSet.IF2C;
+                        I_neighbour2local_ds = interpOpSet.IF2C;
+                        I_local2neighbour_us = interpOpSet.IC2F;
+                        I_local2neighbour_ds = interpOpSet.IC2F;
+                    elseif grid_ratio > 1
+                        I_neighbour2local_us = interpOpSet.IC2F;
+                        I_neighbour2local_ds = interpOpSet.IC2F;
+                        I_local2neighbour_us = interpOpSet.IF2C;
+                        I_local2neighbour_ds = interpOpSet.IF2C;
+                    end
+                case 'AWW'
+                    %String 'C2F' indicates that ICF2 is more accurate.
+                    interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
+                    interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 
+                    if grid_ratio < 1 
+                        % Local is coarser than neighbour
+                        I_neighbour2local_us = interpOpSetC2F.IF2C;
+                        I_neighbour2local_ds = interpOpSetF2C.IF2C;
+                        I_local2neighbour_us = interpOpSetC2F.IC2F;
+                        I_local2neighbour_ds = interpOpSetF2C.IC2F;
+                    elseif grid_ratio > 1
+                        % Local is finer than neighbour 
+                        I_neighbour2local_us = interpOpSetF2C.IC2F;
+                        I_neighbour2local_ds = interpOpSetC2F.IC2F;
+                        I_local2neighbour_us = interpOpSetF2C.IF2C;
+                        I_local2neighbour_ds = interpOpSetC2F.IF2C;
+                    end
+                otherwise
+                    error(['Interpolation type ' obj.interpolation_type ...
+                         ' is not available.' ]);
+                end
+
+             else 
+                % No interpolation required
+                I_neighbour2local_us = speye(m,m);
+                I_neighbour2local_ds = speye(m,m);
+            end    
+             
+             int_damp_us = obj.interpolation_damping{1};
+             int_damp_ds = obj.interpolation_damping{2};
+
+             I = speye(m,m);
+             I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
+             I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
+
+
+             switch boundary
+                 case {'w','W','west','West'}
+                     tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
+                     closure = obj.Hi*tau*obj.e_w';
+                     penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
+
+                     beta = int_damp_ds*obj.a{1}...
+                            *obj.e_w*obj.H_y;
+                     closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w';     
+                 case {'e','E','east','East'}
+                     tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
+                     closure = obj.Hi*tau*obj.e_e';
+                     penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
+
+                     beta = int_damp_us*obj.a{1}...
+                            *obj.e_e*obj.H_y;
+                     closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; 
+                 case {'s','S','south','South'}
+                     tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
+                     closure = obj.Hi*tau*obj.e_s'; 
+                     penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
+
+                     beta = int_damp_ds*obj.a{2}...
+                            *obj.e_s*obj.H_x;
+                     closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s';
+                 case {'n','N','north','North'}
+                     tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
+                     closure = obj.Hi*tau*obj.e_n';
+                     penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
+
+                     beta = int_damp_us*obj.a{2}...
+                            *obj.e_n*obj.H_x;
+                     closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; 
+             end
+             
+                 
+         end
+      
+        function N = size(obj)
+            N = obj.m;
+        end
+
+    end
+
+    methods(Static)
+        % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
+        % and bound_v of scheme schm_v.
+        %   [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
+        function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
+            [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
+            [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
+        end
+    end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/kroneckerDelta.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,6 @@
+function d = kroneckerDelta(i,j)
+
+d = 0;
+if i==j
+	d = 1;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/spdiagVariable.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,17 @@
+function A = spdiagVariable(a,i)
+    default_arg('i',0);
+
+    if isrow(a)
+        a = a';
+    end
+
+    n = length(a)+abs(i);
+
+    if i > 0
+    	a = [sparse(i,1); a];
+    elseif i < 0
+    	a = [a; sparse(abs(i),1)];
+    end
+
+    A = spdiags(a,i,n,n);
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/spdiagsVariablePeriodic.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,42 @@
+function A = spdiagsVariablePeriodic(vals,diags)
+    % Creates an m x m periodic discretization matrix.
+    % vals - m x ndiags matrix of values
+    % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on
+    % vals that are not on main diagonal are going to spill over to 
+    % off-diagonal corners.
+
+    default_arg('diags',0);
+
+    [m, ~] = size(vals); 
+
+    A = sparse(m,m);
+
+    for i = 1:length(diags)
+        
+        d = diags(i);
+        a = vals(:,i);
+
+        % Sub-diagonals
+        if d < 0
+            a_bulk = a(1+abs(d):end);
+            a_corner = a(1:1+abs(d)-1);
+            corner_diag = m-abs(d);
+            A = A + spdiagVariable(a_bulk, d); 
+            A = A + spdiagVariable(a_corner, corner_diag);
+
+        % Super-diagonals
+        elseif d > 0
+            a_bulk = a(1:end-d);
+            a_corner = a(end-d+1:end);
+            corner_diag = -m + d;
+            A = A + spdiagVariable(a_bulk, d); 
+            A = A + spdiagVariable(a_corner, corner_diag);
+
+        % Main diagonal
+        else
+             A = A + spdiagVariable(a, 0);
+        end
+
+    end
+
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stripeMatrixPeriodic.m	Sat Jun 16 14:30:45 2018 -0700
@@ -0,0 +1,8 @@
+% Creates a periodic discretization matrix of size n x n 
+%  with the values of val on the diagonals diag.
+%   A = stripeMatrix(val,diags,n)
+function A = stripeMatrixPeriodic(val,diags,n)
+
+    D = ones(n,1)*val;
+    A = spdiagsVariablePeriodic(D,diags);
+end
\ No newline at end of file