Mercurial > repos > public > sbplib

annotate diracDiscr.m @ 1301:8978521b0f06 default

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix incorrect package name.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 08 Jul 2020 19:11:04 +0200
parents 25efceb0c392
children 60c875c18de3
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1246
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
1 % n-dimensional delta function
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
2 % g: cartesian grid
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
3 % x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
4 % m_order: Number of moment conditions
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
5 % s_order: Number of smoothness conditions
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
6 % H: cell array of 1D norm matrices
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
7 function d = diracDiscr(g, x_s, m_order, s_order, H)
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
8 assertType(g, 'grid.Cartesian');
1246
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1245
diff changeset
9
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
10 dim = g.d;
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 d_1D = cell(dim,1);
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
13 % Allow for non-cell input in 1D
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
14 if dim == 1
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
15 H = {H};
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
16 end
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
17 % Create 1D dirac discr for each coordinate dir.
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
18 for i = 1:dim
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
19 d_1D{i} = diracDiscr1D(x_s(i), g.x{i}, m_order, s_order, H{i});
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
20 end
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
22 d = d_1D{dim};
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
23 for i = dim-1: -1: 1
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
24 % Perform outer product, transpose, and then turn into column vector
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
25 d = (d_1D{i}*d')';
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
26 d = d(:);
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 end
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 end
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 % Helper function for 1D delta functions
1241
d1b201fe328e Turn row vectors into column vectors in diracDiscr comments
Martin Almquist <malmquist@stanford.edu>
parents: 1238
diff changeset
32 function ret = diracDiscr1D(x_s, x, m_order, s_order, H)
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33
1232
52d774e69b1f Clean up diracDiscr, remove obsolete tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1229
diff changeset
34 % Return zeros if x0 is outside grid
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
35 if x_s < x(1) || x_s > x(end)
1232
52d774e69b1f Clean up diracDiscr, remove obsolete tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1229
diff changeset
36 ret = zeros(size(x));
1234
f1806475498b - Pass grids to diracDiscr and adjust tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1232
diff changeset
37 return
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
38 end
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
39
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
40 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
41 S = [];
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
42 M = [];
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
44 % Get interior grid spacing
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
45 middle = floor(length(x)/2);
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
46 h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids.
1241
d1b201fe328e Turn row vectors into column vectors in diracDiscr comments
Martin Almquist <malmquist@stanford.edu>
parents: 1238
diff changeset
47
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
48 index = sourceIndices(x_s, x, tot_order, h);
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
50 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
51 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
53 quadrature = diag(H);
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
54 quadrature_weights = quadrature(index)/h;
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
56 h_polynomial = polynomial(2)-polynomial(1);
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
57 b = zeros(tot_order,1);
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
59 for i = 1:m_order
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
60 b(i,1) = x_0^(i-1);
1232
52d774e69b1f Clean up diracDiscr, remove obsolete tests.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1229
diff changeset
61 end
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
63 for i = 1:tot_order
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
64 for j = 1:m_order
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
65 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
66 end
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
67 end
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
68
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
69 for i = 1:tot_order
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
70 for j = 1:s_order
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
71 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
72 end
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
73 end
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
74
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
75 A = [M;S];
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
76
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
77 d = A\b;
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
78 ret = x*0;
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
79 ret(index) = d/h*h_polynomial;
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 end
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82
1238
dea852e85b77 Merge with refactorization of computing source indices
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1234 1237
diff changeset
83 function I = sourceIndices(x_s, x, tot_order, h)
1236
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
84 % Find the indices that are within range of of the point source location
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
85 I = find(tot_order*h/2 >= abs(x-x_s));
1229
86ee5648e384 Add multi-d dirac discretization with tests
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86
1236
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
87 if length(I) > tot_order
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
88 if length(I) == tot_order + 2
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
89 I = I(2:end-1);
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
90 elseif length(I) == tot_order + 1
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
91 I = I(1:end-1);
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
92 end
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
93 elseif length(I) < tot_order
1238
dea852e85b77 Merge with refactorization of computing source indices
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1234 1237
diff changeset
94 if x_s < x(1) + ceil(tot_order/2)*h
1236
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
95 I = 1:tot_order;
1238
dea852e85b77 Merge with refactorization of computing source indices
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1234 1237
diff changeset
96 elseif x_s > x(end) - ceil(tot_order/2)*h
1236
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
97 I = length(x)-tot_order+1:length(x);
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
98 else
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
99 if I(end) < length(x)
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
100 I = [I; I(end)+1];
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
101 else
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
102 I = [I(1)-1; I];
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
103 end
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
104 end
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
105 end
3722c2579818 Attempt to factor out a function for finding indecies of the source
Jonatan Werpers <jonatan@werpers.com>
parents: 1235
diff changeset
106 end