Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 687:e8fc3aa1faf6 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename elastic scheme.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 09 Feb 2018 13:34:27 -0800
parents
children 60eb7f46d8d9
comparison
equal deleted inserted replaced
686:5ccf6aaf6d6b 687:e8fc3aa1faf6
1 classdef Elastic2dVariable < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
5 % opSet should be cell array of opSets, one per dimension. This
6 % is useful if we have periodic BC in one direction.
7
8 properties
9 m % Number of points in each direction, possibly a vector
10 h % Grid spacing
11
12 grid
13 dim
14
15 order % Order of accuracy for the approximation
16
17 % Diagonal matrices for varible coefficients
18 LAMBDA % Variable coefficient, related to dilation
19 MU % Shear modulus, variable coefficient
20 RHO, RHOi % Density, variable
21
22 D % Total operator
23 D1 % First derivatives
24
25 % Second derivatives
26 D2_lambda
27 D2_mu
28
29 % Traction operators used for BC
30 T_l, T_r
31 tau_l, tau_r
32
33 H, Hi % Inner products
34 phi % Borrowing constant for (d1 - e^T*D1) from R
35 gamma % Borrowing constant for d1 from M
36 H11 % First element of H
37 e_l, e_r
38 d1_l, d1_r % Normal derivatives at the boundary
39 E % E{i}^T picks out component i
40
41 H_boundary % Boundary inner products
42
43 % Kroneckered norms and coefficients
44 RHOi_kron
45 Hi_kron
46 end
47
48 methods
49
50 function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
51 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
52 default_arg('lambda_fun', @(x,y) 0*x+1);
53 default_arg('mu_fun', @(x,y) 0*x+1);
54 default_arg('rho_fun', @(x,y) 0*x+1);
55 dim = 2;
56
57 assert(isa(g, 'grid.Cartesian'))
58
59 lambda = grid.evalOn(g, lambda_fun);
60 mu = grid.evalOn(g, mu_fun);
61 rho = grid.evalOn(g, rho_fun);
62 m = g.size();
63 m_tot = g.N();
64
65 h = g.scaling();
66 lim = g.lim;
67
68 % 1D operators
69 ops = cell(dim,1);
70 for i = 1:dim
71 ops{i} = opSet{i}(m(i), lim{i}, order);
72 end
73
74 % Borrowing constants
75 for i = 1:dim
76 beta = ops{i}.borrowing.R.delta_D;
77 obj.H11{i} = ops{i}.borrowing.H11;
78 obj.phi{i} = beta/obj.H11{i};
79 obj.gamma{i} = ops{i}.borrowing.M.d1;
80 end
81
82 I = cell(dim,1);
83 D1 = cell(dim,1);
84 D2 = cell(dim,1);
85 H = cell(dim,1);
86 Hi = cell(dim,1);
87 e_l = cell(dim,1);
88 e_r = cell(dim,1);
89 d1_l = cell(dim,1);
90 d1_r = cell(dim,1);
91
92 for i = 1:dim
93 I{i} = speye(m(i));
94 D1{i} = ops{i}.D1;
95 D2{i} = ops{i}.D2;
96 H{i} = ops{i}.H;
97 Hi{i} = ops{i}.HI;
98 e_l{i} = ops{i}.e_l;
99 e_r{i} = ops{i}.e_r;
100 d1_l{i} = ops{i}.d1_l;
101 d1_r{i} = ops{i}.d1_r;
102 end
103
104 %====== Assemble full operators ========
105 LAMBDA = spdiag(lambda);
106 obj.LAMBDA = LAMBDA;
107 MU = spdiag(mu);
108 obj.MU = MU;
109 RHO = spdiag(rho);
110 obj.RHO = RHO;
111 obj.RHOi = inv(RHO);
112
113 obj.D1 = cell(dim,1);
114 obj.D2_lambda = cell(dim,1);
115 obj.D2_mu = cell(dim,1);
116 obj.e_l = cell(dim,1);
117 obj.e_r = cell(dim,1);
118 obj.d1_l = cell(dim,1);
119 obj.d1_r = cell(dim,1);
120
121 % D1
122 obj.D1{1} = kron(D1{1},I{2});
123 obj.D1{2} = kron(I{1},D1{2});
124
125 % Boundary operators
126 obj.e_l{1} = kron(e_l{1},I{2});
127 obj.e_l{2} = kron(I{1},e_l{2});
128 obj.e_r{1} = kron(e_r{1},I{2});
129 obj.e_r{2} = kron(I{1},e_r{2});
130
131 obj.d1_l{1} = kron(d1_l{1},I{2});
132 obj.d1_l{2} = kron(I{1},d1_l{2});
133 obj.d1_r{1} = kron(d1_r{1},I{2});
134 obj.d1_r{2} = kron(I{1},d1_r{2});
135
136 % D2
137 for i = 1:dim
138 obj.D2_lambda{i} = sparse(m_tot);
139 obj.D2_mu{i} = sparse(m_tot);
140 end
141 ind = grid.funcToMatrix(g, 1:m_tot);
142
143 for i = 1:m(2)
144 D_lambda = D2{1}(lambda(ind(:,i)));
145 D_mu = D2{1}(mu(ind(:,i)));
146
147 p = ind(:,i);
148 obj.D2_lambda{1}(p,p) = D_lambda;
149 obj.D2_mu{1}(p,p) = D_mu;
150 end
151
152 for i = 1:m(1)
153 D_lambda = D2{2}(lambda(ind(i,:)));
154 D_mu = D2{2}(mu(ind(i,:)));
155
156 p = ind(i,:);
157 obj.D2_lambda{2}(p,p) = D_lambda;
158 obj.D2_mu{2}(p,p) = D_mu;
159 end
160
161 % Quadratures
162 obj.H = kron(H{1},H{2});
163 obj.Hi = inv(obj.H);
164 obj.H_boundary = cell(dim,1);
165 obj.H_boundary{1} = H{2};
166 obj.H_boundary{2} = H{1};
167
168 % E{i}^T picks out component i.
169 E = cell(dim,1);
170 I = speye(m_tot,m_tot);
171 for i = 1:dim
172 e = sparse(dim,1);
173 e(i) = 1;
174 E{i} = kron(I,e);
175 end
176 obj.E = E;
177
178 % Differentiation matrix D (without SAT)
179 D2_lambda = obj.D2_lambda;
180 D2_mu = obj.D2_mu;
181 D1 = obj.D1;
182 D = sparse(dim*m_tot,dim*m_tot);
183 d = @kroneckerDelta; % Kronecker delta
184 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
185 for i = 1:dim
186 for j = 1:dim
187 D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +...
188 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
189 );
190 D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +...
191 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
192 D2_mu{j}*E{i}' ...
193 );
194 end
195 end
196 obj.D = D;
197 %=========================================%
198
199 % Numerical traction operators for BC.
200 % Because d1 =/= e0^T*D1, the numerical tractions are different
201 % at every boundary.
202 T_l = cell(dim,1);
203 T_r = cell(dim,1);
204 tau_l = cell(dim,1);
205 tau_r = cell(dim,1);
206 % tau^{j}_i = sum_k T^{j}_{ik} u_k
207
208 d1_l = obj.d1_l;
209 d1_r = obj.d1_r;
210 e_l = obj.e_l;
211 e_r = obj.e_r;
212 D1 = obj.D1;
213
214 % Loop over boundaries
215 for j = 1:dim
216 T_l{j} = cell(dim,dim);
217 T_r{j} = cell(dim,dim);
218 tau_l{j} = cell(dim,1);
219 tau_r{j} = cell(dim,1);
220
221 % Loop over components
222 for i = 1:dim
223 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
224 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
225 for k = 1:dim
226 T_l{j}{i,k} = ...
227 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
228 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
229 -d(i,k)*MU*e_l{j}*d1_l{j}';
230
231 T_r{j}{i,k} = ...
232 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
233 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
234 +d(i,k)*MU*e_r{j}*d1_r{j}';
235
236 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
237 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
238 end
239
240 end
241 end
242 obj.T_l = T_l;
243 obj.T_r = T_r;
244 obj.tau_l = tau_l;
245 obj.tau_r = tau_r;
246
247 % Kroneckered norms and coefficients
248 I_dim = speye(dim);
249 obj.RHOi_kron = kron(obj.RHOi, I_dim);
250 obj.Hi_kron = kron(obj.Hi, I_dim);
251
252 % Misc.
253 obj.m = m;
254 obj.h = h;
255 obj.order = order;
256 obj.grid = g;
257 obj.dim = dim;
258
259 end
260
261
262 % Closure functions return the operators applied to the own domain to close the boundary
263 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
264 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
265 % type is a cell array of strings specifying the type of boundary condition for each component.
266 % data is a function returning the data that should be applied at the boundary.
267 % neighbour_scheme is an instance of Scheme that should be interfaced to.
268 % neighbour_boundary is a string specifying which boundary to interface to.
269 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
270 default_arg('type',{'free','free'});
271 default_arg('parameter', []);
272
273 % j is the coordinate direction of the boundary
274 % nj: outward unit normal component.
275 % nj = -1 for west, south, bottom boundaries
276 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary);
278 switch nj
279 case 1
280 e = obj.e_r;
281 d = obj.d1_r;
282 tau = obj.tau_r{j};
283 T = obj.T_r{j};
284 case -1
285 e = obj.e_l;
286 d = obj.d1_l;
287 tau = obj.tau_l{j};
288 T = obj.T_l{j};
289 end
290
291 E = obj.E;
292 Hi = obj.Hi;
293 H_gamma = obj.H_boundary{j};
294 LAMBDA = obj.LAMBDA;
295 MU = obj.MU;
296 RHOi = obj.RHOi;
297
298 dim = obj.dim;
299 m_tot = obj.grid.N();
300
301 RHOi_kron = obj.RHOi_kron;
302 Hi_kron = obj.Hi_kron;
303
304 % Preallocate
305 closure = sparse(dim*m_tot, dim*m_tot);
306 penalty = cell(dim,1);
307 for k = 1:dim
308 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
309 end
310
311 % Loop over components that we (potentially) have different BC on
312 for k = 1:dim
313 switch type{k}
314
315 % Dirichlet boundary condition
316 case {'D','d','dirichlet','Dirichlet'}
317
318 tuning = 1.2;
319 phi = obj.phi{j};
320 h = obj.h(j);
321 h11 = obj.H11{j}*h;
322 gamma = obj.gamma{j};
323
324 a_lambda = dim/h11 + 1/(h11*phi);
325 a_mu_i = 2/(gamma*h);
326 a_mu_ij = 2/h11 + 1/(h11*phi);
327
328 d = @kroneckerDelta; % Kronecker delta
329 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
330 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
331 + d(i,j)* a_mu_i*MU ...
332 + db(i,j)*a_mu_ij*MU );
333
334 % Loop over components that Dirichlet penalties end up on
335 for i = 1:dim
336 C = T{k,i};
337 A = -d(i,k)*alpha(i,j);
338 B = A + C;
339 closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' );
340 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
341 end
342
343 % Free boundary condition
344 case {'F','f','Free','free','traction','Traction','t','T'}
345 closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} );
346 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
347
348 % Unknown boundary condition
349 otherwise
350 error('No such boundary condition: type = %s',type);
351 end
352 end
353 end
354
355 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
356 % u denotes the solution in the own domain
357 % v denotes the solution in the neighbour domain
358 tuning = 1.2;
359 % tuning = 20.2;
360 error('Interface not implemented');
361 end
362
363 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
364 function [j, nj] = get_boundary_number(obj, boundary)
365
366 switch boundary
367 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
368 j = 1;
369 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
370 j = 2;
371 otherwise
372 error('No such boundary: boundary = %s',boundary);
373 end
374
375 switch boundary
376 case {'w','W','west','West','s','S','south','South'}
377 nj = -1;
378 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
379 nj = 1;
380 end
381 end
382
383 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
384 function [return_op] = get_boundary_operator(obj, op, boundary)
385
386 switch boundary
387 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
388 j = 1;
389 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
390 j = 2;
391 otherwise
392 error('No such boundary: boundary = %s',boundary);
393 end
394
395 switch op
396 case 'e'
397 switch boundary
398 case {'w','W','west','West','s','S','south','South'}
399 return_op = obj.e_l{j};
400 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
401 return_op = obj.e_r{j};
402 end
403 case 'd'
404 switch boundary
405 case {'w','W','west','West','s','S','south','South'}
406 return_op = obj.d_l{j};
407 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
408 return_op = obj.d_r{j};
409 end
410 otherwise
411 error(['No such operator: operatr = ' op]);
412 end
413
414 end
415
416 function N = size(obj)
417 N = prod(obj.m);
418 end
419 end
420 end