--- a/+anim/setup_time_quantity_plot.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+anim/setup_time_quantity_plot.m	Tue Dec 04 15:24:36 2018 -0800
@@ -16,7 +16,7 @@
         if ishandle(axis_handle)
             % t = [t t_now];
             for j = 1:length(yfun)
-                addpoints(plot_handles(j),t_now,yfun{j}(varargin{:}));
+                addpoints(plot_handles(j),t_now,full(yfun{j}(varargin{:})));
             end
 
             [t,~] = getpoints(plot_handles(1));

--- a/+draw/prompt_point.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+draw/prompt_point.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,22 +1,23 @@
-function [p, button] = prompt_point(s,varargin)
+function [p, button] = prompt_point(s, varargin)
     default_arg('s',[])
 
     set(gcf,'Pointer','crosshair')
 
     if ~isempty(s)
-        fprintf(s,varargin{:});
+        fprintf(s, varargin{:});
     end
 
-    a = gca;
+    fh = gcf();
+    ah = gca();
 
-    function get_point(src,event)
-        cp = a.CurrentPoint;
+    function get_point(src, event)
+        cp = ah.CurrentPoint;
         p = cp(1,1:2)';
-        a.ButtonDownFcn = [];
+        fh.WindowButtonUpFcn = [];
     end
 
-    a.ButtonDownFcn = @get_point;
-    waitfor(a,'ButtonDownFcn', [])
+    fh.WindowButtonUpFcn = @get_point;
+    waitfor(fh,'WindowButtonUpFcn', [])
 
     set(gcf,'Pointer','arrow')
 

--- a/+grid/Cartesian.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+grid/Cartesian.m	Tue Dec 04 15:24:36 2018 -0800
@@ -16,7 +16,7 @@
             obj.d = length(varargin);
 
             for i = 1:obj.d
-                assert(isvector(varargin{i}), 'Coordinate inputs must be a vectors.')
+                assert(isnumeric(varargin{i}), 'Coordinate inputs must be vectors.')
 
                 obj.x{i} = varargin{i};
                 obj.m(i) = length(varargin{i});
@@ -28,7 +28,11 @@
             end
 
             obj.h = [];
-            obj.lim = [];
+
+            obj.lim = cell(1,obj.d);
+            for i = 1:obj.d
+                obj.lim{i} = {obj.x{i}(1), obj.x{i}(end)};
+            end
         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/+grid/Nodes.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,47 @@
+classdef Nodes < grid.Grid
+    properties
+        coords
+    end
+
+    methods
+        % Creates a grid with one point for each row in coords.
+        % The dimension equals the number of columns in coords.
+        function obj = Nodes(coords)
+            obj.coords = coords;
+        end
+
+        function o = N(obj)
+            o = size(obj.coords, 1);
+        end
+
+        % d returns the spatial dimension of the grid
+        function o = D(obj)
+            o = size(obj.coords, 2);
+        end
+
+        % points returns a n x d matrix containing the coordinates for all points.
+        function X = points(obj)
+            X = obj.coords;
+        end
+
+        % Restricts the grid function gf on obj to the subgrid g.
+        function gf = restrictFunc(obj, gf, g)
+            error('Not implemented');
+        end
+
+        % Projects the grid function gf on obj to the grid g.
+        function gf = projectFunc(obj, gf, g)
+            error('Not implemented');
+        end
+
+        % Return the grid.boundaryIdentifiers of all boundaries in a cell array.
+        function bs = getBoundaryNames(obj)
+            error('Not implemented');
+        end
+
+        % Return coordinates for the given boundary
+        function b = getBoundary(obj, name)
+            error('Not implemented');
+        end
+    end
+end

--- a/+grid/evalOn.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+grid/evalOn.m	Tue Dec 04 15:24:36 2018 -0800
@@ -13,18 +13,18 @@
         return
     end
     % func should now be a function_handle
-    assert(g.D == nargin(func),'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.')
+    assert(g.D == nargin(func) || nargin(func) < 0,'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.')
 
     x = num2cell(g.points(),1);
-    k = numberOfComponents(func);
+    k = numberOfComponents(func, g.D);
 
     gf = func(x{:});
     gf = reorderComponents(gf, k);
 end
 
 % Find the number of vector components of func
-function k = numberOfComponents(func)
-    x0 = num2cell(ones(1,nargin(func)));
+function k = numberOfComponents(func, dim)
+    x0 = num2cell(ones(1, dim));
     f0 = func(x0{:});
     assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector');
     k = length(f0);

--- a/+multiblock/+domain/Circle.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+multiblock/+domain/Circle.m	Tue Dec 04 15:24:36 2018 -0800
@@ -65,10 +65,10 @@
             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.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);

--- a/+multiblock/DefCurvilinear.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+multiblock/DefCurvilinear.m	Tue Dec 04 15:24:36 2018 -0800
@@ -48,13 +48,14 @@
             g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
         end
 
-        function show(obj, label, gridLines, varargin)
+        function h = show(obj, label, gridLines, varargin)
             default_arg('label', 'name')
             default_arg('gridLines', false);
 
+            h = [];
             if isempty('label') && ~gridLines
                 for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(2,2);
+                    h = [h, obj.blockMaps{i}.show(2,2)];
                 end
                 axis equal
                 return
@@ -63,7 +64,7 @@
             if gridLines
                 ms = obj.getGridSizes(varargin{:});
                 for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
+                    h = [h, obj.blockMaps{i}.show(ms{i}(1),ms{i}(2))];
                 end
             end
 
@@ -76,7 +77,7 @@
                     for i = 1:obj.nBlocks
                         labels{i} = num2str(i);
                     end
-                otherwise
+                case 'none'
                     axis equal
                     return
             end

--- a/+multiblock/DiffOp.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+multiblock/DiffOp.m	Tue Dec 04 15:24:36 2018 -0800
@@ -10,13 +10,13 @@
     end
 
     methods
-        function obj = DiffOp(doHand, grid, order, doParam, intfTypes)
+        function obj = DiffOp(doHand, g, order, doParam, intfTypes)
             %  doHand -- may either be a function handle or a cell array of
             %            function handles for each grid. The function handle(s)
             %            should be on the form do = doHand(grid, order, ...)
             %            Additional parameters for each doHand may be provided in
             %            the doParam input.
-            %    grid -- a multiblock grid
+            %       g -- a multiblock grid
             %   order -- integer specifying the order of accuracy
             % doParam -- may either be a cell array or a cell array of cell arrays
             %            for each block. If it is a cell array with length equal
@@ -28,12 +28,12 @@
             % intfTypes (optional) -- nBlocks x nBlocks cell array of types for
             %                                 every interface.
             default_arg('doParam', [])
-            default_arg('intfTypes', cell(grid.nBlocks(), grid.nBlocks()) );
+            default_arg('intfTypes', cell(g.nBlocks(), g.nBlocks()) );
 
-            [getHand, getParam] = parseInput(doHand, grid, doParam);
+            [getHand, getParam] = parseInput(doHand, g, doParam);
 
             obj.order = order;
-            nBlocks = grid.nBlocks();
+            nBlocks = g.nBlocks();
 
             % Create the diffOps for each block
             obj.diffOps = cell(1, nBlocks);
@@ -43,7 +43,7 @@
                 if ~iscell(p)
                     p = {p};
                 end
-                obj.diffOps{i} = h(grid.grids{i}, order, p{:});
+                obj.diffOps{i} = h(g.grids{i}, order, p{:});
             end
 
 
@@ -56,7 +56,11 @@
 
 
             % Build the differentiation matrix
-            obj.blockmatrixDiv = {grid.Ns, grid.Ns};
+            Ns = zeros(nBlocks,1);
+            for i = 1:nBlocks
+                Ns(i) = length(obj.diffOps{i}.D);
+            end
+            obj.blockmatrixDiv = {Ns, Ns};
             D = blockmatrix.zero(obj.blockmatrixDiv);
             for i = 1:nBlocks
                 D{i,i} = obj.diffOps{i}.D;
@@ -64,7 +68,7 @@
 
             for i = 1:nBlocks
                 for j = 1:nBlocks
-                    intf = grid.connections{i,j};
+                    intf = g.connections{i,j};
                     if isempty(intf)
                         continue
                     end
@@ -79,14 +83,15 @@
                 end
             end
             obj.D = blockmatrix.toMatrix(D);
+            obj.grid = g;
 
 
-            function [getHand, getParam] = parseInput(doHand, grid, doParam)
-                if ~isa(grid, 'multiblock.Grid')
+            function [getHand, getParam] = parseInput(doHand, g, doParam)
+                if ~isa(g, 'multiblock.Grid')
                     error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
                 end
 
-                if iscell(doHand) && length(doHand) == grid.nBlocks()
+                if iscell(doHand) && length(doHand) == g.nBlocks()
                     getHand = @(i)doHand{i};
                 elseif isa(doHand, 'function_handle')
                     getHand = @(i)doHand;
@@ -106,7 +111,7 @@
 
                 % doParam is a non-empty cell-array
 
-                if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
+                if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
                     % doParam is a cell-array of cell-arrays
                     getParam = @(i)doParam{i};
                     return
@@ -145,6 +150,27 @@
             end
         end
 
+        function op = getBoundaryQuadrature(obj, boundary)
+            opName = 'H';
+            switch class(boundary)
+                case 'cell'
+                    localOpName = [opName '_' boundary{2}];
+                    blockId = boundary{1};
+                    op = obj.diffOps{blockId}.(localOpName);
+
+                    return
+                case 'multiblock.BoundaryGroup'
+                    N = length(boundary);
+                    H_bm = cell(N,N);
+                    for i = 1:N
+                        H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
+                    end
+                    op = blockmatrix.toMatrix(H_bm);
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+        end
+
         % Creates the closure and penalty matrix for a given boundary condition,
         %    boundary -- the name of the boundary on the form {id,name} where
         %                id is the number of a block and name is the name of a
@@ -177,33 +203,8 @@
             [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
 
             % Expand to matrix for full domain.
-            div = obj.blockmatrixDiv;
-            if ~iscell(blockClosure)
-                temp = blockmatrix.zero(div);
-                temp{I,I} = blockClosure;
-                closure = blockmatrix.toMatrix(temp);
-            else
-                for i = 1:length(blockClosure)
-                    temp = blockmatrix.zero(div);
-                    temp{I,I} = blockClosure{i};
-                    closure{i} = blockmatrix.toMatrix(temp);
-                end
-            end
-
-            if ~iscell(blockPenalty)
-                div{2} = size(blockPenalty, 2); % Penalty is a column vector
-                p = blockmatrix.zero(div);
-                p{I} = blockPenalty;
-                penalty = blockmatrix.toMatrix(p);
-            else
-                % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
-                for i = 1:length(blockPenalty)
-                    div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
-                    p = blockmatrix.zero(div);
-                    p{I} = blockPenalty{i};
-                    penalty{i} = blockmatrix.toMatrix(p);
-                end
-            end
+            closure = multiblock.local2globalClosure(blockClosure, obj.blockmatrixDiv, I);
+            penalty = multiblock.local2globalPenalty(blockPenalty, obj.blockmatrixDiv, I);
         end
 
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)

