Mercurial > repos > public > sbplib

annotate diracDiscr.m @ 1336:0666629aa183 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add methods for creating grids with different grid point distributions for each coordinate direction, and also supports constructing periodic grids
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:26:16 +0200
parents 60c875c18de3
children
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;
837
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
11 d_1D = cell(dim,1);
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
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});
837
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
20 end
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
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
837
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
28 end
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
29
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
30
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
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)
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 S = [];
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 M = [];
837
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
diff changeset
43
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 % Get interior grid spacing
1245
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.
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
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);
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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)));
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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;
837
43a1c3ac07b1 Make diracDiscr work for n dimensions. Add tests up to 3D.
Martin Almquist <malmquist@stanford.edu>
parents: 836
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);
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 for i = 1:m_order
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 b(i,1) = x_0^(i-1);
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 for j = 1:m_order
1245
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);
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
69 for i = 1:tot_order
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 for j = 1:s_order
1245
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);
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 A = [M;S];
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 d = A\b;
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 ret = x*0;
1245
0a1c64d3c717 Avoid indentation of whole function
Jonatan Werpers <jonatan@werpers.com>
parents: 1242
diff changeset
79 ret(index) = d/h*h_polynomial;
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 end
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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));
836
049e4c6fa630 Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
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