Mercurial > repos > public > sbplib
changeset 1234:f1806475498b feature/dirac_discr
- Pass grids to diracDiscr and adjust tests.
- Minor edits in diracDiscr
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Tue, 19 Nov 2019 15:16:08 -0800 |
| parents | 57df0bf741dc |
| children | dea852e85b77 |
| files | diracDiscr.m diracDiscrTest.m |
| diffstat | 2 files changed, 88 insertions(+), 98 deletions(-) [+] |
line wrap: on
line diff
--- a/diracDiscr.m Tue Nov 19 13:55:43 2019 -0800 +++ b/diracDiscr.m Tue Nov 19 15:16:08 2019 -0800 @@ -1,30 +1,29 @@ -function d = diracDiscr(x_s, x, m_order, s_order, H) +function d = diracDiscr(g, x_s, m_order, s_order, H) % n-dimensional delta function + % g: cartesian grid % 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); + assertType(g, 'grid.Cartesian'); + dim = g.d; 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); + % Allow for non-cell input in 1D + if dim == 1 + H = {H}; + end + % Create 1D dirac discr for each coordinate dir. + for i = 1:dim + d_1D{i} = diracDiscr1D(x_s(i), g.x{i}, m_order, s_order, H{i}); + end - 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 + 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 @@ -36,10 +35,9 @@ m = length(x); % Return zeros if x0 is outside grid - if(x_s < x(1) || x_s > x(end) ) - + if x_s < x(1) || x_s > x(end) ret = zeros(size(x)); - + return else fnorm = diag(H); @@ -49,7 +47,7 @@ % Get interior grid spacing middle = floor(m/2); - h = x(middle+1) - x(middle); + h = x(middle+1) - x(middle); % Use middle point to allow for staggerd grids. % Find the indices that are within range of of the point source location ind_delta = find(tot_order*h/2 >= abs(x-x_s)); @@ -62,14 +60,14 @@ end % Use first tot_order grid points - if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h; + 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_order grid points - elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h; + 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; @@ -122,9 +120,3 @@ end end - - - - - -
--- 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