Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1234:f1806475498b feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- Pass grids to diracDiscr and adjust tests. - Minor edits in diracDiscr
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 19 Nov 2019 15:16:08 -0800
parents 52d774e69b1f
children dea852e85b77
comparison
equal deleted inserted replaced
1233:57df0bf741dc 1234:f1806475498b
1 1
2 function d = diracDiscr(x_s, x, m_order, s_order, H) 2 function d = diracDiscr(g, x_s, m_order, s_order, H)
3 % n-dimensional delta function 3 % n-dimensional delta function
4 % g: cartesian grid
4 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. 5 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z].
5 % x: cell array of grid point column vectors for each dimension.
6 % m_order: Number of moment conditions 6 % m_order: Number of moment conditions
7 % s_order: Number of smoothness conditions 7 % s_order: Number of smoothness conditions
8 % H: cell array of 1D norm matrices 8 % H: cell array of 1D norm matrices
9 9 assertType(g, 'grid.Cartesian');
10 dim = length(x_s); 10 dim = g.d;
11 d_1D = cell(dim,1); 11 d_1D = cell(dim,1);
12 12
13 % If 1D, non-cell input is accepted 13 % Allow for non-cell input in 1D
14 if dim == 1 && ~iscell(x) 14 if dim == 1
15 d = diracDiscr1D(x_s, x, m_order, s_order, H); 15 H = {H};
16 end
17 % Create 1D dirac discr for each coordinate dir.
18 for i = 1:dim
19 d_1D{i} = diracDiscr1D(x_s(i), g.x{i}, m_order, s_order, H{i});
20 end
16 21
17 else 22 d = d_1D{dim};
18 for i = 1:dim 23 for i = dim-1: -1: 1
19 d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); 24 % Perform outer product, transpose, and then turn into column vector
20 end 25 d = (d_1D{i}*d')';
21 26 d = d(:);
22 d = d_1D{dim};
23 for i = dim-1: -1: 1
24 % Perform outer product, transpose, and then turn into column vector
25 d = (d_1D{i}*d')';
26 d = d(:);
27 end
28 end 27 end
29 28
30 end 29 end
31 30
32 31
34 function ret = diracDiscr1D(x_s , x , m_order, s_order, H) 33 function ret = diracDiscr1D(x_s , x , m_order, s_order, H)
35 34
36 m = length(x); 35 m = length(x);
37 36
38 % Return zeros if x0 is outside grid 37 % Return zeros if x0 is outside grid
39 if(x_s < x(1) || x_s > x(end) ) 38 if x_s < x(1) || x_s > x(end)
40
41 ret = zeros(size(x)); 39 ret = zeros(size(x));
42 40 return
43 else 41 else
44 42
45 fnorm = diag(H); 43 fnorm = diag(H);
46 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for 44 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
47 S = []; 45 S = [];
48 M = []; 46 M = [];
49 47
50 % Get interior grid spacing 48 % Get interior grid spacing
51 middle = floor(m/2); 49 middle = floor(m/2);
52 h = x(middle+1) - x(middle); 50 h = x(middle+1) - x(middle); % Use middle point to allow for staggerd grids.
53 51
54 % Find the indices that are within range of of the point source location 52 % Find the indices that are within range of of the point source location
55 ind_delta = find(tot_order*h/2 >= abs(x-x_s)); 53 ind_delta = find(tot_order*h/2 >= abs(x-x_s));
56 54
57 % Ensure that ind_delta is not too long 55 % Ensure that ind_delta is not too long
60 elseif length(ind_delta) == (tot_order + 1) 58 elseif length(ind_delta) == (tot_order + 1)
61 ind_delta = ind_delta(1:end-1); 59 ind_delta = ind_delta(1:end-1);
62 end 60 end
63 61
64 % Use first tot_order grid points 62 % Use first tot_order grid points
65 if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h; 63 if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h
66 index=1:tot_order; 64 index=1:tot_order;
67 polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1)); 65 polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1));
68 x_0=(x_s-x(1))/(x(tot_order)-x(1)); 66 x_0=(x_s-x(1))/(x(tot_order)-x(1));
69 norm=fnorm(1:tot_order)/h; 67 norm=fnorm(1:tot_order)/h;
70 68
71 % Use last tot_order grid points 69 % Use last tot_order grid points
72 elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h; 70 elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h
73 index = length(x)-tot_order+1:length(x); 71 index = length(x)-tot_order+1:length(x);
74 polynomial = (x(end-tot_order+1:end)-x(end-tot_order+1))/(x(end)-x(end-tot_order+1)); 72 polynomial = (x(end-tot_order+1:end)-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
75 norm = fnorm(end-tot_order+1:end)/h; 73 norm = fnorm(end-tot_order+1:end)/h;
76 x_0 = (x_s-x(end-tot_order+1))/(x(end)-x(end-tot_order+1)); 74 x_0 = (x_s-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
77 75
120 ret = x*0; 118 ret = x*0;
121 ret(index) = d/h*h_polynomial; 119 ret(index) = d/h*h_polynomial;
122 end 120 end
123 121
124 end 122 end
125
126
127
128
129
130