Mercurial > repos > public > sbplib

changeset 1233:57df0bf741dc feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in curvilinear discr of dirac delta
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 19 Nov 2019 13:55:43 -0800
parents 52d774e69b1f (diff) 8a456f6e54cc (current diff)
children f1806475498b 48c9a83260c8
files
diffstat 2 files changed, 89 insertions(+), 267 deletions(-) [+]
line wrap: on
line diff
--- a/diracDiscr.m	Tue Nov 19 13:47:17 2019 -0800
+++ b/diracDiscr.m	Tue Nov 19 13:55:43 2019 -0800
@@ -31,95 +31,95 @@
 
 
 % Helper function for 1D delta functions
-function ret = diracDiscr1D(x_0in , x , m_order, s_order, H)
+function ret = diracDiscr1D(x_s , x , m_order, s_order, H)
 
-m = length(x);
+    m = length(x);
 
-% Return zeros if x0 is outside grid
-if(x_0in < x(1) || x_0in > x(end) )
+    % Return zeros if x0 is outside grid
+    if(x_s < x(1) || x_s > x(end) )
 
-    ret = zeros(size(x));
+        ret = zeros(size(x));
 
-else
+    else
 
-    fnorm = diag(H);
-    eta = abs(x-x_0in);
-    tot = m_order+s_order;
-    S = [];
-    M = [];
+        fnorm = diag(H);
+        tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
+        S = [];
+        M = [];
 
-    % Get interior grid spacing
-    middle = floor(m/2);
-    h = x(middle+1) - x(middle);
+        % Get interior grid spacing
+        middle = floor(m/2);
+        h = x(middle+1) - x(middle);
 
-    poss = find(tot*h/2 >= eta);
+        % Find the indices that are within range of of the point source location
+        ind_delta = find(tot_order*h/2 >= abs(x-x_s));
 
-    % Ensure that poss is not too long
-    if length(poss) == (tot + 2)
-        poss = poss(2:end-1);
-    elseif length(poss) == (tot + 1)
-        poss = poss(1:end-1);
-    end
+        % Ensure that ind_delta is not too long
+        if length(ind_delta) == (tot_order + 2)
+            ind_delta = ind_delta(2:end-1);
+        elseif length(ind_delta) == (tot_order + 1)
+            ind_delta = ind_delta(1:end-1);
+        end
 
-    % Use first tot grid points
-    if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
-        index=1:tot;
-        pol=(x(1:tot)-x(1))/(x(tot)-x(1));
-        x_0=(x_0in-x(1))/(x(tot)-x(1));
-        norm=fnorm(1:tot)/h;
+        % Use first tot_order grid points
+        if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h;
+            index=1:tot_order;
+            polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1));
+            x_0=(x_s-x(1))/(x(tot_order)-x(1));
+            norm=fnorm(1:tot_order)/h;
 
