Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1264:066fdfaa2411 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add staggered Lebedev scheme for Cartesian anisotropic. Traction BC work!
author Martin Almquist <malmquist@stanford.edu>
date Wed, 29 Apr 2020 21:05:01 -0700
parents
children 6696b216b982
comparison
equal deleted inserted replaced
1263:117bc4542b61 1264:066fdfaa2411
1 classdef Elastic2dStaggeredAnisotropic < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
5 % Uses a staggered Lebedev grid
6 % The solution (displacement) is stored on g_u
7 % Stresses (and hance tractions) appear on g_s
8 % Density is evaluated on g_u
9 % The stiffness tensor is evaluated on g_s
10
11 properties
12 m % Number of points in each direction, possibly a vector
13 h % Grid spacing
14
15 grid
16 dim
17 nGrids
18 N % Total number of unknowns stored (2 displacement components on 2 grids)
19
20 order % Order of accuracy for the approximation
21
22 % Diagonal matrices for variable coefficients
23 RHO % Density
24 C % Elastic stiffness tensor
25
26 D % Total operator
27
28 % Boundary operators in cell format, used for BC
29 % T_w, T_e, T_s, T_n
30
31 % Traction operators
32 tau_w, tau_e, tau_s, tau_n % Return scalar field
33
34 % Inner products
35 H, H_u, H_s
36 % , Hi, Hi_kron, H_1D
37
38 % Boundary inner products (for scalar field)
39 H_w_u, H_e_u, H_s_u, H_n_u
40 H_w_s, H_e_s, H_s_s, H_n_s
41
42 % Boundary restriction operators
43 e_w_u, e_e_u, e_s_u, e_n_u % Act on scalar field, return scalar field at boundary
44 e_w_s, e_e_s, e_s_s, e_n_s % Act on scalar field, return scalar field at boundary
45
46 % U{i}^T picks out component i
47 U
48
49 % G{i}^T picks out displacement grid i
50 G
51
52 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
53 h11 % First entry in norm matrix
54
55 end
56
57 methods
58
59 % The coefficients can either be function handles or grid functions
60 function obj = Elastic2dStaggeredAnisotropic(g, order, rho, C)
61 default_arg('rho', @(x,y) 0*x+1);
62 dim = 2;
63 nGrids = 2;
64
65 C_default = cell(dim,dim,dim,dim);
66 for i = 1:dim
67 for j = 1:dim
68 for k = 1:dim
69 for l = 1:dim
70 C_default{i,j,k,l} = @(x,y) 0*x + 1;
71 end
72 end
73 end
74 end
75 default_arg('C', C_default);
76 assert(isa(g, 'grid.Staggered'))
77
78 g_u = g.gridGroups{1};
79 g_s = g.gridGroups{2};
80
81 m_u = {g_u{1}.size(), g_u{2}.size()};
82 m_s = {g_s{1}.size(), g_s{2}.size()};
83
84 if isa(rho, 'function_handle')
85 rho_vec = cell(nGrids, 1);
86 for i = 1:nGrids
87 rho_vec{i} = grid.evalOn(g_u{i}, rho);
88 end
89 rho = rho_vec;
90 end
91 for i = 1:nGrids
92 RHO{i} = spdiag(rho{i});
93 end
94 obj.RHO = RHO;
95
96 C_mat = cell(nGrids, 1);
97 for a = 1:nGrids
98 C_mat{a} = cell(dim,dim,dim,dim);
99 end
100 for a = 1:nGrids
101 for i = 1:dim
102 for j = 1:dim
103 for k = 1:dim
104 for l = 1:dim
105 C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
106 end
107 end
108 end
109 end
110 end
111 C = C_mat;
112 obj.C = C;
113
114 % Reference m for primal grid
115 m = g_u{1}.size();
116 X = g_u{1}.points;
117 lim = cell(dim, 1);
118 for i = 1:dim
119 lim{i} = {min(X(:,i)), max(X(:,i))};
120 end
121
122 % 1D operators
123 ops = cell(dim,1);
124 D1p = cell(dim, 1);
125 D1d = cell(dim, 1);
126 mp = cell(dim, 1);
127 md = cell(dim, 1);
128 Ip = cell(dim, 1);
129 Id = cell(dim, 1);
130 Hp = cell(dim, 1);
131 Hd = cell(dim, 1);
132
133 opSet = @sbp.D1StaggeredUpwind;
134 for i = 1:dim
135 ops{i} = opSet(m(i), lim{i}, order);
136 D1p{i} = ops{i}.D1_dual;
137 D1d{i} = ops{i}.D1_primal;
138 mp{i} = length(ops{i}.x_primal);
139 md{i} = length(ops{i}.x_dual);
140 Ip{i} = speye(mp{i}, mp{i});
141 Id{i} = speye(md{i}, md{i});
142 Hp{i} = ops{i}.H_primal;
143 Hd{i} = ops{i}.H_dual;
144 ep_l{i} = ops{i}.e_primal_l;
145 ep_r{i} = ops{i}.e_primal_r;
146 ed_l{i} = ops{i}.e_dual_l;
147 ed_r{i} = ops{i}.e_dual_r;
148 end
149
150 % Borrowing constants
151 for i = 1:dim
152 obj.h11{i} = ops{i}.H_dual(1,1);
153 end
154
155 %---- Grid layout -------
156 % gu1 = xp o yp;
157 % gu2 = xd o yd;
158 % gs1 = xd o yp;
159 % gs2 = xp o yd;
160 %------------------------
161
162 % Quadratures
163 obj.H_u = cell(nGrids, 1);
164 obj.H_s = cell(nGrids, 1);
165 obj.H_u{1} = kron(Hp{1}, Hp{2});
166 obj.H_u{2} = kron(Hd{1}, Hd{2});
167 obj.H_s{1} = kron(Hd{1}, Hp{2});
168 obj.H_s{2} = kron(Hp{1}, Hd{2});
169
170 obj.H_w_s = cell(nGrids, 1);
171 obj.H_e_s = cell(nGrids, 1);
172 obj.H_s_s = cell(nGrids, 1);
173 obj.H_n_s = cell(nGrids, 1);
174
175 obj.H_w_s{1} = Hp{2};
176 obj.H_w_s{2} = Hd{2};
177 obj.H_e_s{1} = Hp{2};
178 obj.H_e_s{2} = Hd{2};
179
180 obj.H_s_s{1} = Hd{1};
181 obj.H_s_s{2} = Hp{1};
182 obj.H_n_s{1} = Hd{1};
183 obj.H_n_s{2} = Hp{1};
184
185 % Boundary restriction ops
186 e_w_u = cell(nGrids, 1);
187 e_s_u = cell(nGrids, 1);
188 e_e_u = cell(nGrids, 1);
189 e_n_u = cell(nGrids, 1);
190
191 e_w_s = cell(nGrids, 1);
192 e_s_s = cell(nGrids, 1);
193 e_e_s = cell(nGrids, 1);
194 e_n_s = cell(nGrids, 1);
195
196 e_w_u{1} = kron(ep_l{1}, Ip{2});
197 e_e_u{1} = kron(ep_r{1}, Ip{2});
198 e_s_u{1} = kron(Ip{1}, ep_l{2});
199 e_n_u{1} = kron(Ip{1}, ep_r{2});
200
201 e_w_u{2} = kron(ed_l{1}, Id{2});
202 e_e_u{2} = kron(ed_r{1}, Id{2});
203 e_s_u{2} = kron(Id{1}, ed_l{2});
204 e_n_u{2} = kron(Id{1}, ed_r{2});
205
206 e_w_s{1} = kron(ed_l{1}, Ip{2});
207 e_e_s{1} = kron(ed_r{1}, Ip{2});
208 e_s_s{1} = kron(Id{1}, ep_l{2});
209 e_n_s{1} = kron(Id{1}, ep_r{2});
210
211 e_w_s{2} = kron(ep_l{1}, Id{2});
212 e_e_s{2} = kron(ep_r{1}, Id{2});
213 e_s_s{2} = kron(Ip{1}, ed_l{2});
214 e_n_s{2} = kron(Ip{1}, ed_r{2});
215
216 obj.e_w_u = e_w_u;
217 obj.e_e_u = e_e_u;
218 obj.e_s_u = e_s_u;
219 obj.e_n_u = e_n_u;
220
221 obj.e_w_s = e_w_s;
222 obj.e_e_s = e_e_s;
223 obj.e_s_s = e_s_s;
224 obj.e_n_s = e_n_s;
225
226
227 % D1_u2s{a, b}{i} approximates ddi and
228 %takes from u grid number b to s grid number a
229 % Some of D1_x2y{a, b} are 0.
230 D1_u2s = cell(nGrids, nGrids);
231 D1_s2u = cell(nGrids, nGrids);
232
233 N_u = cell(nGrids, 1);
234 N_s = cell(nGrids, 1);
235 for a = 1:nGrids
236 N_u{a} = g_u{a}.N();
237 N_s{a} = g_s{a}.N();
238 end
239
240 %---- Grid layout -------
241 % gu1 = xp o yp;
242 % gu2 = xd o yd;
243 % gs1 = xd o yp;
244 % gs2 = xp o yd;
245 %------------------------
246
247 D1_u2s{1,1}{1} = kron(D1p{1}, Ip{2});
248 D1_s2u{1,1}{1} = kron(D1d{1}, Ip{2});
249
250 D1_u2s{1,2}{2} = kron(Id{1}, D1d{2});
251 D1_u2s{2,1}{2} = kron(Ip{1}, D1p{2});
252
253 D1_s2u{1,2}{2} = kron(Ip{1}, D1d{2});
254 D1_s2u{2,1}{2} = kron(Id{1}, D1p{2});
255
256 D1_u2s{2,2}{1} = kron(D1d{1}, Id{2});
257 D1_s2u{2,2}{1} = kron(D1p{1}, Id{2});
258
259 D1_u2s{1,1}{2} = sparse(N_s{1}, N_u{1});
260 D1_s2u{1,1}{2} = sparse(N_u{1}, N_s{1});
261
262 D1_u2s{2,2}{2} = sparse(N_s{2}, N_u{2});
263 D1_s2u{2,2}{2} = sparse(N_u{2}, N_s{2});
264
265 D1_u2s{1,2}{1} = sparse(N_s{1}, N_u{2});
266 D1_s2u{1,2}{1} = sparse(N_u{1}, N_s{2});
267
268 D1_u2s{2,1}{1} = sparse(N_s{2}, N_u{1});
269 D1_s2u{2,1}{1} = sparse(N_u{2}, N_s{1});
270
271
272 %---- Combine grids and components -----
273
274 % U{a}{i}^T picks out u component i on grid a
275 U = cell(nGrids, 1);
276 for a = 1:2
277 U{a} = cell(dim, 1);
278 I = speye(N_u{a}, N_u{a});
279 for i = 1:dim
280 E = sparse(dim,1);
281 E(i) = 1;
282 U{a}{i} = kron(I, E);
283 end
284 end
285 obj.U = U;
286
287 % Order grids
288 % u1, u2
289 Iu1 = speye(dim*N_u{1}, dim*N_u{1});
290 Iu2 = speye(dim*N_u{2}, dim*N_u{2});
291
292 Gu1 = cell2mat( {Iu1; sparse(dim*N_u{2}, dim*N_u{1})} );
293 Gu2 = cell2mat( {sparse(dim*N_u{1}, dim*N_u{2}); Iu2} );
294 G = {Gu1; Gu2};
295 obj.G = G;
296
297 obj.H = G{1}*(U{1}{1}*obj.H_u{1}*U{1}{1}' + U{1}{2}*obj.H_u{1}*U{1}{2}')*G{1}'...
298 + G{2}*(U{2}{1}*obj.H_u{2}*U{2}{1}' + U{2}{2}*obj.H_u{2}*U{2}{2}')*G{2}';
299
300 % e1_w = (e_scalar_w'*E{1}')';
301 % e1_e = (e_scalar_e'*E{1}')';
302 % e1_s = (e_scalar_s'*E{1}')';
303 % e1_n = (e_scalar_n'*E{1}')';
304
305 % e2_w = (e_scalar_w'*E{2}')';
306 % e2_e = (e_scalar_e'*E{2}')';
307 % e2_s = (e_scalar_s'*E{2}')';
308 % e2_n = (e_scalar_n'*E{2}')';
309
310 % Differentiation matrix D (without SAT)
311 N = dim*(N_u{1} + N_u{2});
312 D = sparse(N, N);
313 for a = 1:nGrids
314 for b = 1:nGrids
315 for c = 1:nGrids
316 for i = 1:dim
317 for j = 1:dim
318 for k = 1:dim
319 for l = 1:dim
320 D = D + (G{a}*U{a}{j})*(RHO{a}\(D1_s2u{a,b}{i}*C{b}{i,j,k,l}*D1_u2s{b,c}{k}*U{c}{l}'*G{c}'));
321 end
322 end
323 end
324 end
325 end
326 end
327 end
328 obj.D = D;
329 clear D;
330 obj.N = N;
331 %=========================================%'
332
333 % Numerical traction operators for BC.
334 %
335 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
336 %
337
338 n_w = obj.getNormal('w');
339 n_e = obj.getNormal('e');
340 n_s = obj.getNormal('s');
341 n_n = obj.getNormal('n');
342
343 tau_w = cell(nGrids, 1);
344 tau_e = cell(nGrids, 1);
345 tau_s = cell(nGrids, 1);
346 tau_n = cell(nGrids, 1);
347
348 for b = 1:nGrids
349 tau_w{b} = cell(dim, 1);
350 tau_e{b} = cell(dim, 1);
351 tau_s{b} = cell(dim, 1);
352 tau_n{b} = cell(dim, 1);
353
354 for j = 1:dim
355 tau_w{b}{j} = sparse(N, m_s{b}(2));
356 tau_e{b}{j} = sparse(N, m_s{b}(2));
357 tau_s{b}{j} = sparse(N, m_s{b}(1));
358 tau_n{b}{j} = sparse(N, m_s{b}(1));
359 end
360
361 for c = 1:nGrids
362 for i = 1:dim
363 for j = 1:dim
364 for k = 1:dim
365 for l = 1:dim
366 sigma_b_ij = C{b}{i,j,k,l}*D1_u2s{b,c}{k}*U{c}{l}'*G{c}';
367
368 tau_w{b}{j} = tau_w{b}{j} + (e_w_s{b}'*n_w(i)*sigma_b_ij)';
369 tau_e{b}{j} = tau_e{b}{j} + (e_e_s{b}'*n_e(i)*sigma_b_ij)';
370 tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)';
371 tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)';
372 end
373 end
374 end
375 end
376 end
377 end
378
379 obj.tau_w = tau_w;
380 obj.tau_e = tau_e;
381 obj.tau_s = tau_s;
382 obj.tau_n = tau_n;
383
384 % D1 = obj.D1;
385
386 % Traction tensors, T_ij
387 % obj.T_w = T_l{1};
388 % obj.T_e = T_r{1};
389 % obj.T_s = T_l{2};
390 % obj.T_n = T_r{2};
391
392 % Restriction operators
393 % obj.e_w = e_w;
394 % obj.e_e = e_e;
395 % obj.e_s = e_s;
396 % obj.e_n = e_n;
397
398 % obj.e1_w = e1_w;
399 % obj.e1_e = e1_e;
400 % obj.e1_s = e1_s;
401 % obj.e1_n = e1_n;
402
403 % obj.e2_w = e2_w;
404 % obj.e2_e = e2_e;
405 % obj.e2_s = e2_s;
406 % obj.e2_n = e2_n;
407
408 % obj.e_scalar_w = e_scalar_w;
409 % obj.e_scalar_e = e_scalar_e;
410 % obj.e_scalar_s = e_scalar_s;
411 % obj.e_scalar_n = e_scalar_n;
412
413 % % First component of traction
414 % obj.tau1_w = tau_l{1}{1};
415 % obj.tau1_e = tau_r{1}{1};
416 % obj.tau1_s = tau_l{2}{1};
417 % obj.tau1_n = tau_r{2}{1};
418
419 % % Second component of traction
420 % obj.tau2_w = tau_l{1}{2};
421 % obj.tau2_e = tau_r{1}{2};
422 % obj.tau2_s = tau_l{2}{2};
423 % obj.tau2_n = tau_r{2}{2};
424
425 % % Traction vectors
426 % obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
427 % obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
428 % obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
429 % obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
430
431 % Misc.
432 obj.m = m;
433 obj.h = [];
434 obj.order = order;
435 obj.grid = g;
436 obj.dim = dim;
437 obj.nGrids = nGrids;
438
439 end
440
441
442 % Closure functions return the operators applied to the own domain to close the boundary
443 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
444 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
445 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
446 % on the first component. Can also be e.g.
447 % {'normal', 'd'} or {'tangential', 't'} for conditions on
448 % tangential/normal component.
449 % data is a function returning the data that should be applied at the boundary.
450 % neighbour_scheme is an instance of Scheme that should be interfaced to.
451 % neighbour_boundary is a string specifying which boundary to interface to.
452
453 % For displacement bc:
454 % bc = {comp, 'd', dComps},
455 % where
456 % dComps = vector of components with displacement BC. Default: 1:dim.
457 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
458 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
459 default_arg('tuning', 1.0);
460
461 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
462 comp = bc{1};
463 type = bc{2};
464 if ischar(comp)
465 comp = obj.getComponent(comp, boundary);
466 end
467
468 e = obj.getBoundaryOperatorForScalarField('e_u', boundary);
469 tau = obj.getBoundaryOperator('tau', boundary);
470 % T = obj.getBoundaryTractionOperator(boundary);
471 % h11 = obj.getBorrowing(boundary);
472 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
473 nu = obj.getNormal(boundary);
474
475 U = obj.U;
476 G = obj.G;
477 H = obj.H_u;
478 RHO = obj.RHO;
479 C = obj.C;
480
481 switch boundary
482 case {'s', 'n'}
483 e = rot90(e,2);
484 G = rot90(G,2);
485 U = rot90(U,2);
486 H = rot90(H,2);
487 RHO = rot90(RHO,2);
488 end
489
490 dim = obj.dim;
491 nGrids = obj.nGrids;
492
493 m_tot = obj.N;
494
495 % Preallocate
496 [~, col] = size(tau{1}{1});
497 closure = sparse(m_tot, m_tot);
498 % penalty = sparse(m_tot, col);
499 penalty = cell(nGrids, 1);
500
501 j = comp;
502 switch type
503
504 % Dirichlet boundary condition
505 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
506 error('not implemented');
507
508 if numel(bc) >= 3
509 dComps = bc{3};
510 else
511 dComps = 1:dim;
512 end
513
514 % Loops over components that Dirichlet penalties end up on
515 % Y: symmetrizing part of penalty
516 % Z: symmetric part of penalty
517 % X = Y + Z.
518
519 % Nonsymmetric part goes on all components to
520 % yield traction in discrete energy rate
521 for i = 1:dim
522 Y = T{j,i}';
523 X = e*Y;
524 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
525 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
526 end
527
528 % Symmetric part only required on components with displacement BC.
529 % (Otherwise it's not symmetric.)
530 for i = dComps
531 Z = sparse(m_tot, m_tot);
532 for l = 1:dim
533 for k = 1:dim
534 Z = Z + nu(l)*C{l,i,k,j}*nu(k);
535 end
536 end
537 Z = -tuning*dim/h11*Z;
538 X = Z;
539 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
540 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
541 end
542
543 % Free boundary condition
544 case {'F','f','Free','free','traction','Traction','t','T'}
545 for a = 1:nGrids
546 closure = closure - G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a}*tau{a}{j}')));
547 penalty{a} = G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a})));
548 end
549
550 % Unknown boundary condition
551 otherwise
552 error('No such boundary condition: type = %s',type);
553 end
554 end
555
556 % type Struct that specifies the interface coupling.
557 % Fields:
558 % -- tuning: penalty strength, defaults to 1.0
559 % -- interpolation: type of interpolation, default 'none'
560 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
561
562 defaultType.tuning = 1.0;
563 defaultType.interpolation = 'none';
564 default_struct('type', defaultType);
565
566 switch type.interpolation
567 case {'none', ''}
568 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
569 case {'op','OP'}
570 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
571 otherwise
572 error('Unknown type of interpolation: %s ', type.interpolation);
573 end
574 end
575
576 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
577 tuning = type.tuning;
578
579 % u denotes the solution in the own domain
580 % v denotes the solution in the neighbour domain
581
582 u = obj;
583 v = neighbour_scheme;
584
585 % Operators, u side
586 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
587 tau_u = u.getBoundaryOperator('tau', boundary);
588 h11_u = u.getBorrowing(boundary);
589 nu_u = u.getNormal(boundary);
590
591 E_u = u.E;
592 C_u = u.C;
593 m_tot_u = u.grid.N();
594
595 % Operators, v side
596 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
597 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
598 h11_v = v.getBorrowing(neighbour_boundary);
599 nu_v = v.getNormal(neighbour_boundary);
600
601 E_v = v.E;
602 C_v = v.C;
603 m_tot_v = v.grid.N();
604
605 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
606 flipFlag = false;
607 e_v_flip = e_v;
608 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
609 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ...
610 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ...
611 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ...
612 (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ...
613 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ...
614 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ...
615 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
616
617 flipFlag = true;
618 e_v_flip = fliplr(e_v);
619
620 t1 = tau_v(:,1:2:end-1);
621 t2 = tau_v(:,2:2:end);
622
623 t1 = fliplr(t1);
624 t2 = fliplr(t2);
625
626 tau_v(:,1:2:end-1) = t1;
627 tau_v(:,2:2:end) = t2;
628 end
629
630 % Operators that are only required for own domain
631 Hi = u.Hi_kron;
632 RHOi = u.RHOi_kron;
633 e_kron = u.getBoundaryOperator('e', boundary);
634 T_u = u.getBoundaryTractionOperator(boundary);
635
636 % Shared operators
637 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
638 H_gamma_kron = u.getBoundaryQuadrature(boundary);
639 dim = u.dim;
640
641 % Preallocate
642 [~, m_int] = size(H_gamma);
643 closure = sparse(dim*m_tot_u, dim*m_tot_u);
644 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
645
646 % ---- Continuity of displacement ------
647
648 % Y: symmetrizing part of penalty
649 % Z: symmetric part of penalty
650 % X = Y + Z.
651
652 % Loop over components to couple across interface
653 for j = 1:dim
654
655 % Loop over components that penalties end up on
656 for i = 1:dim
657 Y = 1/2*T_u{j,i}';
658 Z_u = sparse(m_int, m_int);
659 Z_v = sparse(m_int, m_int);
660 for l = 1:dim
661 for k = 1:dim
662 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
663 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
664 end
665 end
666
667 if flipFlag
668 Z_v = rot90(Z_v,2);
669 end
670
671 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
672 X = Y + Z*e_u';
673 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
674 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
675
676 end
677 end
678
679 % ---- Continuity of traction ------
680 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
681 penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
682
683 % ---- Multiply by inverse of density x quadraure ----
684 closure = RHOi*Hi*closure;
685 penalty = RHOi*Hi*penalty;
686
687 end
688
689 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
690 error('Non-conforming interfaces not implemented yet.');
691 end
692
693 % Returns the component number that is the tangential/normal component
694 % at the specified boundary
695 function comp = getComponent(obj, comp_str, boundary)
696 assertIsMember(comp_str, {'normal', 'tangential'});
697 assertIsMember(boundary, {'w', 'e', 's', 'n'});
698
699 switch boundary
700 case {'w', 'e'}
701 switch comp_str
702 case 'normal'
703 comp = 1;
704 case 'tangential'
705 comp = 2;
706 end
707 case {'s', 'n'}
708 switch comp_str
709 case 'normal'
710 comp = 2;
711 case 'tangential'
712 comp = 1;
713 end
714 end
715 end
716
717 % Returns h11 for the boundary specified by the string boundary.
718 % op -- string
719 function h11 = getBorrowing(obj, boundary)
720 assertIsMember(boundary, {'w', 'e', 's', 'n'})
721
722 switch boundary
723 case {'w','e'}
724 h11 = obj.h11{1};
725 case {'s', 'n'}
726 h11 = obj.h11{2};
727 end
728 end
729
730 % Returns the outward unit normal vector for the boundary specified by the string boundary.
731 function nu = getNormal(obj, boundary)
732 assertIsMember(boundary, {'w', 'e', 's', 'n'})
733
734 switch boundary
735 case 'w'
736 nu = [-1,0];
737 case 'e'
738 nu = [1,0];
739 case 's'
740 nu = [0,-1];
741 case 'n'
742 nu = [0,1];
743 end
744 end
745
746 % Returns the boundary operator op for the boundary specified by the string boundary.
747 % op -- string
748 function o = getBoundaryOperator(obj, op, boundary)
749 assertIsMember(boundary, {'w', 'e', 's', 'n'})
750 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
751
752 switch op
753 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
754 o = obj.([op, '_', boundary]);
755 end
756
757 end
758
759 % Returns the boundary operator op for the boundary specified by the string boundary.
760 % op -- string
761 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
762 assertIsMember(boundary, {'w', 'e', 's', 'n'})
763 assertIsMember(op, {'e_u', 'e_s'})
764
765 switch op
766 case 'e_u'
767 o = obj.(['e_', boundary, '_u']);
768 case 'e_s'
769 o = obj.(['e_', boundary, '_s']);
770 end
771
772 end
773
774 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
775 % Formula: tau_i = T_ij u_j
776 % op -- string
777 function T = getBoundaryTractionOperator(obj, boundary)
778 assertIsMember(boundary, {'w', 'e', 's', 'n'})
779
780 T = obj.(['T', '_', boundary]);
781 end
782
783 % Returns square boundary quadrature matrix, of dimension
784 % corresponding to the number of boundary unknowns
785 %
786 % boundary -- string
787 function H = getBoundaryQuadrature(obj, boundary)
788 assertIsMember(boundary, {'w', 'e', 's', 'n'})
789
790 H = obj.getBoundaryQuadratureForScalarField(boundary);
791 I_dim = speye(obj.dim, obj.dim);
792 H = kron(H, I_dim);
793 end
794
795 % Returns square boundary quadrature matrix, of dimension
796 % corresponding to the number of boundary grid points
797 %
798 % boundary -- string
799 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
800 assertIsMember(boundary, {'w', 'e', 's', 'n'})
801
802 H_b = obj.(['H_', boundary, '_s']);
803 end
804
805 function N = size(obj)
806 N = obj.dim*prod(obj.m);
807 end
808 end
809 end