--- a/diracDiscrTest.m	Sun Nov 04 12:36:30 2018 -0800
+++ b/diracDiscrTest.m	Wed Nov 21 18:29:29 2018 -0800
@@ -10,7 +10,7 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test left boundary grid points
         x0s = xl + [0, h, 2*h];
@@ -33,11 +33,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points near left boundary
         x0s = xl + 2*h*rand(1,10);
@@ -60,11 +60,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test right boundary grid points
         x0s = xr-[0, h, 2*h];
@@ -87,11 +87,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points near right boundary
         x0s = xr - 2*h*rand(1,10);
@@ -114,11 +114,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test interior grid points
         m_half = round(m/2);
@@ -142,11 +142,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points in interior
         x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
@@ -170,11 +170,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test points outisde grid
         x0s = [xl-1.1*h, xr+1.1*h];
@@ -197,11 +197,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test all grid points
         x0s = x;
@@ -224,11 +224,11 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test halfway between all grid points
         x0s = 1/2*( x(2:end)+x(1:end-1) );
@@ -247,81 +247,224 @@
     end
 end
 
-function testAllGPStaggered(testCase)
+% function testAllGPStaggered(testCase)
+
+%     orders = [2, 4, 6];
+%     mom_conds = orders;
+
+%     for o = 1:length(orders)
+%         order = orders(o);
+%         mom_cond = mom_conds(o);
+%         [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+%         % 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(x0, x, mom_cond, 0, H);
+%                 integral = delta'*H*fx;
+%                 err = abs(integral - f(x0));
+%                 testCase.verifyLessThan(err, 1e-12);
+%             end
+%         end
+%     end
+% end
+
+% function testHalfGPStaggered(testCase)
+
+%     orders = [2, 4, 6];
+%     mom_conds = orders;
+
+%     for o = 1:length(orders)
+%         order = orders(o);
+%         mom_cond = mom_conds(o);
+%         [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+%         % 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(x0, x, mom_cond, 0, H);
+%                 integral = delta'*H*fx;
+%                 err = abs(integral - f(x0));
+%                 testCase.verifyLessThan(err, 1e-12);
+%             end
+%         end
+%     end
+% end
+
+% function testRandomStaggered(testCase)
+
+%     orders = [2, 4, 6];
+%     mom_conds = orders;
+
+%     for o = 1:length(orders)
+%         order = orders(o);
+%         mom_cond = mom_conds(o);
+%         [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+%         % Test random points within grid boundaries
+%         x0s = xl + (xr-xl)*rand(1,300);
+
+%         for j = 1:length(fs)
+%                 f = fs{j};
+%                 fx = f(x);
+%             for i = 1:length(x0s)
+%                 x0 = x0s(i);
+%                 delta = diracDiscr(x0, x, 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);
-        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+        [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
 
         % Test all grid points
-        x0s = x;
+        x0s = X;
 
         for j = 1:length(fs)
                 f = fs{j};
-                fx = f(x);
+                fx = f(X(:,1), X(:,2));
             for i = 1:length(x0s)
-                x0 = x0s(i);
+                x0 = x0s(i,:);
                 delta = diracDiscr(x0, x, mom_cond, 0, H);
-                integral = delta'*H*fx;
-                err = abs(integral - f(x0));
+                integral = delta'*H_global*fx;
+                err = abs(integral - f(x0(1), x0(2)));
                 testCase.verifyLessThan(err, 1e-12);
             end
         end
     end
 end
 
-function testHalfGPStaggered(testCase)
+function testAllRandom2D(testCase)
 
     orders = [2, 4, 6];
     mom_conds = orders;
-    
+
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+        [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
 
-        % Test halfway between all grid points
-        x0s = 1/2*( x(2:end)+x(1:end-1) );
+        xl = xlims{1};
+        xr = xlims{2};
+        yl = ylims{1};
+        yr = ylims{2};
+
+        % 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);
+                fx = f(X(:,1), X(:,2));
             for i = 1:length(x0s)
-                x0 = x0s(i);
+                x0 = x0s(i,:);
                 delta = diracDiscr(x0, x, mom_cond, 0, H);
-                integral = delta'*H*fx;
-                err = abs(integral - f(x0));
+                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
 
-function testRandomStaggered(testCase)
+%=============== 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);
-        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+        [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
+        H_global = kron(kron(H{1}, H{2}), H{3});
 
-        % Test random points within grid boundaries
-        x0s = xl + (xr-xl)*rand(1,300);
+        % Test all grid points
+        x0s = X;
 
         for j = 1:length(fs)
                 f = fs{j};
-                fx = f(x);
+                fx = f(X(:,1), X(:,2), X(:,3));
             for i = 1:length(x0s)
-                x0 = x0s(i);
+                x0 = x0s(i,:);
                 delta = diracDiscr(x0, x, mom_cond, 0, H);
-                integral = delta'*H*fx;
-                err = abs(integral - f(x0));
+                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;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
+        H_global = kron(kron(H{1}, H{2}), H{3});
+
+        xl = xlims{1};
+        xr = xlims{2};
+        yl = ylims{1};
+        yr = ylims{2};
+        zl = zlims{1};
+        zr = zlims{2};
+
+        % 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(x0, x, 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
@@ -329,8 +472,10 @@
 end
 
 
+% ======================================================
 % ============== Setup functions =======================
-function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
+% ======================================================
+function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond)
 
     % Grid
     xl = -3;
@@ -353,6 +498,79 @@
 
 end
 
+function [xlims, ylims, m, x, X, h, 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);
+    X = g.points();
+    x = g.x;
+
+    % 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
+
+    % Grid spacing in interior
+    mm = round(m/2);
+    hx = x{1}(mm(1)+1) - x{1}(mm(1));
+    hy = x{2}(mm(2)+1) - x{2}(mm(2));
+    h = {hx, hy};
+
+end
+
+function [xlims, ylims, zlims, m, x, X, h, 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);
+    X = g.points();
+    x = g.x;
+
+    % 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
+
+    % Grid spacing in interior
+    mm = round(m/2);
+    hx = x{1}(mm(1)+1) - x{1}(mm(1));
+    hy = x{2}(mm(2)+1) - x{2}(mm(2));
+    hz = x{3}(mm(3)+1) - x{3}(mm(3));
+    h = {hx, hy, hz};
+
+end
+
 function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond)
 
     % Grid
@@ -375,5 +593,3 @@
     end
 
 end
-
-