Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 837:43a1c3ac07b1 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make diracDiscr work for n dimensions. Add tests up to 3D.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Sep 2018 18:09:18 -0700
parents 049e4c6fa630
children c70131daaa6e 60c875c18de3
comparison
equal deleted inserted replaced
836:049e4c6fa630 837:43a1c3ac07b1
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H, h) 1
2 function d = diracDiscr(x_s, x, m_order, s_order, H)
3 % n-dimensional delta function
4 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z].
5 % x: cell array of grid point column vectors for each dimension.
6 % m_order: Number of moment conditions
7 % s_order: Number of smoothness conditions
8 % H: cell array of 1D norm matrices
9
10 dim = length(x_s);
11 d_1D = cell(dim,1);
12
13 % If 1D, non-cell input is accepted
14 if dim == 1 && ~iscell(x)
15 d = diracDiscr1D(x_s, x, m_order, s_order, H);
16
17 else
18 for i = 1:dim
19 d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i});
20 end
21
22 d = d_1D{dim};
23 for i = dim-1: -1: 1
24 % Perform outer product, transpose, and then turn into column vector
25 d = (d_1D{i}*d')';
26 d = d(:);
27 end
28 end
29
30 end
31
32
33 % Helper function for 1D delta functions
34 function ret = diracDiscr1D(x_0in , x , m_order, s_order, H)
2 35
3 m = length(x); 36 m = length(x);
4 37
5 % Return zeros if x0 is outside grid 38 % Return zeros if x0 is outside grid
6 if(x_0in < x(1) || x_0in > x(end) ) 39 if(x_0in < x(1) || x_0in > x(end) )
12 fnorm = diag(H); 45 fnorm = diag(H);
13 eta = abs(x-x_0in); 46 eta = abs(x-x_0in);
14 tot = m_order+s_order; 47 tot = m_order+s_order;
15 S = []; 48 S = [];
16 M = []; 49 M = [];
17 50
18 % Get interior grid spacing 51 % Get interior grid spacing
19 middle = floor(m/2); 52 middle = floor(m/2);
20 h = x(middle+1) - x(middle); 53 h = x(middle+1) - x(middle);
21 54
22 poss = find(tot*h/2 >= eta); 55 poss = find(tot*h/2 >= eta);
51 end 84 end
52 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); 85 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))); 86 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
54 norm = fnorm(poss)/h; 87 norm = fnorm(poss)/h;
55 index = poss; 88 index = poss;
56 89
57 % Interior 90 % Interior
58 else 91 else
59 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); 92 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))); 93 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
61 norm = fnorm(poss)/h; 94 norm = fnorm(poss)/h;
62 index = poss; 95 index = poss;
63 end 96 end