--- a/diracDiscrTest.m	Tue Nov 19 13:55:43 2019 -0800
+++ b/diracDiscrTest.m	Tue Nov 19 15:16:08 2019 -0800
@@ -2,8 +2,10 @@
 	    tests = functiontests(localfunctions);
 end
 
-%TODO: Test discretizing with smoothness conditions.
-%      Only discretization with moment conditions currently tested.
+%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)
 
@@ -14,7 +16,10 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, ~, h, x, H, fs] = setup1D(order, mom_cond);
+        [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);
@@ -24,7 +29,7 @@
                 fx = f(x);
             for i = 1:length(x0s)
                 x0 = x0s(i);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - f(x0));
                 testCase.verifyLessThan(err, 1e-12);
@@ -42,7 +47,10 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [~, xr, h, x, H, fs] = setup1D(order, mom_cond);
+        [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);
@@ -52,7 +60,7 @@
                 fx = f(x);
             for i = 1:length(x0s)
                 x0 = x0s(i);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - f(x0));
                 testCase.verifyLessThan(err, 1e-12);
@@ -70,7 +78,11 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
+        [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);
@@ -80,7 +92,7 @@
                 fx = f(x);
             for i = 1:length(x0s)
                 x0 = x0s(i);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - f(x0));
                 testCase.verifyLessThan(err, 1e-12);
@@ -98,7 +110,11 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
+        [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];
@@ -108,7 +124,7 @@
                 fx = f(x);
             for i = 1:length(x0s)
                 x0 = x0s(i);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - 0);
                 testCase.verifyLessThan(err, 1e-12);
@@ -125,7 +141,8 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
+       [g, H, fs] = setup1D(order, mom_cond);
+        x = g.x{1};
 
         % Test all grid points
         x0s = x;
@@ -135,7 +152,7 @@
                 fx = f(x);
             for i = 1:length(x0s)
                 x0 = x0s(i);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - f(x0));
                 testCase.verifyLessThan(err, 1e-12);
@@ -152,17 +169,18 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
+        [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) );
+        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);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H*fx;
                 err = abs(integral - f(x0));
                 testCase.verifyLessThan(err, 1e-12);
@@ -180,9 +198,9 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [~, ~, x, X, ~, H, fs] = setup2D(order, mom_cond);
+        [g, H, fs] = setup2D(order, mom_cond);
         H_global = kron(H{1}, H{2});
-
+        X = g.points();
         % Test all grid points
         x0s = X;
 
@@ -191,7 +209,7 @@
                 fx = f(X(:,1), X(:,2));
             for i = 1:length(x0s)
                 x0 = x0s(i,:);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                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);
@@ -209,25 +227,27 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond);
+        [g, H, fs] = setup2D(order, mom_cond);
         H_global = kron(H{1}, H{2});
-
-        xl = xlims{1};
-        xr = xlims{2};
-        yl = ylims{1};
-        yr = ylims{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) ];
+        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(x0, x, mom_cond, 0, H);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H_global*fx;
 
                 % Integral should be 0 if point is outside grid
@@ -251,9 +271,9 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [~, ~, ~, x, X, ~, H, fs] = setup3D(order, mom_cond);
+        [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;
 
@@ -262,7 +282,7 @@
                 fx = f(X(:,1), X(:,2), X(:,3));
             for i = 1:length(x0s)
                 x0 = x0s(i,:);
-                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                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);
@@ -280,28 +300,29 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xlims, ylims, zlims, x, X, h, H, fs] = setup3D(order, mom_cond);
+        [g, 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};
+        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) ];
+        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);
+                delta = diracDiscr(g, x0, mom_cond, 0, H);
                 integral = delta'*H_global*fx;
 
                 % Integral should be 0 if point is outside grid
@@ -320,16 +341,14 @@
 % ======================================================
 % ============== Setup functions =======================
 % ======================================================
-function [xl, xr, h, x, H, fs] = setup1D(order, mom_cond)
+function [g, H, fs] = setup1D(order, mom_cond)
 
     % Grid
     xl = -3;
     xr = 900;
     L = xr-xl;
     m = 101;
-    h = (xr-xl)/(m-1);
     g = grid.equidistant(m, {xl, xr});
-    x = g.points();
 
     % Quadrature
     ops = sbp.D2Standard(m, {xl, xr}, order);
@@ -343,18 +362,15 @@
 
 end
 
-function [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond)
+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);
-    X = g.points();
-    x = g.x;
 
     % Quadrature
     opsx = sbp.D2Standard(m(1), xlims, order);
@@ -368,16 +384,9 @@
     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, x, X, h, H, fs] = setup3D(order, mom_cond)
+function [g, H, fs] = setup3D(order, mom_cond)
 
     % Grid
     xlims = {-3, 20};
@@ -386,11 +395,8 @@
     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);
@@ -406,12 +412,4 @@
     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
\ No newline at end of file