Mercurial > repos > public > sbplib

comparison +sbp/+grid/accurateBoundaryOptimizedGrid.m @ 1297:e53b1e25970a feature/boundary_optimized_grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 +sbp/+util/accurateBoundaryOptimizedGrid.m@e059a43bb675
children 0ffb5bfa65e4
comparison
equal deleted inserted replaced
1296:2853b655c172 1297:e53b1e25970a
1 function [x,h] = accurateBoundaryOptimizedGrid(lim,N,order)
2 assert(iscell(lim) && numel(lim) == 2,'The limits should be cell arrays with 2 elements.');
3 L = lim{2} - lim{1};
4 assert(L>0,'Limits must be given in increasing order.');
5 %%%% Non-equidistant grid points %%%%%
6 xb = boundaryPoints(order);
7 m = length(xb)-1; % Number of non-equidistant points
8 assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1));
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10
11 %%%% Compute h %%%%%%%%%%
12 h = L/(2*xb(end) + N-1-2*m);
13 %%%%%%%%%%%%%%%%%%%%%%%%%
14
15 %%%% Define grid %%%%%%%%
16 x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ];
17 x = x + lim{1};
18 %%%%%%%%%%%%%%%%%%%%%%%%%
19 end
20 function xb = boundaryPoints(order)
21 switch order
22 case 4
23 x0 = 0.0000000000000e+00;
24 x1 = 6.8764546205559e-01;
25 x2 = 1.8022115125776e+00;
26 xb = [x0 x1 x2]';
27 case 6
28 x0 = 0.0000000000000e+00;
29 x1 = 4.4090263368623e-01;
30 x2 = 1.2855984345073e+00;
31 x3 = 2.2638953951239e+00;
32 xb = [x0 x1 x2 x3]';
33 case 8
34 x0 = 0.0000000000000e+00;
35 x1 = 3.8118550247622e-01;
36 x2 = 1.1899550868338e+00;
37 x3 = 2.2476300175641e+00;
38 x4 = 3.3192851303204e+00;
39 xb = [x0 x1 x2 x3 x4]';
40 case 10
41 x0 = 0.0000000000000e+00;
42 x1 = 3.5902433622052e-01;
43 x2 = 1.1436659188355e+00;
44 x3 = 2.2144895894456e+00;
45 x4 = 3.3682742337736e+00;
46 x5 = 4.4309689056870e+00;
47 xb = [x0 x1 x2 x3 x4 x5]';
48 case 12
49 x0 = 0.0000000000000e+00;
50 x1 = 3.6098032343909e-01;
51 x2 = 1.1634317168086e+00;
52 x3 = 2.2975905356987e+00;
53 x4 = 3.6057529790929e+00;
54 x5 = 4.8918275675510e+00;
55 x6 = 6.0000000000000e+00;
56 xb = [x0 x1 x2 x3 x4 x5 x6]';
57 otherwise
58 error('Invalid operator order %d.',order);
59 end
60 end