Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 868:57760d7088ad bcSetupExperiment

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