Mercurial > repos > public > sbplib

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