Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinear.m @ 739:8efc04e97da4 feature/poroelastic

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