Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1250:8ec777fb473e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged in feature/dirac_discr (pull request #17) Add multi-d dirac discretization with tests Approved-by: Vidar Stiernström <vidar.stiernstrom@it.uu.se> Approved-by: Martin Almquist <malmquist@stanford.edu> Approved-by: Jonatan Werpers <jonatan.werpers@it.uu.se>
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 20 Nov 2019 22:24:06 +0000
parents 25efceb0c392
children 60c875c18de3
comparison
equal deleted inserted replaced
1220:4dc295afe473 1250:8ec777fb473e
1 % n-dimensional delta function
2 % g: cartesian grid
3 % x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
4 % m_order: Number of moment conditions
5 % s_order: Number of smoothness conditions
6 % H: cell array of 1D norm matrices
7 function d = diracDiscr(g, x_s, m_order, s_order, H)
8 assertType(g, 'grid.Cartesian');
9
10 dim = g.d;
11 d_1D = cell(dim,1);
12
13 % Allow for non-cell input in 1D
14 if dim == 1
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
21
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
29
30
31 % Helper function for 1D delta functions
32 function ret = diracDiscr1D(x_s, x, m_order, s_order, H)
33
34 % Return zeros if x0 is outside grid
35 if x_s < x(1) || x_s > x(end)
36 ret = zeros(size(x));
37 return
38 end
39
40 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
41 S = [];
42 M = [];
43
44 % Get interior grid spacing
45 middle = floor(length(x)/2);
46 h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids.
47
48 index = sourceIndices(x_s, x, tot_order, h);
49
50 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
51 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
52
53 quadrature = diag(H);
54 quadrature_weights = quadrature(index)/h;
55
56 h_polynomial = polynomial(2)-polynomial(1);
57 b = zeros(tot_order,1);
58
59 for i = 1:m_order
60 b(i,1) = x_0^(i-1);
61 end
62
63 for i = 1:tot_order
64 for j = 1:m_order
65 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
66 end
67 end
68
69 for i = 1:tot_order
70 for j = 1:s_order
71 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
72 end
73 end
74
75 A = [M;S];
76
77 d = A\b;
78 ret = x*0;
79 ret(index) = d/h*h_polynomial;
80 end
81
82
83 function I = sourceIndices(x_s, x, tot_order, h)
84 % Find the indices that are within range of of the point source location
85 I = find(tot_order*h/2 >= abs(x-x_s));
86
87 if length(I) > tot_order
88 if length(I) == tot_order + 2
89 I = I(2:end-1);
90 elseif length(I) == tot_order + 1
91 I = I(1:end-1);
92 end
93 elseif length(I) < tot_order
94 if x_s < x(1) + ceil(tot_order/2)*h
95 I = 1:tot_order;
96 elseif x_s > x(end) - ceil(tot_order/2)*h
97 I = length(x)-tot_order+1:length(x);
98 else
99 if I(end) < length(x)
100 I = [I; I(end)+1];
101 else
102 I = [I(1)-1; I];
103 end
104 end
105 end
106 end