--- a/+multiblock/Grid.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+multiblock/Grid.m	Tue Dec 04 15:24:36 2018 -0800
@@ -77,7 +77,7 @@
             % Collect number of points in each block
             N = zeros(1,nBlocks);
             for i = 1:nBlocks
-                N(i) = obj.grids{i}.N();
+                N(i) = obj.grids{i}.N()*nComponents;
             end
 
             gfs = blockmatrix.fromMatrix(gf, {N,1});

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/Laplace.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,56 @@
+classdef Laplace < scheme.Scheme
+    properties
+        grid
+        order
+        mbDiffOp
+
+        D
+        H
+        J
+    end
+    methods
+        function obj = Laplace(g, order, a, b, opGen)
+            default_arg('order', 4);
+            default_arg('a', 1);
+            default_arg('b', 1);
+            default_arg('opGen', @sbp.D4Variable);
+
+            obj.grid = g;
+            obj.order = order;
+            obj.mbDiffOp = multiblock.DiffOp(@scheme.LaplaceCurvilinear, obj.grid, order, {a,b,opGen});
+
+            obj.D = obj.mbDiffOp.D;
+            obj.J = obj.jacobian();
+            obj.H = obj.mbDiffOp.H;
+        end
+
+        function s = size(obj)
+            s = size(obj.mbDiffOp);
+        end
+
+        function J = jacobian(obj)
+            N = obj.grid.nBlocks;
+            J = cell(N,N);
+
+            for i = 1:N
+                J{i,i} = obj.mbDiffOp.diffOps{i}.J;
+            end
+            J = blockmatrix.toMatrix(J);
+        end
+
+        function op = getBoundaryOperator(obj, opName, boundary)
+            op = getBoundaryOperator(obj.mbDiffOp, opName, boundary);
+        end
+
+        function op = getBoundaryQuadrature(obj, boundary)
+            op = getBoundaryQuadrature(obj.mbDiffOp, boundary);
+        end
+
+        function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
+            [closure, penalty] = boundary_condition(obj.mbDiffOp, boundary, type);
+        end
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+            error('Not implemented')
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/evalOn.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,27 @@
+% Evaluate different function handle for each block in a multiblock.Grid
+% Function handles may optionaly take a time argument
+% f -- cell array of function handles
+%       f{i} = f_i(t,x,y,...)
+% t -- optional time point. If not specified, it is assumed that the functions take only spatial arguments.
+function gf = evalOn(g, f, t)
+    assertType(g, 'multiblock.Grid');
+    assertType(f, 'cell');
+
+    default_arg('t', []);
+
+    grids = g.grids;
+    nBlocks = length(grids);
+    gf = cell(nBlocks, 1);
+
+    if isempty(t)
+        for i = 1:nBlocks
+            gf{i} = grid.evalOn(grids{i}, f{i});
+        end
+    else
+        for i = 1:nBlocks
+            gf{i} = grid.evalOn(grids{i}, @(varargin)f{i}(t,varargin{:}));
+        end
+    end
+
+    gf = blockmatrix.toMatrix(gf);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/local2globalClosure.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,10 @@
+% Takes the block-local closures and turns it into a global closure
+%   local -- The local closure
+%   div   -- block matrix division for the diffOp
+%   I     -- Index of blockmatrix block
+function closure = local2globalClosure(local, div, I)
+    closure_bm = blockmatrix.zero(div);
+    closure_bm{I,I} = local;
+
+    closure = blockmatrix.toMatrix(closure_bm);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/local2globalPenalty.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,11 @@
+% Takes the block-local penalty and turns it into a global penalty
+%   local -- The local penalty
+%   div   -- block matrix division for the diffOp
+%   I     -- Index of blockmatrix block
+function penalty = local2globalPenalty(local, div, I)
+    penaltyDiv = {div{1}, size(local,2)};
+    penalty_bm = blockmatrix.zero(penaltyDiv);
+    penalty_bm{I,1} = local;
+
+    penalty = blockmatrix.toMatrix(penalty_bm);
+end

--- a/+multiblock/multiblockgrid.m	Tue Dec 04 14:54:28 2018 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,70 +0,0 @@
-% Creates a multi block square grid with defined boundary conditions.
-%   x,y defines the grid lines. Rember to think of the indexing as a matrix. Order matters!
-%   bc is a struct defining the boundary conditions on each side of the block.
-%       bc.w = {'dn',[function or value]}
-function [block,conn,bound,ms] = multiblockgrid(x,y,mx,my,bc)
-    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
-    % y = sort(y,'descend');
-
-    % Dimensions of blocks and number of points
-    block = cell(n,m);
-    for i = 1:n
-        for j = 1:m
-            block{i,j} = {
-                {x(j),x(j+1)}, {y(i+1),y(i)};
-            };
-
-            ms{i,j} = [mx(i),my(j)];
-        end
-    end
-
-    % Interface couplings
-    conn = cell(N,N);
-    for i = 1:n
-        for j = 1:m
-            I = flat_index(n,i,j);
-            if i < n
-                J = flat_index(n,i+1,j);
-                conn{I,J} = {'s','n'};
-            end
-
-            if j < m
-                J = flat_index(n,i,j+1);
-                conn{I,J} = {'e','w'};
-            end
-        end
-    end
-
-
-    % Boundary conditions
-    bound = cell(n,m);
-    for i = 1:n
-        if isfield(bc,'w')
-            bound{i,1}.w = bc.w;
-        end
-
-        if isfield(bc,'e')
-            bound{i,n}.e = bc.e;
-        end
-    end
-
-    for j = 1:m
-        if isfield(bc,'n')
-            bound{1,j}.n = bc.n;
-        end
-
-        if isfield(bc,'s')
-            bound{m,j}.s = bc.s;
-        end
-    end
-end
-

--- a/+multiblock/stitchSchemes.m	Tue Dec 04 14:54:28 2018 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,92 +0,0 @@
-% Stitch schemes together given connection matrix and BC vector.
-%     schmHand  - function_handle to a Scheme constructor
-%     order     - order of accuracy
-%     schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays
-%     blocks    - block definitions, On whatever form the scheme expects.
-%     ms        - grid points in each direction for each block. Ex {[10,10], [10, 20]}
-%     conn      - connection matrix
-%     bound     - boundary condition vector, array of structs with fields w,e,s,n
-%                 each field with a parameter array that is sent to schm.boundary_condition
-%
-% Output parameters are cell arrays and cell matrices.
-%
-% Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound)
-function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound)
-    default_arg('schmParam',[]);
-
-    n_blocks = numel(grids);
-
-    % Creating Schemes
-    for i = 1:n_blocks
-        if isempty(schmParam);
-            schms{i} = schmHand(grids{i},order,[]);
-        elseif ~iscell(schmParam)
-            param = schmParam(i);
-            schms{i} = schmHand(grids{i},order,param);
-        else
-            param = schmParam{i};
-            if iscell(param)
-                schms{i} = schmHand(grids{i},order,param{:});
-            else
-                schms{i} = schmHand(grids{i},order,param);
-            end
-        end
-
-        % class(schmParam)
-        % class(ms)
-        % class(blocks)
-        % class(schmParam{i})
-        % class(ms)
-
-
-    end
-
-
-    % Total norm
-    H = cell(n_blocks,n_blocks);
-    for i = 1:n_blocks
-        H{i,i} = schms{i}.H;
-    end
-
-    %% Total system matrix
-
-    % Differentiation terms
-    D = cell(n_blocks,n_blocks);
-    for i = 1:n_blocks
-        D{i,i} = schms{i}.D;
-    end
-
-    % Boundary penalty terms
-    for i = 1:n_blocks
-        if ~isstruct(bound{i})
-            continue
-        end
-
-        fn = fieldnames(bound{i});
-        for j = 1:length(fn);
-            bc = bound{i}.(fn{j});
-            if isempty(bc)
-                continue
-            end
-
-            [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:});
-            D{i,i} = D{i,i}+closure;
-        end
-    end
-
-    % Interface penalty terms
-    for i = 1:n_blocks
-        for j = 1:n_blocks
-            intf = conn{i,j};
-            if isempty(intf)
-                continue
-            end
-
-            [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
-            D{i,i} = D{i,i} + uu;
-            D{i,j} = uv;
-            D{j,j} = D{j,j} + vv;
-            D{j,i} = vu;
-        end
-    end
-end
\ No newline at end of file

--- a/+noname/calculateErrors.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+noname/calculateErrors.m	Tue Dec 04 15:24:36 2018 -0800
@@ -4,27 +4,40 @@
 % m are grid size parameters.
 % N are number of timesteps to use for each gird size
 % timeOpt are options for the timeStepper
+% errorFun is a function_handle taking 2 or 3 arguments, errorFun(trueSolution, approxSolution), errorFun(trueSolution, approxSolution, discr)
 function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt)
+    %TODO: Ability to choose paralell or not
     assertType(schemeFactory, 'function_handle');
     assertNumberOfArguments(schemeFactory, 1);
     assertScalar(T);
     assert(length(m) == length(N), 'Vectors m and N must have the same length');
     assertType(errorFun, 'function_handle');
-    assertNumberOfArguments(errorFun, 2);
-    default_arg('timeOpt');
+
+    if ~ismember(nargin(errorFun), [2,3])
+        error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun));
+    end
 
-    e = [];
-    for i = 1:length(m)
+    default_arg('timeOpt', struct());
+
+
+    e = zeros(1,length(m));
+    parfor i = 1:length(m)
         done = timeTask('m = %3d ', m(i));
 
         [discr, trueSolution] = schemeFactory(m(i));
 
-        timeOpt.k = T/N(i);
-        ts = discr.getTimestepper(timeOpt);
+        timeOptTemp = timeOpt;
+        timeOptTemp.k = T/N(i);
+        ts = discr.getTimestepper(timeOptTemp);
         ts.stepTo(N(i), true);
         approxSolution = discr.getTimeSnapshot(ts);
 
-        e(i) = errorFun(trueSolution, approxSolution);
+        switch nargin(errorFun)
+            case 2
+                e(i) = errorFun(trueSolution, approxSolution);
+            case 3
+                e(i) = errorFun(trueSolution, approxSolution, discr);
+        end
 
         fprintf('e = %.4e', e(i))
         done()

--- a/+parametrization/Curve.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+parametrization/Curve.m	Tue Dec 04 15:24:36 2018 -0800
@@ -181,8 +181,8 @@
         end
 
         function D = mirror(C, a, b)
