Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropicMixedStencil.m @ 1268:af9e52e97154 feature/poroelastic

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