comparison +sbp/+util/minimalBoundaryOptimizedGrid.m @ 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
children 38653d26225c
comparison
equal deleted inserted replaced
1284:5412fd72a7d8 1285:6b68f939d023
1 function [x,h] = minimalBoundaryOptimizedGrid(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
17 function xb = boundaryPoints(order)
18 switch order
19 case 4
20 x0 = 0.0000000000000e+00;
21 x1 = 7.7122987842562e-01;
22 xb = [x0 x1]';
23 case 6
24 x0 = 0.0000000000000e+00;
25 x1 = 4.0842950991998e-01;
26 x2 = 1.1968523189207e+00;
27 xb = [x0 x1 x2]';
28 case 8
29 x0 = 0.0000000000000e+00;
30 x1 = 4.9439570885261e-01;
31 x2 = 1.4051531374839e+00;
32 xb = [x0 x1 x2]';
33 case 10
34 x0 = 0.0000000000000e+00;
35 x1 = 5.8556160757529e-01;
36 x2 = 1.7473267488572e+00;
37 x3 = 3.0000000000000e+00;
38 xb = [x0 x1 x2 x3]';
39 case 12
40 x0 = 0.0000000000000e+00;
41 x1 = 4.6552112904489e-01;
42 x2 = 1.4647984306493e+00;
43 x3 = 2.7620429464763e+00;
44 x4 = 4.0000000000000e+00;
45 xb = [x0 x1 x2 x3 x4]';
46 otherwise
47 error('Invalid operator order %d.',order);
48 end
49 end