-            assert_size(a,[2,1]);
-            assert_size(b,[2,1]);
+            assertSize(a,[2,1]);
+            assertSize(b,[2,1]);
 
             g = C.g;
             gp = C.gp;
@@ -219,8 +219,8 @@
         end
 
         function D = rotate(C,a,rad)
-            assert_size(a, [2,1]);
-            assert_size(rad, [1,1]);
+            assertSize(a, [2,1]);
+            assertSize(rad, [1,1]);
             g = C.g;
             gp = C.gp;
 

--- a/+parametrization/Ti.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+parametrization/Ti.m	Tue Dec 04 15:24:36 2018 -0800
@@ -129,13 +129,13 @@
             S = obj.S;
 
             if(nu>2 || nv>2)
-                h_grid = obj.plot(nu,nv);
-                set(h_grid,'Color',[0 0.4470 0.7410]);
+                h.grid = obj.plot(nu,nv);
+                set(h.grid,'Color',[0 0.4470 0.7410]);
             end
 
-            h_bord = obj.plot(2,2);
-            set(h_bord,'Color',[0.8500 0.3250 0.0980]);
-            set(h_bord,'LineWidth',2);
+            h.border = obj.plot(2,2);
+            set(h.border,'Color',[0.8500 0.3250 0.0980]);
+            set(h.border,'LineWidth',2);
         end
 
 

--- a/+sbp/+implementations/d2_variable_periodic_2.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+sbp/+implementations/d2_variable_periodic_2.m	Tue Dec 04 15:24:36 2018 -0800
@@ -27,7 +27,7 @@
 
     scheme_width = 3;
     scheme_radius = (scheme_width-1)/2;
-    
+
     r = 1:m;
     offset = scheme_width;
     r = r + offset;
@@ -41,7 +41,7 @@
 
         vals = [Mm1,M0,Mp1];
         diags = -scheme_radius : scheme_radius;
-        M = spdiagsVariablePeriodic(vals,diags); 
+        M = spdiagsPeriodic(vals,diags);
 
         M=M/h;
         D2=HI*(-M );

--- a/+sbp/+implementations/d2_variable_periodic_4.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+sbp/+implementations/d2_variable_periodic_4.m	Tue Dec 04 15:24:36 2018 -0800
@@ -30,7 +30,7 @@
 
     scheme_width = 5;
     scheme_radius = (scheme_width-1)/2;
-    
+
     r = 1:m;
     offset = scheme_width;
     r = r + offset;
@@ -47,7 +47,7 @@
 
         vals = -[Mm2,Mm1,M0,Mp1,Mp2];
         diags = -scheme_radius : scheme_radius;
-        M = spdiagsVariablePeriodic(vals,diags); 
+        M = spdiagsPeriodic(vals,diags);
 
         M=M/h;
         D2=HI*(-M );

--- a/+sbp/+implementations/d2_variable_periodic_6.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+sbp/+implementations/d2_variable_periodic_6.m	Tue Dec 04 15:24:36 2018 -0800
@@ -47,12 +47,12 @@
 
         vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3];
         diags = -scheme_radius : scheme_radius;
-        M = spdiagsVariablePeriodic(vals,diags); 
+        M = spdiagsPeriodic(vals,diags);
 
         M=M/h;
         D2=HI*(-M );
     end
     D2 = @D2_fun;
 
