Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1071:92cb03e64ca4 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 20 Sep 2018 12:05:20 +0200
parents b374a8aa9246
children 386ef449df51 21394c78c72e
comparison
equal deleted inserted replaced
1070:f6b3af6febf3 1071:92cb03e64ca4
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
35 phi % Borrowing constant for (d1 - e^T*D1) from R
36 gamma % Borrowing constant for d1 from M
37 H11 % First element of H
38
39 % Borrowing from H, M, and R
40 thH
41 thM
42 thR
43
44 e_l, e_r
45 d1_l, d1_r % Normal derivatives at the boundary
46 E % E{i}^T picks out component i
47
48 H_boundary % Boundary inner products
49
50 % Kroneckered norms and coefficients
51 RHOi_kron
52 Hi_kron
53 end
54
55 methods
56
57 function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
58 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
59 default_arg('lambda_fun', @(x,y) 0*x+1);
60 default_arg('mu_fun', @(x,y) 0*x+1);
61 default_arg('rho_fun', @(x,y) 0*x+1);
62 dim = 2;
63
64 assert(isa(g, 'grid.Cartesian'))
65
66 lambda = grid.evalOn(g, lambda_fun);
67 mu = grid.evalOn(g, mu_fun);
68 rho = grid.evalOn(g, rho_fun);
69 m = g.size();
70 m_tot = g.N();
71
72 h = g.scaling();
73 lim = g.lim;
74 if isempty(lim)
75 x = g.x;
76 lim = cell(length(x),1);
77 for i = 1:length(x)
78 lim{i} = {min(x{i}), max(x{i})};
79 end
80 end
81
82 % 1D operators
83 ops = cell(dim,1);
84 for i = 1:dim
85 ops{i} = opSet{i}(m(i), lim{i}, order);
86 end
87
88 % Borrowing constants
89 for i = 1:dim
90 beta = ops{i}.borrowing.R.delta_D;
91 obj.H11{i} = ops{i}.borrowing.H11;
92 obj.phi{i} = beta/obj.H11{i};
93 obj.gamma{i} = ops{i}.borrowing.M.d1;
94
95 % Better names
96 obj.thR{i} = ops{i}.borrowing.R.delta_D;
97 obj.thM{i} = ops{i}.borrowing.M.d1;
98 obj.thH{i} = ops{i}.borrowing.H11;
99 end
100
101 I = cell(dim,1);
102 D1 = cell(dim,1);
103 D2 = cell(dim,1);
104 H = cell(dim,1);
105 Hi = cell(dim,1);
106 e_l = cell(dim,1);
107 e_r = cell(dim,1);
108 d1_l = cell(dim,1);
109 d1_r = cell(dim,1);
110
111 for i = 1:dim
112 I{i} = speye(m(i));
113 D1{i} = ops{i}.D1;
114 D2{i} = ops{i}.D2;
115 H{i} = ops{i}.H;
116 Hi{i} = ops{i}.HI;
117 e_l{i} = ops{i}.e_l;
118 e_r{i} = ops{i}.e_r;
119 d1_l{i} = ops{i}.d1_l;
120 d1_r{i} = ops{i}.d1_r;
121 end
122
123 %====== Assemble full operators ========
124 LAMBDA = spdiag(lambda);
125 obj.LAMBDA = LAMBDA;
126 MU = spdiag(mu);
127 obj.MU = MU;
128 RHO = spdiag(rho);
129 obj.RHO = RHO;
130 obj.RHOi = inv(RHO);
131
132 obj.D1 = cell(dim,1);
133 obj.D2_lambda = cell(dim,1);
134 obj.D2_mu = cell(dim,1);
135 obj.e_l = cell(dim,1);
136 obj.e_r = cell(dim,1);
137 obj.d1_l = cell(dim,1);
138 obj.d1_r = cell(dim,1);
139
140 % D1
141 obj.D1{1} = kron(D1{1},I{2});
142 obj.D1{2} = kron(I{1},D1{2});
143
144 % Boundary operators
145 obj.e_l{1} = kron(e_l{1},I{2});
146 obj.e_l{2} = kron(I{1},e_l{2});
147 obj.e_r{1} = kron(e_r{1},I{2});
148 obj.e_r{2} = kron(I{1},e_r{2});
149
150 obj.d1_l{1} = kron(d1_l{1},I{2});
151 obj.d1_l{2} = kron(I{1},d1_l{2});
152 obj.d1_r{1} = kron(d1_r{1},I{2});
153 obj.d1_r{2} = kron(I{1},d1_r{2});
154
155 % D2
156 for i = 1:dim
157 obj.D2_lambda{i} = sparse(m_tot);
158 obj.D2_mu{i} = sparse(m_tot);
159 end
160 ind = grid.funcToMatrix(g, 1:m_tot);
161
162 for i = 1:m(2)
163 D_lambda = D2{1}(lambda(ind(:,i)));
164 D_mu = D2{1}(mu(ind(:,i)));
165
166 p = ind(:,i);
167 obj.D2_lambda{1}(p,p) = D_lambda;
168 obj.D2_mu{1}(p,p) = D_mu;
169 end
170
171 for i = 1:m(1)
172 D_lambda = D2{2}(lambda(ind(i,:)));
173 D_mu = D2{2}(mu(ind(i,:)));
174
175 p = ind(i,:);
176 obj.D2_lambda{2}(p,p) = D_lambda;
177 obj.D2_mu{2}(p,p) = D_mu;
178 end
179
180 % Quadratures
181 obj.H = kron(H{1},H{2});
182 obj.Hi = inv(obj.H);
183 obj.H_boundary = cell(dim,1);
184 obj.H_boundary{1} = H{2};
185 obj.H_boundary{2} = H{1};
186
187 % E{i}^T picks out component i.
188 E = cell(dim,1);
189 I = speye(m_tot,m_tot);
190 for i = 1:dim
191 e = sparse(dim,1);
192 e(i) = 1;
193 E{i} = kron(I,e);
194 end
195 obj.E = E;
196
197 % Differentiation matrix D (without SAT)
198 D2_lambda = obj.D2_lambda;
199 D2_mu = obj.D2_mu;
200 D1 = obj.D1;
201 D = sparse(dim*m_tot,dim*m_tot);
202 d = @kroneckerDelta; % Kronecker delta
203 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
204 for i = 1:dim
205 for j = 1:dim
206 D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +...
207 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
208 );
209 D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +...
210 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
211 D2_mu{j}*E{i}' ...
212 );
213 end
214 end
215 obj.D = D;
216 %=========================================%
217
218 % Numerical traction operators for BC.
219 % Because d1 =/= e0^T*D1, the numerical tractions are different
220 % at every boundary.
221 T_l = cell(dim,1);
222 T_r = cell(dim,1);
223 tau_l = cell(dim,1);
224 tau_r = cell(dim,1);
225 % tau^{j}_i = sum_k T^{j}_{ik} u_k
226
227 d1_l = obj.d1_l;
228 d1_r = obj.d1_r;
229 e_l = obj.e_l;
230 e_r = obj.e_r;
231 D1 = obj.D1;
232
233 % Loop over boundaries
234 for j = 1:dim
235 T_l{j} = cell(dim,dim);
236 T_r{j} = cell(dim,dim);
237 tau_l{j} = cell(dim,1);
238 tau_r{j} = cell(dim,1);
239
240 % Loop over components
241 for i = 1:dim
242 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
243 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
244 for k = 1:dim
245 T_l{j}{i,k} = ...
246 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
247 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
248 -d(i,k)*MU*e_l{j}*d1_l{j}';
249
250 T_r{j}{i,k} = ...
251 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
252 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
253 +d(i,k)*MU*e_r{j}*d1_r{j}';
254
255 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
256 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
257 end
258
259 end
260 end
261 obj.T_l = T_l;
262 obj.T_r = T_r;
263 obj.tau_l = tau_l;
264 obj.tau_r = tau_r;
265
266 % Kroneckered norms and coefficients
267 I_dim = speye(dim);
268 obj.RHOi_kron = kron(obj.RHOi, I_dim);
269 obj.Hi_kron = kron(obj.Hi, I_dim);
270
271 % Misc.
272 obj.m = m;
273 obj.h = h;
274 obj.order = order;
275 obj.grid = g;
276 obj.dim = dim;
277
278 end
279
280
281 % Closure functions return the operators applied to the own domain to close the boundary
282 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
283 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
284 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
285 % on the first component.
286 % data is a function returning the data that should be applied at the boundary.
287 % neighbour_scheme is an instance of Scheme that should be interfaced to.
288 % neighbour_boundary is a string specifying which boundary to interface to.
289 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
290 default_arg('tuning', 1.2);
291
292 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
293 comp = bc{1};
294 type = bc{2};
295
296 % j is the coordinate direction of the boundary
297 j = obj.get_boundary_number(boundary);
298 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
299
300 E = obj.E;
301 Hi = obj.Hi;
302 LAMBDA = obj.LAMBDA;
303 MU = obj.MU;
304 RHOi = obj.RHOi;
305
306 dim = obj.dim;
307 m_tot = obj.grid.N();
308
309 % Preallocate
310 closure = sparse(dim*m_tot, dim*m_tot);
311 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
312
313 k = comp;
314 switch type
315
316 % Dirichlet boundary condition
317 case {'D','d','dirichlet','Dirichlet'}
318
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*H_gamma*(e'*E{k}' );
340 penalty = penalty - E{i}*RHOi*Hi*B'*e*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*H_gamma* (e'*tau{k} );
346 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
347
348 % Unknown boundary condition
349 otherwise
350 error('No such boundary condition: type = %s',type);
351 end
352 end
353
354 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
355 % u denotes the solution in the own domain
356 % v denotes the solution in the neighbour domain
357 % Operators without subscripts are from the own domain.
358 tuning = 1.2;
359
360 % j is the coordinate direction of the boundary
361 j = obj.get_boundary_number(boundary);
362 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
363
364 % Get boundary operators
365 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
366 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
367
368 % Operators and quantities that correspond to the own domain only
369 Hi = obj.Hi;
370 RHOi = obj.RHOi;
371 dim = obj.dim;
372
373 %--- Other operators ----
374 m_tot_u = obj.grid.N();
375 E = obj.E;
376 LAMBDA_u = obj.LAMBDA;
377 MU_u = obj.MU;
378 lambda_u = e'*LAMBDA_u*e;
379 mu_u = e'*MU_u*e;
380
381 m_tot_v = neighbour_scheme.grid.N();
382 E_v = neighbour_scheme.E;
383 LAMBDA_v = neighbour_scheme.LAMBDA;
384 MU_v = neighbour_scheme.MU;
385 lambda_v = e_v'*LAMBDA_v*e_v;
386 mu_v = e_v'*MU_v*e_v;
387 %-------------------------
388
389 % Borrowing constants
390 h_u = obj.h(j);
391 thR_u = obj.thR{j}*h_u;
392 thM_u = obj.thM{j}*h_u;
393 thH_u = obj.thH{j}*h_u;
394
395 h_v = neighbour_scheme.h(j_v);
396 thR_v = neighbour_scheme.thR{j_v}*h_v;
397 thH_v = neighbour_scheme.thH{j_v}*h_v;
398 thM_v = neighbour_scheme.thM{j_v}*h_v;
399
400 % alpha = penalty strength for normal component, beta for tangential
401 alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u);
402 alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v);
403 beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u);
404 beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v);
405 alpha = alpha_u + alpha_v;
406 beta = beta_u + beta_v;
407
408 d = @kroneckerDelta; % Kronecker delta
409 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
410 strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta);
411
412 % Preallocate
413 closure = sparse(dim*m_tot_u, dim*m_tot_u);
414 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
415
416 % Loop over components that penalties end up on
417 for i = 1:dim
418 closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}';
419 penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}';
420
421 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
422 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
423
424 % Loop over components that we have interface conditions on
425 for k = 1:dim
426 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
427 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
428 end
429 end
430 end
431
432 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
433 function [j, nj] = get_boundary_number(obj, boundary)
434
435 switch boundary
436 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
437 j = 1;
438 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
439 j = 2;
440 otherwise
441 error('No such boundary: boundary = %s',boundary);
442 end
443
444 switch boundary
445 case {'w','W','west','West','s','S','south','South'}
446 nj = -1;
447 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
448 nj = 1;
449 end
450 end
451
452 % Returns the boundary operator op for the boundary specified by the string boundary.
453 % op: may be a cell array of strings
454 function [varargout] = get_boundary_operator(obj, op, boundary)
455
456 switch boundary
457 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
458 j = 1;
459 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
460 j = 2;
461 otherwise
462 error('No such boundary: boundary = %s',boundary);
463 end
464
465 if ~iscell(op)
466 op = {op};
467 end
468
469 for i = 1:length(op)
470 switch op{i}
471 case 'e'
472 switch boundary
473 case {'w','W','west','West','s','S','south','South'}
474 varargout{i} = obj.e_l{j};
475 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
476 varargout{i} = obj.e_r{j};
477 end
478 case 'd'
479 switch boundary
480 case {'w','W','west','West','s','S','south','South'}
481 varargout{i} = obj.d1_l{j};
482 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
483 varargout{i} = obj.d1_r{j};
484 end
485 case 'H'
486 varargout{i} = obj.H_boundary{j};
487 case 'T'
488 switch boundary
489 case {'w','W','west','West','s','S','south','South'}
490 varargout{i} = obj.T_l{j};
491 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
492 varargout{i} = obj.T_r{j};
493 end
494 case 'tau'
495 switch boundary
496 case {'w','W','west','West','s','S','south','South'}
497 varargout{i} = obj.tau_l{j};
498 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
499 varargout{i} = obj.tau_r{j};
500 end
501 otherwise
502 error(['No such operator: operator = ' op{i}]);
503 end
504 end
505
506 end
507
508 function N = size(obj)
509 N = obj.dim*prod(obj.m);
510 end
511 end
512 end