--- a/diracPrimDiscrTest.m	Tue May 21 18:44:01 2019 -0700
+++ b/diracPrimDiscrTest.m	Wed May 22 15:31:26 2019 -0700
@@ -258,74 +258,88 @@
 
 
 %=============== 2D tests ==============================
-% function testAllGP2D(testCase)
+function testAllGP2D(testCase)
 
-%     orders = [2, 4, 6];
-%     mom_conds = orders;
+    orders = [2, 4, 6];
+    mom_conds = orders;
 
-%     for o = 1:length(orders)
-%         order = orders(o);
-%         mom_cond = mom_conds(o);
-%         [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond);
-%         H_global = kron(H{1}, H{2});
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xlims, ylims, m, x, X, ~, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
 
-%         % Test all grid points
-%         x0s = X;
+        % 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 = diracPrimDiscr(x0, x, 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
+        for j = 1:length(fs)
+                f = fs{j};
+                dfdx = dfdxs{j};
+                dfdy = dfdys{j};
+                fx = f(X(:,1), X(:,2));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+
+                delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1);
+                integral1 = delta_x'*H_global*fx;
 
-% function testAllRandom2D(testCase)
+                delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2);
+                integral2 = delta_y'*H_global*fx;
+
+                err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2)));
 
-%     orders = [2, 4, 6];
-%     mom_conds = orders;
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testAllRandom2D(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
 
-%     for o = 1:length(orders)
-%         order = orders(o);
-%         mom_cond = mom_conds(o);
-%         [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond);
-%         H_global = kron(H{1}, H{2});
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xlims, ylims, m, x, X, h, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond);
+        H_global = kron(H{1}, H{2});
 
-%         xl = xlims{1};
-%         xr = xlims{2};
-%         yl = ylims{1};
-%         yr = ylims{2};
+        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) ];
+        % 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 = diracPrimDiscr(x0, x, mom_cond, 0, H);
-%                 integral = delta'*H_global*fx;
+        for j = 1:length(fs)
+                f = fs{j};
+                dfdx = dfdxs{j};
+                dfdy = dfdys{j};
+                fx = f(X(:,1), X(:,2));
+            for i = 1:length(x0s)
+                x0 = x0s(i,:);
+
+                delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1);
+                integral1 = delta_x'*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
+                delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2);
+                integral2 = delta_y'*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(integral1 - 0) + abs(integral2 - 0);
+                else
+                    err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2)));
+                end
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
 
 %=============== 3D tests ==============================
 % function testAllGP3D(testCase)
@@ -434,39 +448,48 @@
 
 end
 
-% function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond)
+function [xlims, ylims, m, x, X, h, H, fs, fxs, fys] = setup2D(order, mom_cond)
+
+    % Grid
+    xlims = {-3, 20};
+    ylims = {-11,5};
+    Lx = xlims{2} - xlims{1};
+    Ly = ylims{2} - ylims{1};
 
-%     % 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;
 
-%     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};
 
-%     % 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,1);
+    fxs = cell(mom_cond+1,1);
+    fys = cell(mom_cond+1,1);
+    for p = 0:mom_cond
+        fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
+    end
 
-%     % Moment conditions
-%     fs = cell(mom_cond,1);
-%     for p = 0:mom_cond-1
-%         fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
-%     end
+    fxs{1} = @(x,y) 0*x;
+    fys{1} = @(x,y) 0*y;
+    for p = 1:mom_cond
+        fxs{p+1} = @(x,y) p/Lx*(x/Lx + y/Ly).^(p-1);
+        fys{p+1} = @(x,y) p/Ly*(x/Lx + y/Ly).^(p-1);
+    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};
+    % 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
+end
 
 % function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)