Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropicIncompatible.m @ 1345:14f44e81e1e3 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme for Elastic2dVariable, not fully compatible.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 30 Apr 2024 13:33:39 +0200
parents
children 464e6f65c2c6
comparison
equal deleted inserted replaced
1322:412b8ceafbc6 1345:14f44e81e1e3
1 classdef Elastic2dVariableAnisotropicIncompatible < 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 theta_R % Borrowing from R
52
53 nBP % Number of boundary points, related to borrowing.
54 e_shifted_scalar_w, e_shifted_scalar_e, e_shifted_scalar_s, e_shifted_scalar_n; % Act on scalar field, return scalar field
55
56 end
57
58 methods
59
60 % The coefficients can either be function handles or grid functions
61 function obj = Elastic2dVariableAnisotropicIncompatible(g, order, rho, C, opSet, optFlag)
62 default_arg('rho', @(x,y) 0*x+1);
63 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
64 dim = 2;
65
66 C_default = cell(dim,dim,dim,dim);
67 for i = 1:dim
68 for j = 1:dim
69 for k = 1:dim
70 for l = 1:dim
71 C_default{i,j,k,l} = @(x,y) 0*x + 1;
72 end
73 end
74 end
75 end
76 default_arg('C', C_default);
77 assert(isa(g, 'grid.Cartesian'))
78
79 if isa(rho, 'function_handle')
80 rho = grid.evalOn(g, rho);
81 end
82
83 C_mat = cell(dim,dim,dim,dim);
84 for i = 1:dim
85 for j = 1:dim
86 for k = 1:dim
87 for l = 1:dim
88 if isa(C{i,j,k,l}, 'function_handle')
89 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
90 end
91 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
92 end
93 end
94 end
95 end
96 obj.C = C_mat;
97
98 m = g.size();
99 m_tot = g.N();
100 lim = g.lim;
101 if isempty(lim)
102 x = g.x;
103 lim = cell(length(x),1);
104 for i = 1:length(x)
105 lim{i} = {min(x{i}), max(x{i})};
106 end
107 end
108
109 % 1D operators
110 ops = cell(dim,1);
111 h = zeros(dim,1);
112 for i = 1:dim
113 ops{i} = opSet{i}(m(i), lim{i}, order);
114 h(i) = ops{i}.h;
115 end
116
117 % Borrowing constants
118 for i = 1:dim
119 obj.h11{i} = h(i)*ops{i}.borrowing.H11;
120 obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D;
121 end
122
123 switch order
124 case 2
125 width = 3;
126 nBP = 2;
127 case 4
128 width = 5;
129 nBP = 6;
130 case 6
131 width = 7;
132 nBP = 9;
133 end
134
135 I = cell(dim,1);
136 D1 = cell(dim,1);
137 D2 = cell(dim,1);
138 H = cell(dim,1);
139 Hi = cell(dim,1);
140 e_0 = cell(dim,1);
141 e_m = cell(dim,1);
142 d1_0 = cell(dim,1);
143 d1_m = cell(dim,1);
144
145 for i = 1:dim
146 I{i} = speye(m(i));
147 D1{i} = ops{i}.D1;
148 D2{i} = ops{i}.D2;
149 H{i} = ops{i}.H;
150 Hi{i} = ops{i}.HI;
151 e_0{i} = ops{i}.e_l;
152 e_m{i} = ops{i}.e_r;
153 d1_0{i} = ops{i}.d1_l;
154 d1_m{i} = ops{i}.d1_r;
155 end
156
157 %====== Assemble full operators ========
158 I_dim = speye(dim, dim);
159 RHO = spdiag(rho);
160 obj.RHO = RHO;
161 obj.RHOi = inv(RHO);
162 obj.RHOi_kron = kron(obj.RHOi, I_dim);
163
164 obj.D1 = cell(dim,1);
165 D2_temp = cell(dim,dim,dim);
166
167 % D1
168 obj.D1{1} = kron(D1{1},I{2});
169 obj.D1{2} = kron(I{1},D1{2});
170
171 % Boundary restriction operators
172 e_l = cell(dim,1);
173 e_r = cell(dim,1);
174 e_l{1} = kron(e_0{1}, I{2});
175 e_l{2} = kron(I{1}, e_0{2});
176 e_r{1} = kron(e_m{1}, I{2});
177 e_r{2} = kron(I{1}, e_m{2});
178
179 e_scalar_w = e_l{1};
180 e_scalar_e = e_r{1};
181 e_scalar_s = e_l{2};
182 e_scalar_n = e_r{2};
183
184 e_w = kron(e_scalar_w, I_dim);
185 e_e = kron(e_scalar_e, I_dim);
186 e_s = kron(e_scalar_s, I_dim);
187 e_n = kron(e_scalar_n, I_dim);
188
189 % Boundary restriction, shifted up to nBP-1 points
190 e_shifted_scalar_w = cell(nBP, 1);
191 e_shifted_scalar_e = cell(nBP, 1);
192 e_shifted_scalar_s = cell(nBP, 1);
193 e_shifted_scalar_n = cell(nBP, 1);
194 for i = 1:nBP-1
195 e_0_nBP = {0*e_0{1}, 0*e_0{2}};
196 e_m_nBP = {0*e_m{1}, 0*e_m{2}};
197 e_0_nBP{1}(i) = 1;
198 e_0_nBP{2}(i) = 1;
199 e_m_nBP{1}(end+1-i) = 1;
200 e_m_nBP{2}(end+1-i) = 1;
201 e_shifted_scalar_w{i} = kron(e_0_nBP{1}, I{2});
202 e_shifted_scalar_s{i} = kron(I{1}, e_0_nBP{2});
203 e_shifted_scalar_e{i} = kron(e_m_nBP{1}, I{2});
204 e_shifted_scalar_n{i} = kron(I{1}, e_m_nBP{2});
205 end
206
207 % Boundary derivatives
208 d1_l = cell(dim,1);
209 d1_r = cell(dim,1);
210 d1_l{1} = kron(d1_0{1}, I{2});
211 d1_l{2} = kron(I{1}, d1_0{2});
212 d1_r{1} = kron(d1_m{1}, I{2});
213 d1_r{2} = kron(I{1}, d1_m{2});
214
215 % E{i}^T picks out component i.
216 E = cell(dim,1);
217 I = speye(m_tot,m_tot);
218 for i = 1:dim
219 e = sparse(dim,1);
220 e(i) = 1;
221 E{i} = kron(I,e);
222 end
223 obj.E = E;
224
225 e1_w = (e_scalar_w'*E{1}')';
226 e1_e = (e_scalar_e'*E{1}')';
227 e1_s = (e_scalar_s'*E{1}')';
228 e1_n = (e_scalar_n'*E{1}')';
229
230 e2_w = (e_scalar_w'*E{2}')';
231 e2_e = (e_scalar_e'*E{2}')';
232 e2_s = (e_scalar_s'*E{2}')';
233 e2_n = (e_scalar_n'*E{2}')';
234
235 % D2
236 for j = 1:dim
237 for k = 1:dim
238 for l = 1:dim
239 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
240 end
241 end
242 end
243 ind = grid.funcToMatrix(g, 1:m_tot);
244
245 k = 1;
246 for r = 1:m(2)
247 p = ind(:,r);
248 for j = 1:dim
249 for l = 1:dim
250 coeff = C{k,j,k,l};
251 D_kk = D2{1}(coeff(p));
252 D2_temp{j,k,l}(p,p) = D_kk;
253 end
254 end
255 end
256
257 k = 2;
258 for r = 1:m(1)
259 p = ind(r,:);
260 for j = 1:dim
261 for l = 1:dim
262 coeff = C{k,j,k,l};
263 D_kk = D2{2}(coeff(p));
264 D2_temp{j,k,l}(p,p) = D_kk;
265 end
266 end
267 end
268
269 % Quadratures
270 obj.H = kron(H{1},H{2});
271 obj.Hi = inv(obj.H);
272 obj.H_w = H{2};
273 obj.H_e = H{2};
274 obj.H_s = H{1};
275 obj.H_n = H{1};
276 obj.H_1D = {H{1}, H{2}};
277
278 % Differentiation matrix D (without SAT)
279 D1 = obj.D1;
280 D = sparse(dim*m_tot,dim*m_tot);
281 for i = 1:dim
282 for j = 1:dim
283 for k = 1:dim
284 for l = 1:dim
285 if i == k
286 D = D + E{j}*D2_temp{j,k,l}*E{l}';
287 D2_temp{j,k,l} = [];
288 else
289 D = D + E{j}*(D1{i})*C_mat{i,j,k,l}*D1{k}*E{l}';
290 end
291 end
292 end
293 end
294 end
295 clear D2_temp;
296 D = obj.RHOi_kron*D;
297 obj.D = D;
298 clear D;
299 %=========================================%'
300
301 % Numerical traction operators for BC.
302 %
303 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
304 %
305 T_l = cell(dim,1);
306 T_r = cell(dim,1);
307 tau_l = cell(dim,1);
308 tau_r = cell(dim,1);
309
310 D1 = obj.D1;
311
312 d = @kroneckerDelta; % Kronecker delta
313 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
314
315 % Boundary j
316 for j = 1:dim
317 T_l{j} = cell(dim,dim);
318 T_r{j} = cell(dim,dim);
319 tau_l{j} = cell(dim,1);
320 tau_r{j} = cell(dim,1);
321
322 [~, n_l] = size(e_l{j});
323 [~, n_r] = size(e_r{j});
324
325 % Traction component i
326 for i = 1:dim
327 tau_l{j}{i} = sparse(dim*m_tot, n_l);
328 tau_r{j}{i} = sparse(dim*m_tot, n_r);
329
330 % Displacement component l
331 for l = 1:dim
332 T_l{j}{i,l} = sparse(m_tot, n_l);
333 T_r{j}{i,l} = sparse(m_tot, n_r);
334
335 % Derivative direction k
336 for k = 1:dim
337 T_l{j}{i,l} = T_l{j}{i,l} ...
338 - d(j,k)*(e_l{j}'*C_mat{j,i,k,l}*e_l{j}*d1_l{k}')'...
339 - db(j,k)*(e_l{j}'*C_mat{j,i,k,l}*D1{k})';
340 T_r{j}{i,l} = T_r{j}{i,l} ...
341 + d(j,k)*(e_r{j}'*C_mat{j,i,k,l}*e_r{j}*d1_r{k}')'...
342 + db(j,k)*(e_r{j}'*C_mat{j,i,k,l}*D1{k})';
343 end
344 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
345 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
346 end
347 end
348 end
349
350 % Traction tensors, T_ij
351 obj.T_w = T_l{1};
352 obj.T_e = T_r{1};
353 obj.T_s = T_l{2};
354 obj.T_n = T_r{2};
355
356 % Restriction operators
357 obj.e_w = e_w;
358 obj.e_e = e_e;
359 obj.e_s = e_s;
360 obj.e_n = e_n;
361
362 obj.e1_w = e1_w;
363 obj.e1_e = e1_e;
364 obj.e1_s = e1_s;
365 obj.e1_n = e1_n;
366
367 obj.e2_w = e2_w;
368 obj.e2_e = e2_e;
369 obj.e2_s = e2_s;
370 obj.e2_n = e2_n;
371
372 obj.e_scalar_w = e_scalar_w;
373 obj.e_scalar_e = e_scalar_e;
374 obj.e_scalar_s = e_scalar_s;
375 obj.e_scalar_n = e_scalar_n;
376
377
378 obj.e_shifted_scalar_w = e_shifted_scalar_w;
379 obj.e_shifted_scalar_e = e_shifted_scalar_e;
380 obj.e_shifted_scalar_s = e_shifted_scalar_s;
381 obj.e_shifted_scalar_n = e_shifted_scalar_n;
382
383 % First component of traction
384 obj.tau1_w = tau_l{1}{1};
385 obj.tau1_e = tau_r{1}{1};
386 obj.tau1_s = tau_l{2}{1};
387 obj.tau1_n = tau_r{2}{1};
388
389 % Second component of traction
390 obj.tau2_w = tau_l{1}{2};
391 obj.tau2_e = tau_r{1}{2};
392 obj.tau2_s = tau_l{2}{2};
393 obj.tau2_n = tau_r{2}{2};
394
395 % Traction vectors
396 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
397 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
398 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
399 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
400
401 % Kroneckered norms and coefficients
402 obj.Hi_kron = kron(obj.Hi, I_dim);
403
404 % Misc.
405 obj.m = m;
406 obj.h = h;
407 obj.order = order;
408 obj.grid = g;
409 obj.dim = dim;
410 obj.nBP = nBP;
411
412 end
413
414
415 % Closure functions return the operators applied to the own domain to close the boundary
416 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
417 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
418 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
419 % on the first component. Can also be e.g.
420 % {'normal', 'd'} or {'tangential', 't'} for conditions on
421 % tangential/normal component.
422 % data is a function returning the data that should be applied at the boundary.
423 % neighbour_scheme is an instance of Scheme that should be interfaced to.
424 % neighbour_boundary is a string specifying which boundary to interface to.
425
426 % For displacement bc:
427 % bc = {comp, 'd', dComps},
428 % where
429 % dComps = vector of components with displacement BC. Default: 1:dim.
430 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
431 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
432 default_arg('tuning', 1.0);
433
434 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
435 comp = bc{1};
436 type = bc{2};
437 if ischar(comp)
438 comp = obj.getComponent(comp, boundary);
439 end
440
441 e = obj.getBoundaryOperatorForScalarField('e', boundary);
442 e_sh = obj.getBoundaryOperatorForScalarField('e_shifted', boundary);
443 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
444 T = obj.getBoundaryTractionOperator(boundary);
445 [h11, th_R] = obj.getBorrowing(boundary);
446 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
447 nu = obj.getNormal(boundary);
448
449 E = obj.E;
450 Hi = obj.Hi;
451 RHOi = obj.RHOi;
452 C = obj.C;
453
454 dim = obj.dim;
455 m_tot = obj.grid.N();
456
457 % Preallocate
458 [~, col] = size(tau);
459 closure = sparse(dim*m_tot, dim*m_tot);
460 penalty = sparse(dim*m_tot, col);
461
462 j = comp;
463 switch type
464
465 % Dirichlet boundary condition
466 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
467
468 if numel(bc) >= 3
469 dComps = bc{3};
470 else
471 dComps = 1:dim;
472 end
473
474 % Loops over components that Dirichlet penalties end up on
475 % Y: symmetrizing part of penalty
476 % Z: symmetric part of penalty
477 % X = Y + Z.
478
479 % Nonsymmetric part goes on all components to
480 % yield traction in discrete energy rate
481 for i = 1:dim
482 Y = T{j,i}';
483 X = e*Y;
484 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
485 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
486 end
487
488 % Symmetric part only required on components with displacement BC.
489 % (Otherwise it's not symmetric.)
490
491 % 1. Compute beta such that C_bJbL_boundary - \beta*C_bJbL_i >= 0, i = 0, 1, ..., nBP.
492 [~, N] = size(e);
493 beta_vec = ones(N, 1);
494 b = obj.getComponent('normal', boundary);
495 C_boundary = cell(dim, dim, dim, dim);
496 for i = 1:dim
497 for j = 1:dim
498 for k = 1:dim
499 for l = 1:dim
500 C_boundary{i,j,k,l} = e'*C{i,j,k,l}*e ;
501 end
502 end
503 end
504 end
505
506 % Loop over shift inward from boundary
507 for m = 1:obj.nBP-1
508 C_shifted = cell(dim, dim, dim, dim);
509 for i = 1:dim
510 for j = 1:dim
511 for k = 1:dim
512 for l = 1:dim
513 C_shifted{i,j,k,l} = e_sh{m}'*C{i,j,k,l}*e_sh{m};
514 end
515 end
516 end
517 end
518
519 % Loop along boundary
520 for i = 1:N
521 C_boundary_mat = [C_boundary{b,1,b,1}(i,i), C_boundary{b,1,b,2}(i,i); ...
522 C_boundary{b,2,b,1}(i,i), C_boundary{b,2,b,2}(i,i)];
523
524 C_shifted_mat = [C_shifted{b,1,b,1}(i,i), C_shifted{b,1,b,2}(i,i); ...
525 C_shifted{b,2,b,1}(i,i), C_shifted{b,2,b,2}(i,i)];
526
527 beta = obj.computeBorrowingFromC(C_boundary_mat, C_shifted_mat);
528 beta_vec(i) = min(beta_vec(i), beta);
529 end
530 end
531 if min(beta_vec) < 1e-10
532 disp('WARNING! No borrowing from C seems to be possible.')
533 disp(['min(beta) = ' num2str(min(beta_vec))])
534 end
535 beta_inv = e*spdiag(1./beta_vec)*e';
536
537 % Compensate at corners
538 corner_weight = ones(N, 1);
539 corner_weight(1) = 2;
540 corner_weight(end) = 2;
541 corner_weight = e*spdiag(corner_weight)*e';
542
543
544 % 2. Build Z.
545 j = comp;
546 for i = dComps
547 Z = sparse(m_tot, m_tot);
548 for l = 1:dim
549 for k = 1:dim
550 Z = Z + nu(l)*C{l,i,k,j}*nu(k);
551 end
552 end
553 Z = -1.2*tuning*(corner_weight/h11 + beta_inv/th_R)*Z;
554 X = Z;
555 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
556 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
557 end
558
559 % Free boundary condition
560 case {'F','f','Free','free','traction','Traction','t','T'}
561 closure = closure - E{j}*RHOi*Hi*e*H_gamma*tau';
562 penalty = penalty + E{j}*RHOi*Hi*e*H_gamma;
563
564 % Unknown boundary condition
565 otherwise
566 error('No such boundary condition: type = %s',type);
567 end
568 end
569
570 % type Struct that specifies the interface coupling.
571 % Fields:
572 % -- tuning: penalty strength, defaults to 1.0
573 % -- interpolation: type of interpolation, default 'none'
574 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
575
576 defaultType.tuning = 1.0;
577 defaultType.interpolation = 'none';
578 default_struct('type', defaultType);
579
580 switch type.interpolation
581 case {'none', ''}
582 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
583 case {'op','OP'}
584 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
585 otherwise
586 error('Unknown type of interpolation: %s ', type.interpolation);
587 end
588 end
589
590 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
591 tuning = type.tuning;
592
593 % u denotes the solution in the own domain
594 % v denotes the solution in the neighbour domain
595
596 u = obj;
597 v = neighbour_scheme;
598
599 % Operators, u side
600 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
601 tau_u = u.getBoundaryOperator('tau', boundary);
602 h11_u = u.getBorrowing(boundary);
603 nu_u = u.getNormal(boundary);
604
605 E_u = u.E;
606 C_u = u.C;
607 m_tot_u = u.grid.N();
608
609 % Operators, v side
610 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
611 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
612 h11_v = v.getBorrowing(neighbour_boundary);
613 nu_v = v.getNormal(neighbour_boundary);
614
615 E_v = v.E;
616 C_v = v.C;
617 m_tot_v = v.grid.N();
618
619 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
620 flipFlag = false;
621 e_v_flip = e_v;
622 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
623 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ...
624 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ...
625 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ...
626 (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ...
627 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ...
628 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ...
629 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
630
631 flipFlag = true;
632 e_v_flip = fliplr(e_v);
633
634 t1 = tau_v(:,1:2:end-1);
635 t2 = tau_v(:,2:2:end);
636
637 t1 = fliplr(t1);
638 t2 = fliplr(t2);
639
640 tau_v(:,1:2:end-1) = t1;
641 tau_v(:,2:2:end) = t2;
642 end
643
644 % Operators that are only required for own domain
645 Hi = u.Hi_kron;
646 RHOi = u.RHOi_kron;
647 e_kron = u.getBoundaryOperator('e', boundary);
648 T_u = u.getBoundaryTractionOperator(boundary);
649
650 % Shared operators
651 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
652 H_gamma_kron = u.getBoundaryQuadrature(boundary);
653 dim = u.dim;
654
655 % Preallocate
656 [~, m_int] = size(H_gamma);
657 closure = sparse(dim*m_tot_u, dim*m_tot_u);
658 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
659
660 % ---- Continuity of displacement ------
661
662 % Y: symmetrizing part of penalty
663 % Z: symmetric part of penalty
664 % X = Y + Z.
665
666 % Loop over components to couple across interface
667 for j = 1:dim
668
669 % Loop over components that penalties end up on
670 for i = 1:dim
671 Y = 1/2*T_u{j,i}';
672 Z_u = sparse(m_int, m_int);
673 Z_v = sparse(m_int, m_int);
674 for l = 1:dim
675 for k = 1:dim
676 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
677 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
678 end
679 end
680
681 if flipFlag
682 Z_v = rot90(Z_v,2);
683 end
684
685 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
686 X = Y + Z*e_u';
687 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
688 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
689
690 end
691 end
692
693 % ---- Continuity of traction ------
694 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
695 penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
696
697 % ---- Multiply by inverse of density x quadraure ----
698 closure = RHOi*Hi*closure;
699 penalty = RHOi*Hi*penalty;
700
701 end
702
703 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
704 error('Non-conforming interfaces not implemented yet.');
705 end
706
707 % Computes the largest beta such that C_b - beta*C_s >= 0, using bisection.
708 function beta = computeBorrowingFromC(obj, C_b, C_s)
709 tol = 1e-8;
710 beta_min = 0;
711 beta_max = 1;
712 beta = 0.5;
713 err = (beta_max - beta_min)/2;
714 while err > tol
715 if min(eig(C_b - beta*C_s)) < -1e-14
716 % Tried to borrow too much. Reduce beta.
717 beta_max = beta;
718 beta = (beta_min + beta)/2;
719 else
720 % Increase beta
721 beta_min = beta;
722 beta = (beta_max + beta)/2;
723 end
724 err = err/2;
725 end
726 end
727
728 % Returns the component number that is the tangential/normal component
729 % at the specified boundary
730 function comp = getComponent(obj, comp_str, boundary)
731 assertIsMember(comp_str, {'normal', 'tangential'});
732 assertIsMember(boundary, {'w', 'e', 's', 'n'});
733
734 switch boundary
735 case {'w', 'e'}
736 switch comp_str
737 case 'normal'
738 comp = 1;
739 case 'tangential'
740 comp = 2;
741 end
742 case {'s', 'n'}
743 switch comp_str
744 case 'normal'
745 comp = 2;
746 case 'tangential'
747 comp = 1;
748 end
749 end
750 end
751
752 % Returns h11 for the boundary specified by the string boundary.
753 % op -- string
754 function [h11, theta_R] = getBorrowing(obj, boundary)
755 assertIsMember(boundary, {'w', 'e', 's', 'n'})
756
757 switch boundary
758 case {'w','e'}
759 h11 = obj.h11{1};
760 theta_R = obj.theta_R{1};
761 case {'s', 'n'}
762 h11 = obj.h11{2};
763 theta_R = obj.theta_R{2};
764 end
765 end
766
767 % Returns the outward unit normal vector for the boundary specified by the string boundary.
768 function nu = getNormal(obj, boundary)
769 assertIsMember(boundary, {'w', 'e', 's', 'n'})
770
771 switch boundary
772 case 'w'
773 nu = [-1,0];
774 case 'e'
775 nu = [1,0];
776 case 's'
777 nu = [0,-1];
778 case 'n'
779 nu = [0,1];
780 end
781 end
782
783 % Returns the boundary operator op for the boundary specified by the string boundary.
784 % op -- string
785 function o = getBoundaryOperator(obj, op, boundary)
786 assertIsMember(boundary, {'w', 'e', 's', 'n'})
787 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
788
789 switch op
790 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
791 o = obj.([op, '_', boundary]);
792 end
793
794 end
795
796 % Returns the boundary operator op for the boundary specified by the string boundary.
797 % op -- string
798 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
799 assertIsMember(boundary, {'w', 'e', 's', 'n'})
800 assertIsMember(op, {'e', 'e_shifted'})
801
802 switch op
803
804 case 'e'
805 o = obj.(['e_scalar', '_', boundary]);
806 case 'e_shifted'
807 o = obj.(['e_shifted_scalar', '_', boundary]);
808 end
809
810 end
811
812 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
813 % Formula: tau_i = T_ij u_j
814 % op -- string
815 function T = getBoundaryTractionOperator(obj, boundary)
816 assertIsMember(boundary, {'w', 'e', 's', 'n'})
817
818 T = obj.(['T', '_', boundary]);
819 end
820
821 % Returns square boundary quadrature matrix, of dimension
822 % corresponding to the number of boundary unknowns
823 %
824 % boundary -- string
825 function H = getBoundaryQuadrature(obj, boundary)
826 assertIsMember(boundary, {'w', 'e', 's', 'n'})
827
828 H = obj.getBoundaryQuadratureForScalarField(boundary);
829 I_dim = speye(obj.dim, obj.dim);
830 H = kron(H, I_dim);
831 end
832
833 % Returns square boundary quadrature matrix, of dimension
834 % corresponding to the number of boundary grid points
835 %
836 % boundary -- string
837 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
838 assertIsMember(boundary, {'w', 'e', 's', 'n'})
839
840 H_b = obj.(['H_', boundary]);
841 end
842
843 function N = size(obj)
844 N = obj.dim*prod(obj.m);
845 end
846 end
847 end