-    
+
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/+bc/closureSetup.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,25 @@
+% Setup closure and penalty matrices for several boundary conditions at once.
+% Each bc is a struct with the fields
+%  * type     -- Type of boundary condition
+%  * boundary -- Boundary identifier
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
+% Also takes S_sign which modifies the sign of the penalty function, [-1,1]
+% Returns a closure matrix and a penalty matrices for each boundary condition.
+%
+% The boundary data function can either be a function of time or a function of time and space coordinates.
+% In the case where it only depends on time it should return the data as grid function for the boundary.
+% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
+% For example in the 2D case: f(t,x,y).
+function [closure, penalties] = closureSetup(diffOp, bcs)
+    scheme.bc.verifyFormat(bcs, diffOp);
+
+    % Setup storage arrays
+    closure = spzeros(size(diffOp));
+    penalties = cell(1, length(bcs));
+
+    % Collect closures and penalties
+    for i = 1:length(bcs)
+        [localClosure, penalties{i}] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type);
+        closure = closure + localClosure;
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/+bc/forcingSetup.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,86 @@
+% Setup the forcing function for the given boundary conditions and data.
+% Each bc is a struct with the fields
+%  * type     -- Type of boundary condition
+%  * boundary -- Boundary identifier
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
+% S_sign allows changing the sign of the function to put on different sides in the system of ODEs.
+%   default is 1, which the same side as the diffOp.
+% Returns a forcing function S.
+%
+% The boundary data function can either be a function of time or a function of time and space coordinates.
+% In the case where it only depends on time it should return the data as grid function for the boundary.
+% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
+% For example in the 2D case: f(t,x,y).
+
+function S = forcingSetup(diffOp, penalties, bcs, S_sign)
+    default_arg('S_sign', 1);
+
+    assertType(bcs, 'cell');
+    assertIsMember(S_sign, [1, -1]);
+
+    scheme.bc.verifyFormat(bcs, diffOp);
+
+    [gridData, symbolicData] = parseAndSortData(bcs, penalties, diffOp);
+
+    % Setup penalty function
+    O = spzeros(size(diffOp),1);
+    function v = S_fun(t)
+        v = O;
+        for i = 1:length(gridData)
+            v = v + gridData{i}.penalty*gridData{i}.func(t);
+        end
+
+        for i = 1:length(symbolicData)
+            v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:});
+        end
+
+        v = S_sign * v;
+    end
+    S = @S_fun;
+end
+
+% Go through a cell array of boundary condition specifications and return cell arrays
+% of structs for grid and symbolic data.
+function [gridData, symbolicData] = parseAndSortData(bcs, penalties, diffOp)
+    gridData = {};
+    symbolicData = {};
+    for i = 1:length(bcs)
+        [ok, isSymbolic, data] = parseData(bcs{i}, penalties{i}, diffOp.grid);
+
+        if ~ok
+            continue % There was no data
+        end
+
+        if isSymbolic
+            symbolicData{end+1} = data;
+        else
+            gridData{end+1} = data;
+        end
+    end
+end
+
+function [ok, isSymbolic, dataStruct] = parseData(bc, penalty, grid)
+    if ~isfield(bc,'data') || isempty(bc.data)
+        isSymbolic = [];
+        dataStruct = struct();
+        ok = false;
+        return
+    end
+    ok = true;
+
+    nArg = nargin(bc.data);
+
+    if nArg > 1
+        % Symbolic data
+        isSymbolic = true;
+        coord = grid.getBoundary(bc.boundary);
+        dataStruct.penalty = penalty;
+        dataStruct.func = bc.data;
+        dataStruct.coords = num2cell(coord, 1);
+    else
+        % Grid data
+        isSymbolic = false;
+        dataStruct.penalty = penalty;
+        dataStruct.func = bc.data;
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/+bc/verifyFormat.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,31 @@
+% Errors with a more or less detailed error message if there is a problem with the bc specification
+function verifyBcFormat(bcs, diffOp)
+    assertType(bcs, 'cell');
+    for i = 1:length(bcs)
+        assertType(bcs{i}, 'struct');
+        assertStructFields(bcs{i}, {'type', 'boundary'});
+
+        if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data)
+            continue
+        end
+
+        if ~isa(bcs{i}.data, 'function_handle')
+            error('bcs{%d}.data should be a function of time or a function of time and space',i);
+        end
+
+        % Find dimension of boundary
+        b = diffOp.grid.getBoundary(bcs{i}.boundary);
+        dim = size(b,2);
+
+        % Assert that the data function has a valid number of input arguments
+        if ~(nargin(bcs{i}.data) == 1 || nargin(bcs{i}.data) == 1 + dim)
+            error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bcs{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i);
+        end
+
+        if nargin(bcs{i}.data) == 1
+            % Grid data (only function of time)
+            % Assert that the data has the correct dimension
+            assertSize(bcs{i}.data(0), 1, size(b,1));
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Elastic2dCurvilinear.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,621 @@
+classdef Elastic2dCurvilinear < scheme.Scheme
+
+% Discretizes the elastic wave equation in curvilinear coordinates.
+%
+% Untransformed equation:
+% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i 
+%
+% Transformed equation:
+% J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j 
+%                + dk J b_jk mu b_il dl u_j 
+%                + dk J b_jk mu b_jl dl 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
+
+        % Metric coefficients
+        b % Cell matrix of size dim x dim
+        J, Ji
+        beta % Cell array of scale factors
+
+        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_l, H_boundary_r % Boundary inner products
+
+        % Kroneckered norms and coefficients
+        RHOi_kron
+        Ji_kron, J_kron
+        Hi_kron, H_kron
+    end
+
+    methods
+
+        function obj = Elastic2dCurvilinear(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;
+
+            lambda = grid.evalOn(g, lambda_fun);
+            mu = grid.evalOn(g, mu_fun);
+            rho = grid.evalOn(g, rho_fun);
+            m = g.size();
+            obj.m = m;
+            m_tot = g.N();
+
+            % 1D operators
+            ops = cell(dim,1);
+            for i = 1:dim
+                ops{i} = opSet{i}(m(i), {0, 1}, 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 ========
+
+            % Variable coefficients
+            LAMBDA = spdiag(lambda);
+            obj.LAMBDA = LAMBDA;
+            MU = spdiag(mu);
+            obj.MU = MU;
+            RHO = spdiag(rho);
+            obj.RHO = RHO;
+            obj.RHOi = inv(RHO);
+
+            % Allocate
+            obj.D1 = cell(dim,1);
+            obj.D2_lambda = cell(dim,dim,dim);
+            obj.D2_mu = cell(dim,dim,dim);
+            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});
+
+            % -- Metric coefficients ----
+            coords = g.points();
+            x = coords(:,1);
+            y = coords(:,2);
+
+            % Use non-periodic difference operators for metric even if opSet is periodic.
+            xmax = max(ops{1}.x);
+            ymax = max(ops{2}.x);
+            opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
+            opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
+            D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
+            D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); 
+
+            x_xi = D1Metric{1}*x;
+            x_eta = D1Metric{2}*x;
+            y_xi = D1Metric{1}*y;
+            y_eta = D1Metric{2}*y;
+
+            J = x_xi.*y_eta - x_eta.*y_xi;
+
+            b = cell(dim,dim);
+            b{1,1} = y_eta./J;
+            b{1,2} = -x_eta./J;
+            b{2,1} = -y_xi./J;
+            b{2,2} = x_xi./J;
+
+            % Scale factors for boundary integrals
+            beta = cell(dim,1);
+            beta{1} = sqrt(x_eta.^2 + y_eta.^2);
+            beta{2} = sqrt(x_xi.^2 + y_xi.^2);
+
+            J = spdiag(J);
+            Ji = inv(J);
+            for i = 1:dim
+                beta{i} = spdiag(beta{i});
+                for j = 1:dim
+                    b{i,j} = spdiag(b{i,j});
+                end
+            end
+            obj.J = J;
+            obj.Ji = Ji;
+            obj.b = b;
+            obj.beta = beta;
+            %----------------------------
+
+            % 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
+                for j = 1:dim
+                    for k = 1:dim
+                        obj.D2_lambda{i,j,k} = sparse(m_tot);
+                        obj.D2_mu{i,j,k} = sparse(m_tot);
+                    end
+                end
+            end
+            ind = grid.funcToMatrix(g, 1:m_tot);
+
+            % x-dir
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1
+
+                        coeff_lambda = J*b{i,k}*b{j,k}*lambda;
+                        coeff_mu = J*b{j,k}*b{i,k}*mu;
+
+                        for col = 1:m(2)
+                            D_lambda = D2{1}(coeff_lambda(ind(:,col)));
+                            D_mu = D2{1}(coeff_mu(ind(:,col)));
+
+                            p = ind(:,col);
+                            obj.D2_lambda{i,j,k}(p,p) = D_lambda;
+                            obj.D2_mu{i,j,k}(p,p) = D_mu;
+                        end
+
+                    end
+                end
+            end
+
+            % y-dir
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 2
+
+                        coeff_lambda = J*b{i,k}*b{j,k}*lambda;
+                        coeff_mu = J*b{j,k}*b{i,k}*mu;
+
+                        for row = 1:m(1)
+                            D_lambda = D2{2}(coeff_lambda(ind(row,:)));
+                            D_mu = D2{2}(coeff_mu(ind(row,:)));
+
+                            p = ind(row,:);
+                            obj.D2_lambda{i,j,k}(p,p) = D_lambda;
+                            obj.D2_mu{i,j,k}(p,p) = D_mu;
+                        end
+
+                    end
+                end
+            end
+
+            % Quadratures
+            obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
+            obj.H_boundary_l = cell(dim,1);
+            obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2};
+            obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1};
+            obj.H_boundary_r = cell(dim,1);
+            obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2};
+            obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{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
+                    for k = 1:dim
+                        for l = 1:dim
+                            D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_lambda{i,j,k}*E{j}' + ...
+                                                      db(k,l)*D1{k}*J*b{i,k}*b{j,l}*LAMBDA*D1{l}*E{j}' ...
+                                                  );
+
+                            D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{i,j,k}*E{j}' + ...
+                                                      db(k,l)*D1{k}*J*b{j,k}*b{i,l}*MU*D1{l}*E{j}' ...
+                                                  );
+
+                            D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{j,j,k}*E{i}' + ...
+                                                      db(k,l)*D1{k}*J*b{j,k}*b{j,l}*MU*D1{l}*E{i}' ...
+                                                  );
+
+                        end
+                    end
+                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;
+
+            % 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);
+
+                    % Loop over components that T_{ik}^{(j)} acts on
+                    for k = 1:dim
+
+                        T_l{j}{i,k} = sparse(m_tot,m_tot);
+                        T_r{j}{i,k} = sparse(m_tot,m_tot);
+
+                        for m = 1:dim
+                            for l = 1:dim
+                                T_l{j}{i,k} = T_l{j}{i,k} + ... 
+                                -d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
+                                -d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
+                                -d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m});
+
+                                T_r{j}{i,k} = T_r{j}{i,k} + ... 
+                                d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
+                                d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
+                                d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m});
+                            end
+                        end
+
+                        T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k};
+                        T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k}; 
+
+                        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.Ji_kron = kron(obj.Ji, I_dim);
+            obj.Hi_kron = kron(obj.Hi, I_dim);
+            obj.H_kron = kron(obj.H, I_dim);
+            obj.J_kron = kron(obj.J, I_dim);
+
+            % Misc.
+            obj.h = g.scaling();
+            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'.
+        %       bc                  is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
+        %                           on the first 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, bc, tuning)
+            default_arg('tuning', 1.2);
+
+            assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
+            comp = bc{1};
+            type = bc{2};
+
+            % j is the coordinate direction of the boundary
+            j = obj.get_boundary_number(boundary);
+            [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
+
+            E = obj.E;
+            Hi = obj.Hi;
+            LAMBDA = obj.LAMBDA;
+            MU = obj.MU;
+            RHOi = obj.RHOi;
+            Ji = obj.Ji;
+
+            dim = obj.dim;
+            m_tot = obj.grid.N();
+
+            % Preallocate
+            closure = sparse(dim*m_tot, dim*m_tot);
+            penalty = sparse(dim*m_tot, m_tot/obj.m(j));
+
+            % Loop over components that we (potentially) have different BC on
+            k = comp;
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
+
+                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*Ji*B'*e*H_gamma*(e'*E{k}' ); 
+                    penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma;
+                end 
+
+            % Free boundary condition
+            case {'F','f','Free','free','traction','Traction','t','T'}
+                    closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} ); 
+                    penalty = penalty + E{k}*RHOi*Ji*Hi*e*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
+            % Operators without subscripts are from the own domain.
+            error('Not implemented');
+            tuning = 1.2;
+
+            % j is the coordinate direction of the boundary
+            j = obj.get_boundary_number(boundary);
+            j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
+
+            % Get boundary operators
+            [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
+            [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
+
+            % Operators and quantities that correspond to the own domain only
+            Hi = obj.Hi;
+            RHOi = obj.RHOi;
+            dim = obj.dim;
+        
+            %--- Other operators ----
+            m_tot_u = obj.grid.N();
+            E = obj.E;
+            LAMBDA_u = obj.LAMBDA;
+            MU_u = obj.MU;
+            lambda_u = e'*LAMBDA_u*e;
+            mu_u = e'*MU_u*e;
+
+            m_tot_v = neighbour_scheme.grid.N();
+            E_v = neighbour_scheme.E;
+            LAMBDA_v = neighbour_scheme.LAMBDA;
+            MU_v = neighbour_scheme.MU;
+            lambda_v = e_v'*LAMBDA_v*e_v;
+            mu_v = e_v'*MU_v*e_v;
+            %-------------------------
+            
+            % Borrowing constants
+            phi_u = obj.phi{j};
+            h_u = obj.h(j);
+            h11_u = obj.H11{j}*h_u;
+            gamma_u = obj.gamma{j};
+
+            phi_v = neighbour_scheme.phi{j_v};
+            h_v = neighbour_scheme.h(j_v);
+            h11_v = neighbour_scheme.H11{j_v}*h_v;
+            gamma_v = neighbour_scheme.gamma{j_v};
+
+            % E > sum_i 1/(2*alpha_ij)*(tau_i)^2
+            function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) 
+                th1 = h11/(2*dim);
+                th2 = h11*phi/2;
+                th3 = h*gamma;
+                a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3);
+                a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3);
+                alpha_ii = a1 + sqrt(a2 + a1^2);
+
+                alpha_ij = mu*(2/h11 + 1/(phi*h11));
+            end
+
+            [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
+            [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);  
+            sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
+            sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
+
+            d = @kroneckerDelta;  % Kronecker delta
+            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+            sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
+
+            % Preallocate
+            closure = sparse(dim*m_tot_u, dim*m_tot_u);
+            penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+
+            % Loop over components that penalties end up on
+            for i = 1:dim
+                closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
+                penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
+
+                closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
+                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
+
+                % Loop over components that we have interface conditions on
+                for k = 1:dim
+                    closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}'; 
+                    penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}'; 
+                end 
+            end
+        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 boundary operator op for the boundary specified by the string boundary.
+        % op: may be a cell array of strings
+        function [varargout] = 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
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:length(op)
+                switch op{i}
+                    case 'e'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.e_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.e_r{j};
+                        end
+                    case 'd'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.d1_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.d1_r{j};
+                        end
+                    case 'H'
+                        switch boundary 
+                            case {'w','W','west','West','s','S','south','South'}
+                                    varargout{i} = obj.H_boundary_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                    varargout{i} = obj.H_boundary_r{j};
+                        end
+                    case 'T'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.T_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.T_r{j};
+                        end
+                    case 'tau'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.tau_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.tau_r{j};
+                        end                        
+                    otherwise
+                        error(['No such operator: operator = ' op{i}]);
+                end
+            end
+
+        end
+
+        function N = size(obj)
+            N = obj.dim*prod(obj.m);
+        end
+    end
+end

--- a/+scheme/Elastic2dVariable.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+scheme/Elastic2dVariable.m	Tue Dec 04 15:24:36 2018 -0800
@@ -31,9 +31,16 @@
         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
+
+        % Borrowing from H, M, and R
+        thH
+        thM
+        thR
+
         e_l, e_r
         d1_l, d1_r % Normal derivatives at the boundary
         E % E{i}^T picks out component i
@@ -64,6 +71,13 @@
 
             h = g.scaling();
             lim = g.lim;
+            if isempty(lim)
+                x = g.x;
+                lim = cell(length(x),1);
+                for i = 1:length(x)
+                    lim{i} = {min(x{i}), max(x{i})};
+                end
+            end
 
             % 1D operators
             ops = cell(dim,1);
@@ -77,6 +91,11 @@
                 obj.H11{i} = ops{i}.borrowing.H11;
                 obj.phi{i} = beta/obj.H11{i};
                 obj.gamma{i} = ops{i}.borrowing.M.d1;
+
+                % Better names
+                obj.thR{i} = ops{i}.borrowing.R.delta_D;
+                obj.thM{i} = ops{i}.borrowing.M.d1;
+                obj.thH{i} = ops{i}.borrowing.H11;
             end
 
             I = cell(dim,1);
@@ -262,35 +281,24 @@
         % 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.
+        %       bc                  is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
+        %                           on the first 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', []);
+        function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
+            default_arg('tuning', 1.2);
+
+            assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
+            comp = bc{1};
+            type = bc{2};
 
             % 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
+            j = obj.get_boundary_number(boundary);
+            [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
 
             E = obj.E;
             Hi = obj.Hi;
-            H_gamma = obj.H_boundary{j};
             LAMBDA = obj.LAMBDA;
             MU = obj.MU;
             RHOi = obj.RHOi;
@@ -298,66 +306,127 @@
             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
+            penalty = sparse(dim*m_tot, m_tot/obj.m(j));
 
-            % Loop over components that we (potentially) have different BC on
-            for k = 1:dim
-                switch type{k}
+            k = comp;
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
 
-                % Dirichlet boundary condition
-                case {'D','d','dirichlet','Dirichlet'}
+                phi = obj.phi{j};
+                h = obj.h(j);
+                h11 = obj.H11{j}*h;
+                gamma = obj.gamma{j};
 
-                    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);
+                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 );
+                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
+                % 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*H_gamma*(e'*E{k}' );
+                    penalty = penalty - E{i}*RHOi*Hi*B'*e*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;
+            % Free boundary condition
+            case {'F','f','Free','free','traction','Traction','t','T'}
+                    closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
+                    penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
 
