--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/joinGrids.m	Wed Nov 20 14:26:57 2019 -0800
@@ -0,0 +1,90 @@
+% Connects several multiblock grids into one grid
+% gs 		   -- 	Cell array of multiblock grids
+% gConnections -- 	Upper-triangular cell matrix
+% 					gConnections{i,j} specifies all connections between grid i and grid j
+%					Example:
+%					gConnections{i,j} = { {{1,'e'},{2,'w'}}, {{5,'s'},{2,'n'}},... };
+% names 	   --	(Optional) cell array of strings, used for boundary groups in new grid
+%					default: names = {'g1', 'g2', ..., 'gN'};
+% 					Boundary groups from grid i are contained in g.boundaryGroups.(names{i}), etc.
+function g = joinGrids(gs, gConnections, names)
+
+	nGrids = numel(gs);
+
+	% Default names are {'g1', 'g2', ... 'gN'}.
+	defaultNames = cell(nGrids, 1);
+	for i = 1:nGrids
+		defaultNames{i} = sprintf('g%d',i);
+	end
+	default_arg('names', defaultNames);
+
+	nBlocks = 0;
+	for i = 1:nGrids
+		nBlocks = nBlocks + gs{i}.nBlocks();
+	end
+
+	% Create vector of cumulative sum of number of blocks per grid
+	startIndex = zeros(1, nGrids);
+	for i = 2:nGrids
+		startIndex(i) = startIndex(i-1) + gs{i-1}.nBlocks();
+	end
+
+	% Create cell array of all grids
+	grids = cell(nBlocks, 1);
+	for i = 1:nGrids
+		for j = 1:gs{i}.nBlocks();
+			grids{startIndex(i)+j} = gs{i}.grids{j};
+		end
+	end
+
+	% Create cell matrix of connections
+	connections = cell(nBlocks, nBlocks);
+
+	% Connections within grids
+	for i = 1:nGrids
+		for j = 1:gs{i}.nBlocks()
+			for k = 1:gs{i}.nBlocks()
+				connections{startIndex(i)+j,startIndex(i)+k} = gs{i}.connections{j,k};
+			end
+		end
+	end
+
+	% Connections between grids
+	for i = 1:nGrids
+		for j = 1:nGrids
+			for k = 1:numel(gConnections{i,j})
+				b1 = gConnections{i,j}{k}{1};
+				id1 = b1{1};
+				str1 = b1{2};
+
+				b2 = gConnections{i,j}{k}{2};
+				id2 = b2{1};
+				str2 = b2{2};
+
+				connections{startIndex(i)+id1, startIndex(j)+id2} = {str1, str2};
+			end
+		end
+	end
+
+	% Boundary groups
+	boundaryGroups = struct;
+	for i = 1:nGrids
+		bgs = gs{i}.boundaryGroups;
+		bgNames = fieldnames(bgs);
+		for j = 1:numel(bgNames)
+			bg = bgs.(bgNames{j});
+
+			% Shift block id:s in boundary groups
+			for k = 1:length(bg)
+				bg{k}{1} = bg{k}{1} + startIndex(i);
+			end
+
+			bgs.(bgNames{j}) = bg;
+		end
+		boundaryGroups.(names{i}) = bgs;
+	end
+
+	% Create grid object
+	g = multiblock.Grid(grids, connections, boundaryGroups);
+
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D2VariableCompatible.m	Wed Nov 20 14:26:57 2019 -0800
@@ -0,0 +1,83 @@
+classdef D2VariableCompatible < sbp.OpSet
+    properties
+        D1 % SBP operator approximating first derivative
+        H % Norm matrix
+        HI % H^-1
+        Q % Skew-symmetric matrix
+        e_l % Left boundary operator
+        e_r % Right boundary operator
+        D2 % SBP operator for second derivative
+        M % Norm matrix, second derivative
+        d1_l % Left boundary first derivative
+        d1_r % Right boundary first derivative
+        m % Number of grid points.
+        h % Step size
+        x % grid
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D2VariableCompatible(m,lim,order)
+
+            x_l = lim{1};
+            x_r = lim{2};
+            L = x_r-x_l;
+            obj.h = L/(m-1);
+            obj.x = linspace(x_l,x_r,m)';
+
+            switch order
+
+                case 6
+
+                    [obj.H, obj.HI, obj.D1, D2, ...
+                    ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,...
+                     d1_l, d1_r] = ...
+                        sbp.implementations.d4_variable_6(m, obj.h);
+
+                case 4
+                    [obj.H, obj.HI, obj.D1, D2, obj.e_l,...
+                        obj.e_r, d1_l, d1_r] = ...
+                        sbp.implementations.d2_variable_4(m,obj.h);
+                case 2
+                    [obj.H, obj.HI, obj.D1, D2, obj.e_l,...
+                        obj.e_r, d1_l, d1_r] = ...
+                        sbp.implementations.d2_variable_2(m,obj.h);
+
+                otherwise
+                    error('Invalid operator order %d.',order);
+            end
+            obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h,
+            obj.borrowing.M.d1 = obj.H(1,1)/obj.h; % First element in H/h is borrowing also for M
+            obj.borrowing.R.delta_D = inf;
+            obj.m = m;
+            obj.M = [];
+
+
+            D1 = obj.D1;
+            e_r = obj.e_r;
+            e_l = obj.e_l;
+
+            % D2 = Hinv * (-M + br*er*d1r^T - bl*el*d1l^T);
+            % Replace d1' by e'*D1 in D2.
+            correction_l = obj.HI*(e_l*d1_l' - e_l*e_l'*D1);
+            correction_r = - obj.HI*(e_r*d1_r' - e_r*e_r'*D1);
+
+            D2_compatible = @(b) D2(b) + b(1)*correction_l + b(m)*correction_r;
+
+            obj.D2 = D2_compatible;
+            obj.d1_l = (e_l'*D1)';
+            obj.d1_r = (e_r'*D1)';
+
+        end
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+
+
+end
+
+
+
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscr.m	Wed Nov 20 14:26:57 2019 -0800
@@ -0,0 +1,106 @@
+% n-dimensional delta function
+% g: cartesian grid
+% x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
+% m_order: Number of moment conditions
+% s_order: Number of smoothness conditions
+% H: cell array of 1D norm matrices
+function d = diracDiscr(g, x_s, m_order, s_order, H)
+    assertType(g, 'grid.Cartesian');
+    
+    dim = g.d;
+    d_1D = cell(dim,1);
+
+    % Allow for non-cell input in 1D
+    if dim == 1
+        H = {H};
+    end
+    % Create 1D dirac discr for each coordinate dir.
+    for i = 1:dim
+        d_1D{i} = diracDiscr1D(x_s(i), g.x{i}, m_order, s_order, H{i});
+    end
+
+    d = d_1D{dim};
+    for i = dim-1: -1: 1
+        % Perform outer product, transpose, and then turn into column vector
+        d = (d_1D{i}*d')';
+        d = d(:);
+    end
+end
+
+
+% Helper function for 1D delta functions
+function ret = diracDiscr1D(x_s, x, m_order, s_order, H)
+
+    % Return zeros if x0 is outside grid
+    if x_s < x(1) || x_s > x(end)
+        ret = zeros(size(x));
+        return
+    end
+       
+    tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
+    S = [];
+    M = [];
+
+    % Get interior grid spacing
+    middle = floor(length(x)/2);
+    h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids.
+
+    index = sourceIndices(x_s, x, tot_order, h);
+
+    polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
+    x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
+
+    quadrature = diag(H);
+    quadrature_weights = quadrature(index)/h;
+
+    h_polynomial = polynomial(2)-polynomial(1);
+    b = zeros(tot_order,1);
+
+    for i = 1:m_order
+        b(i,1) = x_0^(i-1);
+    end
+
+    for i = 1:tot_order
+        for j = 1:m_order
+            M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
+        end
+    end
+
+    for i = 1:tot_order
+        for j = 1:s_order
+            S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
+        end
+    end
+
+    A = [M;S];
+
+    d = A\b;
+    ret = x*0;
+    ret(index) = d/h*h_polynomial;
+end
+
+
+function I = sourceIndices(x_s, x, tot_order, h)
+    % Find the indices that are within range of of the point source location
+    I = find(tot_order*h/2 >= abs(x-x_s));
+
+    if length(I) > tot_order
+        if length(I) == tot_order + 2
+            I = I(2:end-1);
+        elseif length(I) == tot_order + 1
+            I = I(1:end-1);
+        end
+    elseif length(I) < tot_order
+        if x_s < x(1) + ceil(tot_order/2)*h
+            I = 1:tot_order;
+        elseif x_s > x(end) - ceil(tot_order/2)*h
+            I = length(x)-tot_order+1:length(x);
+        else
+            if I(end) < length(x)
+                I = [I; I(end)+1];
+            else
+                I = [I(1)-1; I];
+            end
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscrCurve.m	Wed Nov 20 14:26:57 2019 -0800
@@ -0,0 +1,63 @@
+% 2-dimensional delta function for single-block curvilinear grid
+% x_s:      source point coordinate vector, e.g. [x; y] or [x; y; z].
+% g:        single-block grid containing the source
+% m_order:  Number of moment conditions
+% s_order:  Number of smoothness conditions
+% order:    Order of SBP derivative approximations
+% opSet:    Cell array of function handle to opSet generator
+function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet)
+    default_arg('order', m_order);
+    default_arg('opSet', {@sbp.D2Variable, @sbp.D2Variable});
+
+    dim = length(x_s);
+    assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.');
+    assertType(g, 'grid.Curvilinear');
+
+    % Compute Jacobian
+    J = jacobian(g, opSet, order);
+
+    % Find approximate logical coordinates of point source
+    X = g.points();
+    U = g.logic.points();
+    U_interp = scatteredInterpolant(X, U(:,1));
+    V_interp = scatteredInterpolant(X, U(:,2));
+    uS = U_interp(x_s);
+    vS = V_interp(x_s);
+
+    % Get quadrature matrices for moment conditions
+    m = g.size();
+    ops_u = opSet{1}(m(1), {0, 1}, order);
+    ops_v = opSet{2}(m(2), {0, 1}, order);
+    H_u =  ops_u.H;
+    H_v =  ops_v.H;
+
+    % Get delta function for logical grid and scale by Jacobian
+    d = (1./J).*diracDiscr(g, [uS; vS], m_order, s_order, {H_u, H_v});
+end
+
+function J = jacobian(g, opSet, order)
+    m = g.size();
+    m_u = m(1);
+    m_v = m(2);
+    ops_u = opSet{1}(m_u, {0, 1}, order);
+    ops_v = opSet{2}(m_v, {0, 1}, order);
+    I_u = speye(m_u);
+    I_v = speye(m_v);
+
+    D1_u = ops_u.D1;
+    D1_v = ops_v.D1;
+
+    Du = kr(D1_u,I_v);
+    Dv = kr(I_u,D1_v);
+
+    coords = g.points();
+    x = coords(:,1);
+    y = coords(:,2);
+
+    x_u = Du*x;
+    x_v = Dv*x;
+    y_u = Du*y;
+    y_v = Dv*y;
+
+    J = x_u.*y_v - x_v.*y_u;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscrTest.m	Wed Nov 20 14:26:57 2019 -0800
@@ -0,0 +1,415 @@
+function tests = diracDiscrTest()
+	    tests = functiontests(localfunctions);
+end
+
+%TODO:
+%   1.  Test discretizing with smoothness conditions.
+%       Only discretization with moment conditions currently tested.
+%   2.  Test using other types of grids. Only equidistant grids currently used.
+
+function testLeftRandom(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    rng(1) % Set seed
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup1D(order, mom_cond);
+        xl = g.lim{1}{1};
+        h = g.scaling();
+        x = g.x{1};
+
+        % Test random points near left boundary
+        x0s = xl + 2*h*rand(1,10);
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testRightRandom(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    rng(1) % Set seed
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup1D(order, mom_cond);
+        xr = g.lim{1}{2};
+        h = g.scaling();
+        x = g.x{1};
+
+        % Test random points near right boundary
+        x0s = xr - 2*h*rand(1,10);
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testInteriorRandom(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    rng(1) % Set seed
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup1D(order, mom_cond);
+        xl = g.lim{1}{1};
+        xr = g.lim{1}{2};
+        h = g.scaling();
+        x = g.x{1};
+
+        % Test random points in interior
+        x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+% x0 outside grid should yield 0 integral!
+function testX0OutsideGrid(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup1D(order, mom_cond);
+        xl = g.lim{1}{1};
+        xr = g.lim{1}{2};
+        h = g.scaling();
+        x = g.x{1};
+
+        % Test points outisde grid
+        x0s = [xl-1.1*h, xr+1.1*h];
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - 0);
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testAllGP(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+       [g, H, fs] = setup1D(order, mom_cond);
+        x = g.x{1};
+
+        % Test all grid points
+        x0s = x;
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testHalfGP(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup1D(order, mom_cond);
+        x = g.x{1};
+
+        % Test halfway between all grid points
+        x0s = 1/2*(x(2:end)+x(1:end-1));
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+%=============== 2D tests ==============================
+function testAllGP2D(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
+        X = g.points();
+        % Test all grid points
+        x0s = X;
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(X(:,1), X(:,2));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H_global*fx;
+                err = abs(integral - f(x0(1), x0(2)));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testAllRandom2D(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    rng(1) % Set seed
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
+        X = g.points();
+        xl = g.lim{1}{1};
+        xr = g.lim{1}{2};
+        yl = g.lim{2}{1};
+        yr = g.lim{2}{2};
+        h = g.scaling();
+
+
+        % Test random points, even outside grid
+        Npoints = 100;
+        x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ...
+               (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1) ];
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(X(:,1), X(:,2));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H_global*fx;
+
+                % Integral should be 0 if point is outside grid
+                if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
+                    err = abs(integral - 0);
+                else
+                    err = abs(integral - f(x0(1), x0(2)));
+                end
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+%=============== 3D tests ==============================
+function testAllGP3D(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup3D(order, mom_cond);
+        H_global = kron(kron(H{1}, H{2}), H{3});
+        X = g.points();
+        % Test all grid points
+        x0s = X;
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(X(:,1), X(:,2), X(:,3));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H_global*fx;
+                err = abs(integral - f(x0(1), x0(2), x0(3)));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testAllRandom3D(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    rng(1) % Set seed
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [g, H, fs] = setup3D(order, mom_cond);
+        H_global = kron(kron(H{1}, H{2}), H{3});
+        X = g.points();
+        xl = g.lim{1}{1};
+        xr = g.lim{1}{2};
+        yl = g.lim{2}{1};
+        yr = g.lim{2}{2};
+        zl = g.lim{3}{1};
+        zr = g.lim{3}{2};
+        h = g.scaling();
+
+        % Test random points, even outside grid
+        Npoints = 200;
+        x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ...
+               (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1), ...
+               (zl-3*h(3)) + (zr-zl+6*h(3))*rand(Npoints,1) ];
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(X(:,1), X(:,2), X(:,3));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
+                integral = delta'*H_global*fx;
+
+                % Integral should be 0 if point is outside grid
+                if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr
+                    err = abs(integral - 0);
+                else
+                    err = abs(integral - f(x0(1), x0(2), x0(3)));
+                end
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+
+% ======================================================
+% ============== Setup functions =======================
+% ======================================================
+function [g, H, fs] = setup1D(order, mom_cond)
+
+    % Grid
+    xl = -3;
+    xr = 900;
+    L = xr-xl;
+    m = 101;
+    g = grid.equidistant(m, {xl, xr});
+
+    % Quadrature
+    ops = sbp.D2Standard(m, {xl, xr}, order);
+    H = ops.H;
+
+    % Moment conditions
+    fs = cell(mom_cond,1);
+    for p = 0:mom_cond-1
+        fs{p+1} = @(x) (x/L).^p;
+    end
+
+end
+
+function [g, H, fs] = setup2D(order, mom_cond)
+
+    % Grid
+    xlims = {-3, 20};
+    ylims = {-11,5};
+    Lx = xlims{2} - xlims{1};
+    Ly = ylims{2} - ylims{1};
+    m = [15, 16];
+    g = grid.equidistant(m, xlims, ylims);
+
+    % Quadrature
+    opsx = sbp.D2Standard(m(1), xlims, order);
+    opsy = sbp.D2Standard(m(2), ylims, order);
+    Hx = opsx.H;
+    Hy = opsy.H;
+    H = {Hx, Hy};
+
+    % Moment conditions
+    fs = cell(mom_cond,1);
+    for p = 0:mom_cond-1
+        fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
+    end
+end
+
+function [g, H, fs] = setup3D(order, mom_cond)
+
+    % Grid
+    xlims = {-3, 20};
+    ylims = {-11,5};
+    zlims = {2,4};
+    Lx = xlims{2} - xlims{1};
+    Ly = ylims{2} - ylims{1};
+    Lz = zlims{2} - zlims{1};
+    m = [13, 14, 15];
+    g = grid.equidistant(m, xlims, ylims, zlims);
+
+    % Quadrature
+    opsx = sbp.D2Standard(m(1), xlims, order);
+    opsy = sbp.D2Standard(m(2), ylims, order);
+    opsz = sbp.D2Standard(m(3), zlims, order);
+    Hx = opsx.H;
+    Hy = opsy.H;
+    Hz = opsz.H;
+    H = {Hx, Hy, Hz};
+
+    % Moment conditions
+    fs = cell(mom_cond,1);
+    for p = 0:mom_cond-1
+        fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p;
+    end
+end