Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1232:52d774e69b1f feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up diracDiscr, remove obsolete tests.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 19 Nov 2019 13:54:41 -0800
parents 86ee5648e384
children f1806475498b 48c9a83260c8
comparison
equal deleted inserted replaced
1229:86ee5648e384 1232:52d774e69b1f
29 29
30 end 30 end
31 31
32 32
33 % Helper function for 1D delta functions 33 % Helper function for 1D delta functions
34 function ret = diracDiscr1D(x_0in , x , m_order, s_order, H) 34 function ret = diracDiscr1D(x_s , x , m_order, s_order, H)
35 35
36 m = length(x); 36 m = length(x);
37 37
38 % Return zeros if x0 is outside grid 38 % Return zeros if x0 is outside grid
39 if(x_0in < x(1) || x_0in > 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 44
45 fnorm = diag(H); 45 fnorm = diag(H);
46 eta = abs(x-x_0in); 46 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
47 tot = m_order+s_order; 47 S = [];
48 S = []; 48 M = [];
49 M = [];
50 49
51 % Get interior grid spacing 50 % Get interior grid spacing
52 middle = floor(m/2); 51 middle = floor(m/2);
53 h = x(middle+1) - x(middle); 52 h = x(middle+1) - x(middle);
54 53
55 poss = find(tot*h/2 >= eta); 54 % Find the indices that are within range of of the point source location
55 ind_delta = find(tot_order*h/2 >= abs(x-x_s));
56 56
57 % Ensure that poss is not too long 57 % Ensure that ind_delta is not too long
58 if length(poss) == (tot + 2) 58 if length(ind_delta) == (tot_order + 2)
59 poss = poss(2:end-1); 59 ind_delta = ind_delta(2:end-1);
60 elseif length(poss) == (tot + 1) 60 elseif length(ind_delta) == (tot_order + 1)
61 poss = poss(1:end-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 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
86 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
87 norm = fnorm(ind_delta)/h;
88 index = ind_delta;
89
90 % Interior
91 else
92 polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
93 x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
94 norm = fnorm(ind_delta)/h;
95 index = ind_delta;
96 end
97
98 h_polynomial = polynomial(2)-polynomial(1);
99 b = zeros(m_order+s_order,1);
100
101 for i = 1:m_order
102 b(i,1) = x_0^(i-1);
103 end
104
105 for i = 1:(m_order+s_order)
106 for j = 1:m_order
107 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i);
108 end
109 end
110
111 for i = 1:(m_order+s_order)
112 for j = 1:s_order
113 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
114 end
115 end
116
117 A = [M;S];
118
119 d = A\b;
120 ret = x*0;
121 ret(index) = d/h*h_polynomial;
62 end 122 end
63
64 % Use first tot grid points
65 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
66 index=1:tot;
67 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
68 x_0=(x_0in-x(1))/(x(tot)-x(1));
69 norm=fnorm(1:tot)/h;
70
71 % Use last tot grid points
72 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
73 index = length(x)-tot+1:length(x);
74 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
75 norm = fnorm(end-tot+1:end)/h;
76 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
77
78 % Interior, compensate for round-off errors.
79 elseif length(poss) < tot
80 if poss(end)<m
81 poss = [poss; poss(end)+1];
82 else
83 poss = [poss(1)-1; poss];
84 end
85 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
86 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
87 norm = fnorm(poss)/h;
88 index = poss;
89
90 % Interior
91 else
92 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
93 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
94 norm = fnorm(poss)/h;
95 index = poss;
96 end
97
98 h_pol = pol(2)-pol(1);
99 b = zeros(m_order+s_order,1);
100
101 for i = 1:m_order
102 b(i,1) = x_0^(i-1);
103 end
104
105 for i = 1:(m_order+s_order)
106 for j = 1:m_order
107 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
108 end
109 end
110
111 for i = 1:(m_order+s_order)
112 for j = 1:s_order
113 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
114 end
115 end
116
117 A = [M;S];
118
119 d = A\b;
120 ret = x*0;
121 ret(index) = d/h*h_pol;
122 end
123 123
124 end 124 end
125 125
126 126
127 127