-                % Unknown boundary condition
-                otherwise
-                    error('No such boundary condition: type = %s',type);
-                end
+            % Unknown boundary condition
+            otherwise
+                error('No such boundary condition: type = %s',type);
             end
         end
 
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
+        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
+            % Operators without subscripts are from the own domain.
             tuning = 1.2;
-            % tuning = 20.2;
-            error('Interface not implemented');
+
+            % j is the coordinate direction of the boundary
+            j = obj.get_boundary_number(boundary);
+            j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
+
+            % Get boundary operators
+            [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
+            [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
+
+            % Operators and quantities that correspond to the own domain only
+            Hi = obj.Hi;
+            RHOi = obj.RHOi;
+            dim = obj.dim;
+
+            %--- Other operators ----
+            m_tot_u = obj.grid.N();
+            E = obj.E;
+            LAMBDA_u = obj.LAMBDA;
+            MU_u = obj.MU;
+            lambda_u = e'*LAMBDA_u*e;
+            mu_u = e'*MU_u*e;
+
+            m_tot_v = neighbour_scheme.grid.N();
+            E_v = neighbour_scheme.E;
+            LAMBDA_v = neighbour_scheme.LAMBDA;
+            MU_v = neighbour_scheme.MU;
+            lambda_v = e_v'*LAMBDA_v*e_v;
+            mu_v = e_v'*MU_v*e_v;
+            %-------------------------
+
+            % Borrowing constants
+            h_u = obj.h(j);
+            thR_u = obj.thR{j}*h_u;
+            thM_u = obj.thM{j}*h_u;
+            thH_u = obj.thH{j}*h_u;
+
+            h_v = neighbour_scheme.h(j_v);
+            thR_v = neighbour_scheme.thR{j_v}*h_v;
+            thH_v = neighbour_scheme.thH{j_v}*h_v;
+            thM_v = neighbour_scheme.thM{j_v}*h_v;
+
+            % alpha = penalty strength for normal component, beta for tangential
+            alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u);
+            alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v);
+            beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u);
+            beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v);
+            alpha = alpha_u + alpha_v;
+            beta = beta_u + beta_v;
+
+            d = @kroneckerDelta;  % Kronecker delta
+            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+            strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta);
+
+            % Preallocate
+            closure = sparse(dim*m_tot_u, dim*m_tot_u);
+            penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+
+            % Loop over components that penalties end up on
+            for i = 1:dim
+                closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}';
+                penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}';
+
+                closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
+                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
+
+                % Loop over components that we have interface conditions on
+                for k = 1:dim
+                    closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
+                    penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
+                end
+            end
         end
 
         % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
@@ -380,8 +449,9 @@
             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)
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op: may be a cell array of strings
+        function [varargout] = get_boundary_operator(obj, op, boundary)
 
             switch boundary
                 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
@@ -392,29 +462,51 @@
                     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]);
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:length(op)
+                switch op{i}
+                    case 'e'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.e_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.e_r{j};
+                        end
+                    case 'd'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.d1_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.d1_r{j};
+                        end
+                    case 'H'
+                        varargout{i} = obj.H_boundary{j};
+                    case 'T'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.T_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.T_r{j};
+                        end
+                    case 'tau'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.tau_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.tau_r{j};
+                        end
+                    otherwise
+                        error(['No such operator: operator = ' op{i}]);
+                end
             end
 
         end
 
         function N = size(obj)
-            N = prod(obj.m);
+            N = obj.dim*prod(obj.m);
         end
     end
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Heat2dCurvilinear.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,385 @@
+classdef Heat2dCurvilinear < scheme.Scheme
+
+% Discretizes the Laplacian with variable coefficent, curvilinear,
+% 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
+        alpha % Vector of borrowing constants
+        
+        % Boundary inner products
+        H_boundary_l, H_boundary_r 
+
+        % Metric coefficients
+        b % Cell matrix of size dim x dim
+        J, Ji
+        beta % Cell array of scale factors
+
+        % Numerical boundary flux operators
+        flux_l, flux_r
+
+    end
+
+    methods
+
+        function obj = Heat2dCurvilinear(g ,order, kappa_fun, opSet)
+            default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
+            default_arg('kappa_fun', @(x,y) 0*x+1);
+            dim = 2;
+
+            kappa = grid.evalOn(g, kappa_fun);
+            m = g.size();
+            m_tot = g.N();
+
+            % 1D operators
+            ops = cell(dim,1);
+            for i = 1:dim
+                ops{i} = opSet{i}(m(i), {0, 1}, 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;
+
+            % Allocate
+            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});
+
+            % -- Metric coefficients ----
+            coords = g.points();
+            x = coords(:,1);
+            y = coords(:,2);
+
+            % Use non-periodic difference operators for metric even if opSet is periodic.
+            xmax = max(ops{1}.x);
+            ymax = max(ops{2}.x);
+            opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
+            opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
+            D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
+            D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); 
+
+            x_xi = D1Metric{1}*x;
+            x_eta = D1Metric{2}*x;
+            y_xi = D1Metric{1}*y;
+            y_eta = D1Metric{2}*y;
+
+            J = x_xi.*y_eta - x_eta.*y_xi;
+
+            b = cell(dim,dim);
+            b{1,1} = y_eta./J;
+            b{1,2} = -x_eta./J;
+            b{2,1} = -y_xi./J;
+            b{2,2} = x_xi./J;
+
+            % Scale factors for boundary integrals
+            beta = cell(dim,1);
+            beta{1} = sqrt(x_eta.^2 + y_eta.^2);
+            beta{2} = sqrt(x_xi.^2 + y_xi.^2);
+
+            J = spdiag(J);
+            Ji = inv(J);
+            for i = 1:dim
+                beta{i} = spdiag(beta{i});
+                for j = 1:dim
+                    b{i,j} = spdiag(b{i,j});
+                end
+            end
+            obj.J = J;
+            obj.Ji = Ji;
+            obj.b = b;
+            obj.beta = beta;
+            %----------------------------
+
+            % 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 coefficients
+            kappa_coeff = cell(dim,dim);
+            for j = 1:dim
+                obj.D2_kappa{j} = sparse(m_tot,m_tot); 
+                kappa_coeff{j} = sparse(m_tot,1);
+                for i = 1:dim
+                    kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa;
+                end
+            end
+            ind = grid.funcToMatrix(g, 1:m_tot);
+
+            % x-dir
+            j = 1;
+            for col = 1:m(2)
+                D_kappa = D2{1}(kappa_coeff{j}(ind(:,col)));
+
+                p = ind(:,col);
+                obj.D2_kappa{j}(p,p) = D_kappa;
+            end
+
+            % y-dir
+            j = 2;
+            for row = 1:m(1)
+                D_kappa = D2{2}(kappa_coeff{j}(ind(row,:)));
+
+                p = ind(row,:);
+                obj.D2_kappa{j}(p,p) = D_kappa;
+            end
+
+            % Quadratures
+            obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
+            obj.H_boundary_l = cell(dim,1);
+            obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2};
+            obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1};
+            obj.H_boundary_r = cell(dim,1);
+            obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2};
+            obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1};
+
+            %=== Differentiation matrix D (without SAT) ===
+            D2_kappa = obj.D2_kappa;
+            D1 = obj.D1;
+            D = sparse(m_tot,m_tot);
+
+            d = @kroneckerDelta;    % Kronecker delta
+            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+
+            % 2nd derivatives
+            for j = 1:dim
+                D = D + Ji*D2_kappa{j};
+            end
+
+            % Mixed terms
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        D = D + db(i,j)*Ji*D1{j}*b{i,j}*J*KAPPA*b{i,k}*D1{k};
+                    end
+                end
+            end
+            obj.D = D;
+            %=========================================%
+
+            % Normal flux operators for BC.
+            flux_l = cell(dim,1);
+            flux_r = cell(dim,1);
+
+            d1_l = obj.d1_l;
+            d1_r = obj.d1_r;
+            e_l = obj.e_l;
+            e_r = obj.e_r;
+
+            % Loop over boundaries
+            for j = 1:dim
+                flux_l{j} = sparse(m_tot,m_tot);
+                flux_r{j} = sparse(m_tot,m_tot);
+
+                % Loop over dummy index
+                for i = 1:dim
+                    % Loop over dummy index
+                    for k = 1:dim
+                        flux_l{j} = flux_l{j} ...
+                                  - beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_l{k}*d1_l{k}' + db(j,k)*D1{k} );
+
+                        flux_r{j} = flux_r{j} ...
+                                  + beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_r{k}*d1_r{k}' + db(j,k)*D1{k} );
+                    end
+
+                end
+            end
+            obj.flux_l = flux_l;
+            obj.flux_r = flux_r;
+
+            % Misc.
+            obj.m = m;
+            obj.h = g.scaling();
+            obj.order = order;
+            obj.grid = g;
+            obj.dim = dim;
+            obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
+
+        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, symmetric, tuning)
+            default_arg('type','Neumann');
+            default_arg('symmetric', false);
+            default_arg('tuning',1.2);
+
+            % 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{j};
+                flux = obj.flux_r{j};
+                H_gamma = obj.H_boundary_r{j};
+            case -1
+                e = obj.e_l{j};
+                flux = obj.flux_l{j};
+                H_gamma = obj.H_boundary_l{j};
+            end
+
+            Hi = obj.Hi;
+            Ji = obj.Ji;
+            KAPPA = obj.KAPPA;
+            kappa_gamma = e'*KAPPA*e; 
+            h = obj.h(j);
+            alpha = h*obj.alpha(j);
+
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
+
+                if ~symmetric
+                    closure = -Ji*Hi*flux'*e*H_gamma*(e' ); 
+                    penalty = Ji*Hi*flux'*e*H_gamma;
+                else
+                    closure = Ji*Hi*flux'*e*H_gamma*(e' )...
+                              -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; 
+                    penalty =  -Ji*Hi*flux'*e*H_gamma ...
+                              +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma;
+                end
+
+            % Normal flux boundary condition
+            case {'N','n','neumann','Neumann'}
+                    closure = -Ji*Hi*e*H_gamma*(e'*flux ); 
+                    penalty =  Ji*Hi*e*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/Heat2dVariable.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+scheme/Heat2dVariable.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,9 +1,9 @@
 classdef Heat2dVariable < scheme.Scheme
 
 % Discretizes the Laplacian with variable coefficent,
-% In the Heat equation way (i.e., the discretization matrix is not necessarily
+% In the Heat equation way (i.e., the discretization matrix is not necessarily 
 % symmetric)
-% u_t = div * (kappa * grad u )
+% 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.
 
@@ -28,7 +28,8 @@
         H, Hi % Inner products
         e_l, e_r
         d1_l, d1_r % Normal derivatives at the boundary
-
+        alpha % Vector of borrowing constants
+        
         H_boundary % Boundary inner products
 
     end
@@ -144,6 +145,7 @@
             obj.order = order;
             obj.grid = g;
             obj.dim = dim;
+            obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
 
         end
 
@@ -155,12 +157,13 @@
         %       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)
+        function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
             default_arg('type','Neumann');
-            default_arg('parameter', []);
+            default_arg('symmetric', false);
+            default_arg('tuning',1.2);
 
             % j is the coordinate direction of the boundary
-            % nj: outward unit normal component.
+            % 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);
@@ -176,19 +179,30 @@
             Hi = obj.Hi;
             H_gamma = obj.H_boundary{j};
             KAPPA = obj.KAPPA;
