Mercurial > repos > public > sbplib
changeset 1285:6b68f939d023 feature/boundary_optimized_grids
Add utility functions for constructing the grids used by the boundary optimized D1Nonequidistant opsets
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 01 Jul 2020 11:15:57 +0200 |
| parents | 5412fd72a7d8 |
| children | 4cb627c7fb90 |
| files | +sbp/+util/accurateBoundaryOptimizedGrid.m +sbp/+util/minimalBoundaryOptimizedGrid.m |
| diffstat | 2 files changed, 105 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+util/accurateBoundaryOptimizedGrid.m Wed Jul 01 11:15:57 2020 +0200 @@ -0,0 +1,56 @@ +function [x,h] = accurateBoundaryOptimizedGrid(L,N,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) ]; + %%%%%%%%%%%%%%%%%%%%%%%%% +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+util/minimalBoundaryOptimizedGrid.m Wed Jul 01 11:15:57 2020 +0200 @@ -0,0 +1,49 @@ +function [x,h] = minimalBoundaryOptimizedGrid(L,N,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) ]; + %%%%%%%%%%%%%%%%%%%%%%%%% +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