Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1236:3722c2579818 feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt to factor out a function for finding indecies of the source
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 20 Nov 2019 00:10:30 +0100
parents 48c9a83260c8
children 6e4cc4b66de0
comparison
equal deleted inserted replaced
1235:48c9a83260c8 1236:3722c2579818
39 if(x_s < x(1) || x_s > x(end) ) 39 if(x_s < x(1) || x_s > x(end) )
40 40
41 ret = zeros(size(x)); 41 ret = zeros(size(x));
42 42
43 else 43 else
44
45 fnorm = diag(H); 44 fnorm = diag(H);
46 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for 45 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
47 S = []; 46 S = [];
48 M = []; 47 M = [];
49 48
50 % Get interior grid spacing 49 % Get interior grid spacing
51 middle = floor(m/2); 50 middle = floor(m/2);
52 h = x(middle+1) - x(middle); 51 h = x(middle+1) - x(middle);
53 52
54 % Find the indices that are within range of of the point source location 53 index = sourceIndecies(x_s, x, tot_order, h)
55 ind_delta = find(tot_order*h/2 >= abs(x-x_s));
56 54
57 % Ensure that ind_delta is not too long 55 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
58 if length(ind_delta) == (tot_order + 2) 56 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
59 ind_delta = ind_delta(2:end-1); 57 norm = fnorm(index)/h;
60 elseif length(ind_delta) == (tot_order + 1)
61 ind_delta = ind_delta(1:end-1);
62 end
63
64 % Use first tot_order grid points
65 if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h;
66 index=1:tot_order;
67 polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1));
68 x_0=(x_s-x(1))/(x(tot_order)-x(1));
69 norm=fnorm(1:tot_order)/h;
70
71 % Use last tot_order grid points
72 elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h;
73 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));
75 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));
77
78 % Interior, compensate for round-off errors.
79 elseif length(ind_delta) < tot_order
80 if ind_delta(end)<m
81 ind_delta = [ind_delta; ind_delta(end)+1];
82 else
83 ind_delta = [ind_delta(1)-1; ind_delta];
84 end
85
86 index = ind_delta;
87 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
88 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
89 norm = fnorm(ind_delta)/h;
90
91 % Interior
92 else
93 index = ind_delta;
94 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
95 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
96 norm = fnorm(ind_delta)/h;
97 end
98 58
99 h_polynomial = polynomial(2)-polynomial(1); 59 h_polynomial = polynomial(2)-polynomial(1);
100 b = zeros(m_order+s_order,1); 60 b = zeros(tot_order,1);
101 61
102 for i = 1:m_order 62 for i = 1:m_order
103 b(i,1) = x_0^(i-1); 63 b(i,1) = x_0^(i-1);
104 end 64 end
105 65
106 for i = 1:(m_order+s_order) 66 for i = 1:tot_order
107 for j = 1:m_order 67 for j = 1:m_order
108 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i); 68 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i);
109 end 69 end
110 end 70 end
111 71
112 for i = 1:(m_order+s_order) 72 for i = 1:tot_order
113 for j = 1:s_order 73 for j = 1:s_order
114 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1); 74 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
115 end 75 end
116 end 76 end
117 77
123 end 83 end
124 84
125 end 85 end
126 86
127 87
88 function I = sourceIndecies(x_s, x, tot_order, h)
89 % Find the indices that are within range of of the point source location
90 I = find(tot_order*h/2 >= abs(x-x_s));
128 91
129 92 if length(I) > tot_order
130 93 if length(I) == tot_order + 2
131 94 I = I(2:end-1);
95 elseif length(I) == tot_order + 1
96 I = I(1:end-1);
97 end
98 elseif length(I) < tot_order
99 if x_s < x(1) + ceil(tot_order/2)*h;
100 I = 1:tot_order;
101 elseif x_s > x(end) - ceil(tot_order/2)*h;
102 I = length(x)-tot_order+1:length(x);
103 else
104 if I(end) < length(x)
105 I = [I; I(end)+1];
106 else
107 I = [I(1)-1; I];
108 end
109 end
110 end
111 end