-            kappa_gamma = e{j}'*KAPPA*e{j};
+            kappa_gamma = e{j}'*KAPPA*e{j}; 
+            h = obj.h(j);
+            alpha = h*obj.alpha(j);
 
             switch type
 
             % Dirichlet boundary condition
             case {'D','d','dirichlet','Dirichlet'}
-                    closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
+
+                if ~symmetric
+                    closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 
                     penalty =  nj*Hi*d{j}*kappa_gamma*H_gamma;
+                else
+                    closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )...
+                              -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ; 
+                    penalty =  -nj*Hi*d{j}*kappa_gamma*H_gamma ...
+                              +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma;
+                end
 
             % 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;
+                    closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 
+                    penalty =  Hi*e{j}*kappa_gamma*H_gamma; 
+                    % penalty is for normal derivative and not for derivative, hence the sign.
 
             % Unknown boundary condition
             otherwise
@@ -196,7 +210,7 @@
             end
         end
 
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
+        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');

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Laplace1D.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,146 @@
+classdef Laplace1D < scheme.Scheme
+    properties
+        grid
+        order % Order accuracy for the approximation
+
+        D % non-stabalized scheme operator
+        H % Discrete norm
+        M % Derivative norm
+        a
+
+        D2
+        Hi
+        e_l
+        e_r
+        d_l
+        d_r
+        gamm
+    end
+
+    methods
+        function obj = Laplace1D(grid, order, a)
+            default_arg('a', 1);
+
+            assertType(grid, 'grid.Cartesian');
+
+            ops = sbp.D2Standard(grid.size(), grid.lim{1}, order);
+
+            obj.D2 = sparse(ops.D2);
+            obj.H =  sparse(ops.H);
+            obj.Hi = sparse(ops.HI);
+            obj.M =  sparse(ops.M);
+            obj.e_l = sparse(ops.e_l);
+            obj.e_r = sparse(ops.e_r);
+            obj.d_l = -sparse(ops.d1_l);
+            obj.d_r = sparse(ops.d1_r);
+
+
+            obj.grid = grid;
+            obj.order = order;
+
+            obj.a = a;
+            obj.D = a*obj.D2;
+
+            obj.gamm = grid.h*ops.borrowing.M.S;
+        end
+
+
+        % Closure functions return the opertors applied to the own doamin to close the boundary
+        % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
+        %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
+        %       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,data)
+            default_arg('type','neumann');
+            default_arg('data',0);
+
+            [e,d,s] = obj.get_boundary_ops(boundary);
+
+            switch type
+                % Dirichlet boundary condition
+                case {'D','dirichlet'}
+                    tuning = 1.1;
+                    tau1 = -tuning/obj.gamm;
+                    tau2 =  1;
+
+                    tau = tau1*e + tau2*d;
+
+                    closure = obj.a*obj.Hi*tau*e';
+                    penalty = obj.a*obj.Hi*tau;
+
+                % Neumann boundary condition
+                case {'N','neumann'}
+                    tau = -e;
+
+                    closure = obj.a*obj.Hi*tau*d';
+                    penalty = -obj.a*obj.Hi*tau;
+
+                % Unknown, boundary condition
+                otherwise
+                    error('No such boundary condition: type = %s',type);
+            end
+        end
+
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+
+            [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
+            [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+
+
+            a_u = obj.a;
+            a_v = neighbour_scheme.a;
+
+            gamm_u = obj.gamm;
+            gamm_v = neighbour_scheme.gamm;
+
+            tuning = 1.1;
+            
+            tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning;
+            tau2 = 1/2*a_u;
+            sig1 = -1/2;
+            sig2 = 0;
+
+            tau = tau1*e_u + tau2*d_u;
+            sig = sig1*e_u + sig2*d_u;
+
+            closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
+            penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
+        end
+
+        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
+        % The right boundary is considered the positive boundary
+        function [e,d,s] = get_boundary_ops(obj,boundary)
+            switch boundary
+                case 'l'
+                    e = obj.e_l;
+                    d = obj.d_l;
+                    s = -1;
+                case 'r'
+                    e = obj.e_r;
+                    d = obj.d_r;
+                    s = 1;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+        end
+
+        function N = size(obj)
+            N = obj.grid.size();
+        end
+
+    end
+
+    methods(Static)
+        % Calculates the matrcis need 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')
+        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

--- a/+scheme/Wave.m	Tue Dec 04 14:54:28 2018 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,175 +0,0 @@
-classdef Wave < scheme.Scheme
-    properties
-        m % Number of points in each direction, possibly a vector
-        h % Grid spacing
-        x % Grid
-        order % Order accuracy for the approximation
-
-        D % non-stabalized scheme operator
-        H % Discrete norm
-        M % Derivative norm
-        alpha
-
-        D2
-        Hi
-        e_l
-        e_r
-        d1_l
-        d1_r
-        gamm
-    end
-
-    methods
-        function obj = Wave(m,xlim,order,alpha)
-            default_arg('a',1);
-            [x, h] = util.get_grid(xlim{:},m);
-
-            ops = sbp.Ordinary(m,h,order);
-
-            obj.D2 = sparse(ops.derivatives.D2);
-            obj.H =  sparse(ops.norms.H);
-            obj.Hi = sparse(ops.norms.HI);
-            obj.M =  sparse(ops.norms.M);
-            obj.e_l = sparse(ops.boundary.e_1);
-            obj.e_r = sparse(ops.boundary.e_m);
-            obj.d1_l = sparse(ops.boundary.S_1);
-            obj.d1_r = sparse(ops.boundary.S_m);
-
-
-            obj.m = m;
-            obj.h = h;
-            obj.order = order;
-
-            obj.alpha = alpha;
-            obj.D = alpha*obj.D2;
-            obj.x = x;
-
-            obj.gamm = h*ops.borrowing.M.S;
-
-        end
-
-
-        % Closure functions return the opertors applied to the own doamin to close the boundary
-        % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
-        %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
-        %       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,data)
-            default_arg('type','neumann');
-            default_arg('data',0);
-
-            [e,d,s] = obj.get_boundary_ops(boundary);
-
-            switch type
-                % Dirichlet boundary condition
-                case {'D','dirichlet'}
-                    alpha = obj.alpha;
-
-                    % tau1 < -alpha^2/gamma
-                    tuning = 1.1;
-                    tau1 = -tuning*alpha/obj.gamm;
-                    tau2 =  s*alpha;
-
-                    p = tau1*e + tau2*d;
-
-                    closure = obj.Hi*p*e';
-
-                    pp = obj.Hi*p;
-                    switch class(data)
-                        case 'double'
-                            penalty = pp*data;
-                        case 'function_handle'
-                            penalty = @(t)pp*data(t);
-                        otherwise
-                            error('Wierd data argument!')
-                    end
-
-
-                % Neumann boundary condition
-                case {'N','neumann'}
-                    alpha = obj.alpha;
-                    tau1 = -s*alpha;
-                    tau2 = 0;
-                    tau = tau1*e + tau2*d;
-
-                    closure = obj.Hi*tau*d';
-
-                    pp = obj.Hi*tau;
-                    switch class(data)
-                        case 'double'
-                            penalty = pp*data;
-                        case 'function_handle'
-                            penalty = @(t)pp*data(t);
-                        otherwise
-                            error('Wierd data argument!')
-                    end
-
-                % Unknown, boundary condition
-                otherwise
-                    error('No such boundary condition: type = %s',type);
-            end
-        end
-
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
-            % u denotes the solution in the own domain
-            % v denotes the solution in the neighbour domain
-            [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
-            [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-
-            tuning = 1.1;
-
-            alpha_u = obj.alpha;
-            alpha_v = neighbour_scheme.alpha;
-
-            gamm_u = obj.gamm;
-            gamm_v = neighbour_scheme.gamm;
-
-            % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v)
-
-            tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning;
-            tau2 = s_u*1/2*alpha_u;
-            sig1 = s_u*(-1/2);
-            sig2 = 0;
-
-            tau = tau1*e_u + tau2*d_u;
-            sig = sig1*e_u + sig2*d_u;
-
-            closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u');
-            penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v');
-        end
-
-        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
-        % The right boundary is considered the positive boundary
-        function [e,d,s] = get_boundary_ops(obj,boundary)
-            switch boundary
-                case 'l'
-                    e = obj.e_l;
-                    d = obj.d1_l;
-                    s = -1;
-                case 'r'
-                    e = obj.e_r;
-                    d = obj.d1_r;
-                    s = 1;
-                otherwise
-                    error('No such boundary: boundary = %s',boundary);
-            end
-        end
-
-        function N = size(obj)
-            N = obj.m;
-        end
-
-    end
-
-    methods(Static)
-        % Calculates the matrcis need 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')
-        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

--- a/+scheme/bcSetup.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+scheme/bcSetup.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,48 +1,20 @@
-% function [closure, S] = bcSetup(diffOp, bc)
 % Takes a diffOp and a cell array of boundary condition definitions.
 % Each bc is a struct with the fields
 %  * type     -- Type of boundary condition
 %  * boundary -- Boundary identifier
-%  * data     -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem
-% Also takes S_sign which modifies the sign of S, [-1,1]
-% Returns a closure matrix and a forcing function S
-function [closure, S] = bcSetup(diffOp, bc, S_sign)
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
+% Also takes S_sign which modifies the sign of the penalty function, [-1,1]
+% Returns a closure matrix and a forcing function S.
+%
+% The boundary data function can either be a function of time or a function of time and space coordinates.
+% In the case where it only depends on time it should return the data as grid function for the boundary.
+% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
+% For example in the 2D case: f(t,x,y).
+function [closure, S] = bcSetup(diffOp, bcs, S_sign)
     default_arg('S_sign', 1);
-    assertType(bc, 'cell');
+    assertType(bcs, 'cell');
     assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
 
-
-    closure = spzeros(size(diffOp));
-    penalties = {};
-    dataFunctions = {};
-    dataParams = {};
-
-    for i = 1:length(bc)
-        assertType(bc{i}, 'struct');
-        [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
-        closure = closure + localClosure;
-
-        if isempty(bc{i}.data)
-            continue
-        end
-        assertType(bc{i}.data, 'function_handle');
-
-        coord = diffOp.grid.getBoundary(bc{i}.boundary);
-        assertNumberOfArguments(bc{i}.data, 1+size(coord,2));
-
-        penalties{end+1} = penalty;
-        dataFunctions{end+1} = bc{i}.data;
-        dataParams{end+1} = num2cell(coord ,1);
-    end
-
-    O = spzeros(size(diffOp),1);
-    function v = S_fun(t)
-        v = O;
-        for i = 1:length(dataFunctions)
-            v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:});
-        end
-
-        v = S_sign * v;
-    end
-    S = @S_fun;
+    [closure, penalties] = scheme.bc.closureSetup(diffOp, bcs);
+    S = scheme.bc.forcingSetup(diffOp, penalties, bcs, S_sign);
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeScaled.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,139 @@
+classdef SBPInTimeScaled < time.Timestepper
+    % The SBP in time method.
+    % Implemented for A*v_t = B*v + f(t), v(0) = v0
+    % The resulting system of equations is
+    %   M*u_next= K*u_prev_end + f
+    properties
+        A,B
+        f
+
+        k % total time step.
+
+        blockSize % number of points in each block
+        N % Number of components
+
+        order
+        nodes
+
+        Mtilde,Ktilde     % System matrices
+        L,U,p,q % LU factorization of M
+        e_T
+
+        scaling
+        S, Sinv % Scaling matrices
+
+        % Time state
+        t
+        vtilde
+        n
+    end
+
+    methods
+        function obj = SBPInTimeScaled(A, B, f, k, t0, v0, scaling, TYPE, order, blockSize)
+            default_arg('TYPE','gauss');
+            default_arg('f',[]);
+
+            if(strcmp(TYPE,'gauss'))
+                default_arg('order',4)
+                default_arg('blockSize',4)
+            else
+                default_arg('order', 8);
+                default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE));
+            end
+
+            obj.A = A;
+            obj.B = B;
+            obj.scaling = scaling;
+
+            if ~isempty(f)
+                obj.f = f;
+            else
+                obj.f = @(t)sparse(length(v0),1);
+            end
+
+            obj.k = k;
+            obj.blockSize = blockSize;
+            obj.N = length(v0);
+
+            obj.n = 0;
+            obj.t = t0;
+
+            %==== Build the time discretization matrix =====%
+            switch TYPE
+                case 'equidistant'
+                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                case 'optimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                case 'minimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                case 'gauss'
+                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+            end
+
+            I = speye(size(A));
+            I_t = speye(blockSize,blockSize);
+
+            D1 = kron(ops.D1, I);
+            HI = kron(ops.HI, I);
+            e_0 = kron(ops.e_l, I);
+            e_T = kron(ops.e_r, I);
+            obj.nodes = ops.x;
+
+            % Convert to form M*w = K*v0 + f(t)
+            tau = kron(I_t, A) * e_0;
+            M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B);
+
+            K = HI*tau;
+
+            obj.S =    kron(I_t, spdiag(scaling));
+            obj.Sinv = kron(I_t, spdiag(1./scaling));
+
+            obj.Mtilde = obj.Sinv*M*obj.S;
+            obj.Ktilde = obj.Sinv*K*spdiag(scaling);
+            obj.e_T = e_T;
+
+
+            % LU factorization
+            [obj.L,obj.U,obj.p,obj.q] = lu(obj.Mtilde, 'vector');
+
+            obj.vtilde = (1./obj.scaling).*v0;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.scaling.*obj.vtilde;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            forcing = zeros(obj.blockSize*obj.N,1);
+
+            for i = 1:obj.blockSize
+                forcing((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i));
+            end
+
+            RHS = obj.Sinv*forcing + obj.Ktilde*obj.vtilde;
+
+            y = obj.L\RHS(obj.p);
+            z = obj.U\y;
+
+            w = zeros(size(z));
+            w(obj.q) = z;
+
+            obj.vtilde = obj.e_T'*w;
+
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+        end
+    end
+
+    methods(Static)
+        function N = smallestBlockSize(order,TYPE)
+            default_arg('TYPE','gauss')
+
+            switch TYPE
+                case 'gauss'
+                    N = 4;
+            end
+        end
+    end
+end

