Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinear.m @ 943:21394c78c72e feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 15:24:36 -0800
parents 1f6b2fb69225
children cb4bfd0d81ea b8bd54332763
comparison
equal deleted inserted replaced
942:35701c85e356 943:21394c78c72e
1 classdef Elastic2dCurvilinear < scheme.Scheme
2
3 % Discretizes the elastic wave equation in curvilinear coordinates.
4 %
5 % Untransformed equation:
6 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
7 %
8 % Transformed equation:
9 % J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j
10 % + dk J b_jk mu b_il dl u_j
11 % + dk J b_jk mu b_jl dl u_i
12 % opSet should be cell array of opSets, one per dimension. This
13 % is useful if we have periodic BC in one direction.
14
15 properties
16 m % Number of points in each direction, possibly a vector
17 h % Grid spacing
18
19 grid
20 dim
21
22 order % Order of accuracy for the approximation
23
24 % Diagonal matrices for varible coefficients
25 LAMBDA % Variable coefficient, related to dilation
26 MU % Shear modulus, variable coefficient
27 RHO, RHOi % Density, variable
28
29 % Metric coefficients
30 b % Cell matrix of size dim x dim
31 J, Ji
32 beta % Cell array of scale factors
33
34 D % Total operator
35 D1 % First derivatives
36
37 % Second derivatives
38 D2_lambda
39 D2_mu
40
41 % Traction operators used for BC
42 T_l, T_r
43 tau_l, tau_r
44
45 H, Hi % Inner products
46 phi % Borrowing constant for (d1 - e^T*D1) from R
47 gamma % Borrowing constant for d1 from M
48 H11 % First element of H
49 e_l, e_r
50 d1_l, d1_r % Normal derivatives at the boundary
51 E % E{i}^T picks out component i
52
53 H_boundary_l, H_boundary_r % Boundary inner products
54
55 % Kroneckered norms and coefficients
56 RHOi_kron
57 Ji_kron, J_kron
58 Hi_kron, H_kron
59 end
60
61 methods
62
63 function obj = Elastic2dCurvilinear(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
64 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
65 default_arg('lambda_fun', @(x,y) 0*x+1);
66 default_arg('mu_fun', @(x,y) 0*x+1);
67 default_arg('rho_fun', @(x,y) 0*x+1);
68 dim = 2;
69
70 lambda = grid.evalOn(g, lambda_fun);
71 mu = grid.evalOn(g, mu_fun);
72 rho = grid.evalOn(g, rho_fun);
73 m = g.size();
74 obj.m = m;
75 m_tot = g.N();
76
77 % 1D operators
78 ops = cell(dim,1);
79 for i = 1:dim
80 ops{i} = opSet{i}(m(i), {0, 1}, order);
81 end
82
83 % Borrowing constants
84 for i = 1:dim
85 beta = ops{i}.borrowing.R.delta_D;
86 obj.H11{i} = ops{i}.borrowing.H11;
87 obj.phi{i} = beta/obj.H11{i};
88 obj.gamma{i} = ops{i}.borrowing.M.d1;
89 end
90
91 I = cell(dim,1);
92 D1 = cell(dim,1);
93 D2 = cell(dim,1);
94 H = cell(dim,1);
95 Hi = cell(dim,1);
96 e_l = cell(dim,1);
97 e_r = cell(dim,1);
98 d1_l = cell(dim,1);
99 d1_r = cell(dim,1);
100
101 for i = 1:dim
102 I{i} = speye(m(i));
103 D1{i} = ops{i}.D1;
104 D2{i} = ops{i}.D2;
105 H{i} = ops{i}.H;
106 Hi{i} = ops{i}.HI;
107 e_l{i} = ops{i}.e_l;
108 e_r{i} = ops{i}.e_r;
109 d1_l{i} = ops{i}.d1_l;
110 d1_r{i} = ops{i}.d1_r;
111 end
112
113 %====== Assemble full operators ========
114
115 % Variable coefficients
116 LAMBDA = spdiag(lambda);
117 obj.LAMBDA = LAMBDA;
118 MU = spdiag(mu);
119 obj.MU = MU;
120 RHO = spdiag(rho);
121 obj.RHO = RHO;
122 obj.RHOi = inv(RHO);
123
124 % Allocate
125 obj.D1 = cell(dim,1);
126 obj.D2_lambda = cell(dim,dim,dim);
127 obj.D2_mu = cell(dim,dim,dim);
128 obj.e_l = cell(dim,1);
129 obj.e_r = cell(dim,1);
130 obj.d1_l = cell(dim,1);
131 obj.d1_r = cell(dim,1);
132
133 % D1
134 obj.D1{1} = kron(D1{1},I{2});
135 obj.D1{2} = kron(I{1},D1{2});
136
137 % -- Metric coefficients ----
138 coords = g.points();
139 x = coords(:,1);
140 y = coords(:,2);
141
142 % Use non-periodic difference operators for metric even if opSet is periodic.
143 xmax = max(ops{1}.x);
144 ymax = max(ops{2}.x);
145 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
146 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
147 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
148 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
149
150 x_xi = D1Metric{1}*x;
151 x_eta = D1Metric{2}*x;
152 y_xi = D1Metric{1}*y;
153 y_eta = D1Metric{2}*y;
154
155 J = x_xi.*y_eta - x_eta.*y_xi;
156
157 b = cell(dim,dim);
158 b{1,1} = y_eta./J;
159 b{1,2} = -x_eta./J;
160 b{2,1} = -y_xi./J;
161 b{2,2} = x_xi./J;
162
163 % Scale factors for boundary integrals
164 beta = cell(dim,1);
165 beta{1} = sqrt(x_eta.^2 + y_eta.^2);
166 beta{2} = sqrt(x_xi.^2 + y_xi.^2);
167
168 J = spdiag(J);
169 Ji = inv(J);
170 for i = 1:dim
171 beta{i} = spdiag(beta{i});
172 for j = 1:dim
173 b{i,j} = spdiag(b{i,j});
174 end
175 end
176 obj.J = J;
177 obj.Ji = Ji;
178 obj.b = b;
179 obj.beta = beta;
180 %----------------------------
181
182 % Boundary operators
183 obj.e_l{1} = kron(e_l{1},I{2});
184 obj.e_l{2} = kron(I{1},e_l{2});
185 obj.e_r{1} = kron(e_r{1},I{2});
186 obj.e_r{2} = kron(I{1},e_r{2});
187
188 obj.d1_l{1} = kron(d1_l{1},I{2});
189 obj.d1_l{2} = kron(I{1},d1_l{2});
190 obj.d1_r{1} = kron(d1_r{1},I{2});
191 obj.d1_r{2} = kron(I{1},d1_r{2});
192
193 % D2
194 for i = 1:dim
195 for j = 1:dim
196 for k = 1:dim
197 obj.D2_lambda{i,j,k} = sparse(m_tot);
198 obj.D2_mu{i,j,k} = sparse(m_tot);
199 end
200 end
201 end
202 ind = grid.funcToMatrix(g, 1:m_tot);
203
204 % x-dir
205 for i = 1:dim
206 for j = 1:dim
207 for k = 1
208
209 coeff_lambda = J*b{i,k}*b{j,k}*lambda;
210 coeff_mu = J*b{j,k}*b{i,k}*mu;
211
212 for col = 1:m(2)
213 D_lambda = D2{1}(coeff_lambda(ind(:,col)));
214 D_mu = D2{1}(coeff_mu(ind(:,col)));
215
216 p = ind(:,col);
217 obj.D2_lambda{i,j,k}(p,p) = D_lambda;
218 obj.D2_mu{i,j,k}(p,p) = D_mu;
219 end
220
221 end
222 end
223 end
224
225 % y-dir
226 for i = 1:dim
227 for j = 1:dim
228 for k = 2
229
230 coeff_lambda = J*b{i,k}*b{j,k}*lambda;
231 coeff_mu = J*b{j,k}*b{i,k}*mu;
232
233 for row = 1:m(1)
234 D_lambda = D2{2}(coeff_lambda(ind(row,:)));
235 D_mu = D2{2}(coeff_mu(ind(row,:)));
236
237 p = ind(row,:);
238 obj.D2_lambda{i,j,k}(p,p) = D_lambda;
239 obj.D2_mu{i,j,k}(p,p) = D_mu;
240 end
241
242 end
243 end
244 end
245
246 % Quadratures
247 obj.H = kron(H{1},H{2});
248 obj.Hi = inv(obj.H);
249 obj.H_boundary_l = cell(dim,1);
250 obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2};
251 obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1};
252 obj.H_boundary_r = cell(dim,1);
253 obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2};
254 obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1};
255
256 % E{i}^T picks out component i.
257 E = cell(dim,1);
258 I = speye(m_tot,m_tot);
259 for i = 1:dim
260 e = sparse(dim,1);
261 e(i) = 1;
262 E{i} = kron(I,e);
263 end
264 obj.E = E;
265
266 % Differentiation matrix D (without SAT)
267 D2_lambda = obj.D2_lambda;
268 D2_mu = obj.D2_mu;
269 D1 = obj.D1;
270 D = sparse(dim*m_tot,dim*m_tot);
271 d = @kroneckerDelta; % Kronecker delta
272 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
273 for i = 1:dim
274 for j = 1:dim
275 for k = 1:dim
276 for l = 1:dim
277 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_lambda{i,j,k}*E{j}' + ...
278 db(k,l)*D1{k}*J*b{i,k}*b{j,l}*LAMBDA*D1{l}*E{j}' ...
279 );
280
281 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{i,j,k}*E{j}' + ...
282 db(k,l)*D1{k}*J*b{j,k}*b{i,l}*MU*D1{l}*E{j}' ...
283 );
284
285 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{j,j,k}*E{i}' + ...
286 db(k,l)*D1{k}*J*b{j,k}*b{j,l}*MU*D1{l}*E{i}' ...
287 );
288
289 end
290 end
291 end
292 end
293 obj.D = D;
294 %=========================================%
295
296 % Numerical traction operators for BC.
297 % Because d1 =/= e0^T*D1, the numerical tractions are different
298 % at every boundary.
299 T_l = cell(dim,1);
300 T_r = cell(dim,1);
301 tau_l = cell(dim,1);
302 tau_r = cell(dim,1);
303 % tau^{j}_i = sum_k T^{j}_{ik} u_k
304
305 d1_l = obj.d1_l;
306 d1_r = obj.d1_r;
307 e_l = obj.e_l;
308 e_r = obj.e_r;
309
310 % Loop over boundaries
311 for j = 1:dim
312 T_l{j} = cell(dim,dim);
313 T_r{j} = cell(dim,dim);
314 tau_l{j} = cell(dim,1);
315 tau_r{j} = cell(dim,1);
316
317 % Loop over components
318 for i = 1:dim
319 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
320 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
321
322 % Loop over components that T_{ik}^{(j)} acts on
323 for k = 1:dim
324
325 T_l{j}{i,k} = sparse(m_tot,m_tot);
326 T_r{j}{i,k} = sparse(m_tot,m_tot);
327
328 for m = 1:dim
329 for l = 1:dim
330 T_l{j}{i,k} = T_l{j}{i,k} + ...
331 -d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
332 -d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
333 -d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m});
334
335 T_r{j}{i,k} = T_r{j}{i,k} + ...
336 d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
337 d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
338 d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m});
339 end
340 end
341
342 T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k};
343 T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k};
344
345 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
346 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
347 end
348
349 end
350 end
351 obj.T_l = T_l;
352 obj.T_r = T_r;
353 obj.tau_l = tau_l;
354 obj.tau_r = tau_r;
355
356 % Kroneckered norms and coefficients
357 I_dim = speye(dim);
358 obj.RHOi_kron = kron(obj.RHOi, I_dim);
359 obj.Ji_kron = kron(obj.Ji, I_dim);
360 obj.Hi_kron = kron(obj.Hi, I_dim);
361 obj.H_kron = kron(obj.H, I_dim);
362 obj.J_kron = kron(obj.J, I_dim);
363
364 % Misc.
365 obj.h = g.scaling();
366 obj.order = order;
367 obj.grid = g;
368 obj.dim = dim;
369
370 end
371
372
373 % Closure functions return the operators applied to the own domain to close the boundary
374 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
375 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
376 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
377 % on the first component.
378 % data is a function returning the data that should be applied at the boundary.
379 % neighbour_scheme is an instance of Scheme that should be interfaced to.
380 % neighbour_boundary is a string specifying which boundary to interface to.
381 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
382 default_arg('tuning', 1.2);
383
384 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
385 comp = bc{1};
386 type = bc{2};
387
388 % j is the coordinate direction of the boundary
389 j = obj.get_boundary_number(boundary);
390 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
391
392 E = obj.E;
393 Hi = obj.Hi;
394 LAMBDA = obj.LAMBDA;
395 MU = obj.MU;
396 RHOi = obj.RHOi;
397 Ji = obj.Ji;
398
399 dim = obj.dim;
400 m_tot = obj.grid.N();
401
402 % Preallocate
403 closure = sparse(dim*m_tot, dim*m_tot);
404 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
405
406 % Loop over components that we (potentially) have different BC on
407 k = comp;
408 switch type
409
410 % Dirichlet boundary condition
411 case {'D','d','dirichlet','Dirichlet'}
412
413 phi = obj.phi{j};
414 h = obj.h(j);
415 h11 = obj.H11{j}*h;
416 gamma = obj.gamma{j};
417
418 a_lambda = dim/h11 + 1/(h11*phi);
419 a_mu_i = 2/(gamma*h);
420 a_mu_ij = 2/h11 + 1/(h11*phi);
421
422 d = @kroneckerDelta; % Kronecker delta
423 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
424 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
425 + d(i,j)* a_mu_i*MU ...
426 + db(i,j)*a_mu_ij*MU );
427
428 % Loop over components that Dirichlet penalties end up on
429 for i = 1:dim
430 C = T{k,i};
431 A = -d(i,k)*alpha(i,j);
432 B = A + C;
433 closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' );
434 penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma;
435 end
436
437 % Free boundary condition
438 case {'F','f','Free','free','traction','Traction','t','T'}
439 closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} );
440 penalty = penalty + E{k}*RHOi*Ji*Hi*e*H_gamma;
441
442 % Unknown boundary condition
443 otherwise
444 error('No such boundary condition: type = %s',type);
445 end
446 end
447
448 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
449 % u denotes the solution in the own domain
450 % v denotes the solution in the neighbour domain
451 % Operators without subscripts are from the own domain.
452 error('Not implemented');
453 tuning = 1.2;
454
455 % j is the coordinate direction of the boundary
456 j = obj.get_boundary_number(boundary);
457 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
458
459 % Get boundary operators
460 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
461 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
462
463 % Operators and quantities that correspond to the own domain only
464 Hi = obj.Hi;
465 RHOi = obj.RHOi;
466 dim = obj.dim;
467
468 %--- Other operators ----
469 m_tot_u = obj.grid.N();
470 E = obj.E;
471 LAMBDA_u = obj.LAMBDA;
472 MU_u = obj.MU;
473 lambda_u = e'*LAMBDA_u*e;
474 mu_u = e'*MU_u*e;
475
476 m_tot_v = neighbour_scheme.grid.N();
477 E_v = neighbour_scheme.E;
478 LAMBDA_v = neighbour_scheme.LAMBDA;
479 MU_v = neighbour_scheme.MU;
480 lambda_v = e_v'*LAMBDA_v*e_v;
481 mu_v = e_v'*MU_v*e_v;
482 %-------------------------
483
484 % Borrowing constants
485 phi_u = obj.phi{j};
486 h_u = obj.h(j);
487 h11_u = obj.H11{j}*h_u;
488 gamma_u = obj.gamma{j};
489
490 phi_v = neighbour_scheme.phi{j_v};
491 h_v = neighbour_scheme.h(j_v);
492 h11_v = neighbour_scheme.H11{j_v}*h_v;
493 gamma_v = neighbour_scheme.gamma{j_v};
494
495 % E > sum_i 1/(2*alpha_ij)*(tau_i)^2
496 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)
497 th1 = h11/(2*dim);
498 th2 = h11*phi/2;
499 th3 = h*gamma;
500 a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3);
501 a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3);
502 alpha_ii = a1 + sqrt(a2 + a1^2);
503
504 alpha_ij = mu*(2/h11 + 1/(phi*h11));
505 end
506
507 [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
508 [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);
509 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
510 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
511
512 d = @kroneckerDelta; % Kronecker delta
513 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
514 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
515
516 % Preallocate
517 closure = sparse(dim*m_tot_u, dim*m_tot_u);
518 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
519
520 % Loop over components that penalties end up on
521 for i = 1:dim
522 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
523 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
524
525 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
526 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
527
528 % Loop over components that we have interface conditions on
529 for k = 1:dim
530 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
531 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
532 end
533 end
534 end
535
536 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
537 function [j, nj] = get_boundary_number(obj, boundary)
538
539 switch boundary
540 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
541 j = 1;
542 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
543 j = 2;
544 otherwise
545 error('No such boundary: boundary = %s',boundary);
546 end
547
548 switch boundary
549 case {'w','W','west','West','s','S','south','South'}
550 nj = -1;
551 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
552 nj = 1;
553 end
554 end
555
556 % Returns the boundary operator op for the boundary specified by the string boundary.
557 % op: may be a cell array of strings
558 function [varargout] = get_boundary_operator(obj, op, boundary)
559
560 switch boundary
561 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
562 j = 1;
563 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
564 j = 2;
565 otherwise
566 error('No such boundary: boundary = %s',boundary);
567 end
568
569 if ~iscell(op)
570 op = {op};
571 end
572
573 for i = 1:length(op)
574 switch op{i}
575 case 'e'
576 switch boundary
577 case {'w','W','west','West','s','S','south','South'}
578 varargout{i} = obj.e_l{j};
579 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
580 varargout{i} = obj.e_r{j};
581 end
582 case 'd'
583 switch boundary
584 case {'w','W','west','West','s','S','south','South'}
585 varargout{i} = obj.d1_l{j};
586 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
587 varargout{i} = obj.d1_r{j};
588 end
589 case 'H'
590 switch boundary
591 case {'w','W','west','West','s','S','south','South'}
592 varargout{i} = obj.H_boundary_l{j};
593 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
594 varargout{i} = obj.H_boundary_r{j};
595 end
596 case 'T'
597 switch boundary
598 case {'w','W','west','West','s','S','south','South'}
599 varargout{i} = obj.T_l{j};
600 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
601 varargout{i} = obj.T_r{j};
602 end
603 case 'tau'
604 switch boundary
605 case {'w','W','west','West','s','S','south','South'}
606 varargout{i} = obj.tau_l{j};
607 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
608 varargout{i} = obj.tau_r{j};
609 end
610 otherwise
611 error(['No such operator: operator = ' op{i}]);
612 end
613 end
614
615 end
616
617 function N = size(obj)
618 N = obj.dim*prod(obj.m);
619 end
620 end
621 end