Mercurial > repos > public > sbplib
diff +sbp/+implementations/d1_noneq_12.m @ 1286:4cb627c7fb90 feature/boundary_optimized_grids
Make D1Nonequidistant use the grid generation functions accurate/minimalBoundaryOptimizedGrid and remove grid generation from +implementations
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Wed, 01 Jul 2020 13:43:32 +0200 |
parents | f7ac3cd6eeaa |
children |
line wrap: on
line diff
--- a/+sbp/+implementations/d1_noneq_12.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_12.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,50 +1,12 @@ -function [D1,H,x,h] = d1_noneq_12(N,L) +function [D1,H] = d1_noneq_12(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<24) error('Operator requires at least 24 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 12; -m = 6; -order = 12; - -%%%% Non-equidistant grid points %%%%% -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; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; -x9 = 9.0000000000000e+00; -x10 = 1.0000000000000e+01; -x11 = 1.1000000000000e+01; -x12 = 1.2000000000000e+01; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% 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) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -73,22 +35,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 12; +d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);