Mercurial > repos > public > sbplib

changeset 836:049e4c6fa630 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Sep 2018 14:46:39 -0700
parents 9d4246ac94c0
children 43a1c3ac07b1
files diracDiscr.m diracDiscrTest.m
diffstat 2 files changed, 476 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscr.m	Wed Sep 05 14:46:39 2018 -0700
@@ -0,0 +1,97 @@
+function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
+
+m = length(x);
+
+% Return zeros if x0 is outside grid
+if(x_0in < x(1) || x_0in > x(end) )
+
+    ret = zeros(size(x));
+
+else
+
+    fnorm = diag(H);
+    eta = abs(x-x_0in);
+    tot = m_order+s_order;
+    S = [];
+    M = [];
+    
+    % Get interior grid spacing
+    middle = floor(m/2);
+    h = x(middle+1) - x(middle);
+
+    poss = find(tot*h/2 >= eta);
+
+    % 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
+
+    % 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 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));
+
+    % Interior, compensate for round-off errors.
+    elseif length(poss) < tot
+        if poss(end)<m
+            poss = [poss; poss(end)+1];
+        else
+            poss = [poss(1)-1; poss];
+        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);
+
+    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);
+        end
+    end
+
+    for i = 1:(m_order+s_order)
+        for j = 1:s_order
+            S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
+        end
+    end
+
+    A = [M;S];
+
+    d = A\b;
+    ret = x*0;
+    ret(index) = d/h*h_pol;
+end
+
+end
+
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscrTest.m	Wed Sep 05 14:46:39 2018 -0700
@@ -0,0 +1,379 @@
+function tests = diracDiscrTest()
+	    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] = setupStuff(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
+
+function testLeftRandom(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] = setupStuff(order, mom_cond);
+
+        % 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(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 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] = setupStuff(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;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+
+        % 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(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 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] = setupStuff(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;
+
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+
+        % 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(x0, x, 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);
+        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+
+        % 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(x0, x, 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);
+        [xl, xr, m, h, x, H, fs] = setupStuff(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 testHalfGP(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] = setupStuff(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 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
+
+
+% ============== Setup functions =======================
+function [xl, xr, m, h, x, H, fs] = setupStuff(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);
+    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 [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
+
+