--- a/+time/SBPInTimeSecondOrderFormImplicit.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Tue Dec 04 15:24:36 2018 -0800
@@ -14,15 +14,16 @@
         % Solves A*u_tt + B*u_t + C*u = f(t)
         % A, B can either both be constants or both be function handles,
         % They can also be omitted by setting them equal to the empty matrix.
-        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize)
+        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, do_scaling, TYPE, order, blockSize)
             default_arg('f', []);
             default_arg('TYPE', []);
             default_arg('order', []);
             default_arg('blockSize',[]);
+            default_arg('do_scaling', false);
 
             m = length(v0);
 
-            default_arg('A', sparse(m, m));
+            default_arg('A', speye(m, m));
             default_arg('B', sparse(m, m));
             default_arg('C', sparse(m, m));
 
@@ -56,7 +57,12 @@
             obj.t = t0;
             obj.n = 0;
 
-            obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+            if do_scaling
+                scaling = [ones(m,1); sqrt(diag(C))];
+                obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize);
+            else
+                obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+            end
         end
 
         function [v,t] = getV(obj)

--- a/+util/ReplaceableString.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+util/ReplaceableString.m	Tue Dec 04 15:24:36 2018 -0800
@@ -58,3 +58,5 @@
 function b = padStr(a, n)
     b = sprintf('%-*s', n, a);
 end
+
+% TODO: Add a debug mode which prints without replacing?

--- a/.hgtags	Tue Dec 04 14:54:28 2018 -0800
+++ b/.hgtags	Tue Dec 04 15:24:36 2018 -0800
@@ -1,1 +1,4 @@
 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1
+0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1
+b723495cdb2f96314d7b3f0aa79723a7dc088c7d v0.2
+08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0

--- a/Color.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/Color.m	Tue Dec 04 15:24:36 2018 -0800
@@ -10,6 +10,10 @@
         black     = [0.000 0.000 0.000];
         white     = [1.000 1.000 1.000];
         colors = { Color.blue, Color.red, Color.yellow, Color.green, Color.purple, Color.lightblue, Color.darkred, Color.black, Color.white};
+        markers = {'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
+        lineStyles = {'-', '--', ':', '-.'};
+
+        solidMarkers = {'o', 'square', 'diamond', 'v', 'pentagram', '^', '>', '<', 'hexagram'};
 
         notabilityYellow     = [100.0   99.0    22.0    ]/100;
         notabilityOrange     = [97.0    61.0    15.0    ]/100;
@@ -34,13 +38,11 @@
 
     methods(Static)
         function sample()
-            markers ={'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
             % Filled and non-filled markers?
-            lineStyles = {'-', '--', ':', '-.'};
 
 
             function showMarkers(x0, y0, lx, ly, color, filled)
-                n = length(markers);
+                n = length(Color.markers);
                 s = ceil(sqrt(n));
 
                 x = linspace(x0, x0 + lx, s);
@@ -50,7 +52,7 @@
 
                 for i = 1:n
                     lh = line(X(i),Y(i));
-                    lh.Marker = markers{i};
+                    lh.Marker = Color.markers{i};
                     lh.MarkerSize = 12;
                     lh.Color = color;
 
@@ -79,13 +81,13 @@
             end
 
             function showLines(y0, ly, A, w)
-                n = length(lineStyles);
+                n = length(Color.lineStyles);
                 x = linspace(0,1,100);
                 y = linspace(y0, y0+ ly, n);
                 for i = 1:n
                     lh = line(x, y(i) + A*sin(pi*x*w));
                     lh.LineWidth = 2;
-                    lh.LineStyle = lineStyles{i};
+                    lh.LineStyle = Color.lineStyles{i};
                 end
             end
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LICENSE.txt	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,25 @@
+MIT License
+
+Copyright (c)
+2015-2018 Jonatan Werpers
+2015-2018 Martin Almquist
+2016-2018 Ylva Rydin
+2018 Vidar Stiernström
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README.md	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,2 @@
+# SBPLIB
+sbplib is a library of primitives and help functions for working with summation-by-parts finite differences in Matlab. To use sbplib download the code and add the sbplib folder to the matlab path.

--- a/TextTable.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/TextTable.m	Tue Dec 04 15:24:36 2018 -0800
@@ -4,28 +4,36 @@
         fmtArray
         vertDiv
         horzDiv
-
-        nCols
-        nRows
     end
 
     methods
-        function obj = TextTable(data, vertDiv, horzDiv);
+        function obj = TextTable(data, vertDiv, horzDiv)
             default_arg('vertDiv', []);
             default_arg('horzDiv', []);
 
-
             obj.data = data;
             obj.vertDiv = vertDiv;
             obj.horzDiv = horzDiv;
 
-            [obj.nRows, obj.nCols] = size(data);
             obj.fmtArray = cell(size(data));
             obj.formatAll('%s');
 
         end
 
+        function n = nRows(obj)
+            n = size(obj.data, 1);
+        end
+
+        function n = nCols(obj)
+            n = size(obj.data, 2);
+        end
+
+        function print(obj)
+            disp(obj.toString());
+        end
+
         function formatAll(obj, fmt)
+            obj.fmtArray = cell(size(obj.data));
             obj.fmtArray(:,:) = {fmt};
         end
 
@@ -33,6 +41,14 @@
             obj.fmtArray{i,j} = fmt;
         end
 
+        function formatGrid(obj, I, J, fmt)
+            for i = I
+                for j = J
+                    obj.fmtArray{i,j} = fmt;
+                end
+            end
+        end
+
         function formatRow(obj, i, fmt)
             obj.fmtArray(i,:) = {fmt};
         end
@@ -58,28 +74,31 @@
 
             str = '';
 
+            N = size(strArray, 2);
+
             % First horzDiv
-            if ismember(0, obj.horzDiv)
+            if isDiv(0, obj.horzDiv, N);
                 str = [str, obj.getHorzDiv(widths)];
             end
 
             for i = 1:obj.nRows
-                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv, obj.horzDiv)];
+                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv)];
 
                 % Interior horzDiv
-                if ismember(i, obj.horzDiv)
+                if isDiv(i, obj.horzDiv, N)
                     str = [str, obj.getHorzDiv(widths)];
                 end
             end
         end
 
         function str = getHorzDiv(obj, widths)
