Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 836:049e4c6fa630 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Sep 2018 14:46:39 -0700
parents
children 43a1c3ac07b1
comparison
equal deleted inserted replaced
811:9d4246ac94c0 836:049e4c6fa630
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
2
3 m = length(x);
4
5 % Return zeros if x0 is outside grid
6 if(x_0in < x(1) || x_0in > x(end) )
7
8 ret = zeros(size(x));
9
10 else
11
12 fnorm = diag(H);
13 eta = abs(x-x_0in);
14 tot = m_order+s_order;
15 S = [];
16 M = [];
17
18 % Get interior grid spacing
19 middle = floor(m/2);
20 h = x(middle+1) - x(middle);
21
22 poss = find(tot*h/2 >= eta);
23
24 % Ensure that poss is not too long
25 if length(poss) == (tot + 2)
26 poss = poss(2:end-1);
27 elseif length(poss) == (tot + 1)
28 poss = poss(1:end-1);
29 end
30
31 % Use first tot grid points
32 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
33 index=1:tot;
34 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
35 x_0=(x_0in-x(1))/(x(tot)-x(1));
36 norm=fnorm(1:tot)/h;
37
38 % Use last tot grid points
39 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
40 index = length(x)-tot+1:length(x);
41 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
42 norm = fnorm(end-tot+1:end)/h;
43 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
44
45 % Interior, compensate for round-off errors.
46 elseif length(poss) < tot
47 if poss(end)<m
48 poss = [poss; poss(end)+1];
49 else
50 poss = [poss(1)-1; poss];
51 end
52 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
53 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
54 norm = fnorm(poss)/h;
55 index = poss;
56
57 % Interior
58 else
59 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
60 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
61 norm = fnorm(poss)/h;
62 index = poss;
63 end
64
65 h_pol = pol(2)-pol(1);
66 b = zeros(m_order+s_order,1);
67
68 for i = 1:m_order
69 b(i,1) = x_0^(i-1);
70 end
71
72 for i = 1:(m_order+s_order)
73 for j = 1:m_order
74 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
75 end
76 end
77
78 for i = 1:(m_order+s_order)
79 for j = 1:s_order
80 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
81 end
82 end
83
84 A = [M;S];
85
86 d = A\b;
87 ret = x*0;
88 ret(index) = d/h*h_pol;
89 end
90
91 end
92
93
94
95
96
97