Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1331:60c875c18de3 feature/D2_boundary_opt

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