Mercurial > repos > public > sbplib
changeset 1297:e53b1e25970a feature/boundary_optimized_grids
Change +sbp/+util/ to +sbp/+grid and change function names to camel case
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Tue, 07 Jul 2020 16:08:08 +0200 |
| parents | 2853b655c172 |
| children | 0ffb5bfa65e4 |
| files | +grid/boundaryOptimized.m +grid/boundaryOptimizedTest.m +grid/boundaryoptimized.m +grid/boundaryoptimizedTest.m +multiblock/+domain/Rectangle.m +sbp/+grid/accurateBoundaryOptimizedGrid.m +sbp/+grid/minimalBoundaryOptimizedGrid.m +sbp/+util/accurateBoundaryOptimizedGrid.m +sbp/+util/minimalBoundaryOptimizedGrid.m |
| diffstat | 9 files changed, 280 insertions(+), 280 deletions(-) [+] |
line wrap: on
line diff
diff -r 2853b655c172 -r e53b1e25970a +grid/boundaryOptimized.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimized.m Tue Jul 07 16:08:08 2020 +0200 @@ -0,0 +1,54 @@ +% Creates a Cartesian grid of dimension length(m) +% over the domain xlim, ylim, ... +% The grid is non-equidistant in the boundary regions, +% with node placement based on boundary-optimized SBP operators. +% Examples: +% g = grid.boundaryOptimized([mx, my], xlim, ylim, order, opt) +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4) - defaults to 'accurate' stencils +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4, 'minimal') +function g = boundaryOptimized(m, varargin) + n = length(m); + + % Check that parameters matches dimensions + matchingParams = false; + if length(varargin) == n+1 % Minimal number of arguments + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]); + elseif length(varargin) == n+2 % Stencil options supplied + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]) && ... + ischar([varargin{n+2}]); + end + assert(matchingParams,'grid:boundaryOptimized:NonMatchingParameters','The number of parameters per dimensions do not match.'); + + % Check that stencil options are passed correctly (if supplied) + if length(varargin) == n+2 % Stencil options supplied + availabe_opts = ["Accurate","accurate","A","Minimal","minimal","M"]; + assert(any(varargin{n+2} == availabe_opts), ... + 'grid:boundaryOptimized:InvalidOption',"The operator option must be 'accurate' or 'minimal.'"); + else %If not passed, populate varargin with default option 'accurate' + varargin(n+2) = {'accurate'}; + end + + % Specify generating function + switch varargin{n+2} + case {'Accurate','accurate','A'} + gridgenerator = @sbp.grid.accurateBoundaryOptimizedGrid; + case {'Minimal','minimal','M'} + gridgenerator = @sbp.grid.minimalBoundaryOptimizedGrid; + end + + X = {}; + h = []; + for i = 1:n + try + [X{i},h(i)] = gridgenerator(varargin{i},m(i),varargin{n+1}); + catch exception % Propagate any errors in the grid generation functions. + msgText = getReport(exception); + error('grid:boundaryOptimized:InvalidParameter',msgText) + end + end + + g = grid.Cartesian(X{:}); + g.h = h; +end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +grid/boundaryOptimizedTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimizedTest.m Tue Jul 07 16:08:08 2020 +0200 @@ -0,0 +1,112 @@ +function tests = boundaryOptimizedTest() + tests = functiontests(localfunctions); +end + +function testErrorInvalidParam(testCase) + in = { + %Invalid order + {[10 10],{0,1},{0,2},3}, + %Invalid grid size + {5, {0,1}, 4}, + {[10 5],{0,1},{0,2},4}, + {[10 5],{0,1},{0,2},6,'M'}, + %Invalid limits + {10,{1},4}, + {[10,10],{0,1},{1},4}, + {[10,10],{1},{1,0},4}, + {10,{1,0},4}, + {[10, 5],{1,0},{0,-1},4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidParameter',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorInvalidOption(testCase) + in = { + {[8 8],{0,1},{0,2},4,'acrurate'}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidOption',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorNonMatchingParam(testCase) + in = { + {[],{1},4}, + {[],{0,1},{0,1},4}, + {[5,5],{0,1},{0,1},{0,1},4}, + {[5,5,4],{0,1},{0,1},4,'accurate'} + {[5,5,4],{0,1},{0,1},{0,1},4,4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +% Tests that the expected grid points are obtained for a boundary optimized grid with a 4th order +% accurate stencil and 8th order minimal stencil. +% The boundary grid point distance weights are taken from the D1Nonequidistant operators and +% grid spacing is calculated according to Mattsson et al 2018. The test uses minimal number of grid +% points required by the operators. +function testCompiles(testCase) + + %% 1D 4th order accurate stencil + % Boundary weights, number of non-equidistantly spaced points for 4th order accurate stencil + bw = [0.0000000000000e+00 6.8764546205559e-01 1.8022115125776e+00]; + n = length(bw)-1; + xi_n = bw(end); + + % Grid points in x-direction. + Lx = 1; + mx = 8; + hx_4 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_4*bw; + bp_r = Lx-flip(hx_4*bw); + interior = [hx_4*(xi_n+1) hx_4*(xi_n+2)]; + x_4 = [bp_l interior bp_r]; + + % Boundary weights, number of non-equidistantly spaced points for 8th order minimal stencil + bw = [0.0000000000000e+00, 4.9439570885261e-01, 1.4051531374839e+00]; + n = length(bw)-1; + xi_n = bw(end); + + %% 2D 8th order minimal stencil + % Grid points in x-direction. + hx_8 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_8*bw; + bp_r = Lx-flip(hx_8*bw); + interior = [hx_8*(xi_n+1) hx_8*(xi_n+2)]; + x_8 = [bp_l interior bp_r]; + + % Grid points in y-direction. + Ly = 2; + my = 9; + hy = Ly/(2*xi_n+(my-2*n-1)); + + bp_l = hy*bw; + bp_r = Ly-flip(hy*bw); + interior = [hy*(xi_n+1) hy*(xi_n+2) hy*(xi_n+3)]; + y = [bp_l interior bp_r]; + + in = { + {mx, {0,Lx},4}, + {[mx, my],{0,Lx},{0,Ly},8,'M'}, + }; + + out = { + {[x_4'],hx_4} + {[kr(x_8',ones(size(y'))),kr(ones(size(x_8')),y')],[hx_8, hy]} + }; + + for i = 1:length(in) + g = grid.boundaryOptimized(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}{1},'AbsTol', 1e-14, 'RelTol', 1e-14); + testCase.verifyEqual(g.scaling(),out{i}{2},'AbsTol', 1e-14, 'RelTol', 1e-14); + end +end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +grid/boundaryoptimized.m --- a/+grid/boundaryoptimized.m Tue Jul 07 16:00:24 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,54 +0,0 @@ -% Creates a Cartesian grid of dimension length(m) -% over the domain xlim, ylim, ... -% The grid is non-equidistant in the boundary regions, -% with node placement based on boundary-optimized SBP operators. -% Examples: -% g = grid.boundaryoptimized([mx, my], xlim, ylim, order, opt) -% g = grid.boundaryoptimized([10, 15], {0,1}, {0,2}, 4) - defaults to 'accurate' stencils -% g = grid.boundaryoptimized([10, 15], {0,1}, {0,2}, 4, 'minimal') -function g = boundaryoptimized(m, varargin) - n = length(m); - - % Check that parameters matches dimensions - matchingParams = false; - if length(varargin) == n+1 % Minimal number of arguments - matchingParams = iscell([varargin{1:n}]) && ... - isfloat([varargin{n+1}]); - elseif length(varargin) == n+2 % Stencil options supplied - matchingParams = iscell([varargin{1:n}]) && ... - isfloat([varargin{n+1}]) && ... - ischar([varargin{n+2}]); - end - assert(matchingParams,'grid:boundaryoptimized:NonMatchingParameters','The number of parameters per dimensions do not match.'); - - % Check that stencil options are passed correctly (if supplied) - if length(varargin) == n+2 % Stencil options supplied - availabe_opts = ["Accurate","accurate","A","Minimal","minimal","M"]; - assert(any(varargin{n+2} == availabe_opts), ... - 'grid:boundaryoptimized:InvalidOption',"The operator option must be 'accurate' or 'minimal.'"); - else %If not passed, populate varargin with default option 'accurate' - varargin(n+2) = {'accurate'}; - end - - % Specify generating function - switch varargin{n+2} - case {'Accurate','accurate','A'} - gridgenerator = @sbp.util.accurateBoundaryOptimizedGrid; - case {'Minimal','minimal','M'} - gridgenerator = @sbp.util.minimalBoundaryOptimizedGrid; - end - - X = {}; - h = []; - for i = 1:n - try - [X{i},h(i)] = gridgenerator(varargin{i},m(i),varargin{n+1}); - catch exception % Propagate any errors in the grid generation functions. - msgText = getReport(exception); - error('grid:boundaryoptimized:InvalidParameter',msgText) - end - end - - g = grid.Cartesian(X{:}); - g.h = h; -end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +grid/boundaryoptimizedTest.m --- a/+grid/boundaryoptimizedTest.m Tue Jul 07 16:00:24 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,112 +0,0 @@ -function tests = boundaryoptimizedTest() - tests = functiontests(localfunctions); -end - -function testErrorInvalidParam(testCase) - in = { - %Invalid order - {[10 10],{0,1},{0,2},3}, - %Invalid grid size - {5, {0,1}, 4}, - {[10 5],{0,1},{0,2},4}, - {[10 5],{0,1},{0,2},6,'M'}, - %Invalid limits - {10,{1},4}, - {[10,10],{0,1},{1},4}, - {[10,10],{1},{1,0},4}, - {10,{1,0},4}, - {[10, 5],{1,0},{0,-1},4}, - }; - - for i = 1:length(in) - testCase.verifyError(@()grid.boundaryoptimized(in{i}{:}),'grid:boundaryoptimized:InvalidParameter',sprintf('in(%d) = %s',i,toString(in{i}))); - end -end - -function testErrorInvalidOption(testCase) - in = { - {[8 8],{0,1},{0,2},4,'acrurate'}, - }; - - for i = 1:length(in) - testCase.verifyError(@()grid.boundaryoptimized(in{i}{:}),'grid:boundaryoptimized:InvalidOption',sprintf('in(%d) = %s',i,toString(in{i}))); - end -end - -function testErrorNonMatchingParam(testCase) - in = { - {[],{1},4}, - {[],{0,1},{0,1},4}, - {[5,5],{0,1},{0,1},{0,1},4}, - {[5,5,4],{0,1},{0,1},4,'accurate'} - {[5,5,4],{0,1},{0,1},{0,1},4,4}, - }; - - for i = 1:length(in) - testCase.verifyError(@()grid.boundaryoptimized(in{i}{:}),'grid:boundaryoptimized:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); - end -end - -% Tests that the expected grid points are obtained for a boundary optimized grid with a 4th order -% accurate stencil and 8th order minimal stencil. -% The boundary grid point distance weights are taken from the D1Nonequidistant operators and -% grid spacing is calculated according to Mattsson et al 2018. The test uses minimal number of grid -% points required by the operators. -function testCompiles(testCase) - - %% 1D 4th order accurate stencil - % Boundary weights, number of non-equidistantly spaced points for 4th order accurate stencil - bw = [0.0000000000000e+00 6.8764546205559e-01 1.8022115125776e+00]; - n = length(bw)-1; - xi_n = bw(end); - - % Grid points in x-direction. - Lx = 1; - mx = 8; - hx_4 = Lx/(2*xi_n+(mx-2*n-1)); - - bp_l = hx_4*bw; - bp_r = Lx-flip(hx_4*bw); - interior = [hx_4*(xi_n+1) hx_4*(xi_n+2)]; - x_4 = [bp_l interior bp_r]; - - % Boundary weights, number of non-equidistantly spaced points for 8th order minimal stencil - bw = [0.0000000000000e+00, 4.9439570885261e-01, 1.4051531374839e+00]; - n = length(bw)-1; - xi_n = bw(end); - - %% 2D 8th order minimal stencil - % Grid points in x-direction. - hx_8 = Lx/(2*xi_n+(mx-2*n-1)); - - bp_l = hx_8*bw; - bp_r = Lx-flip(hx_8*bw); - interior = [hx_8*(xi_n+1) hx_8*(xi_n+2)]; - x_8 = [bp_l interior bp_r]; - - % Grid points in y-direction. - Ly = 2; - my = 9; - hy = Ly/(2*xi_n+(my-2*n-1)); - - bp_l = hy*bw; - bp_r = Ly-flip(hy*bw); - interior = [hy*(xi_n+1) hy*(xi_n+2) hy*(xi_n+3)]; - y = [bp_l interior bp_r]; - - in = { - {mx, {0,Lx},4}, - {[mx, my],{0,Lx},{0,Ly},8,'M'}, - }; - - out = { - {[x_4'],hx_4} - {[kr(x_8',ones(size(y'))),kr(ones(size(x_8')),y')],[hx_8, hy]} - }; - - for i = 1:length(in) - g = grid.boundaryoptimized(in{i}{:}); - testCase.verifyEqual(g.points(),out{i}{1},'AbsTol', 1e-14, 'RelTol', 1e-14); - testCase.verifyEqual(g.scaling(),out{i}{2},'AbsTol', 1e-14, 'RelTol', 1e-14); - end -end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +multiblock/+domain/Rectangle.m --- a/+multiblock/+domain/Rectangle.m Tue Jul 07 16:00:24 2020 +0200 +++ b/+multiblock/+domain/Rectangle.m Tue Jul 07 16:08:08 2020 +0200 @@ -148,7 +148,7 @@ elseif strcmp(varargin{1},'boundaryopt') order = varargin{2}; stenciloption = varargin{3}; - gridgenerator = @(m,xlim,ylim)grid.boundaryoptimized(m,xlim,ylim,... + gridgenerator = @(m,xlim,ylim)grid.boundaryOptimized(m,xlim,ylim,... order,stenciloption); else error('No grid type supplied!');
diff -r 2853b655c172 -r e53b1e25970a +sbp/+grid/accurateBoundaryOptimizedGrid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/accurateBoundaryOptimizedGrid.m Tue Jul 07 16:08:08 2020 +0200 @@ -0,0 +1,60 @@ +function [x,h] = accurateBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limits should be cell arrays with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Compute h %%%%%%%%%% + h = L/(2*xb(end) + N-1-2*m); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Define grid %%%%%%%% + x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 6.8764546205559e-01; + x2 = 1.8022115125776e+00; + xb = [x0 x1 x2]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.4090263368623e-01; + x2 = 1.2855984345073e+00; + x3 = 2.2638953951239e+00; + xb = [x0 x1 x2 x3]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 3.8118550247622e-01; + x2 = 1.1899550868338e+00; + x3 = 2.2476300175641e+00; + x4 = 3.3192851303204e+00; + xb = [x0 x1 x2 x3 x4]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 3.5902433622052e-01; + x2 = 1.1436659188355e+00; + x3 = 2.2144895894456e+00; + x4 = 3.3682742337736e+00; + x5 = 4.4309689056870e+00; + xb = [x0 x1 x2 x3 x4 x5]'; + case 12 + x0 = 0.0000000000000e+00; + x1 = 3.6098032343909e-01; + x2 = 1.1634317168086e+00; + x3 = 2.2975905356987e+00; + x4 = 3.6057529790929e+00; + x5 = 4.8918275675510e+00; + x6 = 6.0000000000000e+00; + xb = [x0 x1 x2 x3 x4 x5 x6]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +sbp/+grid/minimalBoundaryOptimizedGrid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/minimalBoundaryOptimizedGrid.m Tue Jul 07 16:08:08 2020 +0200 @@ -0,0 +1,53 @@ +function [x,h] = minimalBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limits should be cell arrays with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Compute h %%%%%%%%%% + h = L/(2*xb(end) + N-1-2*m); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Define grid %%%%%%%% + x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end + +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 7.7122987842562e-01; + xb = [x0 x1]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.0842950991998e-01; + x2 = 1.1968523189207e+00; + xb = [x0 x1 x2]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 4.9439570885261e-01; + x2 = 1.4051531374839e+00; + xb = [x0 x1 x2]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 5.8556160757529e-01; + x2 = 1.7473267488572e+00; + x3 = 3.0000000000000e+00; + xb = [x0 x1 x2 x3]'; + case 12 + x0 = 0.0000000000000e+00; + x1 = 4.6552112904489e-01; + x2 = 1.4647984306493e+00; + x3 = 2.7620429464763e+00; + x4 = 4.0000000000000e+00; + xb = [x0 x1 x2 x3 x4]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +sbp/+util/accurateBoundaryOptimizedGrid.m --- a/+sbp/+util/accurateBoundaryOptimizedGrid.m Tue Jul 07 16:00:24 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -function [x,h] = accurateBoundaryOptimizedGrid(lim,N,order) - assert(iscell(lim) && numel(lim) == 2,'The limits should be cell arrays with 2 elements.'); - L = lim{2} - lim{1}; - assert(L>0,'Limits must be given in increasing order.'); - %%%% Non-equidistant grid points %%%%% - xb = boundaryPoints(order); - m = length(xb)-1; % Number of non-equidistant points - assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - %%%% Compute h %%%%%%%%%% - h = L/(2*xb(end) + N-1-2*m); - %%%%%%%%%%%%%%%%%%%%%%%%% - - %%%% Define grid %%%%%%%% - x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; - x = x + lim{1}; - %%%%%%%%%%%%%%%%%%%%%%%%% -end -function xb = boundaryPoints(order) - switch order - case 4 - x0 = 0.0000000000000e+00; - x1 = 6.8764546205559e-01; - x2 = 1.8022115125776e+00; - xb = [x0 x1 x2]'; - case 6 - x0 = 0.0000000000000e+00; - x1 = 4.4090263368623e-01; - x2 = 1.2855984345073e+00; - x3 = 2.2638953951239e+00; - xb = [x0 x1 x2 x3]'; - case 8 - x0 = 0.0000000000000e+00; - x1 = 3.8118550247622e-01; - x2 = 1.1899550868338e+00; - x3 = 2.2476300175641e+00; - x4 = 3.3192851303204e+00; - xb = [x0 x1 x2 x3 x4]'; - case 10 - x0 = 0.0000000000000e+00; - x1 = 3.5902433622052e-01; - x2 = 1.1436659188355e+00; - x3 = 2.2144895894456e+00; - x4 = 3.3682742337736e+00; - x5 = 4.4309689056870e+00; - xb = [x0 x1 x2 x3 x4 x5]'; - case 12 - x0 = 0.0000000000000e+00; - x1 = 3.6098032343909e-01; - x2 = 1.1634317168086e+00; - x3 = 2.2975905356987e+00; - x4 = 3.6057529790929e+00; - x5 = 4.8918275675510e+00; - x6 = 6.0000000000000e+00; - xb = [x0 x1 x2 x3 x4 x5 x6]'; - otherwise - error('Invalid operator order %d.',order); - end -end \ No newline at end of file
diff -r 2853b655c172 -r e53b1e25970a +sbp/+util/minimalBoundaryOptimizedGrid.m --- a/+sbp/+util/minimalBoundaryOptimizedGrid.m Tue Jul 07 16:00:24 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -function [x,h] = minimalBoundaryOptimizedGrid(lim,N,order) - assert(iscell(lim) && numel(lim) == 2,'The limits should be cell arrays with 2 elements.'); - L = lim{2} - lim{1}; - assert(L>0,'Limits must be given in increasing order.'); - %%%% Non-equidistant grid points %%%%% - xb = boundaryPoints(order); - m = length(xb)-1; % Number of non-equidistant points - assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - %%%% Compute h %%%%%%%%%% - h = L/(2*xb(end) + N-1-2*m); - %%%%%%%%%%%%%%%%%%%%%%%%% - - %%%% Define grid %%%%%%%% - x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; - x = x + lim{1}; - %%%%%%%%%%%%%%%%%%%%%%%%% -end - -function xb = boundaryPoints(order) - switch order - case 4 - x0 = 0.0000000000000e+00; - x1 = 7.7122987842562e-01; - xb = [x0 x1]'; - case 6 - x0 = 0.0000000000000e+00; - x1 = 4.0842950991998e-01; - x2 = 1.1968523189207e+00; - xb = [x0 x1 x2]'; - case 8 - x0 = 0.0000000000000e+00; - x1 = 4.9439570885261e-01; - x2 = 1.4051531374839e+00; - xb = [x0 x1 x2]'; - case 10 - x0 = 0.0000000000000e+00; - x1 = 5.8556160757529e-01; - x2 = 1.7473267488572e+00; - x3 = 3.0000000000000e+00; - xb = [x0 x1 x2 x3]'; - case 12 - x0 = 0.0000000000000e+00; - x1 = 4.6552112904489e-01; - x2 = 1.4647984306493e+00; - x3 = 2.7620429464763e+00; - x4 = 4.0000000000000e+00; - xb = [x0 x1 x2 x3 x4]'; - otherwise - error('Invalid operator order %d.',order); - end -end \ No newline at end of file