-    % Use last tot grid points
-    elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
-        index = length(x)-tot+1:length(x);
-        pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
-        norm = fnorm(end-tot+1:end)/h;
-        x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
+        % Use last tot_order grid points
+        elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h;
+            index = length(x)-tot_order+1:length(x);
+            polynomial = (x(end-tot_order+1:end)-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
+            norm = fnorm(end-tot_order+1:end)/h;
+            x_0 = (x_s-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
 
-    % Interior, compensate for round-off errors.
-    elseif length(poss) < tot
-        if poss(end)<m
-            poss = [poss; poss(end)+1];
+        % Interior, compensate for round-off errors.
+        elseif length(ind_delta) < tot_order
+            if ind_delta(end)<m
+                ind_delta = [ind_delta; ind_delta(end)+1];
+            else
+                ind_delta = [ind_delta(1)-1; ind_delta];
+            end
+            polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
+            x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
+            norm = fnorm(ind_delta)/h;
+            index = ind_delta;
+
+        % Interior
         else
-            poss = [poss(1)-1; poss];
+            polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
+            x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
+            norm = fnorm(ind_delta)/h;
+            index = ind_delta;
         end
-        pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
-        x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
-        norm = fnorm(poss)/h;
-        index = poss;
 
-    % Interior
-    else
-        pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
-        x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
-        norm = fnorm(poss)/h;
-        index = poss;
-    end
-
-    h_pol = pol(2)-pol(1);
-    b = zeros(m_order+s_order,1);
+        h_polynomial = polynomial(2)-polynomial(1);
+        b = zeros(m_order+s_order,1);
 
-    for i = 1:m_order
-        b(i,1) = x_0^(i-1);
-    end
+        for i = 1:m_order
+            b(i,1) = x_0^(i-1);
+        end
 
-    for i = 1:(m_order+s_order)
-        for j = 1:m_order
-            M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
+        for i = 1:(m_order+s_order)
+            for j = 1:m_order
+                M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i);
+            end
         end
-    end
 
-    for i = 1:(m_order+s_order)
-        for j = 1:s_order
-            S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
+        for i = 1:(m_order+s_order)
+            for j = 1:s_order
+                S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
+            end
         end
-    end
 
-    A = [M;S];
+        A = [M;S];
 
-    d = A\b;
-    ret = x*0;
-    ret(index) = d/h*h_pol;
-end
+        d = A\b;
+        ret = x*0;
+        ret(index) = d/h*h_polynomial;
+    end
 
 end
 
--- a/diracDiscrTest.m	Tue Nov 19 13:47:17 2019 -0800
+++ b/diracDiscrTest.m	Tue Nov 19 13:55:43 2019 -0800
@@ -2,42 +2,19 @@
 	    tests = functiontests(localfunctions);
 end
 
-function testLeftGP(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] = setup1D(order, mom_cond);
-
-        % Test left boundary grid points
-        x0s = xl + [0, h, 2*h];
-
-        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
+%TODO: Test discretizing with smoothness conditions.
+%      Only discretization with moment conditions currently tested.
 
 function testLeftRandom(testCase)
 
     orders = [2, 4, 6];
     mom_conds = orders;
+    rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
 
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [xl, ~, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points near left boundary
         x0s = xl + 2*h*rand(1,10);
@@ -56,42 +33,16 @@
     end
 end
 
-function testRightGP(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] = setup1D(order, mom_cond);
-
-        % Test right boundary grid points
-        x0s = xr-[0, h, 2*h];
-
-        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 testRightRandom(testCase)
 
     orders = [2, 4, 6];
     mom_conds = orders;
+    rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
 
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [~, xr, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points near right boundary
         x0s = xr - 2*h*rand(1,10);
@@ -110,43 +61,16 @@
     end
 end
 
-function testInteriorGP(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] = setup1D(order, mom_cond);
-
-        % Test interior grid points
-        m_half = round(m/2);
-        x0s = xl + (m_half-1:m_half+1)*h;
-
-        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 testInteriorRandom(testCase)
 
     orders = [2, 4, 6];
     mom_conds = orders;
+    rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
 
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test random points in interior
         x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
@@ -174,7 +98,7 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
 
         % Test points outisde grid
         x0s = [xl-1.1*h, xr+1.1*h];
@@ -201,7 +125,7 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
 
         % Test all grid points
         x0s = x;
@@ -228,7 +152,7 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
+        [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
 
         % Test halfway between all grid points
         x0s = 1/2*( x(2:end)+x(1:end-1) );
@@ -247,87 +171,6 @@
     end
 end
 
-% 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)
 
@@ -337,7 +180,7 @@
     for o = 1:length(orders)
         order = orders(o);
         mom_cond = mom_conds(o);
-        [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond);
+        [~, ~, x, X, ~, H, fs] = setup2D(order, mom_cond);
         H_global = kron(H{1}, H{2});
 
         % Test all grid points
@@ -361,11 +204,12 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
+    rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
 
     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);
+        [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond);
         H_global = kron(H{1}, H{2});
 
         xl = xlims{1};
@@ -407,7 +251,7 @@
     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);
+        [~, ~, ~, x, X, ~, H, fs] = setup3D(order, mom_cond);
         H_global = kron(kron(H{1}, H{2}), H{3});
 
         % Test all grid points
@@ -431,11 +275,12 @@
 
     orders = [2, 4, 6];
     mom_conds = orders;
+    rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
 
     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);
+        [xlims, ylims, zlims, x, X, h, H, fs] = setup3D(order, mom_cond);
         H_global = kron(kron(H{1}, H{2}), H{3});
 
         xl = xlims{1};
@@ -475,7 +320,7 @@
 % ======================================================
 % ============== Setup functions =======================
 % ======================================================
-function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond)
+function [xl, xr, h, x, H, fs] = setup1D(order, mom_cond)
 
     % Grid
     xl = -3;
@@ -498,7 +343,7 @@
 
 end
 
-function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond)
+function [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond)
 
     % Grid
     xlims = {-3, 20};
@@ -532,7 +377,7 @@
 
 end
 
-function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)
+function [xlims, ylims, zlims, x, X, h, H, fs] = setup3D(order, mom_cond)
 
     % Grid
     xlims = {-3, 20};
@@ -569,27 +414,4 @@
     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
-    xl = -3;
-    xr = 900;
-    L = xr-xl;
-    m = 101;
-    [~, g_dual] = grid.primalDual1D(m, {xl, xr});
-    x = g_dual.points();
-    h = g_dual.h;
-
-    % Quadrature
-    ops = sbp.D1Staggered(m, {xl, xr}, order);
-    H = ops.H_dual;
-
-    % Moment conditions
-    fs = cell(mom_cond,1);
-    for p = 0:mom_cond-1
-        fs{p+1} = @(x) (x/L).^p;
-    end
-
-end
+end
\ No newline at end of file