Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1201:8f4e79aa32ba feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Thu, 05 Sep 2019 14:13:00 -0700
parents
children 31d7288d0653
comparison
equal deleted inserted replaced
1200:d9da4c1cdaa0 1201:8f4e79aa32ba
1 classdef Elastic2dVariableAnisotropic < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
5 % opSet should be cell array of opSets, one per dimension. This
6 % is useful if we have periodic BC in one direction.
7 % Assumes fully compatible operators
8
9 properties
10 m % Number of points in each direction, possibly a vector
11 h % Grid spacing
12
13 grid
14 dim
15
16 order % Order of accuracy for the approximation
17
18 % Diagonal matrices for variable coefficients
19 RHO, RHOi, RHOi_kron % Density
20 C % Elastic stiffness tensor
21
22 D % Total operator
23 D1 % First derivatives
24 D2 % Second derivatives
25
26 % Boundary operators in cell format, used for BC
27 T_w, T_e, T_s, T_n
28
29 % Traction operators
30 tau_w, tau_e, tau_s, tau_n % Return vector field
31 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
32 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
33
34 % Inner products
35 H, Hi, Hi_kron, H_1D
36
37 % Boundary inner products (for scalar field)
38 H_w, H_e, H_s, H_n
39
40 % Boundary restriction operators
41 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
42 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
43 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
44 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
45
46 % E{i}^T picks out component i
47 E
48
49 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
50 h11 % First entry in norm matrix
51
52 end
53
54 methods
55
56 % The coefficients can either be function handles or grid functions
57 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
58 function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag)
59 default_arg('rho', @(x,y) 0*x+1);
60 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
61 default_arg('optFlag', false);
62 dim = 2;
63
64 C_default = cell(dim,dim,dim,dim);
65 for i = 1:dim
66 for j = 1:dim
67 for k = 1:dim
68 for l = 1:dim
69 C_default{i,j,k,l} = @(x,y) 0*x + 1;
70 end
71 end
72 end
73 end
74 default_arg('C', C_default);
75 assert(isa(g, 'grid.Cartesian'))
76
77 if isa(rho, 'function_handle')
78 rho = grid.evalOn(g, rho);
79 end
80
81 C_mat = cell(dim,dim,dim,dim);
82 for i = 1:dim
83 for j = 1:dim
84 for k = 1:dim
85 for l = 1:dim
86 if isa(C{i,j,k,l}, 'function_handle')
87 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
88 end
89 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
90 end
91 end
92 end
93 end
94 obj.C = C_mat;
95
96 m = g.size();
97 m_tot = g.N();
98 lim = g.lim;
99 if isempty(lim)
100 x = g.x;
101 lim = cell(length(x),1);
102 for i = 1:length(x)
103 lim{i} = {min(x{i}), max(x{i})};
104 end
105 end
106
107 % 1D operators
108 ops = cell(dim,1);
109 h = zeros(dim,1);
110 for i = 1:dim
111 ops{i} = opSet{i}(m(i), lim{i}, order);
112 h(i) = ops{i}.h;
113 end
114
115 % Borrowing constants
116 for i = 1:dim
117 obj.h11{i} = h(i)*ops{i}.borrowing.H11;
118 end
119
120 I = cell(dim,1);
121 D1 = cell(dim,1);
122 D2 = cell(dim,1);
123 H = cell(dim,1);
124 Hi = cell(dim,1);
125 e_0 = cell(dim,1);
126 e_m = cell(dim,1);
127 d1_0 = cell(dim,1);
128 d1_m = cell(dim,1);
129
130 for i = 1:dim
131 I{i} = speye(m(i));
132 D1{i} = ops{i}.D1;
133 D2{i} = ops{i}.D2;
134 H{i} = ops{i}.H;
135 Hi{i} = ops{i}.HI;
136 e_0{i} = ops{i}.e_l;
137 e_m{i} = ops{i}.e_r;
138 d1_0{i} = ops{i}.d1_l;
139 d1_m{i} = ops{i}.d1_r;
140 end
141
142 %====== Assemble full operators ========
143 RHO = spdiag(rho);
144 obj.RHO = RHO;
145 obj.RHOi = inv(RHO);
146
147 obj.D1 = cell(dim,1);
148 obj.D2 = cell(dim,dim,dim);
149
150 % D1
151 obj.D1{1} = kron(D1{1},I{2});
152 obj.D1{2} = kron(I{1},D1{2});
153
154 % Boundary restriction operators
155 e_l = cell(dim,1);
156 e_r = cell(dim,1);
157 e_l{1} = kron(e_0{1}, I{2});
158 e_l{2} = kron(I{1}, e_0{2});
159 e_r{1} = kron(e_m{1}, I{2});
160 e_r{2} = kron(I{1}, e_m{2});
161
162 e_scalar_w = e_l{1};
163 e_scalar_e = e_r{1};
164 e_scalar_s = e_l{2};
165 e_scalar_n = e_r{2};
166
167 I_dim = speye(dim, dim);
168 e_w = kron(e_scalar_w, I_dim);
169 e_e = kron(e_scalar_e, I_dim);
170 e_s = kron(e_scalar_s, I_dim);
171 e_n = kron(e_scalar_n, I_dim);
172
173 % E{i}^T picks out component i.
174 E = cell(dim,1);
175 I = speye(m_tot,m_tot);
176 for i = 1:dim
177 e = sparse(dim,1);
178 e(i) = 1;
179 E{i} = kron(I,e);
180 end
181 obj.E = E;
182
183 e1_w = (e_scalar_w'*E{1}')';
184 e1_e = (e_scalar_e'*E{1}')';
185 e1_s = (e_scalar_s'*E{1}')';
186 e1_n = (e_scalar_n'*E{1}')';
187
188 e2_w = (e_scalar_w'*E{2}')';
189 e2_e = (e_scalar_e'*E{2}')';
190 e2_s = (e_scalar_s'*E{2}')';
191 e2_n = (e_scalar_n'*E{2}')';
192
193
194 % D2
195 for i = 1:dim
196 for k = 2:dim
197 for l = 2:dim
198 obj.D2{i,k,l} = sparse(m_tot);
199 end
200 end
201 end
202 ind = grid.funcToMatrix(g, 1:m_tot);
203
204 k = 1;
205 for r = 1:m(2)
206 p = ind(:,r);
207 for i = 1:dim
208 for l = 1:dim
209 coeff = C{i,k,k,l};
210 D_kk = D2{1}(coeff(p));
211 obj.D2{i,k,l}(p,p) = D_kk;
212 end
213 end
214 end
215
216 k = 2;
217 for r = 1:m(1)
218 p = ind(r,:);
219 for i = 1:dim
220 for l = 1:dim
221 coeff = C{i,k,k,l};
222 D_kk = D2{2}(coeff(p));
223 obj.D2{i,k,l}(p,p) = D_kk;
224 end
225 end
226 end
227
228 % Quadratures
229 obj.H = kron(H{1},H{2});
230 obj.Hi = inv(obj.H);
231 obj.H_w = H{2};
232 obj.H_e = H{2};
233 obj.H_s = H{1};
234 obj.H_n = H{1};
235 obj.H_1D = {H{1}, H{2}};
236
237 % Differentiation matrix D (without SAT)
238 D2 = obj.D2;
239 D1 = obj.D1;
240 D = sparse(dim*m_tot,dim*m_tot);
241 d = @kroneckerDelta; % Kronecker delta
242 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
243 for i = 1:dim
244 for j = 1:dim
245 for k = 1:dim
246 for l = 1:dim
247 D = D + E{i}*inv(RHO)*( d(j,k)*D2{i,k,l}*E{l}' +...
248 db(j,k)*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}' ...
249 );
250 end
251 end
252 end
253 end
254 obj.D = D;
255 %=========================================%'
256
257 % Numerical traction operators for BC.
258 %
259 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
260 %
261 T_l = cell(dim,1);
262 T_r = cell(dim,1);
263 tau_l = cell(dim,1);
264 tau_r = cell(dim,1);
265
266 D1 = obj.D1;
267
268 % Boundary j
269 for j = 1:dim
270 T_l{j} = cell(dim,dim);
271 T_r{j} = cell(dim,dim);
272 tau_l{j} = cell(dim,1);
273 tau_r{j} = cell(dim,1);
274
275 [~, n_l] = size(e_l{j});
276 [~, n_r] = size(e_r{j});
277
278 % Traction component i
279 for i = 1:dim
280 tau_l{j}{i} = sparse(dim*m_tot, n_l);
281 tau_r{j}{i} = sparse(dim*m_tot, n_r);
282
283 % Displacement component l
284 for l = 1:dim
285 T_l{j}{i,l} = sparse(m_tot, n_l);
286 T_r{j}{i,l} = sparse(m_tot, n_r);
287
288 % Derivative direction k
289 for k = 1:dim
290 T_l{j}{i,l} = T_l{j}{i,l} ...
291 - (e_l{j}'*C_mat{i,j,k,l}*D1{k})';
292 T_r{j}{i,l} = T_r{j}{i,l} ...
293 + (e_r{j}'*C_mat{i,j,k,l}*D1{k})';
294 end
295 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
296 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
297 end
298 end
299 end
300
301 % Traction tensors, T_ij
302 obj.T_w = T_l{1};
303 obj.T_e = T_r{1};
304 obj.T_s = T_l{2};
305 obj.T_n = T_r{2};
306
307 % Restriction operators
308 obj.e_w = e_w;
309 obj.e_e = e_e;
310 obj.e_s = e_s;
311 obj.e_n = e_n;
312
313 obj.e1_w = e1_w;
314 obj.e1_e = e1_e;
315 obj.e1_s = e1_s;
316 obj.e1_n = e1_n;
317
318 obj.e2_w = e2_w;
319 obj.e2_e = e2_e;
320 obj.e2_s = e2_s;
321 obj.e2_n = e2_n;
322
323 obj.e_scalar_w = e_scalar_w;
324 obj.e_scalar_e = e_scalar_e;
325 obj.e_scalar_s = e_scalar_s;
326 obj.e_scalar_n = e_scalar_n;
327
328 % First component of traction
329 obj.tau1_w = tau_l{1}{1};
330 obj.tau1_e = tau_r{1}{1};
331 obj.tau1_s = tau_l{2}{1};
332 obj.tau1_n = tau_r{2}{1};
333
334 % Second component of traction
335 obj.tau2_w = tau_l{1}{2};
336 obj.tau2_e = tau_r{1}{2};
337 obj.tau2_s = tau_l{2}{2};
338 obj.tau2_n = tau_r{2}{2};
339
340 % Traction vectors
341 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
342 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
343 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
344 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
345
346 % Kroneckered norms and coefficients
347 obj.RHOi_kron = kron(obj.RHOi, I_dim);
348 obj.Hi_kron = kron(obj.Hi, I_dim);
349
350 % Misc.
351 obj.m = m;
352 obj.h = h;
353 obj.order = order;
354 obj.grid = g;
355 obj.dim = dim;
356
357 end
358
359
360 % Closure functions return the operators applied to the own domain to close the boundary
361 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
362 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
363 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
364 % on the first component. Can also be e.g.
365 % {'normal', 'd'} or {'tangential', 't'} for conditions on
366 % tangential/normal component.
367 % data is a function returning the data that should be applied at the boundary.
368 % neighbour_scheme is an instance of Scheme that should be interfaced to.
369 % neighbour_boundary is a string specifying which boundary to interface to.
370 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
371 default_arg('tuning', 1.0);
372
373 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
374 comp = bc{1};
375 type = bc{2};
376 if ischar(comp)
377 comp = obj.getComponent(comp, boundary);
378 end
379
380 e = obj.getBoundaryOperatorForScalarField('e', boundary);
381 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
382 T = obj.getBoundaryTractionOperator(boundary);
383 h11 = obj.getBorrowing(boundary);
384 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
385 nu = obj.getNormal(boundary);
386
387 E = obj.E;
388 Hi = obj.Hi;
389 RHOi = obj.RHOi;
390 C = obj.C;
391
392 dim = obj.dim;
393 m_tot = obj.grid.N();
394
395 % Preallocate
396 [~, col] = size(tau);
397 closure = sparse(dim*m_tot, dim*m_tot);
398 penalty = sparse(dim*m_tot, col);
399
400 j = comp;
401 switch type
402
403 % Dirichlet boundary condition
404 % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components
405 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
406
407 % Loop over components that Dirichlet penalties end up on
408 % Y: symmetrizing part of penalty
409 % Z: symmetric part of penalty
410 % X = Y + Z.
411 for i = 1:dim
412 Y = T{j,i}';
413
414 Z = sparse(m_tot, m_tot);
415 for l = 1:dim
416 for k = 1:dim
417 Z = Z + nu(l)*C{i,l,k,j}*nu(k);
418 end
419 end
420 Z = -tuning*dim/h11*Z;
421 X = Z + e*Y;
422 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
423 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
424 end
425
426 % Free boundary condition
427 case {'F','f','Free','free','traction','Traction','t','T'}
428 closure = closure - E{j}*RHOi*Hi*e*H_gamma*tau';
429 penalty = penalty + E{j}*RHOi*Hi*e*H_gamma;
430
431 % Unknown boundary condition
432 otherwise
433 error('No such boundary condition: type = %s',type);
434 end
435 end
436
437 % type Struct that specifies the interface coupling.
438 % Fields:
439 % -- tuning: penalty strength, defaults to 1.2
440 % -- interpolation: type of interpolation, default 'none'
441 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
442
443 defaultType.tuning = 1.2;
444 defaultType.interpolation = 'none';
445 default_struct('type', defaultType);
446
447 switch type.interpolation
448 case {'none', ''}
449 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
450 case {'op','OP'}
451 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
452 otherwise
453 error('Unknown type of interpolation: %s ', type.interpolation);
454 end
455 end
456
457 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
458 tuning = type.tuning;
459
460 % u denotes the solution in the own domain
461 % v denotes the solution in the neighbour domain
462 % Operators without subscripts are from the own domain.
463
464 % Get boundary operators
465 e = obj.getBoundaryOperator('e', boundary);
466 tau = obj.getBoundaryOperator('tau', boundary);
467
468 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
469 tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
470
471 H_gamma = obj.getBoundaryQuadrature(boundary);
472
473 % Operators and quantities that correspond to the own domain only
474 Hi = obj.Hi_kron;
475 RHOi = obj.RHOi_kron;
476
477 % Penalty strength operators
478 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
479 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
480
481 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
482 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
483
484 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
485 penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
486
487 closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
488 penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
489
490 end
491
492 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
493 error('Non-conforming interfaces not implemented yet.');
494 end
495
496 % Returns the component number that is the tangential/normal component
497 % at the specified boundary
498 function comp = getComponent(obj, comp_str, boundary)
499 assertIsMember(comp_str, {'normal', 'tangential'});
500 assertIsMember(boundary, {'w', 'e', 's', 'n'});
501
502 switch boundary
503 case {'w', 'e'}
504 switch comp_str
505 case 'normal'
506 comp = 1;
507 case 'tangential'
508 comp = 2;
509 end
510 case {'s', 'n'}
511 switch comp_str
512 case 'normal'
513 comp = 2;
514 case 'tangential'
515 comp = 1;
516 end
517 end
518 end
519
520 % Returns h11 for the boundary specified by the string boundary.
521 % op -- string
522 function h11 = getBorrowing(obj, boundary)
523 assertIsMember(boundary, {'w', 'e', 's', 'n'})
524
525 switch boundary
526 case {'w','e'}
527 h11 = obj.h11{1};
528 case {'s', 'n'}
529 h11 = obj.h11{2};
530 end
531 end
532
533 % Returns the outward unit normal vector for the boundary specified by the string boundary.
534 function nu = getNormal(obj, boundary)
535 assertIsMember(boundary, {'w', 'e', 's', 'n'})
536
537 switch boundary
538 case 'w'
539 nu = [-1,0];
540 case 'e'
541 nu = [1,0];
542 case 's'
543 nu = [0,-1];
544 case 'n'
545 nu = [0,1];
546 end
547 end
548
549 % Returns the boundary operator op for the boundary specified by the string boundary.
550 % op -- string
551 function o = getBoundaryOperator(obj, op, boundary)
552 assertIsMember(boundary, {'w', 'e', 's', 'n'})
553 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
554
555 switch op
556 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
557 o = obj.([op, '_', boundary]);
558 end
559
560 end
561
562 % Returns the boundary operator op for the boundary specified by the string boundary.
563 % op -- string
564 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
565 assertIsMember(boundary, {'w', 'e', 's', 'n'})
566 assertIsMember(op, {'e'})
567
568 switch op
569
570 case 'e'
571 o = obj.(['e_scalar', '_', boundary]);
572 end
573
574 end
575
576 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
577 % Formula: tau_i = T_ij u_j
578 % op -- string
579 function T = getBoundaryTractionOperator(obj, boundary)
580 assertIsMember(boundary, {'w', 'e', 's', 'n'})
581
582 T = obj.(['T', '_', boundary]);
583 end
584
585 % Returns square boundary quadrature matrix, of dimension
586 % corresponding to the number of boundary unknowns
587 %
588 % boundary -- string
589 function H = getBoundaryQuadrature(obj, boundary)
590 assertIsMember(boundary, {'w', 'e', 's', 'n'})
591
592 H = obj.getBoundaryQuadratureForScalarField(boundary);
593 I_dim = speye(obj.dim, obj.dim);
594 H = kron(H, I_dim);
595 end
596
597 % Returns square boundary quadrature matrix, of dimension
598 % corresponding to the number of boundary grid points
599 %
600 % boundary -- string
601 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
602 assertIsMember(boundary, {'w', 'e', 's', 'n'})
603
604 H_b = obj.(['H_', boundary]);
605 end
606
607 function N = size(obj)
608 N = obj.dim*prod(obj.m);
609 end
610 end
611 end