Mercurial > repos > public > sbplib
changeset 837:43a1c3ac07b1 feature/poroelastic
Make diracDiscr work for n dimensions. Add tests up to 3D.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 05 Sep 2018 18:09:18 -0700 |
| parents | 049e4c6fa630 |
| children | bd43256f0c41 |
| files | diracDiscr.m diracDiscrTest.m |
| diffstat | 2 files changed, 265 insertions(+), 16 deletions(-) [+] |
line wrap: on
line diff
--- a/diracDiscr.m Wed Sep 05 14:46:39 2018 -0700 +++ b/diracDiscr.m Wed Sep 05 18:09:18 2018 -0700 @@ -1,4 +1,37 @@ -function ret = diracDiscr(x_0in , x , m_order, s_order, H, h) + +function d = diracDiscr(x_s, x, m_order, s_order, H) + % n-dimensional delta function + % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. + % x: cell array of grid point column vectors for each dimension. + % m_order: Number of moment conditions + % s_order: Number of smoothness conditions + % H: cell array of 1D norm matrices + + dim = length(x_s); + d_1D = cell(dim,1); + + % If 1D, non-cell input is accepted + if dim == 1 && ~iscell(x) + d = diracDiscr1D(x_s, x, m_order, s_order, H); + + else + for i = 1:dim + d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); + end + + d = d_1D{dim}; + for i = dim-1: -1: 1 + % Perform outer product, transpose, and then turn into column vector + d = (d_1D{i}*d')'; + d = d(:); + end + end + +end + + +% Helper function for 1D delta functions +function ret = diracDiscr1D(x_0in , x , m_order, s_order, H) m = length(x); @@ -14,7 +47,7 @@ tot = m_order+s_order; S = []; M = []; - + % Get interior grid spacing middle = floor(m/2); h = x(middle+1) - x(middle); @@ -53,9 +86,9 @@ x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); norm = fnorm(poss)/h; index = poss; - + % Interior - else + 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;
--- a/diracDiscrTest.m Wed Sep 05 14:46:39 2018 -0700 +++ b/diracDiscrTest.m Wed Sep 05 18:09:18 2018 -0700 @@ -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]; @@ -37,7 +37,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 random points near left boundary x0s = xl + 2*h*rand(1,10); @@ -64,7 +64,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 right boundary grid points x0s = xr-[0, h, 2*h]; @@ -91,7 +91,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 random points near right boundary x0s = xr - 2*h*rand(1,10); @@ -118,7 +118,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 interior grid points m_half = round(m/2); @@ -146,7 +146,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 random points in interior x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); @@ -174,7 +174,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 points outisde grid x0s = [xl-1.1*h, xr+1.1*h]; @@ -201,7 +201,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 all grid points x0s = x; @@ -228,7 +228,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 halfway between all grid points x0s = 1/2*( x(2:end)+x(1:end-1) ); @@ -328,9 +328,154 @@ % 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); + [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + + % 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 = diracDiscr(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 + +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}); + + 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(:,1), X(:,2)); + 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 + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +%=============== 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); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + + % Test all grid points + x0s = X; + + 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; + 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 + end +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 - -