Mercurial > repos > public > sbplib

annotate diracDiscr.m @ 1128:3a9262c045d0 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Copy diracDiscr.m from feature/poroelastic
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 May 2019 17:59:30 -0700
parents
children b29892853daf 52d774e69b1f
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1128
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2 function d = diracDiscr(x_s, x, m_order, s_order, H)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % n-dimensional delta function
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z].
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 % x: cell array of grid point column vectors for each dimension.
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % m_order: Number of moment conditions
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 % s_order: Number of smoothness conditions
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 % H: cell array of 1D norm matrices
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 dim = length(x_s);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 d_1D = cell(dim,1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 % If 1D, non-cell input is accepted
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 if dim == 1 && ~iscell(x)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 d = diracDiscr1D(x_s, x, m_order, s_order, H);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 else
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 for i = 1:dim
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i});
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 d = d_1D{dim};
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 for i = dim-1: -1: 1
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 % Perform outer product, transpose, and then turn into column vector
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 d = (d_1D{i}*d')';
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 d = d(:);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 % Helper function for 1D delta functions
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 function ret = diracDiscr1D(x_0in , x , m_order, s_order, H)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 m = length(x);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 % Return zeros if x0 is outside grid
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 if(x_0in < x(1) || x_0in > x(end) )
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 ret = zeros(size(x));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 else
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 fnorm = diag(H);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 eta = abs(x-x_0in);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 tot = m_order+s_order;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 S = [];
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 M = [];
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 % Get interior grid spacing
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 middle = floor(m/2);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 h = x(middle+1) - x(middle);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 poss = find(tot*h/2 >= eta);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 % Ensure that poss is not too long
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 if length(poss) == (tot + 2)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 poss = poss(2:end-1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 elseif length(poss) == (tot + 1)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 poss = poss(1:end-1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 % Use first tot grid points
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 index=1:tot;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 x_0=(x_0in-x(1))/(x(tot)-x(1));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 norm=fnorm(1:tot)/h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 % Use last tot grid points
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 index = length(x)-tot+1:length(x);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 norm = fnorm(end-tot+1:end)/h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 % Interior, compensate for round-off errors.
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 elseif length(poss) < tot
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 if poss(end)<m
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 poss = [poss; poss(end)+1];
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 else
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 poss = [poss(1)-1; poss];
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 norm = fnorm(poss)/h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 index = poss;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 % Interior
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 else
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 norm = fnorm(poss)/h;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 index = poss;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 h_pol = pol(2)-pol(1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 b = zeros(m_order+s_order,1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 for i = 1:m_order
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 b(i,1) = x_0^(i-1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 for i = 1:(m_order+s_order)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 for j = 1:m_order
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 for i = 1:(m_order+s_order)
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 for j = 1:s_order
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 A = [M;S];
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 d = A\b;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 ret = x*0;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 ret(index) = d/h*h_pol;
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 end
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129
3a9262c045d0 Copy diracDiscr.m from feature/poroelastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130