Mercurial > repos > public > sbplib
comparison +scheme/elasticShearVariable.m @ 663:b45ec2b28cc2 feature/poroelastic
First implementation of elastic shear operator with free boundary BC.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Thu, 14 Dec 2017 13:54:20 -0800 |
parents | |
children | 8e6dfd22fc59 |
comparison
equal
deleted
inserted
replaced
662:b189bc409cdb | 663:b45ec2b28cc2 |
---|---|
1 classdef elasticShearVariable < scheme.Scheme | |
2 properties | |
3 m % Number of points in each direction, possibly a vector | |
4 h % Grid spacing | |
5 | |
6 grid | |
7 | |
8 order % Order accuracy for the approximation | |
9 | |
10 a % Variable coefficient lambda of the operator | |
11 rho % Density | |
12 | |
13 D % Total operator | |
14 D1 % First derivatives | |
15 D2 % Second derivatives | |
16 H, Hi % Inner products | |
17 e_l, e_r | |
18 d_l, d_r % Normal derivatives at the boundary | |
19 H_boundary % Boundary inner products | |
20 end | |
21 | |
22 methods | |
23 % Implements the shear part of the elastic wave equation, i.e. | |
24 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i | |
25 % where a = lambda. | |
26 | |
27 function obj = elasticShearVariable(g ,order, a_fun, rho_fun, opSet) | |
28 default_arg('opSet',@sbp.D2Variable); | |
29 default_arg('a_fun', @(x,y) 0*x+1); | |
30 default_arg('rho_fun', @(x,y) 0*x+1); | |
31 dim = 2; | |
32 | |
33 assert(isa(g, 'grid.Cartesian')) | |
34 | |
35 a = grid.evalOn(g, a_fun); | |
36 rho = grid.evalOn(g, rho_fun); | |
37 m = g.size(); | |
38 m_tot = g.N(); | |
39 | |
40 h = g.scaling(); | |
41 | |
42 % 1D operators | |
43 ops = cell(dim,1); | |
44 for i = 1:dim | |
45 ops{i} = opSet(m(i), {0, 1}, order); | |
46 end | |
47 | |
48 I = cell(dim,1); | |
49 D1 = cell(dim,1); | |
50 D2 = cell(dim,1); | |
51 H = cell(dim,1); | |
52 Hi = cell(dim,1); | |
53 e_l = cell(dim,1); | |
54 e_r = cell(dim,1); | |
55 d1_l = cell(dim,1); | |
56 d1_r = cell(dim,1); | |
57 | |
58 for i = 1:dim | |
59 I{i} = speye{m(i)); | |
60 D1{i} = ops{i}.D1; | |
61 D2{i} = ops{i}.D2; | |
62 H{i} = ops{i}.H; | |
63 Hi{i} = ops{i}.HI; | |
64 e_l{i} = ops{i}.e_l; | |
65 e_r{i} = ops{i}.e_r; | |
66 d1_l{i} = ops{i}.d1_l; | |
67 d1_r{i} = ops{i}.d1_r; | |
68 end | |
69 | |
70 %====== Assemble full operators ======== | |
71 A = spdiag(a); | |
72 RHO = spdiag(rho); | |
73 | |
74 obj.D1 = cell(dim,1); | |
75 obj.D2 = cell(dim,1); | |
76 obj.e_l = cell(dim,1); | |
77 obj.e_r = cell(dim,1); | |
78 obj.d1_l = cell(dim,1); | |
79 obj.d1_r = cell(dim,1); | |
80 | |
81 % D1 | |
82 obj.D1{1} = kron{D1{1},I(2)}; | |
83 obj.D1{2} = kron{I{1},D1(2)}; | |
84 | |
85 % Boundary operators | |
86 obj.e_l{1} = kron{e_l{1},I(2)}; | |
87 obj.e_l{2} = kron{I{1},e_l(2)}; | |
88 obj.e_r{1} = kron{e_r{1},I(2)}; | |
89 obj.e_r{2} = kron{I{1},e_r(2)}; | |
90 | |
91 obj.d1_l{1} = kron{d1_l{1},I(2)}; | |
92 obj.d1_l{2} = kron{I{1},d1_l(2)}; | |
93 obj.d1_r{1} = kron{d1_r{1},I(2)}; | |
94 obj.d1_r{2} = kron{I{1},d1_r(2)}; | |
95 | |
96 % D2 | |
97 for i = 1:dim | |
98 obj.D2{i} = sparse(m_tot); | |
99 end | |
100 ind = grid.funcToMatrix(g, 1:m_tot); | |
101 | |
102 for i = 1:m(2) | |
103 D = D2{1}(a(ind(:,i))); | |
104 p = ind(:,i); | |
105 obj.D2{1}(p,p) = D; | |
106 end | |
107 | |
108 for i = 1:m(1) | |
109 D = D2{2}(a(ind(i,:))); | |
110 p = ind(i,:); | |
111 obj.D2{2}(p,p) = D; | |
112 end | |
113 | |
114 % Quadratures | |
115 obj.H = kron(H{1},H{2}); | |
116 obj.H_boundary = cell(dim,1); | |
117 obj.H_boundary{1} = H{2}; | |
118 obj.H_boundary{2} = H{1}; | |
119 | |
120 % Boundary coefficient matrices and quadratures | |
121 obj.A_boundary_l = cell(dim,1); | |
122 obj.A_boundary_r = cell(dim,1); | |
123 for i = 1:dim | |
124 obj.A_boundary_l{i} = e_l{i}'*A*e_l{i}; | |
125 obj.A_boundary_r{i} = e_r{i}'*A*e_r{i}; | |
126 end | |
127 | |
128 % E{i}^T picks out component i. | |
129 E = cell(dim,1); | |
130 I = speye{mtot,mtot}; | |
131 for i = 1:dim | |
132 e = sparse(dim,1); | |
133 e(i) = 1; | |
134 E{i} = kron(e,I) | |
135 end | |
136 obj.E = E; | |
137 | |
138 % Differentiation matrix D (without SAT) | |
139 D = 0; | |
140 d = @kroneckerDelta; % Kronecker delta | |
141 db = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta | |
142 for i = 1:dim | |
143 for j = 1:dim | |
144 D = D + E{i}*inv(rho)*( d(i,j)*D2{i}*E{j}' +... | |
145 db(i,j)*D1{j}*A*D1{i}*E{j}' + ... | |
146 D2{j}*E{i}' ... | |
147 ); | |
148 end | |
149 end | |
150 obj.D = D; | |
151 %=========================================% | |
152 | |
153 | |
154 | |
155 % Misc. | |
156 obj.m = m; | |
157 obj.h = h; | |
158 obj.order = order; | |
159 obj.grid = g; | |
160 | |
161 obj.a = a; | |
162 obj.b = b; | |
163 | |
164 % obj.gamm_u = h_u*ops_u.borrowing.M.d1; | |
165 % obj.gamm_v = h_v*ops_v.borrowing.M.d1; | |
166 end | |
167 | |
168 | |
169 % Closure functions return the operators applied to the own domain to close the boundary | |
170 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
171 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
172 % type is a string specifying the type of boundary condition if there are several. | |
173 % data is a function returning the data that should be applied at the boundary. | |
174 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
175 % neighbour_boundary is a string specifying which boundary to interface to. | |
176 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) | |
177 default_arg('type','free'); | |
178 default_arg('parameter', []); | |
179 | |
180 delta = @kroneckerDelta; % Kronecker delta | |
181 delta_b = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta | |
182 | |
183 % j is the coordinate direction of the boundary | |
184 % nj: outward unit normal component. | |
185 % nj = -1 for west, south, bottom boundaries | |
186 % nj = 1 for east, north, top boundaries | |
187 [j, nj] = obj.get_boundary_number(boundary); | |
188 switch nj | |
189 case 1 | |
190 e = obj.e_r; | |
191 d = obj.d_r; | |
192 case -1 | |
193 e = obj.e_l; | |
194 d = obj.d_l; | |
195 end | |
196 | |
197 E = obj.E; | |
198 Hi = obj.Hi; | |
199 H_gamma = obj.H_boundary{j}; | |
200 A = obj.A; | |
201 | |
202 switch type | |
203 % Dirichlet boundary condition | |
204 case {'D','d','dirichlet'} | |
205 error('Dirichlet not implemented') | |
206 tuning = 1.2; | |
207 % tuning = 20.2; | |
208 | |
209 b1 = gamm*obj.lambda./obj.a11.^2; | |
210 b2 = gamm*obj.lambda./obj.a22.^2; | |
211 | |
212 tau1 = tuning * spdiag(-1./b1 - 1./b2); | |
213 tau2 = 1; | |
214 | |
215 tau = (tau1*e + tau2*d)*H_b; | |
216 | |
217 closure = obj.a*obj.Hi*tau*e'; | |
218 penalty = -obj.a*obj.Hi*tau; | |
219 | |
220 | |
221 % Free boundary condition | |
222 case {'F','f','Free','free'} | |
223 closure = 0; | |
224 penalty = 0; | |
225 % Loop over components | |
226 for i = 1:3 | |
227 closure = closure + E{i}*(-sign)*Hi*e{j}*H_gamma*(... | |
228 e{j}'*A*e{j}*d{j}'*E{i}' + ... | |
229 delta(i,j)*e{j}'*A*e{i}*d{i}*E{j}' + ... | |
230 delta_b(i,j)*A*D1{i}*E{j}' ... | |
231 ); | |
232 penalty = penalty - E{i}*(-sign)*Hi*e{j}*H_gamma; | |
233 end | |
234 | |
235 % Unknown boundary condition | |
236 otherwise | |
237 error('No such boundary condition: type = %s',type); | |
238 end | |
239 end | |
240 | |
241 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
242 % u denotes the solution in the own domain | |
243 % v denotes the solution in the neighbour domain | |
244 tuning = 1.2; | |
245 % tuning = 20.2; | |
246 error('Interface not implemented'); | |
247 end | |
248 | |
249 % Ruturns the coordinate number and outward normal component for the boundary specified by the string boundary. | |
250 function [j, nj] = get_boundary_number(obj, boundary) | |
251 | |
252 switch boundary | |
253 case {'w','W','west','West', 'e', 'E', 'east', 'East'} | |
254 j = 1; | |
255 case {'s','S','south','South', 'n', 'N', 'north', 'North'} | |
256 j = 2; | |
257 otherwise | |
258 error('No such boundary: boundary = %s',boundary); | |
259 end | |
260 | |
261 switch boundary | |
262 case {'w','W','west','West','s','S','south','South'} | |
263 nj = -1; | |
264 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} | |
265 nj = 1; | |
266 end | |
267 end | |
268 | |
269 function N = size(obj) | |
270 N = prod(obj.m); | |
271 end | |
272 end | |
273 end |