Mercurial > repos > public > sbplib

comparison +sbp/+grid/minimalBoundaryOptimizedGrid.m @ 1300:196123459178

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in feature/boundary_optimized_grids
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 08 Jul 2020 18:22:54 +0200
parents 0ffb5bfa65e4
children
comparison
equal deleted inserted replaced
1250:8ec777fb473e 1300:196123459178
1 % Computes the grid points x and grid spacing h used by the boundary optimized SBP operators
2 % with minimal number of non-equidistant boundary points, presented in
3 % 'Boundary optimized diagonal-norm SBP operators - Mattsson, Almquist, van der Weide 2018'.
4 %
5 % lim - cell array with domain limits
6 % N - Number of grid points
7 % order - order of accuracy of sbp operator.
8 function [x,h] = minimalBoundaryOptimizedGrid(lim,N,order)
9 assert(iscell(lim) && numel(lim) == 2,'The limit should be a cell array with 2 elements.');
10 L = lim{2} - lim{1};
11 assert(L>0,'Limits must be given in increasing order.');
12 %%%% Non-equidistant grid points %%%%%
13 xb = boundaryPoints(order);
14 m = length(xb)-1; % Number of non-equidistant points
15 assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1));
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17
18 %%%% Compute h %%%%%%%%%%
19 h = L/(2*xb(end) + N-1-2*m);
20 %%%%%%%%%%%%%%%%%%%%%%%%%
21
22 %%%% Define grid %%%%%%%%
23 x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ];
24 x = x + lim{1};
25 %%%%%%%%%%%%%%%%%%%%%%%%%
26 end
27
28 function xb = boundaryPoints(order)
29 switch order
30 case 4
31 x0 = 0.0000000000000e+00;
32 x1 = 7.7122987842562e-01;
33 xb = [x0 x1]';
34 case 6
35 x0 = 0.0000000000000e+00;
36 x1 = 4.0842950991998e-01;
37 x2 = 1.1968523189207e+00;
38 xb = [x0 x1 x2]';
39 case 8
40 x0 = 0.0000000000000e+00;
41 x1 = 4.9439570885261e-01;
42 x2 = 1.4051531374839e+00;
43 xb = [x0 x1 x2]';
44 case 10
45 x0 = 0.0000000000000e+00;
46 x1 = 5.8556160757529e-01;
47 x2 = 1.7473267488572e+00;
48 x3 = 3.0000000000000e+00;
49 xb = [x0 x1 x2 x3]';
50 case 12
51 x0 = 0.0000000000000e+00;
52 x1 = 4.6552112904489e-01;
53 x2 = 1.4647984306493e+00;
54 x3 = 2.7620429464763e+00;
55 x4 = 4.0000000000000e+00;
56 xb = [x0 x1 x2 x3 x4]';
57 otherwise
58 error('Invalid operator order %d.',order);
59 end
60 end