Mercurial > repos > public > sbplib
diff diracDiscr.m @ 1129:b29892853daf feature/laplace_curvilinear_test
Refactor diracDiscr.m by moving the helper function diracDiscr1D to a separate file.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Tue, 21 May 2019 18:10:06 -0700 |
parents | 3a9262c045d0 |
children |
line wrap: on
line diff
--- a/diracDiscr.m Tue May 21 17:59:30 2019 -0700 +++ b/diracDiscr.m Tue May 21 18:10:06 2019 -0700 @@ -26,101 +26,6 @@ d = d(:); end end - -end - - -% Helper function for 1D delta functions -function ret = diracDiscr1D(x_0in , x , m_order, s_order, H) - -m = length(x); - -% Return zeros if x0 is outside grid -if(x_0in < x(1) || x_0in > x(end) ) - - ret = zeros(size(x)); - -else - - fnorm = diag(H); - eta = abs(x-x_0in); - tot = m_order+s_order; - S = []; - M = []; - - % Get interior grid spacing - middle = floor(m/2); - h = x(middle+1) - x(middle); - - poss = find(tot*h/2 >= eta); - - % Ensure that poss is not too long - if length(poss) == (tot + 2) - poss = poss(2:end-1); - elseif length(poss) == (tot + 1) - poss = poss(1:end-1); - end - - % Use first tot grid points - if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h; - index=1:tot; - pol=(x(1:tot)-x(1))/(x(tot)-x(1)); - x_0=(x_0in-x(1))/(x(tot)-x(1)); - norm=fnorm(1:tot)/h; - - % Use last tot grid points - elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h; - index = length(x)-tot+1:length(x); - pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1)); - norm = fnorm(end-tot+1:end)/h; - x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1)); - - % Interior, compensate for round-off errors. - elseif length(poss) < tot - if poss(end)<m - poss = [poss; poss(end)+1]; - else - poss = [poss(1)-1; poss]; - end - pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); - x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); - norm = fnorm(poss)/h; - index = poss; - - % Interior - else - pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); - x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); - norm = fnorm(poss)/h; - index = poss; - end - - h_pol = pol(2)-pol(1); - b = zeros(m_order+s_order,1); - - for i = 1:m_order - b(i,1) = x_0^(i-1); - end - - for i = 1:(m_order+s_order) - for j = 1:m_order - M(j,i) = pol(i)^(j-1)*h_pol*norm(i); - end - end - - for i = 1:(m_order+s_order) - for j = 1:s_order - S(j,i) = (-1)^(i-1)*pol(i)^(j-1); - end - end - - A = [M;S]; - - d = A\b; - ret = x*0; - ret(index) = d/h*h_pol; -end - end