Mercurial > repos > public > sbplib
changeset 971:e54c2f54dbfe feature/getBoundaryOperator
Merge with feature/poroelastic. Use only the changes made to multiblock.DiffOp and scheme.Elastic2dVariable. DiffOp.getBoundaryOperator/Quadrature now use scheme methods instead of propeties.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 25 Dec 2018 07:50:07 +0100 |
| parents | 2412f407749a (diff) 23d9ca6755be (current diff) |
| children | 47b48a97c675 |
| files | +multiblock/DiffOp.m +multiblock/Grid.m +sbp/+implementations/d2_variable_2.m +sbp/+implementations/d2_variable_4.m +sbp/+implementations/d4_variable_6.m +time/+rkparameters/rk4.m +time/ExplicitRungeKuttaDiscreteData.m +time/ExplicitRungeKuttaSecondOrderDiscreteData.m .hgtags diracDiscr.m diracDiscrTest.m |
| diffstat | 11 files changed, 5 insertions(+), 1018 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Tue Dec 25 07:23:38 2018 +0100 +++ b/+multiblock/DiffOp.m Tue Dec 25 07:50:07 2018 +0100 @@ -151,6 +151,7 @@ end end + % Get square matrix that integrates the solution restricted to a boundary function op = getBoundaryQuadrature(obj, boundary) switch class(boundary) case 'cell'
--- a/+multiblock/Grid.m Tue Dec 25 07:23:38 2018 +0100 +++ b/+multiblock/Grid.m Tue Dec 25 07:50:07 2018 +0100 @@ -132,27 +132,6 @@ end - % Pads a grid function that lives on a subgrid with - % zeros and gives it the size that mathces obj. - function gf = expandFunc(obj, gfSub, subGridId) - nComponents = length(gfSub)/obj.grids{subGridId}.N(); - nBlocks = numel(obj.grids); - - % Create sparse block matrix - gfs = cell(nBlocks,1); - for i = 1:nBlocks - N = obj.grids{i}.N()*nComponents; - gfs{i} = sparse(N, 1); - end - - % Insert local gf - gfs{subGridId} = gfSub; - - % Convert cell to vector - gf = blockmatrix.toMatrix(gfs); - - end - % Find all non interface boundaries of all blocks. % Return their grid.boundaryIdentifiers in a cell array. function bs = getBoundaryNames(obj)
--- a/+sbp/+implementations/d2_variable_2.m Tue Dec 25 07:23:38 2018 +0100 +++ b/+sbp/+implementations/d2_variable_2.m Tue Dec 25 07:50:07 2018 +0100 @@ -27,7 +27,7 @@ diags = -1:1; stencil = [-1/2 0 1/2]; D1 = stripeMatrix(stencil, diags, m); - + D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; D1(m,m-1)=-1;D1(m,m)=1; D1=D1/h; @@ -40,7 +40,7 @@ scheme_radius = (scheme_width-1)/2; r = (1+scheme_radius):(m-scheme_radius); - function [D2, B] = D2_fun(c) + function D2 = D2_fun(c) Mm1 = -c(r-1)/2 - c(r)/2; M0 = c(r-1)/2 + c(r) + c(r+1)/2; @@ -54,8 +54,6 @@ M=M/h; D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); - B = HI*M; end D2 = @D2_fun; - end \ No newline at end of file
--- a/+sbp/+implementations/d2_variable_4.m Tue Dec 25 07:23:38 2018 +0100 +++ b/+sbp/+implementations/d2_variable_4.m Tue Dec 25 07:50:07 2018 +0100 @@ -49,7 +49,7 @@ N = m; - function [D2, B] = D2_fun(c) + function D2 = D2_fun(c) M = 78+(N-12)*5; %h = 1/(N-1); @@ -131,8 +131,6 @@ cols(40+(i-7)*5:44+(i-7)*5) = [i-2;i-1;i;i+1;i+2]; end D2 = sparse(rows,cols,D2); - - B = HI*( c(end)*e_r*d1_r' - c(1)*e_l*d1_l') - D2; end D2 = @D2_fun; end \ No newline at end of file
--- a/+sbp/+implementations/d4_variable_6.m Tue Dec 25 07:23:38 2018 +0100 +++ b/+sbp/+implementations/d4_variable_6.m Tue Dec 25 07:50:07 2018 +0100 @@ -85,7 +85,7 @@ scheme_radius = (scheme_width-1)/2; r = (1+scheme_radius):(m-scheme_radius); - function [D2, B] = D2_fun(c) + function D2 = D2_fun(c) Mm3 = c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r); Mm2 = c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2); @@ -128,7 +128,6 @@ M=M/h; D2 = HI*(-M - c(1)*e_l*d1_l' + c(m)*e_r*d1_r'); - B = HI*M; end D2 = @D2_fun;
--- a/+time/+rkparameters/rk4.m Tue Dec 25 07:23:38 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function [a,b,c,s] = rk4() - -% Butcher tableau for classical RK$ -s = 4; -a = sparse(s,s); -a(2,1) = 1/2; -a(3,2) = 1/2; -a(4,3) = 1; -b = 1/6*[1; 2; 2; 1]; -c = [0; 1/2; 1/2; 1]; - -end \ No newline at end of file
--- a/+time/ExplicitRungeKuttaDiscreteData.m Tue Dec 25 07:23:38 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,121 +0,0 @@ -classdef ExplicitRungeKuttaDiscreteData < time.Timestepper - properties - D - S % Function handle for time-dependent data - data % Matrix of data vectors, one column per stage - k - t - v - m - n - order - a, b, c, s % Butcher tableau - K % Stage rates - U % Stage approximations - T % Stage times - end - - - methods - function obj = ExplicitRungeKuttaDiscreteData(D, S, data, k, t0, v0, order) - default_arg('order', 4); - default_arg('S', []); - default_arg('data', []); - - obj.D = D; - obj.S = S; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.m = length(v0); - obj.n = 0; - obj.order = order; - obj.data = data; - - switch order - case 4 - [obj.a, obj.b, obj.c, obj.s] = time.rkparameters.rk4(); - otherwise - error('That RK method is not available'); - end - - obj.K = sparse(obj.m, obj.s); - obj.U = sparse(obj.m, obj.s); - - end - - function [v,t,U,T,K] = getV(obj) - v = obj.v; - t = obj.t; - U = obj.U; % Stage approximations in previous time step. - T = obj.T; % Stage times in previous time step. - K = obj.K; % Stage rates in previous time step. - end - - function [a,b,c,s] = getTableau(obj) - a = obj.a; - b = obj.b; - c = obj.c; - s = obj.s; - end - - % Returns quadrature weights for stages in one time step - function quadWeights = getTimeStepQuadrature(obj) - [~, b] = obj.getTableau(); - quadWeights = obj.k*b; - end - - function obj = step(obj) - v = obj.v; - a = obj.a; - b = obj.b; - c = obj.c; - s = obj.s; - S = obj.S; - dt = obj.k; - K = obj.K; - U = obj.U; - D = obj.D; - data = obj.data; - - for i = 1:s - U(:,i) = v; - for j = 1:i-1 - U(:,i) = U(:,i) + dt*a(i,j)*K(:,j); - end - - K(:,i) = D*U(:,i); - obj.T(i) = obj.t + c(i)*dt; - - % Data from continuous function and discrete time-points. - if ~isempty(S) - K(:,i) = K(:,i) + S(obj.T(i)); - end - if ~isempty(data) - K(:,i) = K(:,i) + data(:,obj.n*s + i); - end - - end - - obj.v = v + dt*K*b; - obj.t = obj.t + dt; - obj.n = obj.n + 1; - obj.U = U; - obj.K = K; - end - end - - - methods (Static) - function k = getTimeStep(lambda, order) - default_arg('order', 4); - switch order - case 4 - k = time.rk4.get_rk4_time_step(lambda); - otherwise - error('Time-step function not available for this order'); - end - end - end - -end \ No newline at end of file
--- a/+time/ExplicitRungeKuttaSecondOrderDiscreteData.m Tue Dec 25 07:23:38 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,129 +0,0 @@ -classdef ExplicitRungeKuttaSecondOrderDiscreteData < time.Timestepper - properties - k - t - w - m - D - E - M - C_cont % Continuous part (function handle) of forcing on first order form. - C_discr% Discrete part (matrix) of forcing on first order form. - n - order - tsImplementation % Time stepper object, RK first order form, - % which this wraps around. - end - - - methods - % Solves u_tt = Du + Eu_t + S by - % Rewriting on first order form: - % w_t = M*w + C(t) - % where - % M = [ - % 0, I; - % D, E; - % ] - % and - % C(t) = [ - % 0; - % S(t) - % ] - % D, E, should be matrices (or empty for zero) - % They can also be omitted by setting them equal to the empty matrix. - % S = S_cont + S_discr, where S_cont is a function handle - % S_discr a matrix of data vectors, one column per stage. - function obj = ExplicitRungeKuttaSecondOrderDiscreteData(D, E, S_cont, S_discr, k, t0, v0, v0t, order) - default_arg('order', 4); - default_arg('S_cont', []); - default_arg('S_discr', []); - obj.D = D; - obj.E = E; - obj.m = length(v0); - obj.n = 0; - - default_arg('D', sparse(obj.m, obj.m) ); - default_arg('E', sparse(obj.m, obj.m) ); - - obj.k = k; - obj.t = t0; - obj.w = [v0; v0t]; - - I = speye(obj.m); - O = sparse(obj.m,obj.m); - - obj.M = [ - O, I; - D, E; - ]; - - % Build C_cont - if ~isempty(S_cont) - obj.C_cont = @(t)[ - sparse(obj.m,1); - S_cont(t) - ]; - else - obj.C_cont = []; - end - - % Build C_discr - if ~isempty(S_discr) - [~, nt] = size(S_discr); - obj.C_discr = [sparse(obj.m, nt); - S_discr - ]; - else - obj.C_discr = []; - end - obj.tsImplementation = time.ExplicitRungeKuttaDiscreteData(obj.M, obj.C_cont, obj.C_discr,... - k, obj.t, obj.w, order); - end - - function [v,t,U,T,K] = getV(obj) - [w,t,U,T,K] = obj.tsImplementation.getV(); - - v = w(1:end/2); - U = U(1:end/2, :); % Stage approximations in previous time step. - K = K(1:end/2, :); % Stage rates in previous time step. - % T: Stage times in previous time step. - end - - function [vt,t,U,T,K] = getVt(obj) - [w,t,U,T,K] = obj.tsImplementation.getV(); - - vt = w(end/2+1:end); - U = U(end/2+1:end, :); % Stage approximations in previous time step. - K = K(end/2+1:end, :); % Stage rates in previous time step. - % T: Stage times in previous time step. - end - - function [a,b,c,s] = getTableau(obj) - [a,b,c,s] = obj.tsImplementation.getTableau(); - end - - % Returns quadrature weights for stages in one time step - function quadWeights = getTimeStepQuadrature(obj) - [~, b] = obj.getTableau(); - quadWeights = obj.k*b; - end - - % Use RK for first order form to step - function obj = step(obj) - obj.tsImplementation.step(); - [v, t] = obj.tsImplementation.getV(); - obj.w = v; - obj.t = t; - obj.n = obj.n + 1; - end - end - - methods (Static) - function k = getTimeStep(lambda, order) - default_arg('order', 4); - k = obj.tsImplementation.getTimeStep(lambda, order); - end - end - -end \ No newline at end of file
--- a/.hgtags Tue Dec 25 07:23:38 2018 +0100 +++ b/.hgtags Tue Dec 25 07:50:07 2018 +0100 @@ -2,4 +2,3 @@ 0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 b723495cdb2f96314d7b3f0aa79723a7dc088c7d v0.2 08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0 -08f3ffe63f484d02abce8df4df61e826f568193f Heimisson2018
--- a/diracDiscr.m Tue Dec 25 07:23:38 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,130 +0,0 @@ - -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); - -% 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 - - - - - -
--- a/diracDiscrTest.m Tue Dec 25 07:23:38 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,595 +0,0 @@ -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] = 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 - -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] = setup1D(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] = 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; - - 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 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] = 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; - - 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 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] = setup1D(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] = setup1D(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] = setup1D(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 - -%=============== 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] = 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); - 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 [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 - 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