Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1238:dea852e85b77 feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with refactorization of computing source indices
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 19 Nov 2019 16:06:03 -0800
parents f1806475498b 6e4cc4b66de0
children a0940c1db455 d1b201fe328e
comparison
equal deleted inserted replaced
1234:f1806475498b 1238:dea852e85b77
37 % Return zeros if x0 is outside grid 37 % Return zeros if x0 is outside grid
38 if x_s < x(1) || x_s > x(end) 38 if x_s < x(1) || x_s > x(end)
39 ret = zeros(size(x)); 39 ret = zeros(size(x));
40 return 40 return
41 else 41 else
42
43 fnorm = diag(H);
44 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for 42 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
45 S = []; 43 S = [];
46 M = []; 44 M = [];
47 45
48 % Get interior grid spacing 46 % Get interior grid spacing
49 middle = floor(m/2); 47 middle = floor(m/2);
50 h = x(middle+1) - x(middle); % Use middle point to allow for staggerd grids. 48 h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids.
51 49
52 % Find the indices that are within range of of the point source location 50 index = sourceIndices(x_s, x, tot_order, h);
53 ind_delta = find(tot_order*h/2 >= abs(x-x_s));
54 51
55 % Ensure that ind_delta is not too long 52 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
56 if length(ind_delta) == (tot_order + 2) 53 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
57 ind_delta = ind_delta(2:end-1); 54
58 elseif length(ind_delta) == (tot_order + 1) 55 quadrature = diag(H);
59 ind_delta = ind_delta(1:end-1); 56 quadrature_weights = quadrature(index)/h;
60 end
61
62 % Use first tot_order grid points
63 if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h
64 index=1:tot_order;
65 polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1));
66 x_0=(x_s-x(1))/(x(tot_order)-x(1));
67 norm=fnorm(1:tot_order)/h;
68
69 % Use last tot_order grid points
70 elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h
71 index = length(x)-tot_order+1:length(x);
72 polynomial = (x(end-tot_order+1:end)-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
73 norm = fnorm(end-tot_order+1:end)/h;
74 x_0 = (x_s-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
75
76 % Interior, compensate for round-off errors.
77 elseif length(ind_delta) < tot_order
78 if ind_delta(end)<m
79 ind_delta = [ind_delta; ind_delta(end)+1];
80 else
81 ind_delta = [ind_delta(1)-1; ind_delta];
82 end
83 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
84 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
85 norm = fnorm(ind_delta)/h;
86 index = ind_delta;
87
88 % Interior
89 else
90 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
91 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
92 norm = fnorm(ind_delta)/h;
93 index = ind_delta;
94 end
95 57
96 h_polynomial = polynomial(2)-polynomial(1); 58 h_polynomial = polynomial(2)-polynomial(1);
97 b = zeros(m_order+s_order,1); 59 b = zeros(tot_order,1);
98 60
99 for i = 1:m_order 61 for i = 1:m_order
100 b(i,1) = x_0^(i-1); 62 b(i,1) = x_0^(i-1);
101 end 63 end
102 64
103 for i = 1:(m_order+s_order) 65 for i = 1:tot_order
104 for j = 1:m_order 66 for j = 1:m_order
105 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i); 67 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
106 end 68 end
107 end 69 end
108 70
109 for i = 1:(m_order+s_order) 71 for i = 1:tot_order
110 for j = 1:s_order 72 for j = 1:s_order
111 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1); 73 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
112 end 74 end
113 end 75 end
114 76
118 ret = x*0; 80 ret = x*0;
119 ret(index) = d/h*h_polynomial; 81 ret(index) = d/h*h_polynomial;
120 end 82 end
121 83
122 end 84 end
85
86
87 function I = sourceIndices(x_s, x, tot_order, h)
88 % Find the indices that are within range of of the point source location
89 I = find(tot_order*h/2 >= abs(x-x_s));
90
91 if length(I) > tot_order
92 if length(I) == tot_order + 2
93 I = I(2:end-1);
94 elseif length(I) == tot_order + 1
95 I = I(1:end-1);
96 end
97 elseif length(I) < tot_order
98 if x_s < x(1) + ceil(tot_order/2)*h
99 I = 1:tot_order;
100 elseif x_s > x(end) - ceil(tot_order/2)*h
101 I = length(x)-tot_order+1:length(x);
102 else
103 if I(end) < length(x)
104 I = [I; I(end)+1];
105 else
106 I = [I(1)-1; I];
107 end
108 end
109 end
110 end