-            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv, obj.horzDiv);
+            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv);
             str(find(' ' == str)) = '-';
             str(find('|' == str)) = '+';
         end
 
         function strArray = getStringArray(obj)
+            assert(all(size(obj.data) == size(obj.fmtArray)), 'Sizes of format matrix and data matrix do not match.')
             strArray = cell(size(obj.data));
 
             for i = 1:obj.nRows
@@ -91,32 +110,42 @@
     end
 
     methods (Static)
-        function str = rowToString(strs, widths, vertDiv, horzDiv)
+        function str = rowToString(strs, widths, vertDiv)
+            N = length(strs);
+
             % First vertDiv
-            if ismember(0, vertDiv)
-                str = '| ';
+            if isDiv(0, vertDiv, N)
+                prefix = '| ';
             else
-                str = ' ';
+                prefix = ' ';
             end
 
-            % Interior cols
-            for j = 1:length(strs) - 1
-                str = [str, sprintf('%*s ', widths(j), strs{j})];
+            % Pad strings
+            for i = 1:N
+                strs{i} = sprintf('%*s', widths(i), strs{i});
+            end
 
-                % Interior vertDiv
-                if ismember(j, vertDiv)
-                    str = [str, '| '];
+            % Column delimiters
+            delims = cell(1,N-1);
+            for i = 1:length(delims)
+                if isDiv(i, vertDiv, N);
+                    delims{i} = '| ';
+                else
+                    delims{i} = ' ';
                 end
             end
 
-            % Last col
-            str = [str, sprintf('%*s ', widths(end), strs{end})];
-
-            if ismember(length(strs), vertDiv)
-                str = [str, '|'];
+            if isDiv(N, vertDiv, N);
+                suffix = '|';
+            else
+                suffix = '';
             end
 
-            str = [str, sprintf('\n')];
+            str = [prefix, strjoin(strs, delims), suffix, sprintf('\n')];
         end
     end
+end
+
+function b = isDiv(i, div, N)
+    b = ismember(i, div) || ismember(i, N+div+1);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertIsMember.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,3 @@
+function assertIsMember(v, allowed)
+    assert(ismember(v, allowed), 'Expected ''%s'' to be in the set %s', inputname(1), toString(allowed));
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertSize.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,16 @@
+% Assert that array A has the size s.
+function assertSize(A,varargin)
+    if length(varargin) == 1
+        s = varargin{1};
+        errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), toString(s), toString(size(A)));
+        assert(all(size(A) == s), errmsg);
+    elseif length(varargin) == 2
+        dim = varargin{1};
+        s = varargin{2};
+
+        errmsg = sprintf('Expected %s to have size %d along dimension %d, got: %d',inputname(1), s, dim, size(A,dim));
+        assert(size(A,dim) == s, errmsg);
+    else
+        error('Expected 2 or 3 arguments to assertSize()');
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertStructFields.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,12 @@
+% Assert that the struct s has the all the field names in the cell array fns.
+function assertStructFields(s, fns)
+    assertType(s, 'struct');
+    assertType(fns, 'cell');
+
+    ok = ismember(fns, fieldnames(s));
+    if ~all(ok)
+        str1 = sprintf("'%s' must have the fields %s\n", inputname(1), toString(fns));
+        str2 = sprintf("The following fields are missing: %s", toString(fns(~ok)));
+        error(str1 + str2);
+    end
+end

--- a/assert_size.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/assert_size.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,16 +1,5 @@
 % Assert that array A has the size s.
 function assert_size(A,s)
-    errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), format_vector(s), format_vector(size(A)));
-    assert(all(size(A) == s),errmsg);
-end
-
-function str = format_vector(a)
-    l = length(a);
-    str = sprintf('[%d',a(1));
-
-    for i = 2:l
-        str = [str sprintf(', %d',a(i))];
-    end
-
-    str = [str ']'];
+    warning('Use assertSize() instead!')
+    assertSize(A,s);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/convergencePlot.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,55 @@
+function hand = convergencePlot(orders, h, e)
+    N = length(orders);
+
+    fh = figure();
+    ah = axes();
+    ah.XScale = 'log';
+    ah.YScale = 'log';
+    hold on
+    ph = {};
+    phc = {};
+    legends = {};
+    for i = 1:N
+        ph{i} = loglog(h{i}, e{i});
+        phc{i} = plotConvergenceFit(orders{i}, h{i}, e{i});
+
+        ph{i}.LineStyle = 'none';
+        ph{i}.Marker = Color.solidMarkers{i};
+        ph{i}.MarkerSize = 12;
+        ph{i}.Color = Color.colors{i};
+        ph{i}.MarkerFaceColor = Color.colors{i};
+
+        legends{i} = sprintf('$o = %d$', orders{i});
+    end
+    hold off
+
+    lh = legend([ph{:}], legends);
+    lh.Interpreter = 'latex';
+    lh.Location = 'SouthEast';
+
+    for i = 1:N
+        uistack(phc{i}, 'bottom');
+    end
+
+    xlabel('$h$', 'interpreter', 'latex')
+    ylabel('Error', 'interpreter', 'latex')
+
+    % xlim([0.7e-2, 1e-1])
+    % ylim([3e-5, 4])
+
+    grid on
+
+    ah = gca();
+    ah.TickLabelInterpreter = 'latex';
+    setFontSize(fh);
+
+    % if savePngs
+    %     savepng(fh, 'fig/conv/conv',600)
+    % end
+
+    hand = struct();
+    hand.fig = fh;
+    hand.data = ph;
+    hand.fits = phc;
+    hand.legend = lh;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dealStruct.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,18 @@
+function varargout = dealStruct(s, fields)
+    default_arg('fields', []);
+
+    if isempty(fields)
+        out = dealFields(s, fieldnames(s));
+        varargout = out(1:nargout);
+    else
+        assert(nargout == length(fields), 'Number of output arguements must match the number of fieldnames provided');
+        varargout = dealFields(s, fields);
+    end
+end
+
+function out = dealFields(s, fields)
+    out = cell(1, length(fields));
+    for i = 1:length(fields)
+        out{i} = s.(fields{i});
+    end
+end

--- a/gaussian.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/gaussian.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,3 +1,3 @@
 function z = gaussian(x,x0,d)
-    z = exp(-norm(x-x0).^2/d^2);
+    z = exp(-sum((x-x0).^2,2)/d^2);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/hgRevision.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,8 @@
+% Returns the short mercurial revision Id.
+%  ok is false if there are uncommited changes.
+function [revId, ok] = hgRevision()
+    [~, s] = system('hg id -i');
+    revId = strtrim(s);
+
+    ok = s(end) ~= '+';
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mononomial.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,17 @@
+% calculate a N-D mononomial with powers k in points x:
+%  z = x(:,1).^k(1) * x(:,2).^k(2) * ...
+function z = mononomial(x, k)
+    assert(size(x,2) == length(k), 'k must have the same length as the width of x');
+
+    if any(k < 0)
+        z = x(:,1)*0;
+        return
+    end
+
+    denom = prod(factorial(k));
+
+    for i = 1:length(k)
+        x(:,i) = x(:,i).^k(i);
+    end
+    z = prod(x,2)/denom;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/nextColor.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,5 @@
+function c = nextColor(ah)
+    default_arg('ah', gca);
+
+    c = ah.ColorOrder(ah.ColorOrderIndex, :);
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pointIndex.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,4 @@
+% Get the index of the points p within the tall array of points ps
+function [I, ok] = pointIndex(p, ps)
+    [ok, I] = ismember(p,  ps, 'rows');
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rickerWavelet.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,3 @@
+function y = rickerWavelet(x, x0, A)
+    y = (1-2*pi^2*A^2*(x-x0).^2).*exp(-pi^2*A^2*(x-x0).^2);
+end

--- a/spdiag.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/spdiag.m	Tue Dec 04 15:24:36 2018 -0800
@@ -5,6 +5,6 @@
         a = a';
     end
 
-    n = length(a)-abs(i);
+    n = length(a)+abs(i);
     A = spdiags(a,i,n,n);
 end
\ No newline at end of file

--- a/spdiagVariable.m	Tue Dec 04 14:54:28 2018 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-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/spdiagsPeriodic.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,60 @@
+function A = spdiagsPeriodic(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
+
+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

--- a/spdiagsVariablePeriodic.m	Tue Dec 04 14:54:28 2018 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,42 +0,0 @@
-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/stencilEquation.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,13 @@
+% Find the equation for the stencil for d^k/dx^k
+function [A,b] = stencilEquation(k, offsets, order)
+    q = sym('q', [1, length(offsets)]);
+
+    p = 0:(order-1+k);
+
+    v     = vandermonde(offsets, p);
+    vdiff = vandermonde(      0, p-k);
+
+    eq = q*v == vdiff;
+
+    [A,b] = equationsToMatrix(eq, q);
+end

--- a/stripeMatrixPeriodic.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/stripeMatrixPeriodic.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,8 +1,8 @@
-% Creates a periodic discretization matrix of size n x n 
+% 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);
+    A = spdiagsPeriodic(D,diags);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/structArray.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,16 @@
+% % Usage example:
+% c = structArray({'a','b'}, {
+%     1, 2;
+%     3, 4;
+% });
+
+function c = structArray(fields, values)
+    assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal');
+    c = struct();
+
+    for i = 1:size(values, 1)
+        for j = 1:length(fields)
+            c(i).(fields{j}) = values{i,j};
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/structCellArray.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,16 @@
+% % Usage example:
+% c = structCellArray({'a','b'}, {
+%     1, 2;
+%     3, 4;
+% });
+
+function c = structCellArray(fields, values)
+    assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal');
+    c = cell(1, size(values, 1));
+
+    for i = 1:size(values, 1)
+        for j = 1:length(fields)
+            c{i}.(fields{j}) = values{i,j};
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stuffStruct.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,8 @@
+function s = stuffStruct(varargin)
+    s = struct();
+
+    for i = 1:nargin
+        assert(~isempty(inputname(i)), 'All inputs must be variables.');
+        s.(inputname(i)) = varargin{i};
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vandermonde.m	Tue Dec 04 15:24:36 2018 -0800
@@ -0,0 +1,15 @@
+% Create vandermonde matrix for points x and polynomials of order p
+% x is a list of N points of size [N,dim],
+% p is a list of polynomial orders of size [M, dim].
+% the given mononomials are evaluated and the NxM matrix V is returned.
+function V = vandermonde(x, p)
+    assert(size(x,2) == size(p,2), 'x and p must have the same number of columns')
+    n = size(x,1);
+    m = size(p,1);
+
+    for i = 1:m
+        V(:,i) = mononomial(x, p(i,:));
+    end
+
+    assertSize(V,[n,m]);
+end