Mercurial > repos > public > sbplib

comparison +sbp/+util/accurateBoundaryOptimizedGrid.m @ 1285:6b68f939d023 feature/boundary_optimized_grids

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